Changeset 392 in svn

Timestamp:

May 18, 2009, 6:52:13 PM (15 years ago)

Author:

Xavier Rouby

Message:

new CaloIsolation for muons : EtRatio up-to-date

Location:

trunk

Files:

• interface/SmearUtil.h (modified) (1 diff)
• src/SmearUtil.cc (modified) (3 diffs)

trunk/interface/SmearUtil.h

@@ -260 +260 @@
 
   /// Lepton isolation based on calorimetry (optional. Default: off)
-  float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers);
+  float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta);
 
   //********************* returns a segmented value for eta and phi, for calo towers *****

trunk/src/SmearUtil.cc

@@ -121 +121 @@
   ISOL_Cone       = 0.5;      //Cone  for isolation criteria
   ISOL_Calo_ET    = 1E99;     //minimal tower energy for isolation criteria. Default off = 1E99
-  ISOL_Calo_Grid  = 3;        //Grid size (N x N) for calorimetric isolation
+  ISOL_Calo_Grid  = 3;        //Grid size (N x N) for calorimetric isolation -- should be odd
 
   // General jet variable
@@ -592 +592 @@
     else if(strstr(temp_string.c_str(),"PdgTableFilename"))     {curstring >> varname >> svalue; PdgTableFilename = svalue;}
   }
-  
+ 
+  if(ISOL_Calo_Grid%2 ==0) {
+    ISOL_Calo_Grid++;
+    cout <<"**        WARNING: ISOL_Calo_Grid is not odd. Set it to  "<< ISOL_Calo_Grid << "        **" << endl;
+  }
+ 
   //jet stuffs not defined in the input datacard
   JET_overlap     = 0.75;
@@ -1226 +1231 @@
 
 // ******* Calorimetric isolation
-float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers) {
+float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta) {
  // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy 
  // in a 3×3 grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
    if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
-   float etrat=0;
+   float et_sum=0;
    // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
-/*   for(unsigned int i=0; i < towers.size(); i++) {
-     if(towers[i].E > ISOL_Calo_ET) {
-        float genDeltaR = DeltaR(part.Phi(),part.Eta(),towers[i].getPhi(),towers[i].getEta());
-        if(genDeltaR < ISOL_Calo_Cone) {
-           ptiso += towers[i].getET();
-        }           
-     }
-   } // loop on towers
-   ptiso -=
-*/
-  etrat = 100*etrat/part.Pt();
+   // Get the EtaCalo/PhiCalo of the muon ; 
+   // transform it into iEta/iPhi to get the towers, and their neighbourh (i-1, i-2, etc)
+
+   unsigned int N = ISOL_Calo_Grid;
+   int index= iUNDEFINED; // index of the central tower of the grid in TOWER_eta_edges[.]; 
+                          // !! TOWER_eta_edges is only with eta>0
+   // finds the index of the central tower of the NxN grid
+   for (unsigned int i=1; i< TOWER_number+1; i++) {
+       if(fabs(iEta) >= TOWER_eta_edges[i-1] && fabs(iEta) < TOWER_eta_edges[i]) {
+           index = i-1;
+           break;
+       }
+   }
+   if(index != iUNDEFINED) {
+           // finds the size in phi of the cells for this eta
+           float dphi = TOWER_dphi[index]*pi/180.; // in rad
+
+           //cout << "Grid " << " ----------\n";
+           for (unsigned int i_eta=0; i_eta<N; i_eta++) {
+              unsigned int real_index = (iEta>0) ? index+i_eta-(N-1)/2 : index+1+i_eta-(N-1)/2 ;
+              float eta_ith_tower = TOWER_eta_edges[real_index];
+              if(iEta<0) eta_ith_tower *= -1;
+                                        
+              for (unsigned int i_phi=0; i_phi<N; i_phi++) {
+                 float phi_ith_tower = iPhi + (float)(i_phi - (N-1.)/2.)*dphi;
+                 D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
+                 if(calMuon.getEta() != UNDEFINED && calMuon.getE() > ISOL_Calo_ET) { 
+                    et_sum += calMuon.getE();
+                    //cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getE() << " GeV" << endl;
+                 } 
+                 //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
+              }
+    
+           } // NxN grid
+        //cout << "-----------" << etrat << endl;
+   } 
+   else if (CEN_max_mu < CEN_max_calo_fwd) 
+        cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl; 
+        // should never happen ! this would be a bug
+  //cout << "etrat = " << et_sum << "\t Pt=" << part.Pt() << endl;
+  
+  // should return a number between 0 and 99 (due to LHCO definitions)
+  // which is Pt(muon) / sum(ET)
+  float etrat = 0.;
+  if(et_sum==0) etrat = 99.;
+  else if(et_sum>0) etrat = 100*part.Pt()/et_sum;
   if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
-  else if(etrat>99) cout << "Error: etrat shoud be in [0;99] in CaloIsolation (" << etrat <<")\n";
+  //else if(etrat>99) cout << "Error: etrat should be in [0;99] in CaloIsolation (" << etrat <<")\n";
+  if(etrat>99) etrat = 99;
   return etrat; 
 }