Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 39

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

ok for ETmis
File size: 15.6 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 "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
32#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
33
34// get info on how fastjet was configured
35#include "Utilities/Fastjet/include/fastjet/config.h"
36
37// make sure we have what is needed
38#ifdef ENABLE_PLUGIN_SISCONE
39# include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
40#endif
41#ifdef ENABLE_PLUGIN_CDFCONES
42# include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
43# include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
44#endif
45
46#include<vector>
47#include<iostream>
48
49#include "interface/TreeClasses.h"
50using namespace std;
51
52//------------------------------------------------------------------------------
53
54// //********************************** PYTHIA INFORMATION*********************************
55
56TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
57 {
58 TIter it((TCollection*)GEN);
59 it.Reset();
60 TRootGenParticle *gen1;
61 TSimpleArray<TRootGenParticle> array,array2;
62
63 while((gen1 = (TRootGenParticle*) it.Next()))
64 {
65 array.Add(gen1);
66 }
67 it.Reset();
68 bool tauhad;
69 while((gen1 = (TRootGenParticle*) it.Next()))
70 {
71 tauhad=false;
72 if(abs(gen1->PID)==15)
73 {
74 int d1=gen1->D1;
75 int d2=gen1->D2;
76
77 if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
78 {
79 tauhad=true;
80 for(int d=d1; d < d2+1; d++)
81 {
82 if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
83 }
84 }
85 }
86 if(tauhad)array2.Add(gen1);
87 }
88 return array2;
89 }
90
91
92
93void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, vector<fastjet::PseudoJet> sorted_jetsS)
94{
95 JETSm.SetPxPyPzE(0,0,0,0);
96 float deltaRtest=5000;
97 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
98 TLorentzVector Att;
99 Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
100 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
101 {
102 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
103 if(deltaRtest < 0.25)
104 {
105 JETSm = Att;
106 }
107 }
108 }
109}
110
111
112int main(int argc, char *argv[])
113{
114 int appargc = 2;
115 char *appName = "Smearing";
116 char *appargv[] = {appName, "-b"};
117 TApplication app(appName, &appargc, appargv);
118
119 if(argc != 4 && argc != 3) {
120 cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
121 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
122 cout << " output_file - output file." << endl;
123 cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
124 exit(1);
125 }
126
127 srand (time (NULL)); /* Initialisation du générateur */
128
129 //read the input TROOT file
130 string inputFileList(argv[1]), outputfilename(argv[2]);
131 if(outputfilename.find(".root") > outputfilename.length() ) {
132 cout << "output_file should be a .root file!\n";
133 return -1;
134 }
135 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
136 outputFile->Close();
137
138 string line;
139 ifstream infile(inputFileList.c_str());
140 infile >> line; // the first line determines the type of input files
141
142 DataConverter *converter=0;
143
144 if(strstr(line.c_str(),".hep"))
145 {
146 cout<<"*************************************************************************"<<endl;
147 cout<<"************ StdHEP file format detected **************"<<endl;
148 cout<<"************ Starting convertion to TRoot format **************"<<endl;
149 cout<<"*************************************************************************"<<endl;
150 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
151 }
152 else if(strstr(line.c_str(),".lhe"))
153 {
154 cout<<"*************************************************************************"<<endl;
155 cout<<"************ LHEF file format detected **************"<<endl;
156 cout<<"************ Starting convertion to TRoot format **************"<<endl;
157 cout<<"*************************************************************************"<<endl;
158 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
159 }
160 else if(strstr(line.c_str(),".root"))
161 {
162 cout<<"*************************************************************************"<<endl;
163 cout<<"************ h2root file format detected **************"<<endl;
164 cout<<"************ Starting convertion to TRoot format **************"<<endl;
165 cout<<"*************************************************************************"<<endl;
166 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
167 }
168 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
169
170 TChain chain("GEN");
171 chain.Add(outputfilename.c_str());
172 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
173 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
174 TIter itGen((TCollection*)branchGen);
175
176 //write the output root file
177 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
178 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
179 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
180 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
181 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
182 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
183
184 TRootGenParticle *particle;
185
186 RESOLELEC * elementElec;
187 RESOLMUON *elementMuon;
188 RESOLJET *elementJet;
189 TAUHAD *elementTaujet;
190 ETMIS *elementEtmis;
191
192
193 //read the datacard input file
194 string DetDatacard("");
195 if(argc==4) DetDatacard =argv[3];
196 RESOLution *DET = new RESOLution();
197 DET->ReadDataCard(DetDatacard);
198
199 TLorentzVector genMomentum(0,0,0,0);
200 LorentzVector jetMomentum;
201 vector<TLorentzVector> TrackCentral;
202
203 vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
204 vector<fastjet::PseudoJet> inclusive_jets;
205 vector<fastjet::PseudoJet> sorted_jets;
206
207 vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
208 vector<fastjet::PseudoJet> inclusive_jetsS;
209 vector<fastjet::PseudoJet> sorted_jetsS;
210
211 vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
212 vector<fastjet::PseudoJet> inclusive_jetsT;
213 vector<fastjet::PseudoJet> sorted_jetsT;
214
215 vector<PhysicsTower> towers;
216
217 fastjet::JetDefinition jet_def;
218 fastjet::JetDefinition jet_defS;
219 fastjet::JetDefinition jet_defT;
220 fastjet::JetDefinition::Plugin * plugins;
221 fastjet::JetDefinition::Plugin * pluginsS;
222 fastjet::JetDefinition::Plugin * pluginsT;
223
224 // set up a CDF midpoint jet definition
225#ifdef ENABLE_PLUGIN_CDFCONES
226 plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
227 jet_def = fastjet::JetDefinition(plugins);
228#else
229 plugins = NULL;
230#endif
231
232 // set up a CDF midpoint jet definition
233#ifdef ENABLE_PLUGIN_CDFCONES
234 pluginsS = new fastjet::CDFJetCluPlugin(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
235 jet_defS = fastjet::JetDefinition(pluginsS);
236#else
237 pluginsS = NULL;
238#endif
239
240 // set up a CDF midpoint jet definition
241 #ifdef ENABLE_PLUGIN_CDFCONES
242 pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
243 jet_defT = fastjet::JetDefinition(pluginsT);
244 #else
245 pluginsT = NULL;
246 #endif
247
248
249 // Loop over all events
250 Long64_t entry, allEntries = treeReader->GetEntries();
251 cout << "** Chain contains " << allEntries << " events" << endl;
252 for(entry = 0; entry < allEntries; ++entry)
253 {
254 TLorentzVector PTmisS(0,0,0,0);
255 TLorentzVector PTmis(0,0,0,0);
256 treeReader->ReadEntry(entry);
257 treeWriter->Clear();
258
259 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
260
261 TSimpleArray<TRootGenParticle> bGen;
262 itGen.Reset();
263 TrackCentral.clear();
264 TSimpleArray<TRootGenParticle> NFCentralQ;
265
266 sorted_jetsS.clear();
267 input_particlesS.clear();
268 inclusive_jetsS.clear();
269
270 sorted_jetsT.clear();
271 input_particlesT.clear();
272 inclusive_jetsT.clear();
273
274 sorted_jets.clear();
275 input_particles.clear();
276 inclusive_jets.clear();
277 towers.clear();
278
279 // Loop over all particles in event
280 while( (particle = (TRootGenParticle*) itGen.Next()) )
281 {
282 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
283
284 int pid = abs(particle->PID);
285 float eta = fabs(particle->Eta);
286
287 if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
288 {
289 input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
290 }
291
292 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
293 // keeps only final particles, visible by the central detector, including the fiducial volume
294 // the ordering of conditions have been optimised for speed : put first the STATUS condition
295 if( (particle->Status == 1) &&
296 (
297 (pid == pMU && eta < DET->MAX_MU) ||
298 (pid != pMU && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
299 )
300 {
301 //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
302
303 Float_t nonS=0;
304 switch(pid) {
305 case pE: // all electrons with eta < DET->MAX_CALO_FWD
306 nonS = genMomentum.E();
307 DET->SmearElectron(genMomentum);
308 if(eta < DET->MAX_TRACKER){
309 elementElec=(RESOLELEC*) branchelec->NewEntry();
310 elementElec->NonSmeareE = nonS;
311 elementElec->SmeareE = genMomentum.E();}
312 break; // case pE
313 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
314 DET->SmearElectron(genMomentum);
315 break; // case pGAMMA
316 case pMU: // all muons with eta < DET->MAX_MU
317 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
318 elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
319 DET->SmearMu(genMomentum);
320 elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
321 break; // case pMU
322 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
323 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
324 DET->SmearHadron(genMomentum, 0.7);
325 break; // case hadron
326 default: // all other final particles with eta < DET->MAX_CALO_FWD
327 DET->SmearHadron(genMomentum, 1.0);
328 break;
329 } // switch (pid)
330
331 if(pid != pMU)
332 {
333 towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
334 input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
335 }
336
337 // all final charged particles
338 if(
339 ((rand()%100) < DET->TRACKING_EFF) &&
340 (genMomentum.E()!=0) &&
341 (fabs(particle->Eta) < DET->MAX_TRACKER) &&
342 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account
343 (pid != pGAMMA) &&
344 (pid != pPI0) &&
345 (pid != pK0L) &&
346 (pid != pN) &&
347 (pid != pSIGMA0) &&
348 (pid != pDELTA0) &&
349 (pid != pK0S) // not charged particles : invisible by tracker
350 )
351 TrackCentral.push_back(genMomentum);
352 } // switch
353 } // while
354
355 TLorentzVector Att(0.,0.,0.,0.);
356 float ScalarEt=0;
357 for(unsigned int i=0; i < towers.size(); i++)
358 {
359 Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
360 ScalarEt = ScalarEt + Att.Et();
361 PTmisS = PTmisS + Att;
362 }
363//cout<<"non smeare "<<PTmis.Pt()<<" "<<PTmisS.Pt()<<endl;
364//cout<<"smeare "<<PTmis.Px()<<" "<<PTmisS.Px()<<endl;
365//cout<<"**********"<<endl;
366
367 elementEtmis= (ETMIS*) branchetmis->NewEntry();
368 elementEtmis->Et = (PTmis).Pt();
369 elementEtmis->Ex = (-PTmis).Px();
370 elementEtmis->SEt = ScalarEt;
371
372 elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
373 elementEtmis->ExSmeare = (-PTmisS).Px()-(PTmis).Px();
374
375 //*****************************
376
377 double ptmin=1;
378 if(input_particles.size()!=0)
379 {
380 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
381 inclusive_jets = clust_seq.inclusive_jets(ptmin);
382 sorted_jets = sorted_by_pt(inclusive_jets);
383 }
384
385 if(input_particlesS.size()!=0)
386 {
387 fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
388 inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
389 sorted_jetsS = sorted_by_pt(inclusive_jetsS);
390 }
391
392 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
393
394 TLorentzVector JETSm(0,0,0,0);
395 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
396 TLorentzVector JET(0,0,0,0);
397 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
398 PairingJet(JETSm,JET,sorted_jetsS);
399 if(JETSm.Pt()>1)
400 {
401 elementJet= (RESOLJET*) branchjet->NewEntry();
402 elementJet->NonSmearePT = JET.Et();
403 elementJet->SmearePT = JETSm.Et()/JET.Et();
404 }
405 }
406 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
407 TLorentzVector JETT(0,0,0,0);
408 JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
409 if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS))
410 {
411 for(Int_t i=0; i<TausHadr.GetEntries();i++)
412 {
413 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
414 {
415 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
416 elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
417 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JETT.Eta(),JETT.Phi());
418 }
419 }
420 }
421
422
423 } // for itJet : loop on all jets
424
425 treeWriter->Fill();
426 } // Loop over all events
427 treeWriter->Write();
428
429 cout << "** Exiting..." << endl;
430
431 delete treeWriter;
432 delete treeReader;
433 delete DET;
434 if(converter) delete converter;
435}
436
Note: See TracBrowser for help on using the repository browser.