Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 261

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

add header
File size: 15.0 KB
RevLine 
[260]1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \----------------------------------------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
[9]31
32#include "TChain.h"
33#include "TApplication.h"
[227]34#include "TFile.h"
[9]35
[227]36#include "ExRootTreeReader.h"
37#include "ExRootTreeWriter.h"
38#include "ExRootTreeBranch.h"
39#include "TreeClasses.h"
[9]40
[227]41#include "DataConverter.h"
42#include "HEPEVTConverter.h"
43#include "LHEFConverter.h"
44#include "STDHEPConverter.h"
[9]45
[227]46#include "SmearUtil.h"
47#include "JetsUtil.h"
48#include "BFieldProp.h"
[19]49
50#include<vector>
51#include<iostream>
52
[9]53using namespace std;
54
55//------------------------------------------------------------------------------
56
[39]57// //********************************** PYTHIA INFORMATION*********************************
[26]58
[39]59TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
60 {
61 TIter it((TCollection*)GEN);
62 it.Reset();
63 TRootGenParticle *gen1;
64 TSimpleArray<TRootGenParticle> array,array2;
[26]65
[39]66 while((gen1 = (TRootGenParticle*) it.Next()))
67 {
68 array.Add(gen1);
69 }
70 it.Reset();
71 bool tauhad;
72 while((gen1 = (TRootGenParticle*) it.Next()))
73 {
74 tauhad=false;
75 if(abs(gen1->PID)==15)
76 {
[258]77cout<<"au moins on a un tau "<<endl;
[39]78 int d1=gen1->D1;
79 int d2=gen1->D2;
[258]80cout<<"il a des filles? "<<endl;
[39]81 if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
82 {
83 tauhad=true;
84 for(int d=d1; d < d2+1; d++)
[258]85 {cout<<abs(array[d]->PID)<<" "<<endl;
[39]86 if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
87 }
88 }
89 }
90 if(tauhad)array2.Add(gen1);
91 }
92 return array2;
93 }
94
[258]95double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
96 double Energie=0;
97 for(unsigned int i=0; i < towers.size(); i++) {
98 if(towers[i].Pt() < JET_seed) continue;
99 if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
100 Energie += towers[i].E();
101 }
102 }
103 return Energie;
104}
[39]105
106
[258]107void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
[9]108{
109 JETSm.SetPxPyPzE(0,0,0,0);
110 float deltaRtest=5000;
[258]111 TIter itJet((TCollection*)branchJet);
112 TRootJet *jet;
113 itJet.Reset();
114 while( (jet = (TRootJet*) itJet.Next()) )
115 {
[19]116 TLorentzVector Att;
[258]117 Att.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
[19]118 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
119 {
120 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
121 if(deltaRtest < 0.25)
122 {
123 JETSm = Att;
124 }
125 }
126 }
[9]127}
128
[258]129void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
130{
131 ELECSm.SetPxPyPzE(0,0,0,0);
132 float deltaRtest=5000;
133 TIter itElec((TCollection*)branchElec);
134 TRootElectron *elec;
135 itElec.Reset();
136 while( (elec = (TRootElectron*) itElec.Next()) )
137 {
138 TLorentzVector Att;
139 Att.SetPxPyPzE(elec->Px,elec->Py,elec->Pz,elec->E);
140 if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
141 {
142 deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
143 if(deltaRtest < 0.025)
144 {
145 ELECSm = Att;
146 }
147 }
148 }
149}
150
151void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
152{
153 MUONSm.SetPxPyPzE(0,0,0,0);
154 float deltaRtest=5000;
155 TIter itMuon((TCollection*)branchMuon);
156 TRootMuon *muon;
157 itMuon.Reset();
158 while( (muon = (TRootMuon*) itMuon.Next()) )
159 {
160 TLorentzVector Att;
161 Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
162 if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
163 {
164 deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
165 if(deltaRtest < 0.025)
166 {
167 MUONSm = Att;
168 }
169 }
170 }
171}
172
173unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
174 unsigned int numtrack=0;
175 TIter itTrack((TCollection*)branchTracks);
176 TRootTracks *track;
177 itTrack.Reset();
178 while( (track = (TRootTracks*) itTrack.Next()) )
179 {
180 if((track->PT < pt_track )||
181 (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
182 )continue;
183 numtrack++;
184 }
185 return numtrack;
186}
187
188
189
[9]190int main(int argc, char *argv[])
191{
192 int appargc = 2;
[258]193 char *appName = "Resolution";
[9]194 char *appargv[] = {appName, "-b"};
195 TApplication app(appName, &appargc, appargv);
196
[258]197 if(argc != 3) {
198 cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
199 cout << " input_list - list of files in root format," << endl;
[9]200 cout << " output_file - output file." << endl;
201 exit(1);
202 }
203
204 srand (time (NULL)); /* Initialisation du générateur */
205
206 //read the input TROOT file
[258]207 string inputfilename(argv[1]), outputfilename(argv[2]);
208
[9]209 if(outputfilename.find(".root") > outputfilename.length() ) {
210 cout << "output_file should be a .root file!\n";
211 return -1;
212 }
[88]213
[258]214
215
[9]216 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
217 outputFile->Close();
[26]218
[258]219 TChain chainGEN("GEN");
220 chainGEN.Add(inputfilename.c_str());
221 ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
222 TChain chain("Analysis");
223 chain.Add(inputfilename.c_str());
[9]224 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
[258]225 const TClonesArray *branchJet = treeReader->UseBranch("Jet");
226 const TClonesArray *branchElec = treeReader->UseBranch("Electron");
227 const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
228 const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
229 const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
230 const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
[9]231 TIter itGen((TCollection*)branchGen);
232
233 //write the output root file
234 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
235 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
236 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
237 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
238 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
239 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
[26]240
[9]241 TRootGenParticle *particle;
[26]242
[39]243 RESOLELEC * elementElec;
[9]244 RESOLMUON *elementMuon;
245 RESOLJET *elementJet;
246 TAUHAD *elementTaujet;
247 ETMIS *elementEtmis;
[26]248
[258]249 int numTau=0;
250 int numTauRec=0;
[26]251
[258]252 RESOLution *DET = new RESOLution();
253 DET->ReadDataCard("data/Datacard_CMS.dat");
[88]254
[258]255
[88]256 //Jet information
[258]257 JetsUtil *JETRUN = new JetsUtil("data/Datacard_CMS.dat");
[88]258
259 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
[258]260 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
[9]261 LorentzVector jetMomentum;
[26]262
[88]263 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
264 vector<fastjet::PseudoJet> sorted_jetsGEN;
[26]265
[258]266 vector<TLorentzVector> towers;
[152]267
[9]268 // Loop over all events
269 Long64_t entry, allEntries = treeReader->GetEntries();
270 cout << "** Chain contains " << allEntries << " events" << endl;
271 for(entry = 0; entry < allEntries; ++entry)
272 {
[88]273 TLorentzVector PTmisReco(0,0,0,0);
274 TLorentzVector PTmisGEN(0,0,0,0);
[9]275 treeReader->ReadEntry(entry);
[258]276 treeReaderGEN->ReadEntry(entry);
[9]277 treeWriter->Clear();
278 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
279
280 TSimpleArray<TRootGenParticle> bGen;
281 itGen.Reset();
282 TSimpleArray<TRootGenParticle> NFCentralQ;
[89]283
[88]284 input_particlesGEN.clear();
[23]285 towers.clear();
[26]286
[9]287 // Loop over all particles in event
288 while( (particle = (TRootGenParticle*) itGen.Next()) )
289 {
290 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
[26]291
[9]292 int pid = abs(particle->PID);
293 float eta = fabs(particle->Eta);
[89]294
[88]295 //input generator level particle for jet algorithm
[94]296 if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
[26]297 {
[88]298 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
[26]299 }
[89]300
[88]301 //Calculate ETMIS from generated particles
302 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
[258]303
304 if( (particle->Status == 1) &&
305 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
306 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
307 )
[88]308 {
[258]309 eta=fabs(genMomentum.Eta());
[152]310
[258]311 switch(pid) {
312
313 case pE: // all electrons with eta < DET->MAX_CALO_FWD
314 PairingElec(recoMomentum,genMomentum,branchElec);
315 if(recoMomentum.Pt()>1){
316 elementElec=(RESOLELEC*) branchelec->NewEntry();
317 elementElec->E = genMomentum.E();
318 elementElec->SmearedE = recoMomentum.E();}
319 break; // case pE
320 case pMU: // all muons with eta < DET->MAX_MU
321 PairingMuon(recoMomentum,genMomentum,branchMuon);
322 if(recoMomentum.E() !=0){
323 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
324 elementMuon->OverPT = (1/genMomentum.Pt());
325 elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
326 break; // case pMU
327 default:
328 break;
329 } // switch (pid)
330 }
331
[9]332 } // while
[89]333
334 //compute missing transverse energy from calo towers
[258]335 TIter itCalo((TCollection*)branchTowers);
336 TRootCalo *calo;
337 itCalo.Reset();
[89]338 TLorentzVector Att(0.,0.,0.,0.);
339 float ScalarEt=0;
[258]340 while( (calo = (TRootCalo*) itCalo.Next()) )
341 {
342 //Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
343 Att.SetPtEtaPhiE(calo->PT,calo->Eta,calo->Phi,calo->E);
344 towers.push_back(Att);
[152]345 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
346 {
347 ScalarEt = ScalarEt + Att.Et();
348 PTmisReco = PTmisReco + Att;
349 }
350 }
[89]351
352 elementEtmis= (ETMIS*) branchetmis->NewEntry();
353 elementEtmis->Et = (PTmisGEN).Pt();
354 elementEtmis->Ex = (-PTmisGEN).Px();
355 elementEtmis->SEt = ScalarEt;
[88]356
[89]357 elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
358 elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
[9]359
360 //*****************************
[26]361
[88]362 sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
363
[39]364 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
[258]365cout<<"nombre de tau-jets "<<TausHadr.GetEntries()<<endl;
[39]366
[88]367 TLorentzVector JETreco(0,0,0,0);
368 for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
369 TLorentzVector JETgen(0,0,0,0);
370 JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
[258]371 PairingJet(JETreco,JETgen,branchJet);
[88]372 if(JETreco.Pt()>1)
[19]373 {
374 elementJet= (RESOLJET*) branchjet->NewEntry();
[88]375 elementJet->PT = JETgen.Et();
376 elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
[19]377 }
[39]378 }
[152]379 numTau = numTau+TausHadr.GetEntries();
[258]380
381 TIter itJet((TCollection*)branchJet);
382 TRootJet *jet;
383 itJet.Reset();
384 while( (jet = (TRootJet*) itJet.Next()) )
385 {
386 TLorentzVector JETT(0,0,0,0);
387 JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
388 if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
[39]389 {
390 for(Int_t i=0; i<TausHadr.GetEntries();i++)
391 {
392 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
393 {
394 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
[258]395 elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
396 elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
397 if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95
398 && (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
[152]399
[39]400 }
401 }
402 }
403
[26]404
405 } // for itJet : loop on all jets
406
[9]407 treeWriter->Fill();
408 } // Loop over all events
409 treeWriter->Write();
[260]410 float frac = numTauRec/numTau;
411 cout<<numTauRec<<endl;
412 cout<<numTau<<endl;
[152]413
[9]414 cout << "** Exiting..." << endl;
[152]415 cout<<frac<<endl;
416
[9]417
418 delete treeWriter;
419 delete treeReader;
420 delete DET;
421}
422
Note: See TracBrowser for help on using the repository browser.