Changeset 310 in svn for trunk/Delphes.cpp

Timestamp:

Mar 10, 2009, 11:52:16 AM (16 years ago)

Author:

severine ovyn

Message:

add Ehad/Eem

File:

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

trunk/Delphes.cpp

@@ -23 +23 @@
 **         Center for Particle Physics and Phenomenology (CP3)        ** 
 **                Universite catholique de Louvain (UCL)              **       
-**                     Louvain-la-Neuve, Belgium                      **            
+**                     Louvain-la-Neuve, Belgium                      **       
 **                                                                    **
 **                      Copyright (C) 2008-2009,                      **
@@ -326 +326 @@
   vector<int> NTrackJet;
 
-  bool FLAG_lhco = true;
-  
   TSimpleArray<TRootGenParticle> NFCentralQ;
   
@@ -380 +378 @@
       input_particles.clear();
       NTrackJet.clear();     
- 
-      {
-        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();
-        GenParticle *particleG; 
-        while( (particleG = (GenParticle*) itGen.Next()) ) {
-          
+      
+      
+      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();
+      GenParticle *particleG; 
+      while( (particleG = (GenParticle*) itGen.Next()) ) 
+        {
           TRootGenParticle *particle = new TRootGenParticle(particleG);
           int pid = abs(particle->PID);
@@ -402 +401 @@
               fabs(particle->Eta) < DET->CEN_max_tracker &&
               particle->Status != 1 &&
-              particle->PT > DET->PT_QUARKS_MIN ) {
-            NFCentralQ.Add(particle);
-          }
+              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) ){
+          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
+              // 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->Charge==0 ||
+                  (!DET->FLAG_bfield && particle->Charge!=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.1 Central solenoidal magnetic field 
-            TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo
-            // 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->Charge==0 ||
-                (!DET->FLAG_bfield && particle->Charge!=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->Charge !=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) 
+                    {
+                      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 = (TRootTracks*) branchTrack->NewEntry();
+                      elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
+                      TrackCentral.push_back(*elementTrack); // tracks at vertex!
+                      // 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_RP==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);
+                    }
+                  // 2.1a.2.6: Zero degree calorimeter
+                  if(DET->FLAG_vfd==1) 
+                    {
+                      VFD->ZDC(treeWriter,branchZDC,particle);
+                    }
+                }
             
-            // 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->Charge !=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) {
-                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 = (TRootTracks*) branchTrack->NewEntry();
-                elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
-                //TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex!
-                TrackCentral.push_back(*elementTrack); // tracks at vertex!
-                // 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_RP==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);
-                }
-                // 2.1a.2.6: Zero degree calorimeter
-                if(DET->FLAG_vfd==1) {
-                  VFD->ZDC(treeWriter,branchZDC,particle);
-                }
-              }
-            
-          } // 2.1a.2 : if visible particle
+            } // 2.1a.2 : if visible particle
           delete 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++) {
+        }// 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.sortElements();
+      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);
-          }
+          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++) {
+      
+          // 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())); }
-          }
+          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); 
-      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,DET->ISOL_PT,DET->ISOL_Cone);//isolation based on tracks
-      }                                                                                 /////////////// HARDCODING
+      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,DET->ISOL_PT,DET->ISOL_Cone);//isolation based on tracks
+          
+          D_CaloTower calElec(list_of_active_towers.getElement(electron[i].EtaCalo(),electron[i].PhiCalo()));
+          elementElec->EHoverEE = calElec.getEhad()/calElec.getEem();
+        }                                                                                       /////////////// HARDCODING
+      
       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,DET->ISOL_PT,DET->ISOL_Cone); /////////////// HARDCODING
-      }
+      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,DET->ISOL_PT,DET->ISOL_Cone); /////////////// HARDCODING
+          D_CaloTower calMuon(list_of_active_towers.getElement(muon[i].EtaCalo(),muon[i].PhiCalo()));
+          elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
+        }
+
       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());
-      }
+      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());
+          D_CaloTower calGamma(list_of_active_towers.getElement(gamma[i].EtaCalo(),gamma[i].PhiCalo()));
+          elementPhoton->EHoverEE = calGamma.getEhad()/calGamma.getEem();
+        }
       
       // 2.2.2 ************* computes the Missing Transverse Momentum
@@ -617 +643 @@
       
       // 2.2.3 ************* jets, B-tag, tau jets
-      vector<int> NTrackJet;
-      sorted_jets=JETRUN->RunJets(input_particles, TrackCentral,NTrackJet);
-      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ,NTrackJet);
-      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral,NTrackJet);
+      vector<int> NTrackJet; //for number of tracks
+      vector<float> EHADEEM; //for energyHad over energyEm
+      sorted_jets=JETRUN->RunJets(input_particles, TrackCentral,NTrackJet,EHADEEM,list_of_active_towers);
+      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ,NTrackJet,EHADEEM);
+      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral,NTrackJet,EHADEEM);
       
       treeWriter->Fill();
     } // 2. Loop over all events ('for' loop)
   
-  
-  cout <<"**                                                                 **"<< endl;
   cout <<"**                 Exiting detector simulation...                  **"<< endl;