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"
|
---|
50 | using namespace std;
|
---|
51 |
|
---|
52 | //------------------------------------------------------------------------------
|
---|
53 |
|
---|
54 | void 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 |
|
---|
74 | int 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 |
|
---|
310 | cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl;
|
---|
311 | cout<<"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 |
|
---|