- Timestamp:
- Jan 7, 2009, 12:23:20 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Resolutions.cpp
r100 r152 181 181 ETMIS *elementEtmis; 182 182 183 int numTau=0; 184 int numTauRec=0; 183 185 184 186 //read the datacard input file … … 195 197 196 198 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information 199 TLorentzVector recoMomentumCalo(0,0,0,0); 197 200 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information 198 201 LorentzVector jetMomentum; … … 206 209 207 210 vector<PhysicsTower> towers; 208 211 212 float iPhi=0,iEta=0; 213 209 214 // Loop over all events 210 215 Long64_t entry, allEntries = treeReader->GetEntries(); … … 251 256 recoMomentum = genMomentum; 252 257 //use of the magnetic field propagation 253 TRACP->Propagation(particle,recoMomentum);258 //TRACP->Propagation(particle,recoMomentum); 254 259 // cout<<"eta "<<eta<<endl; 255 260 eta=fabs(recoMomentum.Eta()); … … 285 290 286 291 //information to reconstruct jets from reconstructed towers 287 int charge=Charge(pid); 288 if(recoMomentum.E() !=0 && pid != pMU) { 289 if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->TRACK_ptmin)){ 290 towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()))); 291 input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())); 292 } 293 } 292 int charge=Charge(pid); 293 if(recoMomentum.E() !=0 && pid != pMU) { 294 if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->TRACK_ptmin)){ 295 DET->BinEtaPhi(recoMomentum.Phi(), recoMomentum.Eta(), iPhi, iEta); 296 if(iEta != -100 && iPhi != -100) 297 { 298 recoMomentumCalo.SetPtEtaPhiE(recoMomentum.Pt(),iEta,iPhi,recoMomentum.E()); 299 300 PhysicsTower Tower(LorentzVector(recoMomentumCalo.Px(),recoMomentumCalo.Py(),recoMomentumCalo.Pz(), recoMomentumCalo.E())); 301 towers.push_back(Tower); 302 } 303 } 304 } 294 305 295 306 // all final charged particles 296 307 if( 297 ( genMomentum.E()!=0) &&298 (fabs( genMomentum.Eta()) < DET->CEN_max_tracker) &&299 ( genMomentum.Pt() > DET->TRACK_ptmin ) && // pt too small to be taken into account308 (recoMomentum.E()!=0) && 309 (fabs(recoMomentum.Eta()) < DET->CEN_max_tracker) && 310 (recoMomentum.Pt() > DET->TRACK_ptmin ) && // pt too small to be taken into account 300 311 ((rand()%100) < DET->TRACK_eff) && 301 312 (charge!=0) 302 313 ) 303 {TrackCentral.push_back( genMomentum);}314 {TrackCentral.push_back(recoMomentum);} 304 315 305 316 } // switch … … 309 320 TLorentzVector Att(0.,0.,0.,0.); 310 321 float ScalarEt=0; 311 for(unsigned int i=0; i < towers.size(); i++) 312 { 313 Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E); 314 ScalarEt = ScalarEt + Att.Et(); 315 PTmisReco = PTmisReco + Att; 316 } 322 for(unsigned int i=0; i < towers.size(); i++) 323 { 324 Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E); 325 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd) 326 { 327 ScalarEt = ScalarEt + Att.Et(); 328 PTmisReco = PTmisReco + Att; 329 // create a fastjet::PseudoJet with these components and put it onto 330 // back of the input_particles vector 331 input_particlesReco.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E)); 332 } 333 } 317 334 318 335 elementEtmis= (ETMIS*) branchetmis->NewEntry(); … … 343 360 } 344 361 } 362 numTau = numTau+TausHadr.GetEntries(); 345 363 for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) { 346 364 TLorentzVector JETT(0,0,0,0); … … 355 373 elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E()); 356 374 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi()); 375 if( (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E()) > 0.95 376 && (DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi()))==1)numTauRec++; 377 357 378 } 358 379 } … … 365 386 } // Loop over all events 366 387 treeWriter->Write(); 367 388 float frac = numTauRec/numTau; 389 cout<<numTauRec<<endl; 390 cout<<numTau<<endl; 391 368 392 cout << "** Exiting..." << endl; 393 cout<<frac<<endl; 394 369 395 370 396 delete treeWriter;
Note:
See TracChangeset
for help on using the changeset viewer.