Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 151

Last change on this file since 151 was 100, checked in by severine ovyn, 16 years ago

Remove datacard bug + CaloTowers OK
File size: 13.7 KB
Line 
1/*
2 ---- FastSim ----
3 A Fast Simulator for general purpose LHC detector
4 S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5
6 Center for Particle Physics and Phenomenology (CP3)
7 Universite Catholique de Louvain (UCL)
8 Louvain-la-Neuve, Belgium
9*/
10
11/// \file Smearing.cpp
12/// \brief executable for the FastSim
13
14#include "TChain.h"
15#include "TApplication.h"
16
17#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
18#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
19#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
20
21#include "H_BeamParticle.h"
22#include "H_BeamLine.h"
23#include "H_RomanPot.h"
24
25#include "interface/DataConverter.h"
26#include "interface/HEPEVTConverter.h"
27#include "interface/LHEFConverter.h"
28#include "interface/STDHEPConverter.h"
29
30#include "interface/SmearUtil.h"
31#include "interface/JetUtils.h"
32#include "interface/BFieldProp.h"
33
34#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
35#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
36
37#include<vector>
38#include<iostream>
39
40#include "interface/TreeClasses.h"
41using namespace std;
42
43//------------------------------------------------------------------------------
44
45// //********************************** PYTHIA INFORMATION*********************************
46
47TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
48 {
49 TIter it((TCollection*)GEN);
50 it.Reset();
51 TRootGenParticle *gen1;
52 TSimpleArray<TRootGenParticle> array,array2;
53
54 while((gen1 = (TRootGenParticle*) it.Next()))
55 {
56 array.Add(gen1);
57 }
58 it.Reset();
59 bool tauhad;
60 while((gen1 = (TRootGenParticle*) it.Next()))
61 {
62 tauhad=false;
63 if(abs(gen1->PID)==15)
64 {
65 int d1=gen1->D1;
66 int d2=gen1->D2;
67
68 if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
69 {
70 tauhad=true;
71 for(int d=d1; d < d2+1; d++)
72 {
73 if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
74 }
75 }
76 }
77 if(tauhad)array2.Add(gen1);
78 }
79 return array2;
80 }
81
82
83
84void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, vector<fastjet::PseudoJet> sorted_jetsS)
85{
86 JETSm.SetPxPyPzE(0,0,0,0);
87 float deltaRtest=5000;
88 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
89 TLorentzVector Att;
90 Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
91 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
92 {
93 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
94 if(deltaRtest < 0.25)
95 {
96 JETSm = Att;
97 }
98 }
99 }
100}
101
102int main(int argc, char *argv[])
103{
104 int appargc = 2;
105 char *appName = "Smearing";
106 char *appargv[] = {appName, "-b"};
107 TApplication app(appName, &appargc, appargv);
108
109 if(argc != 4 && argc != 3) {
110 cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
111 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
112 cout << " output_file - output file." << endl;
113 cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
114 exit(1);
115 }
116
117 srand (time (NULL)); /* Initialisation du générateur */
118
119 //read the input TROOT file
120 string inputFileList(argv[1]), outputfilename(argv[2]);
121 if(outputfilename.find(".root") > outputfilename.length() ) {
122 cout << "output_file should be a .root file!\n";
123 return -1;
124 }
125
126 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
127 outputFile->Close();
128
129 string line;
130 ifstream infile(inputFileList.c_str());
131 infile >> line; // the first line determines the type of input files
132
133 DataConverter *converter=0;
134
135 if(strstr(line.c_str(),".hep"))
136 {
137 cout<<"*************************************************************************"<<endl;
138 cout<<"************ StdHEP file format detected **************"<<endl;
139 cout<<"************ Starting convertion to TRoot format **************"<<endl;
140 cout<<"*************************************************************************"<<endl;
141 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
142 }
143 else if(strstr(line.c_str(),".lhe"))
144 {
145 cout<<"*************************************************************************"<<endl;
146 cout<<"************ LHEF file format detected **************"<<endl;
147 cout<<"************ Starting convertion to TRoot format **************"<<endl;
148 cout<<"*************************************************************************"<<endl;
149 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
150 }
151 else if(strstr(line.c_str(),".root"))
152 {
153 cout<<"*************************************************************************"<<endl;
154 cout<<"************ h2root file format detected **************"<<endl;
155 cout<<"************ Starting convertion to TRoot format **************"<<endl;
156 cout<<"*************************************************************************"<<endl;
157 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
158 }
159 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
160
161 TChain chain("GEN");
162 chain.Add(outputfilename.c_str());
163 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
164 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
165 TIter itGen((TCollection*)branchGen);
166
167 //write the output root file
168 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
169 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
170 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
171 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
172 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
173 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
174
175 TRootGenParticle *particle;
176
177 RESOLELEC * elementElec;
178 RESOLMUON *elementMuon;
179 RESOLJET *elementJet;
180 TAUHAD *elementTaujet;
181 ETMIS *elementEtmis;
182
183
184 //read the datacard input file
185 string DetDatacard("data/DataCardDet.dat");
186 if(argc==4) DetDatacard =argv[3];
187 RESOLution *DET = new RESOLution();
188 DET->ReadDataCard(DetDatacard);
189
190 //Jet information
191 JetsUtil *JETRUN = new JetsUtil(DetDatacard);
192
193 //Propagation of tracks in the B field
194 TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
195
196 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
197 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
198 LorentzVector jetMomentum;
199 vector<TLorentzVector> TrackCentral;
200
201 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
202 vector<fastjet::PseudoJet> sorted_jetsGEN;
203
204 vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
205 vector<fastjet::PseudoJet> sorted_jetsReco;
206
207 vector<PhysicsTower> towers;
208
209 // Loop over all events
210 Long64_t entry, allEntries = treeReader->GetEntries();
211 cout << "** Chain contains " << allEntries << " events" << endl;
212 for(entry = 0; entry < allEntries; ++entry)
213 {
214 TLorentzVector PTmisReco(0,0,0,0);
215 TLorentzVector PTmisGEN(0,0,0,0);
216 treeReader->ReadEntry(entry);
217 treeWriter->Clear();
218 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
219
220 TSimpleArray<TRootGenParticle> bGen;
221 itGen.Reset();
222 TrackCentral.clear();
223 TSimpleArray<TRootGenParticle> NFCentralQ;
224
225 input_particlesReco.clear();
226 input_particlesGEN.clear();
227 towers.clear();
228
229 // Loop over all particles in event
230 while( (particle = (TRootGenParticle*) itGen.Next()) )
231 {
232 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
233
234 int pid = abs(particle->PID);
235 float eta = fabs(particle->Eta);
236
237 //input generator level particle for jet algorithm
238 if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
239 {
240 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
241 }
242
243 //Calculate ETMIS from generated particles
244 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
245
246 if( (particle->Status == 1) &&
247 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
248 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
249 )
250 {
251 recoMomentum = genMomentum;
252 //use of the magnetic field propagation
253 TRACP->Propagation(particle,recoMomentum);
254 // cout<<"eta "<<eta<<endl;
255 eta=fabs(recoMomentum.Eta());
256 // cout<<"eta apres"<<eta<<endl;
257
258 switch(pid) {
259
260 case pE: // all electrons with eta < DET->MAX_CALO_FWD
261 DET->SmearElectron(recoMomentum);
262 if(recoMomentum.E() !=0 && eta < DET->CEN_max_tracker){
263 elementElec=(RESOLELEC*) branchelec->NewEntry();
264 elementElec->E = genMomentum.E();
265 elementElec->SmearedE = recoMomentum.E();}
266 break; // case pE
267 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
268 DET->SmearElectron(recoMomentum);
269 break; // case pGAMMA
270 case pMU: // all muons with eta < DET->MAX_MU
271 DET->SmearMu(recoMomentum);
272 if(recoMomentum.E() !=0){
273 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
274 elementMuon->OverPT = (1/genMomentum.Pt());
275 elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
276 break; // case pMU
277 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
278 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
279 DET->SmearHadron(recoMomentum, 0.7);
280 break; // case hadron
281 default: // all other final particles with eta < DET->MAX_CALO_FWD
282 DET->SmearHadron(recoMomentum, 1.0);
283 break;
284 } // switch (pid)
285
286 //information to reconstruct jets from reconstructed towers
287 int charge=Charge(pid);
288 if(recoMomentum.E() !=0 && pid != pMU) {
289 if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->TRACK_ptmin)){
290 towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())));
291 input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()));
292 }
293 }
294
295 // all final charged particles
296 if(
297 (genMomentum.E()!=0) &&
298 (fabs(genMomentum.Eta()) < DET->CEN_max_tracker) &&
299 (genMomentum.Pt() > DET->TRACK_ptmin ) && // pt too small to be taken into account
300 ((rand()%100) < DET->TRACK_eff) &&
301 (charge!=0)
302 )
303 {TrackCentral.push_back(genMomentum);}
304
305 } // switch
306 } // while
307
308 //compute missing transverse energy from calo towers
309 TLorentzVector Att(0.,0.,0.,0.);
310 float ScalarEt=0;
311 for(unsigned int i=0; i < towers.size(); i++)
312 {
313 Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
314 ScalarEt = ScalarEt + Att.Et();
315 PTmisReco = PTmisReco + Att;
316 }
317
318 elementEtmis= (ETMIS*) branchetmis->NewEntry();
319 elementEtmis->Et = (PTmisGEN).Pt();
320 elementEtmis->Ex = (-PTmisGEN).Px();
321 elementEtmis->SEt = ScalarEt;
322
323 elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
324 elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
325
326 //*****************************
327
328 sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
329 sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
330
331 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
332
333 TLorentzVector JETreco(0,0,0,0);
334 for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
335 TLorentzVector JETgen(0,0,0,0);
336 JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
337 PairingJet(JETreco,JETgen,sorted_jetsReco);
338 if(JETreco.Pt()>1)
339 {
340 elementJet= (RESOLJET*) branchjet->NewEntry();
341 elementJet->PT = JETgen.Et();
342 elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
343 }
344 }
345 for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
346 TLorentzVector JETT(0,0,0,0);
347 JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
348 if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
349 {
350 for(Int_t i=0; i<TausHadr.GetEntries();i++)
351 {
352 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
353 {
354 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
355 elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
356 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi());
357 }
358 }
359 }
360
361
362 } // for itJet : loop on all jets
363
364 treeWriter->Fill();
365 } // Loop over all events
366 treeWriter->Write();
367
368 cout << "** Exiting..." << endl;
369
370 delete treeWriter;
371 delete treeReader;
372 delete DET;
373 if(converter) delete converter;
374}
375
Note: See TracBrowser for help on using the repository browser.