Changeset 307 in svn for trunk/Delphes.cpp
- Timestamp:
- Mar 9, 2009, 9:47:21 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Delphes.cpp
r306 r307 47 47 #include "LHEFConverter.h" 48 48 #include "STDHEPConverter.h" 49 #include "LHCOConverter.h" 49 50 50 51 #include "SmearUtil.h" … … 77 78 int main(int argc, char *argv[]) 78 79 { 79 80 80 81 int appargc = 2; 81 82 char *appName= new char[20]; … … 87 88 delete [] appName; 88 89 delete [] appOpt; 89 90 90 91 if(argc != 3 && argc != 4 && argc != 5) { 91 92 cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl; … … 96 97 exit(1); 97 98 } 98 99 100 101 102 cout << endl << endl;103 104 cout <<"*********************************************************************"<< endl;105 cout <<"*********************************************************************"<< endl;106 cout <<"** **"<< endl;107 cout <<"** Welcome to **"<< endl;108 cout <<"** **"<< endl;109 cout <<"** **"<< endl;110 cout <<"** .ddddddd- lL hH **"<< endl;111 cout <<"** -Dd` `dD: Ll hH` **"<< endl;112 cout <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;113 cout <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;114 cout <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;115 cout <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;116 cout <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;117 cout <<"** Pp **"<< endl;118 cout <<"** **"<< endl;119 cout <<"** Delphes, a framework for the fast simulation **"<< endl;120 cout <<"** of a generic collider experiment **"<< endl;121 cout <<"** **"<< endl;122 cout <<"** --- Version 1.4beta of Delphes --- **"<< endl;123 cout <<"** Last date of change: 9 February 2009 **"<< endl;124 cout <<"** **"<< endl;125 cout <<"** **"<< endl;126 cout <<"** This package uses: **"<< endl;127 cout <<"** ------------------ **"<< endl;128 cout <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;129 cout <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;130 cout <<"** FROG: [hep-ex/0901.2718v1] **"<< endl;131 cout <<"** **"<< endl;132 cout <<"**-----------------------------------------------------------------**"<< endl;133 cout <<"** **"<< endl;134 cout <<"** Main authors: **"<< endl;135 cout <<"** ------------- **"<< endl;136 cout <<"** **"<< endl;137 cout <<"** Séverine Ovyn Xavier Rouby **"<< endl;138 cout <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;139 cout <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;140 cout <<"** Universite Catholique de Louvain (UCL) **"<< endl;141 cout <<"** Louvain-la-Neuve, Belgium **"<< endl;142 cout <<"** **"<< endl;143 cout <<"**-----------------------------------------------------------------**"<< endl;144 cout <<"** **"<< endl;145 cout <<"** Former Delphes versions and documentation can be found on : **"<< endl;146 cout <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;147 cout <<"** **"<< endl;148 cout <<"** **"<< endl;149 cout <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl;150 cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;151 cout <<"** give us your feedbacks about potential bugs **"<< endl;152 cout <<"** **"<< endl;153 cout <<"*********************************************************************"<< endl;154 cout <<"*********************************************************************"<< endl;155 156 // 1. ********** initialisation ***********157 99 100 101 102 103 cout << endl << endl; 104 105 cout <<"*********************************************************************"<< endl; 106 cout <<"*********************************************************************"<< endl; 107 cout <<"** **"<< endl; 108 cout <<"** Welcome to **"<< endl; 109 cout <<"** **"<< endl; 110 cout <<"** **"<< endl; 111 cout <<"** .ddddddd- lL hH **"<< endl; 112 cout <<"** -Dd` `dD: Ll hH` **"<< endl; 113 cout <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl; 114 cout <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl; 115 cout <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl; 116 cout <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl; 117 cout <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl; 118 cout <<"** Pp **"<< endl; 119 cout <<"** **"<< endl; 120 cout <<"** Delphes, a framework for the fast simulation **"<< endl; 121 cout <<"** of a generic collider experiment **"<< endl; 122 cout <<"** **"<< endl; 123 cout <<"** --- Version 1.4beta of Delphes --- **"<< endl; 124 cout <<"** Last date of change: 9 February 2009 **"<< endl; 125 cout <<"** **"<< endl; 126 cout <<"** **"<< endl; 127 cout <<"** This package uses: **"<< endl; 128 cout <<"** ------------------ **"<< endl; 129 cout <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl; 130 cout <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl; 131 cout <<"** FROG: [hep-ex/0901.2718v1] **"<< endl; 132 cout <<"** **"<< endl; 133 cout <<"**-----------------------------------------------------------------**"<< endl; 134 cout <<"** **"<< endl; 135 cout <<"** Main authors: **"<< endl; 136 cout <<"** ------------- **"<< endl; 137 cout <<"** **"<< endl; 138 cout <<"** Séverine Ovyn Xavier Rouby **"<< endl; 139 cout <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl; 140 cout <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl; 141 cout <<"** Universite Catholique de Louvain (UCL) **"<< endl; 142 cout <<"** Louvain-la-Neuve, Belgium **"<< endl; 143 cout <<"** **"<< endl; 144 cout <<"**-----------------------------------------------------------------**"<< endl; 145 cout <<"** **"<< endl; 146 cout <<"** Former Delphes versions and documentation can be found on : **"<< endl; 147 cout <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl; 148 cout <<"** **"<< endl; 149 cout <<"** **"<< endl; 150 cout <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl; 151 cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl; 152 cout <<"** give us your feedbacks about potential bugs **"<< endl; 153 cout <<"** **"<< endl; 154 cout <<"*********************************************************************"<< endl; 155 cout <<"*********************************************************************"<< endl; 156 157 // 1. ********** initialisation *********** 158 158 159 srand (time (NULL)); /* Initialisation du générateur */ 159 160 TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch; 160 161 globalwatch.Start(); 161 162 162 163 163 164 //read the output TROOT file … … 182 183 string DetDatacard("data/DetectorCard.dat"); //for detector smearing parameters 183 184 string TrigDatacard("data/TriggerCard.dat"); //for trigger selection 184 185 185 186 string lineCard1,lineCard2; 186 187 bool detecCard=false,trigCard=false; … … 190 191 infile1 >> lineCard1; // the first line determines the type of input files 191 192 if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==true) 192 193 cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl; 193 194 else if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[3]; detecCard=true;} 194 195 else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==true) 195 196 cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl; 196 197 else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[3]; trigCard=true;} 197 198 } … … 201 202 infile2 >> lineCard2; // the first line determines the type of input files 202 203 if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==true) 203 204 cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl; 204 205 else if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[4]; detecCard=true;} 205 206 else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==true) 206 207 cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl; 207 208 else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[4]; trigCard=true;} 208 209 } … … 221 222 << right << setw(2) <<"**"<<""<<endl; 222 223 cout <<"** **"<< endl; 223 224 224 225 //Trigger information 225 226 cout <<"** ########### Start reading TRIGGER card ########## **"<< endl; 226 227 if(trigCard==false) 227 228 229 230 228 { 229 cout <<"** WARNING: Datacard not found, use default card **" << endl; 230 TrigDatacard="data/TriggerCard.dat"; 231 } 231 232 TriggerTable *TRIGT = new TriggerTable(); 232 233 TRIGT->TriggerCardReader(TrigDatacard.c_str()); … … 248 249 //VFD information 249 250 VeryForward * VFD = new VeryForward(DET); 250 251 251 252 // data converters 252 253 DataConverter *converter=NULL; 253 254 cout <<"** **"<<endl; 254 255 cout <<"** ####### Start convertion to TRoot format ######## **"<< endl; 255 256 256 257 if(strstr(line.c_str(),".hep")) 257 258 { … … 273 274 } 274 275 else { 275 276 277 278 279 280 281 276 cerr << left << setw(4) <<"** "<<"" 277 << left << setw(63) << line.c_str() <<"" 278 << right << setw(2) <<"**"<<endl; 279 cerr <<"** ERROR: File format not identified -- Exiting... **"<< endl; 280 cout <<"** **"<< endl; 281 cout <<"*********************************************************************"<< endl; 282 return -1;}; 282 283 cout <<"** Exiting conversion... **"<< endl; 283 284 284 285 TChain chain("GEN"); 285 286 chain.Add(outputfilename.c_str()); 286 287 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain); 287 288 const TClonesArray *branchGen = treeReader->UseBranch("Particle"); 288 289 289 290 TIter itGen((TCollection*)branchGen); 290 291 … … 322 323 vector<D_Particle> muon; 323 324 vector<D_Particle> gamma; 324 325 bool FLAG_lhco = true; 326 325 326 vector<int> NTrackJet; 327 328 bool FLAG_lhco = true; 329 327 330 TSimpleArray<TRootGenParticle> NFCentralQ; 328 329 D_CaloList list_of_calorimeters;330 D_CaloElement CentralCalo("centralcalo",331 -DET->CEN_max_calo_cen, DET->CEN_max_calo_cen,332 DET->ELG_Ccen, DET->ELG_Ncen, DET->ELG_Scen,333 334 D_CaloElement ForwardCalo("forwardcalo",335 DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd,336 337 338 D_CaloElement BackwardCalo("backwardcalo",339 -DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen,340 341 342 //D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0);343 list_of_calorimeters.addElement(CentralCalo);344 list_of_calorimeters.addElement(ForwardCalo);345 list_of_calorimeters.addElement(BackwardCalo);346 //list_of_calorimeters.addElement(CastorCalo);347 list_of_calorimeters.sortElements();348 349 350 // 2. ********** Loop over all events ***********331 332 D_CaloList list_of_calorimeters; 333 D_CaloElement CentralCalo("centralcalo", 334 -DET->CEN_max_calo_cen, DET->CEN_max_calo_cen, 335 DET->ELG_Ccen, DET->ELG_Ncen, DET->ELG_Scen, 336 DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal); 337 D_CaloElement ForwardCalo("forwardcalo", 338 DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd, 339 DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd, 340 DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf ); 341 D_CaloElement BackwardCalo("backwardcalo", 342 -DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen, 343 DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd, 344 DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf ); 345 //D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0); 346 list_of_calorimeters.addElement(CentralCalo); 347 list_of_calorimeters.addElement(ForwardCalo); 348 list_of_calorimeters.addElement(BackwardCalo); 349 //list_of_calorimeters.addElement(CastorCalo); 350 list_of_calorimeters.sortElements(); 351 352 353 // 2. ********** Loop over all events *********** 351 354 Long64_t entry, allEntries = treeReader->GetEntries(); 352 355 cout <<"** **"<<endl; … … 355 358 << left << setw(15) << allEntries <<"" 356 359 << right << setw(2) <<"**"<<endl; 357 360 358 361 ExRootProgressBar *Progress = new ExRootProgressBar(allEntries); 359 362 360 363 loopwatch.Start(); 361 364 362 365 // loop on all events 363 366 for(entry = 0; entry < allEntries; ++entry) … … 376 379 towers.clear(); 377 380 input_particles.clear(); 378 381 NTrackJet.clear(); 382 379 383 { 380 381 382 383 384 385 386 387 388 389 390 391 392 393 384 D_CaloTowerList list_of_active_towers; 385 D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates 386 387 // 2.1a Loop over all particles in event, to fill the towers 388 itGen.Reset(); 389 GenParticle *particleG; 390 while( (particleG = (GenParticle*) itGen.Next()) ) { 391 392 TRootGenParticle *particle = new TRootGenParticle(particleG); 393 int pid = abs(particle->PID); 394 particle->Charge=ChargeVal(particle->PID); 395 particle->setFractions(); // init 396 397 394 398 // 2.1a.1********************* preparation for the b-tagging 395 399 //// This subarray is needed for the B-jet algorithm … … 401 405 NFCentralQ.Add(particle); 402 406 } 403 407 404 408 // 2.1a.2********************* visible particles only 405 409 if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) ){ 406 410 407 411 // 2.1a.2.1 Central solenoidal magnetic field 408 412 TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo … … 412 416 particle->Charge==0 || 413 417 (!DET->FLAG_bfield && particle->Charge!=0 && particle->PT > DET->TRACK_ptmin)) 414 if(415 (particle->EtaCalo > list_of_calorimeters.getEtamin() ) &&416 (particle->EtaCalo < list_of_calorimeters.getEtamax() )417 418 if( 419 (particle->EtaCalo > list_of_calorimeters.getEtamin() ) && 420 (particle->EtaCalo < list_of_calorimeters.getEtamax() ) 421 ) { 418 422 float iEta=UNDEFINED, iPhi=UNDEFINED; 419 423 DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta 420 424 if (iEta != UNDEFINED && iPhi != UNDEFINED) 421 422 D_CaloTower tower(iEta,iPhi); // new tower423 tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() );424 425 list_of_active_towers.addTower(tower);426 // this list may contain several times the same calotower, as several particles427 // may leave some energy in the same calotower428 // After the loop on particles, identical cells in the list should be merged429 } // iEta and iPhi must be defined430 }431 425 { 426 D_CaloTower tower(iEta,iPhi); // new tower 427 tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() ); 428 429 list_of_active_towers.addTower(tower); 430 // this list may contain several times the same calotower, as several particles 431 // may leave some energy in the same calotower 432 // After the loop on particles, identical cells in the list should be merged 433 } // iEta and iPhi must be defined 434 } 435 432 436 // 2.1a.2.3 charged particles in tracker: energy flow 433 437 // if bfield not simulated, pt should be high enough to be taken into account … … 435 439 if( particle->Charge !=0 && 436 440 fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available 437 441 ( DET->FLAG_bfield || 438 442 (!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin) 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 443 ) 444 ) { 445 // 2.1a.2.3.1 Filling the particle properties + smearing 446 // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction 447 // This is the EnergyFlow hypothesis 448 particle->SetEtaPhi(particle->Eta,particle->Phi); 449 float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks 450 451 // 2.1a.2.3.2 Muons 452 if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu) { 453 TLorentzVector p; 454 float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT ); 455 if (sPT > 0 && sPT > DET->PTCUT_muon) { 456 p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta)); 457 muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo)); 458 } 459 sET = (sPT >0)? sPT : 0; 460 } 461 // 2.1a.2.3.3 Electrons 462 else if (pid == pE) { 463 // Finds in which calorimeter the particle has gone, to know its resolution 464 465 D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo); 466 if(currentCalo.getName() == dummyCalo.getName()) { 467 cout << "** Warning: the calo coverage behind the tracker is not complete! **" << endl; } 468 469 // final smeared EM energy // electromagnetic fraction F_em =1 for electrons; 470 float sE = currentCalo.getElectromagneticResolution().Smear(particle->E); 471 if (sE>0) { 472 sET = sE/cosh(particle->Eta); 473 // NB: ET is found via the calorimetry and not via the track curvature 474 475 TLorentzVector p; 476 p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE); 477 if (sET > DET->PTCUT_elec) 478 electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo)); 479 } else { sET=0;} // if negative smeared energy -- needed for the tracks 480 } 481 // 2.1a.2.3.4 Other charged particles : smear them for the tracks! 482 else { //other particles 483 D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo); 484 float sEem = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem()); 485 float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad()); 486 float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 ); 487 sET = sE/cosh(particle->EtaCalo); 488 } 489 490 // 2.1a.2.3.5 Tracks 491 if( (rand()%100) < DET->TRACK_eff && sET!=0) { 492 elementTrack = (TRootTracks*) branchTrack->NewEntry(); 493 elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge); 494 //TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex! 495 TrackCentral.push_back(*elementTrack); // tracks at vertex! 496 // TODO!!! apply a smearing on the position of the origin of the track 497 // TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout); 498 } 495 499 } // 2.1a.2.3 : if tracker/energy-flow 496 500 // 2.1a.2.4 Photons 497 501 // stays in the tracker -> track available -> gamma ID 498 502 else if( (pid == pGAMMA) && fabs(particle->EtaCalo)< DET->CEN_max_tracker ) { 499 500 501 502 503 503 float iEta=UNDEFINED, iPhi=UNDEFINED; 504 DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta 505 D_CaloTower tower(iEta,iPhi); 506 // stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search 507 list_of_towers_with_photon.addTower(tower); 504 508 } 505 509 // 2.1a.2.5 : very forward detectors … … 507 511 { 508 512 if (DET->FLAG_RP==1) { 509 510 511 512 513 // for the moment, only protons are transported 514 // BUT !!! could be a beam of other particles! (heavy ions?) 515 // BUT ALSO !!! if very forward muons, or others! 516 VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle); 513 517 } 514 518 // 2.1a.2.6: Zero degree calorimeter 515 519 if(DET->FLAG_vfd==1) { 516 520 VFD->ZDC(treeWriter,branchZDC,particle); 517 521 } 518 522 } 519 523 520 524 } // 2.1a.2 : if visible particle 521 522 525 delete particle; 526 } // loop on all particles 2.1a 523 527 524 525 526 527 528 529 530 531 532 533 float iEta = list_of_active_towers[i].getEta();534 float iPhi = list_of_active_towers[i].getPhi();535 float e = list_of_active_towers[i].getE();536 if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) {537 538 elementCalo->set(list_of_active_towers[i]);539 540 TLorentzVector p;541 p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e );542 543 544 }545 546 547 548 549 550 551 552 553 528 529 // 2.1b loop on all (activated) towers 530 // at this stage, list_of_active_towers may contain several times the same tower 531 // first step is to merge identical towers, by matching their (iEta,iPhi) 532 list_of_active_towers.mergeDuplicates(); 533 // Calotower smearing 534 list_of_active_towers.smearTowers(list_of_calorimeters); 535 536 for(unsigned int i=0; i<list_of_active_towers.size(); i++) { 537 float iEta = list_of_active_towers[i].getEta(); 538 float iPhi = list_of_active_towers[i].getPhi(); 539 float e = list_of_active_towers[i].getE(); 540 if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) { 541 elementCalo = (TRootCalo*) branchCalo->NewEntry(); 542 elementCalo->set(list_of_active_towers[i]); 543 // not beautiful : should be improved! 544 TLorentzVector p; 545 p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e ); 546 PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E())); 547 towers.push_back(Tower); 548 } 549 } // loop on towers 550 551 // 2.1c photon ID 552 // list_of_towers_with_photon is the list of towers with photon candidates 553 // already smeared ! 554 // sorts the vector and smears duplicates 555 list_of_towers_with_photon.mergeDuplicates(); 556 557 for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) { 554 558 float eta = list_of_towers_with_photon[i].getEta(); 555 559 float phi = list_of_towers_with_photon[i].getPhi(); … … 560 564 if (cal.getET() > DET->PTCUT_gamma) { gamma.push_back(D_Particle(p,pGAMMA,p.Eta(),p.Phi())); } 561 565 } 562 563 564 } // IF NEW ALGORITHM with energy flow565 566 567 568 566 } // for -- list of photons 567 568 } // IF NEW ALGORITHM with energy flow 569 570 571 572 569 573 // 2.2 ********** Output preparation & complex objects *********** 570 574 // 2.2.1 ********************* sorting collections by decreasing pt … … 599 603 Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E); 600 604 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd) 601 {602 603 604 605 606 }605 { 606 PTmis = PTmis + Att; 607 // create a fastjet::PseudoJet with these components and put it onto 608 // back of the input_particles vector 609 input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E)); 610 } 607 611 } 608 612 elementEtmis = (TRootETmis*) branchETmis->NewEntry(); … … 613 617 614 618 // 2.2.3 ************* jets, B-tag, tau jets 615 sorted_jets=JETRUN->RunJets(input_particles); 616 JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ); 617 JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral); 619 vector<int> NTrackJet; 620 sorted_jets=JETRUN->RunJets(input_particles, TrackCentral,NTrackJet); 621 JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ,NTrackJet); 622 JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral,NTrackJet); 618 623 619 624 treeWriter->Fill(); 620 625 } // 2. Loop over all events ('for' loop) 621 622 626 627 623 628 cout <<"** **"<< endl; 624 629 cout <<"** Exiting detector simulation... **"<< endl; 625 626 630 631 627 632 treeWriter->Write(); 628 633 delete treeWriter; 629 634 loopwatch.Stop(); 630 631 632 635 636 637 633 638 // 3. ********** Trigger & Frog *********** 634 639 // 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard … … 636 641 if(DET->FLAG_trigger == 1) 637 642 { 638 639 640 643 cout <<"** **"<<endl; 644 cout <<"** ########### Start Trigger selection ########### **"<< endl; 645 641 646 // input 642 647 TChain chainT("Analysis"); … … 654 659 ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger"); 655 660 ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class()); 656 657 661 662 658 663 Long64_t entryT, allEntriesT = treeReaderT->GetEntries(); 659 664 // loop on all entries 660 665 for(entryT = 0; entryT < allEntriesT; ++entryT) { 661 662 663 664 666 treeWriterT->Clear(); 667 treeReaderT->ReadEntry(entryT); 668 TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger); 669 treeWriterT->Fill(); 665 670 } // loop on all entries 666 671 cout <<"** Exiting trigger simulation... **"<< endl; … … 670 675 delete treeReaderT; 671 676 } // trigger 672 673 674 677 triggerwatch.Stop(); 678 679 675 680 // 3.2 ************** FROG display 676 681 frogwatch.Start(); 677 682 if(DET->FLAG_frog == 1) { 678 679 680 681 682 683 684 683 cout <<"** **"<<endl; 684 cout <<"** ################## Start FROG ################# **"<< endl; 685 686 FrogDisplay *FROG = new FrogDisplay(DET); 687 FROG->BuildEvents(outputfilename); 688 FROG->BuildGeom(); 689 delete FROG; 685 690 } 686 691 frogwatch.Stop(); 687 692 688 689 690 691 // 4. ********** End & Exit *********** 692 693 // 3.3 *************** LHCO output 694 if(DET->FLAG_lhco == 1){ 695 cout <<"** **"<<endl; 696 cout <<"** ############ Start LHCO convertion ############ **"<< endl; 697 /* 698 //create output log-file name 699 string forLogLHCO = outputfilename; 700 string LogNameLHCO = forLogLHCO.erase(forLogLHCO.find(".root")); 701 LogNameLHCO = LogNameLHCO+"_events.lhco"; 702 */ 703 704 //LHCOConverter *LHCO = new LHCOConverter(outputfilename,LogNameLHCO); 705 LHCOConverter *LHCO = new LHCOConverter(outputfilename,""); 706 LHCO->CopyRunLogFile(); 707 LHCO->ConvertExRootAnalysisToLHCO(); 708 delete LHCO; 709 } 710 711 712 713 // 4. ********** End & Exit *********** 714 693 715 globalwatch.Stop(); 694 716 cout <<"** **"<< endl; 695 717 cout <<"** ################## Time report ################# **"<< endl; 696 718 cout << left << setw(32) <<"** Time report for "<<"" 697 698 719 << left << setw(15) << allEntries <<"" 720 << right << setw(22) <<"events **"<<endl; 699 721 cout <<"** **"<< endl; 700 722 cout << left << setw(10) <<"**"<<"" … … 715 737 if(DET->FLAG_trigger == 1) 716 738 { 717 cout << left << setw(10) <<"**"<<""718 719 720 721 739 cout << left << setw(10) <<"**"<<"" 740 << left << setw(15) <<" + Trigger:"<<"" 741 << right << setw(15) <<triggerwatch.CpuTime()<<"" 742 << right << setw(15) <<triggerwatch.RealTime()<<"" 743 << right << setw(14) <<"**"<<endl; 722 744 } 723 745 if(DET->FLAG_frog == 1) 724 746 { 725 cout << left << setw(10) <<"**"<<""726 727 728 729 730 731 } 732 733 734 735 736 737 738 739 747 cout << left << setw(10) <<"**"<<"" 748 << left << setw(15) <<" + Frog:"<<"" 749 << right << setw(15) <<frogwatch.CpuTime()<<"" 750 << right << setw(15) <<frogwatch.RealTime()<<"" 751 << right << setw(14) <<"**"<<endl; 752 753 } 754 755 cout <<"** **"<< endl; 756 cout <<"** Exiting Delphes ... **"<< endl; 757 cout <<"** **"<< endl; 758 cout <<"*********************************************************************"<< endl; 759 cout <<"*********************************************************************"<< endl; 760 761 740 762 delete treeReader; 741 763 delete DET; … … 746 768 delete converter; 747 769 748 // todo("TODO");770 // todo("TODO"); 749 771 }
Note:
See TracChangeset
for help on using the changeset viewer.