Changeset 384 in svn for trunk/Delphes.cpp

Timestamp:

May 15, 2009, 7:28:55 PM (16 years ago)

Author:

Xavier Rouby

Message:

energy flow

File:

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

trunk/Delphes.cpp

@@ -99 +99 @@
   }
   
-  
-  
-  
   cout << endl << endl;
   
@@ -180 +177 @@
   string line;
   ifstream infile(inputFileList.c_str());
+  if(!infile.good()) {
+    cout << "**    ERROR: Input list (" << left << setw(13) << inputFileList << ") not found. Exiting...      **"<< endl; 
+    cout <<"*********************************************************************"<< endl;
+    exit(1);
+  }
   infile >> line; // the first line determines the type of input files
   
@@ -390 +392 @@
       NTrackJet.clear();     
       
-      
+      // 'list_of_active_towers' contains the exact list of calorimetric towers which have some deposits inside (E>0).
+      // The towers of this list will be smeared according to the calo resolution, afterwards
       D_CaloTowerList list_of_active_towers; 
+ 
+      // 'list_of_towers_with_photon' and 'list_of_centowers_with_neutrals' are list of towers, whose energy is **not** computed.
+      // They are only used to store the eta/phi of some towers, in order to search later inside 'list_of_active_towers'.
+      // 'list_of_towers_with_photon' contains the towers hit by photons only
+      // 'list_of_centowers_with_neutrals' is used to the jet-E-flow calculation: contains the towers with eta < CEN_max_tracker,
+      //     i.e. towers behind the tracker.
       D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates
-       
-      
+
+      D_CaloTowerList list_of_centowers_with_neutrals;                // list of towers with neutral particles : for jet E-flow
+      float etamax_calocoverage_behindtracker = DET->CEN_max_tracker; // finds the extension in eta of the furthest 
+      for (unsigned int i=1; i< DET->TOWER_number+1; i++) {           // cell (at least) partially behind the tracker
+        if(DET->TOWER_eta_edges[i] > DET->CEN_max_tracker) break;
+        etamax_calocoverage_behindtracker = DET->TOWER_eta_edges[i];
+      }
       // 2.1a Loop over all particles in event, to fill the towers
       itGen.Reset();
@@ -472 +486 @@
                         {
                           p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta));
-                          muon.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
+                          muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo));
                         }
                       sET = (sPT >0)? sPT : 0;
@@ -498 +512 @@
                           if (sET > DET->PTCUT_elec)
                             electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
+                          //if(DET->JET_Eflow) input_particles.push_back(fastjet::PseudoJet(p.Px(),p.Py(),p.Pz(),p.E()));
                         } 
                       else { sET=0;} // if negative smeared energy -- needed for the tracks
@@ -517 +532 @@
                       elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
                       TrackCentral.push_back(*elementTrack); // tracks at vertex!
+                      if(DET->JET_Eflow) 
+                         input_particles.push_back(fastjet::PseudoJet(particle->Px,particle->Py,particle->Pz,particle->E));
                       // TODO!!! apply a smearing on the position of the origin of the track
                       // TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout);
@@ -531 +548 @@
                   list_of_towers_with_photon.addTower(tower); 
                 }
-              // 2.1a.2.5 : very forward detectors
+              // 2.1a.2.5 Neutrals within tracker -- for jet energy flow
+              else if(  particle->Charge ==0 && fabs(particle->EtaCalo)< etamax_calocoverage_behindtracker) 
+                {
+                  float iEta=UNDEFINED, iPhi=UNDEFINED;
+                  DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
+                  D_CaloTower tower(iEta,iPhi); 
+                  list_of_centowers_with_neutrals.addTower(tower);
+                }
+              // 2.1a.2.6 : very forward detectors
               else
                 {
@@ -551 +576 @@
           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
@@ -580 +604 @@
         } // 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 
+      // 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));
+          D_CaloTower cal(list_of_active_towers.getElement(eta,phi)); //// <---------- buh???????
           if(cal.getEta() != UNDEFINED && cal.getPhi() != UNDEFINED &&  cal.getE() > 0) 
             {
@@ -596 +620 @@
             }
         } // for -- list of photons
-      
+
+      // 2.1d jet-E-flow -- taking into account the neutrals within tracker
+      if(DET->JET_Eflow) {
+        list_of_centowers_with_neutrals.mergeDuplicates();
+        for(unsigned int i=0; i<list_of_centowers_with_neutrals.size(); i++) {
+             float eta = list_of_centowers_with_neutrals[i].getEta();
+             float phi = list_of_centowers_with_neutrals[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() );
+                 //cout << "**************list: " << p.Px() << " " << p.Py() << " " << p.Pz() << " " << p.E() << endl;
+                 input_particles.push_back(fastjet::PseudoJet(p.Px(),p.Py(),p.Pz(),p.E()));
+               }
+         } // for - list_of_centowers
+      } // JET_Eflow
+
       // 2.2 ********** Output preparation & complex objects  ***********
       // 2.2.1 ********************* sorting collections by decreasing pt
@@ -652 +693 @@
               // 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));
+              if(!DET->JET_Eflow) 
+                input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
+              else { if(fabs(Att.Eta()) > DET->CEN_max_tracker)
+                input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
+              }
             }
         }