Changeset 264 in svn for trunk/Delphes.cpp

Timestamp:

Feb 11, 2009, 10:22:30 AM (16 years ago)

Author:

Xavier Rouby

Message:

first test 2.0

File:

• trunk/Delphes.cpp (modified) (22 diffs)

trunk/Delphes.cpp

@@ -31 +31 @@
 
 /// \file Delphes.cpp
-/// \brief executable for the Delphes
+/// \brief Executable for Delphes
 
 #include "TChain.h"
@@ -41 +41 @@
 #include "ExRootTreeWriter.h"
 #include "ExRootTreeBranch.h"
+#include "ExRootProgressBar.h"
 
 #include "DataConverter.h"
@@ -48 +49 @@
 
 #include "SmearUtil.h"
+#include "CaloUtil.h"
 #include "BFieldProp.h"
 #include "TriggerUtil.h"
@@ -56 +58 @@
 #include <vector>
 #include <iostream>
-
-#include "Utilities/ExRootAnalysis/interface/ExRootProgressBar.h"
 
 using namespace std;
@@ -77 +77 @@
 int main(int argc, char *argv[])
 {
+unsigned int nnnn=0, mmmm=0;
+
   int appargc = 2;
   char *appName= new char[20];
@@ -86 +88 @@
   delete [] appName;
   delete [] appOpt;
-
 
   if(argc != 3 && argc != 4 && argc != 5) {
@@ -242 +243 @@
   //Propagation of tracks in the B field
   TrackPropagation *TRACP = new TrackPropagation(DET); 
-  //TrackPropagation *TRACP = new TrackPropagation(DetDatacard); 
   
   //Jet information
   JetsUtil *JETRUN = new JetsUtil(DET);
-  //JetsUtil *JETRUN = new JetsUtil(DetDatacard);
   
   //VFD information
   VeryForward * VFD = new VeryForward(DET);
-  //VeryForward * VFD = new VeryForward(DetDatacard);
  
   // data converters
@@ -299 +297 @@
   ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
   ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
-  ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::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());
@@ -306 +305 @@
   ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
   
-  TRootGenParticle *particle; 
   TRootETmis *elementEtmis;
   TRootElectron *elementElec;
   TRootMuon *elementMu;
   TRootPhoton *elementPhoton;
-  TRootTracks *elementTracks;
+  D_Track * elementTrack;
   TRootCalo *elementCalo;
   
@@ -323 +321 @@
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-  vector<ParticleUtil> electron;
-  vector<ParticleUtil> muon;
-  vector<ParticleUtil> gamma;
+  //vector<ParticleUtil> electron;
+  vector<D_Particle> electron;
+  //vector<ParticleUtil> muon;
+  vector<D_Particle> muon;
+  //vector<ParticleUtil> gamma;
+  vector<D_Particle> gamma;
 
   TSimpleArray<TRootGenParticle> NFCentralQ;
-  float iPhi=0,iEta=0;
-
+
+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();
 
 
@@ -360 +379 @@
       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;
 
         
@@ -381 +404 @@
           // the ordering of conditions have been optimised for speed : put first the STATUS condition
           if( (particle->Status == 1)    &&
-              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+              (pid != pNU1) && (pid != pNU2) && (pid != pNU3) &&
               (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
               ) 
@@ -431 +454 @@
               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 */ }
+                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 != -100 && iPhi != -100) { 
+                  if(iEta != UNDEFINED && iPhi != UNDEFINED) { 
                       momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
                       elementCalo = (TRootCalo*) branchCalo->NewEntry();
@@ -441 +464 @@
                       PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
                       towers.push_back(Tower);
-                  } // if iEta != -100
+                  } // if iEta != UNDEFINED
                 } // else : when particles reach the calos
               } // 2.1.2.3 calotowers
-              
+
 
               // 2.1.2.4 ********************* central detector: tracks
@@ -451 +474 @@
                  (genMomentumBfield.E()!=0) &&
                  (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) && 
-                 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&     
+                 (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)  &&
@@ -475 +498 @@
              VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
           }
-          
-       } // 2.1 while : loop on all particles of the event.
-     
+} // 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; i<list_of_active_towers.size(); i++) {
+        float iEta = list_of_active_towers[i].getEta();
+        float iPhi = list_of_active_towers[i].getPhi();
+        float e = list_of_active_towers[i].getE(); 
+        if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) {
+           elementCalo = (TRootCalo*) branchCalo->NewEntry();
+           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<list_of_towers_with_photon.size(); i++) {
+          float eta = list_of_towers_with_photon[i].getEta(); 
+          float phi = list_of_towers_with_photon[i].getPhi();
+          D_CaloTower cal(list_of_active_towers.getElement(eta,phi));
+          if(cal.getEta() != UNDEFINED && cal.getPhi() != UNDEFINED &&  cal.getE() > 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); 
@@ -487 +698 @@
         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
@@ -495 +708 @@
         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);
       }
@@ -501 +716 @@
         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();
       }
       
@@ -522 +739 @@
       elementEtmis->Py = (-PTmis).Py();
       
-      // 2.2.3 ************* B-tag, tau jets
+      // 2.2.3 ************* jets, B-tag, tau jets
       sorted_jets=JETRUN->RunJets(input_particles);
       JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
@@ -529 +746 @@
       treeWriter->Fill();
     } // 2. Loop over all events ('for' loop)
- 
+
+
   cout <<"**                                                                 **"<< endl;
   cout <<"**                 Exiting detector simulation...                  **"<< endl;
@@ -649 +867 @@
   delete treeReader;
   delete DET;
-//  delete TRIGT;
+  delete TRIGT;
   delete TRACP;
   delete JETRUN;