/*********************************************************************** ** ** ** /----------------------------------------------\ ** ** | Delphes, a framework for the fast simulation | ** ** | of a generic collider experiment | ** ** \----------------------------------------------/ ** ** ** ** ** ** This package uses: ** ** ------------------ ** ** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] ** ** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] ** ** FROG: [hep-ex/0901.2718v1] ** ** ** ** ------------------------------------------------------------------ ** ** ** ** Main authors: ** ** ------------- ** ** ** ** Severine Ovyn Xavier Rouby ** ** severine.ovyn@uclouvain.be xavier.rouby@cern ** ** ** ** Center for Particle Physics and Phenomenology (CP3) ** ** Universite catholique de Louvain (UCL) ** ** Louvain-la-Neuve, Belgium ** ** ** ** Copyright (C) 2008-2009, ** ** All rights reserved. ** ** ** ***********************************************************************/ /// \file Delphes.cpp /// \brief Executable for Delphes #include "TChain.h" #include "TApplication.h" #include "TStopwatch.h" #include "TFile.h" #include "ExRootTreeReader.h" #include "ExRootTreeWriter.h" #include "ExRootTreeBranch.h" #include "ExRootProgressBar.h" #include "DataConverter.h" #include "HEPEVTConverter.h" #include "LHEFConverter.h" #include "STDHEPConverter.h" #include "SmearUtil.h" #include "CaloUtil.h" #include "BFieldProp.h" #include "TriggerUtil.h" #include "VeryForward.h" #include "JetsUtil.h" #include "FrogUtil.h" #include #include using namespace std; //------------------------------------------------------------------------------ void todo(string filename) { ifstream infile(filename.c_str()); cout << "** TODO list ..." << endl; while(infile.good()) { string temp; getline(infile,temp); cout << "*" << temp << endl; } cout << "** done...\n"; } //------------------------------------------------------------------------------ int main(int argc, char *argv[]) { unsigned int nnnn=0, mmmm=0; int appargc = 2; char *appName= new char[20]; char *appOpt= new char[20]; sprintf(appName,"Delphes"); sprintf(appOpt,"-b"); char *appargv[] = {appName,appOpt}; TApplication app(appName, &appargc, appargv); delete [] appName; delete [] appOpt; if(argc != 3 && argc != 4 && argc != 5) { cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl; cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl; cout << " output_file - output file." << endl; cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "< outputfilename.length()) { cout <<"** ERROR: 'output_file' should be a .root file. Exiting... **"<< endl; exit(1); } //create output log-file name string forLog = outputfilename; string LogName = forLog.erase(forLog.find(".root")); LogName = LogName+"_run.log"; TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after outputFile->Close(); string line; ifstream infile(inputFileList.c_str()); infile >> line; // the first line determines the type of input files //read the datacard input file string DetDatacard(""); //for detector smearing parameters string TrigDatacard("data/trigger.dat"); //for trigger selection string lineCard1,lineCard2; bool detecCard=false,trigCard=false; if(argv[3]) { ifstream infile1(argv[3]); infile1 >> lineCard1; // the first line determines the type of input files if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==true) cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl; else if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[3]; detecCard=true;} else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==true) cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl; else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[3]; trigCard=true;} } if(argv[4]) { ifstream infile2(argv[4]); infile2 >> lineCard2; // the first line determines the type of input files if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==true) cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl; else if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[4]; detecCard=true;} else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==true) cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl; else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[4]; trigCard=true;} } //Smearing information RESOLution *DET = new RESOLution(); cout <<"** **"<< endl; cout <<"** ####### Start reading DETECTOR parameters ####### **"<< endl; cout << left << setw(40) <<"** Opening configuration card: "<<"" << left << setw(27) << DetDatacard <<"" << right << setw(2) <<"**"<<""<ReadDataCard(DetDatacard); DET->Logfile(LogName); cout << left << setw(42) <<"** Parameters summarised in file: "<<"" << left << setw(27) << LogName <<"" << right << setw(2) <<"**"<<""<TriggerCardReader(TrigDatacard.c_str()); TRIGT->PrintTriggerTable(LogName); if(DET->FLAG_trigger == 1) { cout << left << setw(40) <<"** Opening configuration card: "<<"" << left << setw(27) << TrigDatacard <<"" << right << setw(2) <<"**"<<""<UseBranch("Particle"); TIter itGen((TCollection*)branchGen); //Output file : contents of the analysis object data ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis"); ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class()); ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class()); ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class()); ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class()); ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class()); //ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class()); ExRootTreeBranch *branchTrack = treeWriter->NewBranch("Tracks", D_Track::Class()); ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class()); ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class()); ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class()); ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class()); ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class()); TRootETmis *elementEtmis; TRootElectron *elementElec; TRootMuon *elementMu; TRootPhoton *elementPhoton; D_Track * elementTrack; TRootCalo *elementCalo; TLorentzVector genMomentum(0,0,0,0); // four-momentum at the vertex TLorentzVector genMomentumBfield(0,0,0,0); // four-momentum at the exit of the tracks TLorentzVector momentumCaloSegmentation(0,0,0,0); // four-momentum in the calo, after applying the calo segmentation LorentzVector jetMomentum; vector input_particles;//for FastJet algorithm vector sorted_jets; vector TrackCentral; vector towers; //vector electron; vector electron; //vector muon; vector muon; //vector gamma; vector gamma; TSimpleArray NFCentralQ; D_CaloList list_of_calorimeters; D_CaloElement CentralCalo("centralcalo", -DET->CEN_max_calo_cen, DET->CEN_max_calo_cen, DET->ELG_Ccen, DET->ELG_Ncen, DET->ELG_Scen, DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal); D_CaloElement ForwardCalo("forwardcalo", DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd, DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd, DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf ); D_CaloElement BackwardCalo("backwardcalo", -DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen, DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd, DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf ); //D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0); list_of_calorimeters.addElement(CentralCalo); list_of_calorimeters.addElement(ForwardCalo); list_of_calorimeters.addElement(BackwardCalo); //list_of_calorimeters.addElement(CastorCalo); list_of_calorimeters.sortElements(); // 2. ********** Loop over all events *********** Long64_t entry, allEntries = treeReader->GetEntries(); cout <<"** **"<Update(entry); TLorentzVector PTmis(0,0,0,0); treeReader->ReadEntry(entry); treeWriter->Clear(); electron.clear(); muon.clear(); gamma.clear(); NFCentralQ.Clear(); TrackCentral.clear(); towers.clear(); input_particles.clear(); /* if(0) { // OLD SMEARING ALGORITHM without energy flow // 2.1 Loop over all particles in event itGen.Reset(); while( (particle = (TRootGenParticle*) itGen.Next()) ) { int pid = abs(particle->PID); float iPhi=0,iEta=0; // 2.1.1********************* preparation for the b-tagging //// This subarray is needed for the B-jet algorithm // optimization for speed : put first PID condition, then ETA condition, then either pt or status if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ? fabs(particle->Eta) < DET->CEN_max_tracker && particle->Status != 1 && particle->PT > DET->PT_QUARKS_MIN ) { NFCentralQ.Add(particle); } // 2.1.2 ********************* central detector: keep only visible particles // keeps only final particles, visible by the central detector, including the fiducial volume // the ordering of conditions have been optimised for speed : put first the STATUS condition if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && (fabs(particle->Eta) < DET->CEN_max_calo_fwd) ) { genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E); genMomentumBfield = genMomentum; // 2.1.2.1 ********************* central detector: magnetic field // genMomentumBfield is then changed with respect to the magnetic field if(DET->FLAG_bfield==1) TRACP->Propagation(particle,genMomentumBfield); float eta=fabs(genMomentumBfield.Eta()); // 2.1.2.2 ********************* central detector: smearing (calorimeters, muon chambers) switch(pid) { case pE: // all electrons with eta < DET->MAX_CALO_FWD DET->SmearElectron(genMomentumBfield); if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_elec){ electron.push_back(ParticleUtil(genMomentumBfield,particle->PID)); } break; // case pE case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD DET->SmearElectron(genMomentumBfield); if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_gamma) { gamma.push_back(ParticleUtil(genMomentumBfield,particle->PID)); } break; // case pGAMMA case pMU: // all muons with eta < DET->MAX_MU DET->SmearMu(genMomentumBfield); if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_mu && genMomentumBfield.Pt() > DET->PTCUT_muon){ muon.push_back(ParticleUtil(genMomentumBfield,particle->PID)); } break; // case pMU case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD case pK0S: // all K0s with eta < DET->MAX_CALO_FWD DET->SmearHadron(genMomentumBfield, 0.7); break; // case hadron default: // all other final particles with eta < DET->MAX_CALO_FWD DET->SmearHadron(genMomentumBfield, 1.0); break; } // 2.1.2.2 switch (pid) // 2.1.2.3 ********************* central detector: calotowers // all final particles but muons and neutrinos // for calorimetric towers and missing PT int charge=Charge(pid); if(genMomentumBfield.E() !=0 && pid != pMU) { // in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) {}// particules do not reach calos else { // particles reach calos // applies the calo segmentation and returns iEta & iPhi DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta); if(iEta != UNDEFINED && iPhi != UNDEFINED) { momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E()); elementCalo = (TRootCalo*) branchCalo->NewEntry(); elementCalo->Set(momentumCaloSegmentation); PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E())); towers.push_back(Tower); } // if iEta != UNDEFINED } // else : when particles reach the calos } // 2.1.2.3 calotowers // 2.1.2.4 ********************* central detector: tracks // all final charged particles if( (genMomentumBfield.E()!=0) && (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) && (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentum.Pt() > DET->TRACK_ptmin )) && // if bfield not simulated, pt should be high enough to be taken into account ((rand()%100) < DET->TRACK_eff) && (charge!=0) ) { elementTracks = (TRootTracks*) branchTracks->NewEntry(); elementTracks->Set(genMomentum); // fills px,py,pz,pt,e,eta,phi at vertex elementTracks->Etaout = genMomentumBfield.Eta(); elementTracks->Phiout = genMomentumBfield.Phi(); // TODO!!! apply a smearing on the position of the origin of the track // elementTracks->SetPosition(particle->X,particle->Y,particle->Z); // uses the output of the bfield computation : Xout=... Yout=... Zout... // TODO!!! elementTrakcs->SetPositionOut(Xout,Yout,Zout); TrackCentral.push_back(genMomentum); // tracks at vertex! } // 2.1.2.4 tracks } // 2.1.2 central detector // 2.1.3 ********************* forward detectors: zdc if(DET->FLAG_vfd==1) { VFD->ZDC(treeWriter,branchZDC,particle); VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle); } } // IF OLD ALGORITHM (= no energy flow ) else // IF NEW ALGORITHM with energy flow */ { D_CaloTowerList list_of_active_towers; D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates // 2.1a Loop over all particles in event, to fill the towers itGen.Reset(); TRootGenParticle *particle; while( (particle = (TRootGenParticle*) itGen.Next()) ) { int pid = abs(particle->PID); particle->Charge=Charge(pid); particle->setFractions(); // init // 2.1a.1********************* preparation for the b-tagging //// This subarray is needed for the B-jet algorithm // optimization for speed : put first PID condition, then ETA condition, then either pt or status if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ? fabs(particle->Eta) < DET->CEN_max_tracker && particle->Status != 1 && particle->PT > DET->PT_QUARKS_MIN ) { NFCentralQ.Add(particle); } // 2.1a.2********************* visible particles only if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) ){ // 2.1a.2.1 Central solenoidal magnetic field TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo particle->SetEtaPhi(particle->EtaCalo,particle->PhiCalo); // ???? check this line // 2.1a.2.2 Filling the calorimetric towers -- includes also forward detectors ? // first checks if the charged particles reach the calo! if( DET->FLAG_bfield || particle->getCharge()==0 || (!DET->FLAG_bfield && particle->getCharge()!=0 && particle->PT > DET->TRACK_ptmin)) if( (particle->EtaCalo > list_of_calorimeters.getEtamin() ) && (particle->EtaCalo < list_of_calorimeters.getEtamax() ) ) { float iEta=UNDEFINED, iPhi=UNDEFINED; DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta if (iEta != UNDEFINED && iPhi != UNDEFINED) { D_CaloTower tower(iEta,iPhi); // new tower tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() ); list_of_active_towers.addTower(tower); // this list may contain several times the same calotower, as several particles // may leave some energy in the same calotower // After the loop on particles, identical cells in the list should be merged } // iEta and iPhi must be defined } // 2.1a.2.3 charged particles in tracker: energy flow // if bfield not simulated, pt should be high enough to be taken into account // it is supposed here that DET->MAX_calo > DET->CEN_max_tracker > DET->CEN_max_mu > 0 if( particle->getCharge() !=0 && fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available ( DET->FLAG_bfield || (!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin) ) ) { // 2.1a.2.3.1 Filling the particle properties + smearing // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction // This is the EnergyFlow hypothesis particle->SetEtaPhi(particle->Eta,particle->Phi); float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks // 2.1a.2.3.2 Muons if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu && particle->PT > DET->PTCUT_muon ) { TLorentzVector p; float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT ); if (sPT > 0 && sPT > DET->PTCUT_muon) { p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta)); muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo)); } sET = (sPT >0)? sPT : 0; } // 2.1a.2.3.3 Electrons else if (pid == pE) { // Finds in which calorimeter the particle has gone, to know its resolution D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo); if(currentCalo.getName() == dummyCalo.getName()) { cout << "** Warning: the calo coverage behind the tracker is not complete! **" << endl; } // final smeared EM energy // electromagnetic fraction F_em =1 for electrons; float sE = currentCalo.getElectromagneticResolution().Smear(particle->E); if (sE>0) { sET = sE/cosh(particle->Eta); // NB: ET is found via the calorimetry and not via the track curvature TLorentzVector p; p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE); if (sET > DET->PTCUT_elec) electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo)); } else { sET=0;} // if negative smeared energy -- needed for the tracks } // 2.1a.2.3.4 Other charged particles : smear them for the tracks! else { //other particles D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo); float sEem = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem()); float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad()); float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 ); sET = sE/cosh(particle->EtaCalo); } // 2.1a.2.3.5 Tracks if( (rand()%100) < DET->TRACK_eff && sET!=0) { elementTrack = (D_Track*) branchTrack->NewEntry(); elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET); TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex! // TODO!!! associates the tracks to the calo where it points to // TODO!!! apply a smearing on the position of the origin of the track // TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout); } } // 2.1a.2.3 : if tracker/energy-flow // 2.1a.2.4 Photons // stays in the tracker -> track available -> gamma ID else if( (pid == pGAMMA) && fabs(particle->EtaCalo)< DET->CEN_max_tracker ) { float iEta=UNDEFINED, iPhi=UNDEFINED; DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta D_CaloTower tower(iEta,iPhi); // stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search list_of_towers_with_photon.addTower(tower); } // 2.1a.2.5 : very forward detectors else if (DET->FLAG_vfd==1) { // for the moment, only protons are transported // BUT !!! could be a beam of other particles! (heavy ions?) // BUT ALSO !!! if very forward muons, or others! VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle); VFD->ZDC(treeWriter,branchZDC,particle); } // 2.1a.2.6: Zero degree calorimeter //else if(DET->FLAG_zdc==1) { //VFD->ZDC(treeWriter,branchZDC,particle); //} } // 2.1a.2 : if visible particle } // loop on all particles 2.1a // 2.1b loop on all (activated) towers // at this stage, list_of_active_towers may contain several times the same tower // first step is to merge identical towers, by matching their (iEta,iPhi) list_of_active_towers.mergeDuplicates(); // Calotower smearing list_of_active_towers.smearTowers(list_of_calorimeters); for(unsigned int i=0; iNewEntry(); elementCalo->set(list_of_active_towers[i]); // not beautiful : should be improved! TLorentzVector p; p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e ); PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E())); towers.push_back(Tower); } } // loop on towers // 2.1c photon ID // list_of_towers_with_photon is the list of towers with photon candidates // already smeared ! // sorts the vector and smears duplicates list_of_towers_with_photon.mergeDuplicates(); for(unsigned int i=0; i 0) { TLorentzVector p; p.SetPtEtaPhiE(cal.getET(), eta,phi,cal.getE() ); if (cal.getET() > DET->PTCUT_gamma) { gamma.push_back(D_Particle(p,pGAMMA,p.Eta(),p.Phi())); } } } // for -- list of photons } // IF NEW ALGORITHM with energy flow // 2.2 ********** Output preparation & complex objects *********** // 2.2.1 ********************* sorting collections by decreasing pt DET->SortedVector(electron); for(unsigned int i=0; i < electron.size(); i++) { elementElec = (TRootElectron*) branchElectron->NewEntry(); elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E()); elementElec->EtaCalo = electron[i].EtaCalo(); elementElec->PhiCalo = electron[i].PhiCalo(); elementElec->Charge = sign(electron[i].PID()); elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks } DET->SortedVector(muon); for(unsigned int i=0; i < muon.size(); i++) { elementMu = (TRootMuon*) branchMuon->NewEntry(); elementMu->Charge = sign(muon[i].PID()); elementMu->Set(muon[i].Px(),muon[i].Py(),muon[i].Pz(),muon[i].E()); elementMu->EtaCalo = muon[i].EtaCalo(); elementMu->PhiCalo = muon[i].PhiCalo(); elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0); } DET->SortedVector(gamma); for(unsigned int i=0; i < gamma.size(); i++) { elementPhoton = (TRootPhoton*) branchPhoton->NewEntry(); elementPhoton->Set(gamma[i].Px(),gamma[i].Py(),gamma[i].Pz(),gamma[i].E()); elementPhoton->EtaCalo = gamma[i].EtaCalo(); elementPhoton->PhiCalo = gamma[i].PhiCalo(); } // 2.2.2 ************* computes the Missing Transverse Momentum TLorentzVector Att(0.,0.,0.,0.); for(unsigned int i=0; i < towers.size(); i++) { Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E); if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd) { PTmis = PTmis + Att; // create a fastjet::PseudoJet with these components and put it onto // back of the input_particles vector input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E)); } } elementEtmis = (TRootETmis*) branchETmis->NewEntry(); elementEtmis->ET = (PTmis).Pt(); elementEtmis->Phi = (-PTmis).Phi(); elementEtmis->Px = (-PTmis).Px(); elementEtmis->Py = (-PTmis).Py(); // 2.2.3 ************* jets, B-tag, tau jets sorted_jets=JETRUN->RunJets(input_particles); JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ); JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral); treeWriter->Fill(); } // 2. Loop over all events ('for' loop) cout <<"** **"<< endl; cout <<"** Exiting detector simulation... **"<< endl; treeWriter->Write(); delete treeWriter; loopwatch.Stop(); // 3. ********** Trigger & Frog *********** // 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard triggerwatch.Start(); if(DET->FLAG_trigger == 1) { cout <<"** **"<UseBranch("Electron"); TClonesArray *branchMuonTrig = treeReaderT->UseBranch("Muon"); TClonesArray *branchJetTrig = treeReaderT->UseBranch("Jet"); TClonesArray *branchTauJetTrig = treeReaderT->UseBranch("TauJet"); TClonesArray *branchPhotonTrig = treeReaderT->UseBranch("Photon"); TClonesArray *branchETmisTrig = treeReaderT->UseBranch("ETmis"); ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger"); ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class()); Long64_t entryT, allEntriesT = treeReaderT->GetEntries(); // loop on all entries for(entryT = 0; entryT < allEntriesT; ++entryT) { treeWriterT->Clear(); treeReaderT->ReadEntry(entryT); TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger); treeWriterT->Fill(); } // loop on all entries cout <<"** Exiting trigger simulation... **"<< endl; treeWriterT->Write(); delete treeWriterT; delete treeReaderT; } // trigger triggerwatch.Stop(); // 3.2 ************** FROG display frogwatch.Start(); if(DET->FLAG_frog == 1) { cout <<"** **"<BuildEvents(outputfilename); FROG->BuildGeom(); delete FROG; } frogwatch.Stop(); // 4. ********** End & Exit *********** globalwatch.Stop(); cout <<"** **"<< endl; cout <<"** ################## Time report ################# **"<< endl; cout << left << setw(32) <<"** Time report for "<<"" << left << setw(15) << allEntries <<"" << right << setw(22) <<"events **"<FLAG_trigger == 1) { cout << left << setw(10) <<"**"<<"" << left << setw(15) <<" + Trigger:"<<"" << right << setw(15) <FLAG_frog == 1) { cout << left << setw(10) <<"**"<<"" << left << setw(15) <<" + Frog:"<<"" << right << setw(15) <