Fork me on GitHub

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

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

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