Fork me on GitHub

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

Last change on this file since 857 was 557, checked in by severine ovyn, 15 years ago

SumEt and SumPt + change ISOL names
File size: 12.1 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34// local includes
35#include "LHCOConverter.h"
36#include "ExRootTreeReader.h"
37#include "BlockClasses.h"
38#include "SmearUtil.h"
39
40// root includes
41#include "TFile.h"
42#include "TChain.h"
43#include "TClonesArray.h"
44
45// c++ includes
46#include <iostream>
47#include <string>
48#include <fstream>
49#include <sstream>
50#include <iomanip>
51using namespace std;
52
53
54//------------------------------------------------------------------------------
55extern const float UNDEFINED;
56
57LHCOConverter::LHCOConverter(const string& inputroot, const string& inputlog):
58 inputfilename_(inputroot), inputlogname_(inputlog), outputfilename_(inputroot), valid_(true) {
59 if (inputfilename_ == "") {valid_ = false;}
60 else {
61 TFile ftest(inputfilename_.c_str(),"READ");
62 if (!ftest.IsOpen()) { cout << "ERROR: file " << inputfilename_ << " could not be opened\n"; valid_=false;
63 }
64 else {
65 ftest.Close();
66 const unsigned int r_find = outputfilename_.rfind(".root");
67 if(r_find > outputfilename_.size()) {
68 cerr << "ERROR: the extension of the input file is not '.root'. Exiting...\n"; valid_=false;
69 }
70 else {
71 outputfilename_.replace(r_find,5,"_events.lhco");
72 ofstream outfile( outputfilename_.c_str());
73 outfile.close();
74 cout << left << setw(20) <<"** INFO: "<<""
75 << left << setw(29) << outputfilename_ <<""
76 << right << setw(20) <<"created **"<<endl;
77
78 }
79 }
80 }
81}
82
83
84
85void LHCOConverter::CopyRunLogFile() {
86 if(!valid_) {/*cout << "Nothing done\n";*/ return; }
87 if(inputlogname_ == "") {
88 // re-creates the logrun file name from the rootfile name
89 inputlogname_ = inputfilename_;
90 const unsigned int r_find = inputlogname_.rfind(".root");
91 inputlogname_.replace(r_find,5,"_run.log");
92 }
93 ifstream infile( inputlogname_.c_str() );
94 if (!infile.good()) { cout << "Warning: no input logfile found\n"; return; }
95
96 ofstream outfile( outputfilename_.c_str(),ios_base::app);
97
98 // else, if find is found
99 string linereader;
100 while ( getline(infile,linereader) ) {
101 outfile << " # " << linereader << endl;
102 }
103 outfile << " #" << endl;
104 infile.close();
105 outfile.close();
106}
107
108
109
110void LHCOConverter::ConvertExRootAnalysisToLHCO() {
111 if(!valid_) {/*cout << "Nothing done\n";*/ return; }
112
113 //output file
114 ofstream outfile( outputfilename_.c_str(),ios_base::app);
115 //# typ eta phi pt jmas ntrk btag had/em dum1 dum2
116 outfile << " ## More info on LHCO files: http://cp3wks05.fynu.ucl.ac.be/Manual/lhco.html\n #" << endl;
117 outfile << " # typ eta phi pt jmas ntrk btag had/em dum1 dum2" << endl;
118
119 // input files
120 TChain analysisChain("Analysis");
121 TChain triggerChain("Trigger");
122 analysisChain.Add(inputfilename_.c_str());
123 triggerChain.Add(inputfilename_.c_str());
124 ExRootTreeReader *analysisTree = new ExRootTreeReader(&analysisChain);
125 ExRootTreeReader *triggerTree = new ExRootTreeReader(&triggerChain);
126
127 TClonesArray *branchPhoton = analysisTree->UseBranch("Photon");
128 TClonesArray *branchElectron = analysisTree->UseBranch("Electron");
129 TClonesArray *branchMuon = analysisTree->UseBranch("Muon");
130 TClonesArray *branchTauJet = analysisTree->UseBranch("TauJet");
131 TClonesArray *branchJet = analysisTree->UseBranch("Jet");
132 TClonesArray *branchETmis = analysisTree->UseBranch("ETmis");
133
134 Long64_t Nevents = analysisTree->GetEntries();
135 Nevents = min(Nevents,triggerTree->GetEntries());
136 if (Nevents != analysisTree->GetEntries())
137 cout << "** WARNING: not the same number of entries in 'Analysis' and **\n** 'Trigger' trees\n";
138 for(Long64_t event = 0; event < Nevents; ++event) {
139 analysisTree->ReadEntry(event);
140 triggerTree->ReadEntry(event);
141 unsigned int triginfo = 0;
142 unsigned int line =0;
143
144 outfile << setw(3) << line++ << setw(4) << " " << setw(7) << event << setw(7) << triginfo << endl;
145
146 // 0 photon data
147 BranchReader(branchPhoton,line,lhcoPhotonID);
148
149 // 1 electron/positron data
150 BranchReader(branchElectron,line,lhcoElectronID);
151
152 // 2 muon data
153 BranchReaderMuon(branchMuon,branchJet,line);
154
155 // 3 tau-jets
156 BranchReader(branchTauJet,line,lhcoTauJetID);
157
158 // 4 jets
159 BranchReader(branchJet,line,lhcoJetID);
160
161 // 6 MET
162 BranchReaderETmis(branchETmis,line);
163
164 } // event loop
165
166 outfile.close();
167 delete triggerTree;
168 delete analysisTree;
169}
170
171
172//ostringstream LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
173void LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
174 //ostringstream outfile;
175 ofstream outfile( outputfilename_.c_str(),ios_base::app);
176 unsigned int branch_size = branch->GetEntries();
177 TRootParticle * particle = 0;
178
179 // parsing the branch
180 for (unsigned int i=0; i< branch_size; i++) {
181 particle = (TRootParticle*) branch->At(i);
182 TLorentzVector pmu;
183 pmu.SetPtEtaPhiE(particle->PT, particle->Eta, particle->Phi, particle->E);
184 float jmass = pmu.M();
185//if(jmass < 0) cout << "jmass negatif" << particle->PT << " " << particle->Eta << " " << particle->Phi << " " << particle->E << " " << jmass << endl;
186 float ntrk = 0.0;
187 float btag =0;
188 double ratioE = 0;
189
190 switch(lhcoID) {
191 case lhcoPhotonID: {
192 TRootPhoton gam(*((TRootPhoton*) branch->At(i)));
193 jmass=0; ratioE = gam.EHoverEE;
194 } break;
195 case lhcoElectronID: {
196 TRootElectron elec(*((TRootElectron*) branch->At(i)));
197 ntrk = elec.Charge; ratioE = elec.EHoverEE;
198 } break;
199 case lhcoTauJetID: {
200 TRootTauJet taujet(*((TRootTauJet*) branch->At(i)));
201 ntrk = taujet.Charge; ratioE = taujet.EHoverEE;
202 } break;
203 case lhcoJetID: {
204 TRootJet jet(*((TRootJet*) branch->At(i)));
205 ntrk = jet.NTracks; btag = jet.Btag; ratioE = jet.EHoverEE;
206 } break;
207 case lhcoMuonID: {
208 cout << "ERROR: LHCOConverter: for muons use LHCOConverter::BranchReaderMuon!\n";
209 } break;
210 case lhcoETmisID: {
211 cout << "ERROR: LHCOConverter: for ETmis use LHCOConverter::BranchReaderETmis!\n";
212 } break;
213
214 default: break;
215 } // switch
216
217 if(lhcoID != lhcoETmisID) {
218 outfile << fixed << setprecision(3)
219 << setw(3) << line++ // line counter
220 << setw(4) << lhcoID << " " // object ID in lhco format
221 << setw(7) << particle->Eta <<" "
222 << setw(7) << particle->Phi <<" "
223 << setw(7) << particle->PT <<" "
224 << setw(7) << jmass <<" " // invariant mass
225 << setw(7) << ntrk <<" " // number of tracks
226 << setw(7) << btag <<" "
227 << setw(7) << ratioE <<" " // E_had/E_em
228 << setw(7) << 0.0 << setw(7) << 0.0 // dummy/dummy
229 << endl;
230 }
231 } // branch parsing
232
233 outfile.close();
234 // return outfile;
235}
236
237//ostringstream LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
238void LHCOConverter::BranchReaderETmis(const TClonesArray* branch, unsigned int& line) const {
239 //ostringstream outfile;
240 ofstream outfile( outputfilename_.c_str(),ios_base::app);
241 TRootETmis * particle = 0;
242 particle = (TRootETmis*) branch->At(0);
243 outfile << fixed << setprecision(3)
244 << setw(3) << line++ // line counter
245 << setw(4) << 6 << " " // object ID in lhco format
246 << setw(7) << 0.000 <<" "
247 << setw(7) << particle->Phi <<" "
248 << setw(7) << particle->ET <<" "
249 << setw(7) << 0. <<" " // invariant mass
250 << setw(7) << 0. <<" " // number of tracks
251 << setw(7) << 0. <<" "
252 << setw(7) << 0. <<" " // E_had/E_em
253 << setw(7) << 0.0 << setw(7) << 0.0 // dummy/dummy
254 << endl;
255
256 outfile.close();
257 // return outfile;
258}
259
260
261void LHCOConverter::BranchReaderMuon(const TClonesArray* branchMuon, const TClonesArray* branchJet, unsigned int& line) const {
262
263 ofstream outfile( outputfilename_.c_str(),ios_base::app);
264 unsigned int branch_size = branchMuon->GetEntries();
265 TRootMuon * muon = 0;
266 for (unsigned int i=0; i< branch_size; i++) {
267 muon = (TRootMuon*) branchMuon->At(i);
268
269 TLorentzVector pmu; pmu.SetPtEtaPhiE(muon->PT, muon->Eta, muon->Phi, muon->E);
270 float jmass = pmu.M(); if(jmass <0) jmass =0; // small negative values are possible due to smearing. < 0.1 GeV
271
272 // loop on jets, for muon isolation
273 unsigned int closest_jet_ID = 0;
274 float closest_jet_dR = 1E99; // initialisation value, should be high. No further impact
275 for (unsigned int j=0; j< (unsigned int) branchJet->GetEntries(); j++) {
276 TRootJet jet(*((TRootJet*)branchJet->At(j)));
277 float dr_ij = DeltaR(muon->Phi, muon->Eta, jet.Phi, jet.Eta);
278 if( dr_ij < closest_jet_dR ) {
279 closest_jet_ID = j;
280 closest_jet_dR = dr_ij;
281 }
282 } // loop on jets
283
284 int ratio = 0; char ratio_string[3];
285 if(muon->EtRatio != UNDEFINED) ratio = static_cast<int>(muon->EtRatio);
286 if(ratio<0) sprintf(ratio_string,"00");
287 else if(ratio<10) sprintf(ratio_string,"0%d",ratio);
288 else if(ratio<100)sprintf(ratio_string,"%d",ratio);
289 else { cout << "Error: muon EtRatio > 99" << endl; sprintf(ratio_string,"99"); }
290 // output
291 outfile << fixed << setprecision(3)
292 << setw(3) << line++ // line counter
293 << setw(4) << lhcoMuonID << " " // object ID in lhco format
294 << setw(7) << muon->Eta <<" "
295 << setw(7) << muon->Phi <<" "
296 << setw(7) << muon->PT <<" "
297 << setw(7) << jmass <<" " // invariant mass
298 << setw(7) << (float)muon->Charge <<" " // -1 for mu-, +1 for mu+
299 << setw(7) << (float)closest_jet_ID << " ";
300 outfile << fixed << setprecision(0)
301 << setw(4) << floor(muon->SumPt) << "." << ratio_string<<" "; // E_had/E_em
302 outfile << fixed << setprecision(3)
303 << setw(7) << 0.0 << setw(7) << 0.0 // dummy/dummy
304 << endl;
305 } // muon branch loop
306
307 outfile.close();
308}
Note: See TracBrowser for help on using the repository browser.