Fork me on GitHub

source: svn/trunk/src/LHCOConverter.cc@ 302

Last change on this file since 302 was 288, checked in by Xavier Rouby, 16 years ago

charge added
File size: 8.8 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \----------------------------------------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
31
32// local includes
33#include "LHCOConverter.h"
34#include "ExRootTreeReader.h"
35#include "BlockClasses.h"
36
37// root includes
38#include "TFile.h"
39#include "TChain.h"
40#include "TClonesArray.h"
41
42// c++ includes
43#include <iostream>
44#include <string>
45#include <fstream>
46#include <sstream>
47#include <iomanip>
48using namespace std;
49
50
51//------------------------------------------------------------------------------
52extern const float UNDEFINED;
53
54LHCOConverter::LHCOConverter(const string& inputroot, const string& inputlog):
55 inputfilename_(inputroot), inputlogname_(inputlog), outputfilename_(inputroot), valid_(true) {
56 if (inputfilename_ == "") {valid_ = false;}
57 else {
58 TFile ftest(inputfilename_.c_str(),"READ");
59 if (!ftest.IsOpen()) { cout << "ERROR: file " << inputfilename_ << " could not be opened\n"; valid_=false;
60 }
61 else {
62 ftest.Close();
63 const unsigned int r_find = outputfilename_.rfind(".root");
64 if(r_find > outputfilename_.size()) {
65 cout << "ERROR: the extension of the input file is not '.root'. Exiting...\n"; valid_=false;
66 }
67 else {
68 outputfilename_.replace(r_find,5,".lhco");
69 ofstream outfile( outputfilename_.c_str());
70 outfile.close();
71 cout << "INFO: " << outputfilename_ << " has been created\n";
72 }
73 }
74 }
75}
76
77
78
79void LHCOConverter::CopyRunLogFile() {
80 if(!valid_) {/*cout << "Nothing done\n";*/ return; }
81 if(inputlogname_ == "") {
82 // re-creates the logrun file name from the rootfile name
83 inputlogname_ = inputfilename_;
84 const unsigned int r_find = inputlogname_.rfind(".root");
85 inputlogname_.replace(r_find,5,"_run.log");
86 }
87 ifstream infile( inputlogname_.c_str() );
88 if (!infile.good()) { cout << "Warning: no input logfile found\n"; return; }
89
90 ofstream outfile( outputfilename_.c_str(),ios_base::app);
91
92 // else, if find is found
93 string linereader;
94 while ( getline(infile,linereader) ) {
95 outfile << " # " << linereader << endl;
96 }
97 outfile << " #" << endl;
98 infile.close();
99 outfile.close();
100}
101
102
103
104void LHCOConverter::ConvertExRootAnalysisToLHCO() {
105 if(!valid_) {/*cout << "Nothing done\n";*/ return; }
106
107 //output file
108 ofstream outfile( outputfilename_.c_str(),ios_base::app);
109 //# typ eta phi pt jmas ntrk btag had/em dum1 dum2
110 outfile << " ## More info on LHCO files: http://cp3wks05.fynu.ucl.ac.be/Manual/lhco.html\n #" << endl;
111 outfile << " # typ eta phi pt jmas ntrk btag had/em dum1 dum2" << endl;
112
113 // input files
114 TChain analysisChain("Analysis");
115 TChain triggerChain("Trigger");
116 analysisChain.Add(inputfilename_.c_str());
117 triggerChain.Add(inputfilename_.c_str());
118 ExRootTreeReader *analysisTree = new ExRootTreeReader(&analysisChain);
119 ExRootTreeReader *triggerTree = new ExRootTreeReader(&triggerChain);
120
121 TClonesArray *branchPhoton = analysisTree->UseBranch("Photon");
122 TClonesArray *branchElectron = analysisTree->UseBranch("Electron");
123 TClonesArray *branchMuon = analysisTree->UseBranch("Muon");
124 TClonesArray *branchTauJet = analysisTree->UseBranch("TauJet");
125 TClonesArray *branchJet = analysisTree->UseBranch("Jet");
126 TClonesArray *branchETmis = analysisTree->UseBranch("ETmis");
127
128 Long64_t Nevents = analysisTree->GetEntries();
129 cout << "** TTree 'Analysis' contains " << Nevents << " events" << endl;
130 Nevents = min(Nevents,triggerTree->GetEntries());
131 if (Nevents != analysisTree->GetEntries())
132 cout << "** WARNING: not the same number of entries in 'Analysis' and **\n** 'Trigger' trees\n";
133 for(Long64_t event = 0; event < Nevents; ++event) {
134 analysisTree->ReadEntry(event);
135 triggerTree->ReadEntry(event);
136 unsigned int triginfo = 0;
137 unsigned int line =0;
138
139 outfile << setw(3) << line++ << setw(4) << " " << setw(7) << event << setw(7) << triginfo << endl;
140 outfile.close();
141
142 // 0 photon data
143 //outfile << BranchReader(branchPhoton,line,lhcoPhotonID);
144 BranchReader(branchPhoton,line,lhcoPhotonID);
145
146 // 1 electron/positron data
147 BranchReader(branchElectron,line,lhcoElectronID);
148
149 // 2 muon data
150 BranchReader(branchMuon,line,lhcoMuonID);
151
152 // 3 tau-jets
153 BranchReader(branchTauJet,line,lhcoTauJetID);
154
155 // 4 jets
156 BranchReader(branchJet,line,lhcoJetID);
157
158 // 6 MET
159 BranchReader(branchETmis,line,lhcoETmisID);
160
161 } // event loop
162
163 delete triggerTree;
164 delete analysisTree;
165}
166
167
168//ostringstream LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
169void LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
170 //ostringstream outfile;
171 ofstream outfile( outputfilename_.c_str(),ios_base::app);
172 unsigned int branch_size = branch->GetEntries();
173 TRootParticle * particle = 0;
174 for (unsigned int i=0; i< branch_size; i++) {
175 particle = (TRootParticle*) branch->At(i);
176 double temp = sqrt(particle->PT*particle->PT + particle->Pz * particle->Pz);
177 double jmass = (particle->E - temp) * (particle->E + temp); // prevents some numerical sensitivity
178 TLorentzVector pmu;
179 pmu.SetPtEtaPhiE(particle->PT, particle->Eta, particle->Phi, particle->E);
180 jmass = pmu.M();
181 float ntrk = 0.0;
182 if(lhcoID == lhcoElectronID) { TRootElectron elec(*((TRootElectron*) branch->At(i))); ntrk = elec.Charge; }
183 else if (lhcoID == lhcoMuonID) { TRootMuon muon(*((TRootMuon*) branch->At(i))); ntrk = muon.Charge; }
184 else if (lhcoID == lhcoTauJetID) { TRootTauJet taujet(*((TRootTauJet*) branch->At(i))); ntrk = taujet.Charge; }
185 float btag =0;
186 double ratioE = 0;
187 if(lhcoID != lhcoETmisID) {
188 outfile << fixed << setprecision(3)
189 << setw(3) << line++ // line counter
190 << setw(4) << lhcoID << " " // object ID in lhco format
191 << setw(7) << particle->Eta <<" "
192 << setw(7) << particle->Phi <<" "
193 << setw(7) << particle->PT <<" "
194 << setw(7) << jmass <<" " // invariant mass
195 << setw(7) << ntrk <<" " // number of tracks
196 << setw(7) << btag <<" "
197 << setw(7) << ratioE <<" " // E_had/E_em
198 << setw(7) << 0.0 << setw(7) << 0.0 // dummy/dummy
199 << endl;
200 } else {
201 outfile << fixed << setprecision(3)
202 << setw(3) << line++ // line counter
203 << setw(4) << lhcoID << " " // object ID in lhco format
204 << setw(7) << 0.000 <<" "
205 << setw(7) << particle->Phi <<" "
206 << setw(7) << particle->PT <<" "
207 << setw(7) << 0. <<" " // invariant mass
208 << setw(7) << 0 <<" " // number of tracks
209 << setw(7) << 0. <<" "
210 << setw(7) << 0. <<" " // E_had/E_em
211 << setw(7) << 0.0 << setw(7) << 0.0 // dummy/dummy
212 << endl;
213 }
214 }
215 outfile.close();
216 // return outfile;
217}
Note: See TracBrowser for help on using the repository browser.