Changeset 550 in svn for trunk/src/SmearUtil.cc

Timestamp:

Feb 22, 2010, 7:38:15 PM (15 years ago)

Author:

severine ovyn

Message:

changements pour EtRatio muons et electrons. Bugge

File:

• trunk/src/SmearUtil.cc (modified) (5 diffs)

trunk/src/SmearUtil.cc

@@ -1440 +1440 @@
 
 // ******* Calorimetric isolation
-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.
+float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
+  // iPhi and iEta are the coordinates of the center of the tower crossed by the particle
+  // iPhi is in [-pi ; pi]    and   iEta is in [-etamax ; etamax]
+  // iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range
+
+  // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy 
+  // in a N×N grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
+  // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
+
    if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
    float et_sum=0;
-   // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
-   // 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(N%2 ==0 || N==0) { cout << "Warning: ISOL_Calo_Grid must be odd. ISOL_Calo_Grid = 3 will be used\n"; N=3;}
+   int index= iTower; // index of the central tower of the grid in TOWER_eta_edges[.]; 
+                      // !! TOWER_eta_edges is only with eta>0
+
+cout << "iEta= " << iEta << "\tiPhi= " << iPhi << endl;
+   //cout << "la particule est dans la " << index << "eme cellule, entre " << TOWER_eta_edges[index] << " et " << TOWER_eta_edges[index+1] << endl;
+   //float centre = (iEta>0) ? ((TOWER_eta_edges[index]+TOWER_eta_edges[index+1])/2.0) : -((TOWER_eta_edges[index]+TOWER_eta_edges[index+1])/2.0);
+   //cout << "le centre de cette tour est " << centre << endl; 
+
+   if(index != iUNDEFINED) {
+      int range = (int) (N-1)/2;           // the (int) conversion is needed, as N is "unsigned int"
+      for(int i= -range; i<= range; i++) { // shift of index in eta: e.g.: i = -1, 0, 1 if N=3;
+      // goal : to identify the center of each cell in the NxN grid.
+      // the eta/phi coordinates of each center will be in (eta_ith_tower,phi_ith_tower)
+
+         unsigned int i_eta = (index+i>=0) ? index + i : (int) -1*(index + i +1); // i_eta is the index in [0; index_eta_max]
+         //cout << i_eta << "\t";
+         float dphi = TOWER_dphi[i_eta]*pi/180.; // in rad // size in phi of the cells for this eta
+         float eta_ith_tower = (TOWER_eta_edges[i_eta]+TOWER_eta_edges[i_eta+1])/2.0;
+         if(iEta<0) eta_ith_tower *= -1;     // needed if the central tower is in negative region
+         if(index+i<0) eta_ith_tower *= -1;  // needed if the central part is crossed by the grid
+         //cout << "pour eta_ith_tower=" << eta_ith_tower << ", on va dphi = " << dphi << endl;
+         // !!! at this point, eta_ith_tower is the value in eta of the center of the current calo cell
+
+         for(int j= -range; j<= range; j++) { // shift of index in phi. e.g.: j = -1, 0, 1 if N=3;
+            float phi_ith_tower = iPhi + j*dphi; // in radian
+            //cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower << "\tj= " << j << "\tdphi=" << dphi << endl;
+            if(phi_ith_tower > pi) phi_ith_tower -= 2.*pi;
+            else if (phi_ith_tower < -pi) phi_ith_tower += 2.*pi;
+            
+            D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
+            cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower <<"\t"
+                 << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";
+
+            if(calMuon.getEta() != UNDEFINED )
+             if(calMuon.getET() > ISOL_Calo_ET) {
+              et_sum += calMuon.getET();
+              //cout << "eta/phi/Et = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getET() << " GeV" << endl;
+              cout << calMuon.getET() << " GeV";
+            }
+            //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
+            else cout << phi_ith_tower << "\tnot active";
+            cout << endl;
+         }
+      }
+   }  // undefined index
+   else {
+     cout << "out of range" << endl;
+     if (CEN_max_mu < CEN_max_calo_fwd) 
+        cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl; 
    }
-   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;
@@ -1490 +1510 @@
   // 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(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 should be in [0;99] in CaloIsolation (" << etrat <<")\n";
   if(etrat>99) etrat = 99;
+*/
+  etrat = et_sum/part.Pt();
   return etrat; 
 }
@@ -1500 +1522 @@
 
   //********** returns a segmented value for eta and phi, for calo towers *****
-void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
+int RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
    iEta = UNDEFINED;
    int index= iUNDEFINED;
@@ -1510 +1532 @@
         }
    }
-   if(index==UNDEFINED) return; 
    iPhi = UNDEFINED;
+   if(index==iUNDEFINED) return index; 
+
    float dphi = TOWER_dphi[index]*pi/180.; 
    for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
@@ -1522 +1545 @@
    }
    if (phi > pi-dphi) iPhi = pi-dphi;
+   return index;
 }