Fork me on GitHub

Changeset 550 in svn


Ignore:
Timestamp:
Feb 22, 2010, 7:38:15 PM (15 years ago)
Author:
severine ovyn
Message:

changements pour EtRatio muons et electrons. Bugge

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/Delphes.cpp

    r547 r550  
    628628          elementElec->IsolFlag = DET->Isolation(electron[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoEl);
    629629          elementElec->IsolPt = ptisoEl;
    630           DET->BinEtaPhi(elementElec->PhiCalo,elementElec->EtaCalo,iPhiEl,iEtaEl);
     630
     631cout << "\nC'est un electron\n";
     632          int electron_tower_index = DET->BinEtaPhi(elementElec->PhiCalo,elementElec->EtaCalo,iPhiEl,iEtaEl);
    631633          D_CaloTower calElec(list_of_active_towers.getElement(iEtaEl,iPhiEl));
    632634          elementElec->EHoverEE = calElec.getEhad()/calElec.getEem();
     635          elementElec->EtRatio =  DET->CaloIsolation(muon[i], list_of_active_towers,iPhiEl,iEtaEl,electron_tower_index);
    633636        }                       
    634637     
     
    644647          elementMu->IsolFlag = DET->Isolation(muon[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoMu);
    645648          elementMu->IsolPt = ptisoMu;
    646           DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
     649
     650cout << "\nC'est un muon\n";
     651          // binning of the muon position in calo
     652          int muon_tower_index = DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
    647653          D_CaloTower calMuon(list_of_active_towers.getElement(iEtaMu,iPhiMu));
    648654          if( calMuon.getEem() !=0 ) elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
    649655          else elementMu->EHoverEE = UNDEFINED;
    650           elementMu->EtRatio = DET->CaloIsolation(muon[i], list_of_active_towers,iPhiMu,iEtaMu);
     656          elementMu->EtRatio = DET->CaloIsolation(muon[i], list_of_active_towers,iPhiMu,iEtaMu,muon_tower_index);
    651657        }
    652658
  • trunk/Utilities/ExRootAnalysis/interface/BlockClasses.h

    r547 r550  
    262262  float PhiCalo; // particle azimuthal angle in rad when entering the calo
    263263  float EHoverEE;
     264  float EtRatio;
    264265
    265266  void SetEtaPhiCalo(const float eta, const float phi) {EtaCalo=eta; PhiCalo=phi;};
  • trunk/Utilities/ExRootAnalysis/src/BlockClasses.cc

    r539 r550  
    2929**                                                                    **
    3030***********************************************************************/
     31 
    3132
    3233#include "BlockClasses.h"
  • trunk/interface/SmearUtil.h

    r531 r550  
    278278
    279279  /// Lepton isolation based on calorimetry (optional. Default: off)
    280   float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta);
     280  float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower);
    281281
    282282  //********************* returns a segmented value for eta and phi, for calo towers *****
    283   void BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta);
     283  int BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta);
    284284
    285285  void setNames(const string& list, const string& det, const string& trig);
  • trunk/src/CaloUtil.cc

    r443 r550  
    237237                return _calotowers[i];
    238238  }
     239cout << "Xav!!!!! not found in getElement" << endl;
    239240  return dummyTower; // if not found
    240241}
     
    247248                return _calotowers[i];
    248249  }
     250cout << "Xav!!!!! not found in getElement const" << endl;
    249251  return dummyTower; // if not found
    250252}
  • trunk/src/SmearUtil.cc

    r537 r550  
    14401440
    14411441// ******* Calorimetric isolation
    1442 float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta) {
    1443  // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy
    1444  // 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.
     1442float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
     1443  // iPhi and iEta are the coordinates of the center of the tower crossed by the particle
     1444  // iPhi is in [-pi ; pi]    and   iEta is in [-etamax ; etamax]
     1445  // iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range
     1446
     1447  // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy
     1448  // 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.
     1449  // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
     1450
    14451451   if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
    14461452   float et_sum=0;
    1447    // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
    1448    // Get the EtaCalo/PhiCalo of the muon ;
    1449    // transform it into iEta/iPhi to get the towers, and their neighbourh (i-1, i-2, etc)
    14501453
    14511454   unsigned int N = ISOL_Calo_Grid;
    1452    int index= iUNDEFINED; // index of the central tower of the grid in TOWER_eta_edges[.];
    1453                           // !! TOWER_eta_edges is only with eta>0
    1454    // finds the index of the central tower of the NxN grid
    1455    for (unsigned int i=1; i< TOWER_number+1; i++) {
    1456        if(fabs(iEta) >= TOWER_eta_edges[i-1] && fabs(iEta) < TOWER_eta_edges[i]) {
    1457            index = i-1;
    1458            break;
    1459        }
     1455   if(N%2 ==0 || N==0) { cout << "Warning: ISOL_Calo_Grid must be odd. ISOL_Calo_Grid = 3 will be used\n"; N=3;}
     1456   int index= iTower; // index of the central tower of the grid in TOWER_eta_edges[.];
     1457                      // !! TOWER_eta_edges is only with eta>0
     1458
     1459cout << "iEta= " << iEta << "\tiPhi= " << iPhi << endl;
     1460   //cout << "la particule est dans la " << index << "eme cellule, entre " << TOWER_eta_edges[index] << " et " << TOWER_eta_edges[index+1] << endl;
     1461   //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);
     1462   //cout << "le centre de cette tour est " << centre << endl;
     1463
     1464   if(index != iUNDEFINED) {
     1465      int range = (int) (N-1)/2;           // the (int) conversion is needed, as N is "unsigned int"
     1466      for(int i= -range; i<= range; i++) { // shift of index in eta: e.g.: i = -1, 0, 1 if N=3;
     1467      // goal : to identify the center of each cell in the NxN grid.
     1468      // the eta/phi coordinates of each center will be in (eta_ith_tower,phi_ith_tower)
     1469
     1470         unsigned int i_eta = (index+i>=0) ? index + i : (int) -1*(index + i +1); // i_eta is the index in [0; index_eta_max]
     1471         //cout << i_eta << "\t";
     1472         float dphi = TOWER_dphi[i_eta]*pi/180.; // in rad // size in phi of the cells for this eta
     1473         float eta_ith_tower = (TOWER_eta_edges[i_eta]+TOWER_eta_edges[i_eta+1])/2.0;
     1474         if(iEta<0) eta_ith_tower *= -1;     // needed if the central tower is in negative region
     1475         if(index+i<0) eta_ith_tower *= -1;  // needed if the central part is crossed by the grid
     1476         //cout << "pour eta_ith_tower=" << eta_ith_tower << ", on va dphi = " << dphi << endl;
     1477         // !!! at this point, eta_ith_tower is the value in eta of the center of the current calo cell
     1478
     1479         for(int j= -range; j<= range; j++) { // shift of index in phi. e.g.: j = -1, 0, 1 if N=3;
     1480            float phi_ith_tower = iPhi + j*dphi; // in radian
     1481            //cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower << "\tj= " << j << "\tdphi=" << dphi << endl;
     1482            if(phi_ith_tower > pi) phi_ith_tower -= 2.*pi;
     1483            else if (phi_ith_tower < -pi) phi_ith_tower += 2.*pi;
     1484           
     1485            D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
     1486            cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower <<"\t"
     1487                 << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";
     1488
     1489            if(calMuon.getEta() != UNDEFINED )
     1490             if(calMuon.getET() > ISOL_Calo_ET) {
     1491              et_sum += calMuon.getET();
     1492              //cout << "eta/phi/Et = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getET() << " GeV" << endl;
     1493              cout << calMuon.getET() << " GeV";
     1494            }
     1495            //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
     1496            else cout << phi_ith_tower << "\tnot active";
     1497            cout << endl;
     1498         }
     1499      }
     1500   }  // undefined index
     1501   else {
     1502     cout << "out of range" << endl;
     1503     if (CEN_max_mu < CEN_max_calo_fwd)
     1504        cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl;
    14601505   }
    1461    if(index != iUNDEFINED) {
    1462            // finds the size in phi of the cells for this eta
    1463            float dphi = TOWER_dphi[index]*pi/180.; // in rad
    1464 
    1465            //cout << "Grid " << " ----------\n";
    1466            for (unsigned int i_eta=0; i_eta<N; i_eta++) {
    1467               unsigned int real_index = (iEta>0) ? index+i_eta-(N-1)/2 : index+1+i_eta-(N-1)/2 ;
    1468               float eta_ith_tower = TOWER_eta_edges[real_index];
    1469               if(iEta<0) eta_ith_tower *= -1;
    1470                                        
    1471               for (unsigned int i_phi=0; i_phi<N; i_phi++) {
    1472                  float phi_ith_tower = iPhi + (float)(i_phi - (N-1.)/2.)*dphi;
    1473                  D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
    1474                  if(calMuon.getEta() != UNDEFINED && calMuon.getE() > ISOL_Calo_ET) {
    1475                     et_sum += calMuon.getE();
    1476                     //cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getE() << " GeV" << endl;
    1477                  }
    1478                  //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
    1479               }
    1480    
    1481            } // NxN grid
    1482         //cout << "-----------" << etrat << endl;
    1483    }
    1484    else if (CEN_max_mu < CEN_max_calo_fwd)
    1485         cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl;
    14861506        // should never happen ! this would be a bug
    14871507  //cout << "etrat = " << et_sum << "\t Pt=" << part.Pt() << endl;
     
    14901510  // which is Pt(muon) / sum(ET)
    14911511  float etrat = 0.;
    1492   if(et_sum==0) etrat = 99.;
    1493   else if(et_sum>0) etrat = 100*part.Pt()/et_sum;
     1512/*  if(et_sum==0) etrat = 99.;
     1513  //else if(et_sum>0) etrat = 100*part.Pt()/et_sum;
    14941514  if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
    14951515  //else if(etrat>99) cout << "Error: etrat should be in [0;99] in CaloIsolation (" << etrat <<")\n";
    14961516  if(etrat>99) etrat = 99;
     1517*/
     1518  etrat = et_sum/part.Pt();
    14971519  return etrat;
    14981520}
     
    15001522
    15011523  //********** returns a segmented value for eta and phi, for calo towers *****
    1502 void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
     1524int RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
    15031525   iEta = UNDEFINED;
    15041526   int index= iUNDEFINED;
     
    15101532        }
    15111533   }
    1512    if(index==UNDEFINED) return;
    15131534   iPhi = UNDEFINED;
     1535   if(index==iUNDEFINED) return index;
     1536
    15141537   float dphi = TOWER_dphi[index]*pi/180.;
    15151538   for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
     
    15221545   }
    15231546   if (phi > pi-dphi) iPhi = pi-dphi;
     1547   return index;
    15241548}
    15251549
Note: See TracChangeset for help on using the changeset viewer.