Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 24

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

minor change: appName
File size: 13.1 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
54void PairingJet(TLorentzVector &JETSm, const TLorentzVector& JET, vector<fastjet::PseudoJet> sorted_jetsS)
55{
56 JETSm.SetPxPyPzE(0,0,0,0);
57 float deltaRtest=5000;
58 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
59 TLorentzVector Att;
60 Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
61 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
62 {
63 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
64 if(deltaRtest < 0.25)
65 {
66 JETSm = Att;
67 }
68 }
69 }
70}
71
72//------------------------------------------------------------------------------
73
74int main(int argc, char *argv[])
75{
76 int appargc = 2;
77 char appName[100]; sprintf(appName,argv[0]);
78 char *appargv[] = {appName, "-b"};
79 TApplication app(appName, &appargc, appargv);
80
81 if(argc != 4 && argc != 3) {
82 cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
83 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
84 cout << " output_file - output file." << endl;
85 cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
86 exit(1);
87 }
88
89 srand (time (NULL)); /* Initialisation du générateur */
90
91 //read the input TROOT file
92 string inputFileList(argv[1]), outputfilename(argv[2]);
93 if(outputfilename.find(".root") > outputfilename.length() ) {
94 cout << "output_file should be a .root file!\n";
95 return -1;
96 }
97 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
98 outputFile->Close();
99
100 string line;
101 ifstream infile(inputFileList.c_str());
102 infile >> line; // the first line determines the type of input files
103
104 DataConverter *converter=0;
105
106 if(strstr(line.c_str(),".hep"))
107 {
108 cout<<"*************************************************************************"<<endl;
109 cout<<"************ StdHEP file format detected **************"<<endl;
110 cout<<"************ Starting convertion to TRoot format **************"<<endl;
111 cout<<"*************************************************************************"<<endl;
112 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
113 }
114 else if(strstr(line.c_str(),".lhe"))
115 {
116 cout<<"*************************************************************************"<<endl;
117 cout<<"************ LHEF file format detected **************"<<endl;
118 cout<<"************ Starting convertion to TRoot format **************"<<endl;
119 cout<<"*************************************************************************"<<endl;
120 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
121 }
122 else if(strstr(line.c_str(),".root"))
123 {
124 cout<<"*************************************************************************"<<endl;
125 cout<<"************ h2root file format detected **************"<<endl;
126 cout<<"************ Starting convertion to TRoot format **************"<<endl;
127 cout<<"*************************************************************************"<<endl;
128 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
129 }
130 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
131
132 TChain chain("GEN");
133 chain.Add(outputfilename.c_str());
134 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
135 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
136 TIter itGen((TCollection*)branchGen);
137
138 //write the output root file
139 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
140 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
141 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
142 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
143 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
144 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
145
146 TRootGenParticle *particle;
147 TRootETmis *etmisc;
148
149 RESOLELEC *elementElec;
150 RESOLMUON *elementMuon;
151 RESOLJET *elementJet;
152 TAUHAD *elementTaujet;
153 ETMIS *elementEtmis;
154
155
156 //read the datacard input file
157 string DetDatacard("");
158 if(argc==4) DetDatacard =argv[3];
159 RESOLution *DET = new RESOLution();
160 DET->ReadDataCard(DetDatacard);
161
162 TLorentzVector genMomentum(0,0,0,0);
163 LorentzVector jetMomentum;
164 vector<TLorentzVector> TrackCentral;
165
166 vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
167 vector<fastjet::PseudoJet> inclusive_jets;
168 vector<fastjet::PseudoJet> sorted_jets;
169
170 vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
171 vector<fastjet::PseudoJet> inclusive_jetsS;
172 vector<fastjet::PseudoJet> sorted_jetsS;
173 vector<PhysicsTower> towers;
174
175 fastjet::JetDefinition jet_def;
176 fastjet::JetDefinition jet_defS;
177 fastjet::JetDefinition::Plugin * plugins;
178 fastjet::JetDefinition::Plugin * pluginsS;
179
180 // set up a CDF midpoint jet definition
181 #ifdef ENABLE_PLUGIN_CDFCONES
182 plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
183 jet_def = fastjet::JetDefinition(plugins);
184 #else
185 plugins = NULL;
186 #endif
187
188 // set up a CDF midpoint jet definition
189 #ifdef ENABLE_PLUGIN_CDFCONES
190 pluginsS = new fastjet::CDFJetCluPlugin(2,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
191 jet_defS = fastjet::JetDefinition(pluginsS);
192 #else
193 pluginsS = NULL;
194 #endif
195
196
197 // Loop over all events
198 Long64_t entry, allEntries = treeReader->GetEntries();
199 cout << "** Chain contains " << allEntries << " events" << endl;
200 for(entry = 0; entry < allEntries; ++entry)
201 {
202 TLorentzVector PTmisS(0,0,0,0);
203 TLorentzVector PTmis(0,0,0,0);
204 treeReader->ReadEntry(entry);
205 treeWriter->Clear();
206
207 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
208
209 TSimpleArray<TRootGenParticle> bGen;
210 itGen.Reset();
211 TrackCentral.clear();
212 TSimpleArray<TRootGenParticle> NFCentralQ;
213 input_particles.clear();
214 inclusive_jets.clear();
215 sorted_jets.clear();
216 input_particlesS.clear();
217 inclusive_jetsS.clear();
218 sorted_jetsS.clear();
219 towers.clear();
220
221 // Loop over all particles in event
222 while( (particle = (TRootGenParticle*) itGen.Next()) )
223 {
224 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
225
226 int pid = abs(particle->PID);
227 float eta = fabs(particle->Eta);
228
229 if(particle->Status == 1)
230 {
231 input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
232 }
233
234 // keeps only final particles, visible by the central detector, including the fiducial volume
235 // the ordering of conditions have been optimised for speed : put first the STATUS condition
236 if( (particle->Status == 1) &&
237 (
238 (pid == pMU && eta < DET->MAX_MU) ||
239 (pid != pMU && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
240 )
241 ) {
242 if(pid != pMU)PTmis = PTmis + genMomentum;//ptmis
243 switch(pid) {
244
245 case pE: // all electrons with eta < DET->MAX_CALO_FWD
246 DET->SmearElectron(genMomentum);
247 break; // case pE
248
249 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
250 DET->SmearElectron(genMomentum);
251 break; // case pGAMMA
252
253 case pMU: // all muons with eta < DET->MAX_MU
254 DET->SmearMu(genMomentum);
255 break; // case pMU
256
257 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
258 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
259 DET->SmearHadron(genMomentum, 0.7);
260 break; // case hadron
261
262 default: // all other final particles with eta < DET->MAX_CALO_FWD
263 DET->SmearHadron(genMomentum, 1.0);
264 break;
265 } // switch (pid)
266
267 // all final particles but muons and neutrinos
268 // for calorimetric towers and mission PT
269 //if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
270
271 if(pid != pMU)
272 {
273 PTmisS = PTmisS + genMomentum;
274 towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
275 input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
276 }
277
278 // all final charged particles
279 if(
280 ((rand()%100) < DET->TRACKING_EFF) &&
281 (genMomentum.E()!=0) &&
282 (fabs(particle->Eta) < DET->MAX_TRACKER) &&
283 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account
284 (pid != pGAMMA) &&
285 (pid != pPI0) &&
286 (pid != pK0L) &&
287 (pid != pN) &&
288 (pid != pSIGMA0) &&
289 (pid != pDELTA0) &&
290 (pid != pK0S) // not charged particles : invisible by tracker
291 )
292 TrackCentral.push_back(genMomentum);
293 } // switch
294 } // while
295
296 TLorentzVector toWerS(0,0,0,0);
297 //calcul de ETmis au depart des tours calorimetriques..
298/* for(unsigned int i=0; i < towers.size(); i++) {
299 if(towers[i].fourVector.pt() < 0.5) continue;
300 toWerS.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
301 PTmisS = PTmisS + toWerS;
302 }
303*/
304//PTmis = (0,0,0,0)-PTmis;
305//PTmisS = (0,0,0,0)-PTmisS;
306
307//float ET=PTmis.Et()<<endl;
308//float ETS=PTmisS.Et()<<endl;
309
310cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl;
311cout<<"Ptmis "<<(-PTmis).Pt()<<" PTmis smeare "<<(-PTmisS).Pt()<<endl;
312
313 elementEtmis= (ETMIS*) branchetmis->NewEntry();
314 elementEtmis->Et = (PTmis).Et();
315 elementEtmis->EtSmeare = (PTmisS).Et();
316
317
318
319 //*****************************
320
321 double ptmin=1;
322 if(input_particles.size()!=0)
323 {
324 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
325 inclusive_jets = clust_seq.inclusive_jets(ptmin);
326 sorted_jets = sorted_by_pt(inclusive_jets);
327 }
328
329 if(input_particlesS.size()!=0)
330 {
331 fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
332 inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
333 sorted_jetsS = sorted_by_pt(inclusive_jetsS);
334 }
335
336 TLorentzVector JETSm(0,0,0,0);
337 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
338 TLorentzVector JET(0,0,0,0);
339 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
340 PairingJet(JETSm,JET,sorted_jetsS);
341 if(JETSm.Pt()>1)
342 {
343 float NonSmeareEt=(JET.E()*sin(JET.Theta()));
344 float SmeareEt=(JETSm.E()*sin(JETSm.Theta()));
345
346 elementJet= (RESOLJET*) branchjet->NewEntry();
347 elementJet->NonSmearePT = JET.Et();
348 elementJet->SmearePT = JETSm.Et()/JET.Et();
349
350 }
351
352 } // for itJet : loop on all jets
353
354 treeWriter->Fill();
355 } // Loop over all events
356 treeWriter->Write();
357
358 cout << "** Exiting..." << endl;
359
360 delete treeWriter;
361 delete treeReader;
362 delete DET;
363 if(converter) delete converter;
364}
365
Note: See TracBrowser for help on using the repository browser.