Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 16

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

first commit
File size: 10.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 "Utilities/CDFCones/interface/JetCluAlgorithm.h"
22#include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
23#include "Utilities/CDFCones/interface/PhysicsTower.h"
24#include "Utilities/CDFCones/interface/Cluster.h"
25
26#include "interface/DataConverter.h"
27#include "interface/HEPEVTConverter.h"
28#include "interface/LHEFConverter.h"
29#include "interface/STDHEPConverter.h"
30#include "interface/SmearUtil.h"
31#include "interface/TreeClasses.h"
32
33using namespace std;
34
35//------------------------------------------------------------------------------
36
37
38
39void PairingJet(TLorentzVector &JETSm, TLorentzVector JET, vector<Cluster> jetsS)
40{
41 JETSm.SetPxPyPzE(0,0,0,0);
42 vector<Cluster>::iterator itJetS;
43 float deltaRtest=5000;
44 LorentzVector jetMomentumS;
45 for(itJetS = jetsS.begin(); itJetS != jetsS.end(); ++itJetS)
46 {
47 jetMomentumS = itJetS->fourVector;
48 TLorentzVector Att;
49 Att.SetPxPyPzE(jetMomentumS.px,jetMomentumS.py,jetMomentumS.pz,jetMomentumS.E);
50 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
51 {
52 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
53 if(deltaRtest < 0.25)
54 {
55 JETSm = Att;
56 }
57 }
58 }
59}
60
61
62int main(int argc, char *argv[])
63{
64 int appargc = 2;
65 char *appName = "Smearing";
66 char *appargv[] = {appName, "-b"};
67 TApplication app(appName, &appargc, appargv);
68
69 if(argc != 4 && argc != 3) {
70 cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
71 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
72 cout << " output_file - output file." << endl;
73 cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
74 exit(1);
75 }
76
77 srand (time (NULL)); /* Initialisation du générateur */
78
79 //read the input TROOT file
80 string inputFileList(argv[1]), outputfilename(argv[2]);
81 if(outputfilename.find(".root") > outputfilename.length() ) {
82 cout << "output_file should be a .root file!\n";
83 return -1;
84 }
85 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
86 outputFile->Close();
87
88 string line;
89 ifstream infile(inputFileList.c_str());
90 infile >> line; // the first line determines the type of input files
91
92 DataConverter *converter=0;
93
94 if(strstr(line.c_str(),".hep"))
95 {
96 cout<<"*************************************************************************"<<endl;
97 cout<<"************ StdHEP file format detected **************"<<endl;
98 cout<<"************ Starting convertion to TRoot format **************"<<endl;
99 cout<<"*************************************************************************"<<endl;
100 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
101 }
102 else if(strstr(line.c_str(),".lhe"))
103 {
104 cout<<"*************************************************************************"<<endl;
105 cout<<"************ LHEF file format detected **************"<<endl;
106 cout<<"************ Starting convertion to TRoot format **************"<<endl;
107 cout<<"*************************************************************************"<<endl;
108 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
109 }
110 else if(strstr(line.c_str(),".root"))
111 {
112 cout<<"*************************************************************************"<<endl;
113 cout<<"************ h2root file format detected **************"<<endl;
114 cout<<"************ Starting convertion to TRoot format **************"<<endl;
115 cout<<"*************************************************************************"<<endl;
116 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
117 }
118 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
119
120 TChain chain("GEN");
121 chain.Add(outputfilename.c_str());
122 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
123 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
124 TIter itGen((TCollection*)branchGen);
125
126 //write the output root file
127 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
128 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
129 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
130 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
131 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
132 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
133
134 TRootGenParticle *particle;
135 TRootETmis *etmisc;
136
137 RESOLELEC *elementElec;
138 RESOLMUON *elementMuon;
139 RESOLJET *elementJet;
140 TAUHAD *elementTaujet;
141 ETMIS *elementEtmis;
142
143
144 //read the datacard input file
145 string DetDatacard("");
146 if(argc==4) DetDatacard =argv[3];
147 RESOLution *DET = new RESOLution();
148 DET->ReadDataCard(DetDatacard);
149
150
151 TLorentzVector genMomentum(0,0,0,0);
152 LorentzVector jetMomentum;
153 vector<TLorentzVector> TrackCentral;
154 vector<PhysicsTower> towersS;
155 vector<Cluster> jetsS;
156 vector<PhysicsTower> towers;
157 vector<Cluster> jets;
158
159
160 // Loop over all events
161 Long64_t entry, allEntries = treeReader->GetEntries();
162 cout << "** Chain contains " << allEntries << " events" << endl;
163 for(entry = 0; entry < allEntries; ++entry)
164 {
165 TLorentzVector PTmis(0,0,0,0);
166 treeReader->ReadEntry(entry);
167 treeWriter->Clear();
168
169 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
170
171 TSimpleArray<TRootGenParticle> bGen;
172 itGen.Reset();
173 TrackCentral.clear();
174 towers.clear();
175 towersS.clear();
176 TSimpleArray<TRootGenParticle> NFCentralQ;
177
178
179 // Loop over all particles in event
180 while( (particle = (TRootGenParticle*) itGen.Next()) )
181 {
182 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
183
184 int pid = abs(particle->PID);
185 float eta = fabs(particle->Eta);
186
187if(particle->Status == 1)towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
188
189 // keeps only final particles, visible by the central detector, including the fiducial volume
190 // the ordering of conditions have been optimised for speed : put first the STATUS condition
191 if( (particle->Status == 1) &&
192 (
193 (pid == pMU && eta < DET->MAX_MU) ||
194 (pid != pMU && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
195 )
196 ) {
197 switch(pid) {
198
199 case pE: // all electrons with eta < DET->MAX_CALO_FWD
200 DET->SmearElectron(genMomentum);
201 break; // case pE
202
203 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
204 DET->SmearElectron(genMomentum);
205 break; // case pGAMMA
206
207 case pMU: // all muons with eta < DET->MAX_MU
208 DET->SmearMu(genMomentum);
209 break; // case pMU
210
211 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
212 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
213 DET->SmearHadron(genMomentum, 0.7);
214 break; // case hadron
215
216 default: // all other final particles with eta < DET->MAX_CALO_FWD
217 DET->SmearHadron(genMomentum, 1.0);
218 break;
219 } // switch (pid)
220
221 // all final particles but muons and neutrinos
222 // for calorimetric towers and mission PT
223 if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
224
225 if(pid != pMU && genMomentum.Pt() > DET->PT_TRACKS_MIN)
226 towersS.push_back
227 (
228 PhysicsTower
229 (
230 LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())
231 )
232 );
233 // all final charged particles
234 if(
235 ((rand()%100) < DET->TRACKING_EFF) &&
236 (genMomentum.E()!=0) &&
237 (fabs(particle->Eta) < DET->MAX_TRACKER) &&
238 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account
239 (pid != pGAMMA) &&
240 (pid != pPI0) &&
241 (pid != pK0L) &&
242 (pid != pN) &&
243 (pid != pSIGMA0) &&
244 (pid != pDELTA0) &&
245 (pid != pK0S) // not charged particles : invisible by tracker
246 )
247 TrackCentral.push_back(genMomentum);
248 } // switch
249 } // while
250
251
252 //*****************************
253 JetCluAlgorithm jetAlgoS(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
254 jetAlgoS.run(towersS, jetsS);
255
256 JetCluAlgorithm jetAlgo(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
257 jetAlgo.run(towers, jets);
258
259
260 vector<Cluster>::iterator itJet;
261 TLorentzVector JETSm(0,0,0,0);
262 for(itJet = jets.begin(); itJet != jets.end(); ++itJet) {
263 jetMomentum = itJet->fourVector;
264 TLorentzVector JET(0,0,0,0);
265 JET.SetPxPyPzE(jetMomentum.px,jetMomentum.py,jetMomentum.pz,jetMomentum.E);
266 PairingJet(JETSm,JET,jetsS);
267 if(JETSm.Pt()>3)
268 {
269 elementJet= (RESOLJET*) branchjet->NewEntry();
270 elementJet->NonSmearePT = JET.Pt();
271 elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
272 /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
273 cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
274 cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
275 cout<<"************"<<endl;
276*/
277 }
278
279 } // for itJet : loop on all jets
280
281 treeWriter->Fill();
282 } // Loop over all events
283 treeWriter->Write();
284
285 cout << "** Exiting..." << endl;
286
287 delete treeWriter;
288 delete treeReader;
289 delete DET;
290 if(converter) delete converter;
291}
292
Note: See TracBrowser for help on using the repository browser.