Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 57

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

general jet variable
File size: 15.3 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->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->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->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 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
259
260 TSimpleArray<TRootGenParticle> bGen;
261 itGen.Reset();
262 TrackCentral.clear();
263 TSimpleArray<TRootGenParticle> NFCentralQ;
264
265 sorted_jetsS.clear();
266 input_particlesS.clear();
267 inclusive_jetsS.clear();
268
269 sorted_jetsT.clear();
270 input_particlesT.clear();
271 inclusive_jetsT.clear();
272
273 sorted_jets.clear();
274 input_particles.clear();
275 inclusive_jets.clear();
276 towers.clear();
277
278 // Loop over all particles in event
279 while( (particle = (TRootGenParticle*) itGen.Next()) )
280 {
281 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
282
283 int pid = abs(particle->PID);
284 float eta = fabs(particle->Eta);
285
286 if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
287 {
288 input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
289 }
290
291 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
292 // keeps only final particles, visible by the central detector, including the fiducial volume
293 // the ordering of conditions have been optimised for speed : put first the STATUS condition
294 if( (particle->Status == 1) &&
295 (
296 (pid == pMU && eta < DET->MAX_MU) ||
297 (pid != pMU && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
298 )
299 {
300 //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
301
302 Float_t nonS=0;
303 switch(pid) {
304 case pE: // all electrons with eta < DET->MAX_CALO_FWD
305 nonS = genMomentum.E();
306 DET->SmearElectron(genMomentum);
307 if(eta < DET->MAX_TRACKER){
308 elementElec=(RESOLELEC*) branchelec->NewEntry();
309 elementElec->NonSmeareE = nonS;
310 elementElec->SmeareE = genMomentum.E();}
311 break; // case pE
312 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
313 DET->SmearElectron(genMomentum);
314 break; // case pGAMMA
315 case pMU: // all muons with eta < DET->MAX_MU
316 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
317 elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
318 DET->SmearMu(genMomentum);
319 elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
320 break; // case pMU
321 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
322 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
323 DET->SmearHadron(genMomentum, 0.7);
324 break; // case hadron
325 default: // all other final particles with eta < DET->MAX_CALO_FWD
326 DET->SmearHadron(genMomentum, 1.0);
327 break;
328 } // switch (pid)
329
330 if(pid != pMU && genMomentum.Et() !=0)
331 {
332 towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
333 input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
334 }
335
336 // all final charged particles
337 if(
338 ((rand()%100) < DET->TRACKING_EFF) &&
339 (genMomentum.E()!=0) &&
340 (fabs(particle->Eta) < DET->MAX_TRACKER) &&
341 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account
342 (pid != pGAMMA) &&
343 (pid != pPI0) &&
344 (pid != pK0L) &&
345 (pid != pN) &&
346 (pid != pSIGMA0) &&
347 (pid != pDELTA0) &&
348 (pid != pK0S) // not charged particles : invisible by tracker
349 )
350 TrackCentral.push_back(genMomentum);
351 } // switch
352 } // while
353
354 TLorentzVector Att(0.,0.,0.,0.);
355 for(unsigned int i=0; i < towers.size(); i++)
356 {
357 Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
358 PTmisS = PTmisS + Att;
359 }
360
361 elementEtmis= (ETMIS*) branchetmis->NewEntry();
362 elementEtmis->Et = (PTmis).Pt();
363 elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
364
365 //*****************************
366
367 double ptmin=1;
368 if(input_particles.size()!=0)
369 {
370 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
371 inclusive_jets = clust_seq.inclusive_jets(ptmin);
372 sorted_jets = sorted_by_pt(inclusive_jets);
373 }
374
375 if(input_particlesS.size()!=0)
376 {
377 fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
378 inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
379 sorted_jetsS = sorted_by_pt(inclusive_jetsS);
380 }
381
382 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
383
384 TLorentzVector JETSm(0,0,0,0);
385 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
386 TLorentzVector JET(0,0,0,0);
387 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
388 PairingJet(JETSm,JET,sorted_jetsS);
389 if(JETSm.Pt()>1)
390 {
391 elementJet= (RESOLJET*) branchjet->NewEntry();
392 elementJet->NonSmearePT = JET.Et();
393 elementJet->SmearePT = JETSm.Et()/JET.Et();
394 }
395 }
396 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
397 TLorentzVector JETT(0,0,0,0);
398 JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
399 if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS))
400 {
401 for(Int_t i=0; i<TausHadr.GetEntries();i++)
402 {
403 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
404 {
405 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
406 elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
407 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JETT.Eta(),JETT.Phi());
408 }
409 }
410 }
411
412
413 } // for itJet : loop on all jets
414
415 treeWriter->Fill();
416 } // Loop over all events
417 treeWriter->Write();
418
419 cout << "** Exiting..." << endl;
420
421 delete treeWriter;
422 delete treeReader;
423 delete DET;
424 if(converter) delete converter;
425}
426
Note: See TracBrowser for help on using the repository browser.