- Timestamp:
- Dec 4, 2008, 3:23:03 PM (16 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Resolutions.cpp
r43 r88 29 29 30 30 #include "interface/SmearUtil.h" 31 #include "interface/JetUtils.h" 32 #include "interface/BFieldProp.h" 33 31 34 #include "Utilities/Fastjet/include/fastjet/PseudoJet.hh" 32 35 #include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh" 33 34 // get info on how fastjet was configured35 #include "Utilities/Fastjet/include/fastjet/config.h"36 37 // make sure we have what is needed38 #ifdef ENABLE_PLUGIN_SISCONE39 # include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"40 #endif41 #ifdef ENABLE_PLUGIN_CDFCONES42 # include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"43 # include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"44 #endif45 36 46 37 #include<vector> … … 109 100 } 110 101 111 112 102 int main(int argc, char *argv[]) 113 103 { … … 133 123 return -1; 134 124 } 125 135 126 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after 136 127 outputFile->Close(); … … 196 187 RESOLution *DET = new RESOLution(); 197 188 DET->ReadDataCard(DetDatacard); 198 199 TLorentzVector genMomentum(0,0,0,0); 189 190 //Jet information 191 JetsUtil *JETRUN = new JetsUtil(); 192 193 //Propagation of tracks in the B field 194 TrackPropagation *TRACP = new TrackPropagation(); 195 196 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information 197 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information 200 198 LorentzVector jetMomentum; 201 199 vector<TLorentzVector> TrackCentral; 202 200 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; 201 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm 202 vector<fastjet::PseudoJet> sorted_jetsGEN; 203 204 vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm 205 vector<fastjet::PseudoJet> sorted_jetsReco; 214 206 215 207 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 definition225 #ifdef ENABLE_PLUGIN_CDFCONES226 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 #else229 plugins = NULL;230 #endif231 232 // set up a CDF midpoint jet definition233 #ifdef ENABLE_PLUGIN_CDFCONES234 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 #else237 pluginsS = NULL;238 #endif239 240 // set up a CDF midpoint jet definition241 #ifdef ENABLE_PLUGIN_CDFCONES242 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 #else245 pluginsT = NULL;246 #endif247 248 208 249 209 // Loop over all events … … 252 212 for(entry = 0; entry < allEntries; ++entry) 253 213 { 254 TLorentzVector PTmis S(0,0,0,0);255 TLorentzVector PTmis (0,0,0,0);214 TLorentzVector PTmisReco(0,0,0,0); 215 TLorentzVector PTmisGEN(0,0,0,0); 256 216 treeReader->ReadEntry(entry); 257 217 treeWriter->Clear(); … … 263 223 TSimpleArray<TRootGenParticle> NFCentralQ; 264 224 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(); 225 input_particlesReco.clear(); 226 input_particlesGEN.clear(); 276 227 towers.clear(); 277 228 … … 283 234 int pid = abs(particle->PID); 284 235 float eta = fabs(particle->Eta); 285 236 237 //input generator level particle for jet algorithm 286 238 if(particle->Status == 1 && eta < DET->MAX_CALO_FWD) 287 239 { 288 input_particles .push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));240 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 289 241 } 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; 242 243 //Calculate ETMIS from generated particles 244 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum; 245 246 if( (particle->Status == 1) && 247 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) && 248 (fabs(particle->Eta) < DET->MAX_CALO_FWD) 249 ) 250 { 251 recoMomentum = genMomentum; 252 //use of the magnetic field propagation 253 TRACP->Propagation(particle,recoMomentum); 254 float eta=fabs(recoMomentum.Eta()); 255 303 256 switch(pid) { 257 304 258 case pE: // all electrons with eta < DET->MAX_CALO_FWD 305 nonS = genMomentum.E(); 306 DET->SmearElectron(genMomentum); 259 DET->SmearElectron(recoMomentum); 307 260 if(eta < DET->MAX_TRACKER){ 308 261 elementElec=(RESOLELEC*) branchelec->NewEntry(); 309 elementElec-> NonSmeareE = nonS;310 elementElec->Smeare E = genMomentum.E();}262 elementElec->E = genMomentum.E(); 263 elementElec->SmearedE = recoMomentum.E();} 311 264 break; // case pE 312 265 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD 313 DET->SmearElectron( genMomentum);266 DET->SmearElectron(recoMomentum); 314 267 break; // case pGAMMA 315 268 case pMU: // all muons with eta < DET->MAX_MU 269 DET->SmearMu(recoMomentum); 316 270 elementMuon = (RESOLMUON*) branchmuon->NewEntry(); 317 elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt()); 318 DET->SmearMu(genMomentum); 319 elementMuon->OneoverSmearePT = (1/genMomentum.Pt()); 271 elementMuon->OverPT = (1/genMomentum.Pt()); 272 elementMuon->OverSmearedPT = (1/recoMomentum.Pt()); 320 273 break; // case pMU 321 274 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD 322 275 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD 323 DET->SmearHadron( genMomentum, 0.7);276 DET->SmearHadron(recoMomentum, 0.7); 324 277 break; // case hadron 325 278 default: // all other final particles with eta < DET->MAX_CALO_FWD 326 DET->SmearHadron( genMomentum, 1.0);279 DET->SmearHadron(recoMomentum, 1.0); 327 280 break; 328 281 } // 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())); 282 283 //information to reconstruct jets from reconstructed towers 284 int charge=Charge(pid); 285 if(recoMomentum.E() !=0 && pid != pMU) { 286 if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->PT_TRACKS_MIN)){ 287 towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()))); 288 input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())); 334 289 } 290 } 335 291 336 292 // 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); 293 if( 294 (genMomentum.E()!=0) && 295 (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) && 296 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account 297 ((rand()%100) < DET->TRACKING_EFF) && 298 (charge!=0) 299 ) 300 {TrackCentral.push_back(genMomentum);} 301 351 302 } // switch 352 303 } // while 353 304 305 //compute missing transverse energy from calo towers 354 306 TLorentzVector Att(0.,0.,0.,0.); 355 307 for(unsigned int i=0; i < towers.size(); i++) 356 308 { 357 309 Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E); 358 PTmis S = PTmisS+ Att;310 PTmisReco = PTmisReco + Att; 359 311 } 360 312 361 313 elementEtmis= (ETMIS*) branchetmis->NewEntry(); 362 elementEtmis->Et = (PTmis ).Pt();363 elementEtmis-> EtSmeare = (PTmisS).Pt()-(PTmis).Pt();314 elementEtmis->Et = (PTmisGEN).Pt(); 315 elementEtmis->SmearedEt = (PTmisReco).Pt()-(PTmisGEN).Pt(); 364 316 365 317 //***************************** 366 318 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 319 sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN); 320 sorted_jetsReco=JETRUN->RunJets(input_particlesReco); 321 382 322 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen); 383 323 384 TLorentzVector JET Sm(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(JET Sm,JET,sorted_jetsS);389 if(JET Sm.Pt()>1)324 TLorentzVector JETreco(0,0,0,0); 325 for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) { 326 TLorentzVector JETgen(0,0,0,0); 327 JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E()); 328 PairingJet(JETreco,JETgen,sorted_jetsReco); 329 if(JETreco.Pt()>1) 390 330 { 391 331 elementJet= (RESOLJET*) branchjet->NewEntry(); 392 elementJet-> NonSmearePT = JET.Et();393 elementJet->Smeare PT = JETSm.Et()/JET.Et();332 elementJet->PT = JETgen.Et(); 333 elementJet->SmearedPT = JETreco.Et()/JETgen.Et(); 394 334 } 395 335 } 396 for (unsigned int i = 0; i < sorted_jets S.size(); i++) {336 for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) { 397 337 TLorentzVector JETT(0,0,0,0); 398 JETT.SetPxPyPzE(sorted_jets S[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());338 JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E()); 399 339 if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) 400 340 { -
trunk/interface/TreeClasses.h
r17 r88 13 13 { 14 14 public: 15 Float_t Smeare PT;16 Float_t NonSmearePT;15 Float_t SmearedPT; 16 Float_t PT; 17 17 18 18 static TCompare *fgCompare; //! … … 30 30 { 31 31 public: 32 Float_t SmeareE;33 Float_t NonSmeareE;32 Float_t E; 33 Float_t SmearedE; 34 34 35 35 static TCompare *fgCompare; //! … … 46 46 { 47 47 public: 48 Float_t O neoverSmearePT;49 Float_t O neoverNonSmearePT;48 Float_t OverSmearedPT; 49 Float_t OverPT; 50 50 51 51 static TCompare *fgCompare; //! … … 77 77 public: 78 78 Float_t Et; 79 Float_t EtSmeare;79 Float_t SmearedEt; 80 80 81 81 static TCompare *fgCompare; //!
Note:
See TracChangeset
for help on using the changeset viewer.