Fork me on GitHub

Changeset 321 in svn for trunk


Ignore:
Timestamp:
Mar 11, 2009, 12:53:50 PM (16 years ago)
Author:
Xavier Rouby
Message:

Isolation updated. ptiso implemented. etrat prepared but not finished

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/Delphes.cpp

    r319 r321  
    571571          // sorts the vector and smears duplicates
    572572      list_of_towers_with_photon.mergeDuplicates();
    573      
    574       for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++)
    575         {
     573      for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) {
    576574          float eta = list_of_towers_with_photon[i].getEta();
    577575          float phi = list_of_towers_with_photon[i].getPhi();
    … …  
    588586      // 2.2.1 ********************* sorting collections by decreasing pt
    589587      DET->SortedVector(electron);
    590       float iPhiEl=0,iEtaEl=0;
     588      float iPhiEl=0,iEtaEl=0,ptisoEl=0;
    591589      for(unsigned int i=0; i < electron.size(); i++)
    592590        {
    … …  
    596594          elementElec->PhiCalo = electron[i].PhiCalo();
    597595          elementElec->Charge = sign(electron[i].PID());
    598           elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,DET->ISOL_PT,DET->ISOL_Cone);     
     596          elementElec->IsolFlag = DET->Isolation(electron[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoEl);
     597          elementElec->IsolPt = ptisoEl;
    599598          DET->BinEtaPhi(elementElec->PhiCalo,elementElec->EtaCalo,iPhiEl,iEtaEl);
    600599          D_CaloTower calElec(list_of_active_towers.getElement(iEtaEl,iPhiEl));
    … …  
    603602     
    604603      DET->SortedVector(muon);
    605       float iPhiMu=0,iEtaMu=0;
     604      float iPhiMu=0,iEtaMu=0,ptisoMu=0;
    606605      for(unsigned int i=0; i < muon.size(); i++)
    607606        {
    … …  
    611610          elementMu->EtaCalo = muon[i].EtaCalo();
    612611          elementMu->PhiCalo = muon[i].PhiCalo();
    613           elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,DET->ISOL_PT,DET->ISOL_Cone);
    614           DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
     612          elementMu->IsolFlag = DET->Isolation(muon[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoMu);
     613          elementMu->IsolPt = ptisoMu;
     614          DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
    615615          D_CaloTower calMuon(list_of_active_towers.getElement(iEtaMu,iPhiMu));
    616           elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
     616          if( calMuon.getEem() !=0 ) elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
     617          else elementMu->EHoverEE = UNDEFINED;
     618          elementMu->EtRatio = DET->CaloIsolation(muon[i], list_of_active_towers);
    617619        }
    618620

  • trunk/interface/SmearUtil.h

    r319 r321  
    155155  float ISOL_PT;      //minimal pt of tracks for isolation criteria
    156156  float ISOL_Cone;    //Cone  for isolation criteria
     157  float ISOL_Calo_ET;  //minimal tower energy for isolation criteria
     158  float ISOL_Calo_Cone;//Cone for calorimetric isolation
     159  unsigned int ISOL_Calo_Grid; //Grid size (N x N) for calorimetric isolation
    157160 
    158161  //General jet variable
    … …  
    230233  bool Btaggedjet(const TLorentzVector &JET, const TSimpleArray<GenParticle> &subarray);
    231234
    232   /// Lepton isolation
    233   bool Isolation(const float phi, const float eta,const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone);
     235  /// Lepton isolation based on tracking
     236  bool Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso);
     237
     238  /// Lepton isolation based on calorimetry (optional. Default: off)
     239  float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers);
    234240
    235241  //********************* returns a segmented value for eta and phi, for calo towers *****

  • trunk/src/LHCOConverter.cc

    r316 r321  
    205205          case lhcoMuonID: {
    206206           cout << "ERROR: LHCOConverter: for muons use LHCOConverter::BranchReaderMuon!\n";
    207                 TRootMuon muon(*((TRootMuon*) branch->At(i)));
    208                 ntrk = muon.Charge;  ratioE = muon.EHoverEE;
    209207          } break;
    210208          case lhcoETmisID: {
    211209           cout << "ERROR: LHCOConverter: for ETmis use LHCOConverter::BranchReaderETmis!\n";
    212                 TRootMuon muon(*((TRootMuon*) branch->At(i)));
    213                 ntrk = muon.Charge;  ratioE = muon.EHoverEE;
    214210          } break;
    215211
    … …  
    284280        } // loop on jets
    285281
     282        int ratio = 0; char ratio_string[3];
     283        if(muon->EtRatio != UNDEFINED) ratio =  static_cast<int>(muon->EtRatio);
     284        if(ratio<0) sprintf(ratio_string,"00");
     285        else if(ratio<10) sprintf(ratio_string,"0%d",ratio);
     286        else if(ratio<100)sprintf(ratio_string,"%d",ratio);
     287        else { cout << "Error: muon EtRatio > 99" << endl; sprintf(ratio_string,"99"); }
    286288        // output
    287289        outfile << fixed << setprecision(3)
    … …  
    293295                << setw(7) << jmass <<" "       // invariant mass
    294296                << setw(7) << (float)muon->Charge  <<" "        // -1 for mu-, +1 for mu+
    295                 << setw(7) << (float)closest_jet_ID << " "
    296                 << setw(7) << muon->EHoverEE <<" "      // E_had/E_em
     297                << setw(7) << (float)closest_jet_ID << " ";
     298        outfile << fixed << setprecision(0)
     299                << setw(4) << floor(muon->IsolPt) << "." << ratio_string<<" ";  // E_had/E_em
     300        outfile << fixed << setprecision(3)
    297301                << setw(7) << 0.0 << setw(7) << 0.0  // dummy/dummy
    298302                << endl;

  • trunk/src/SmearUtil.cc

    r319 r321  
    102102  PTCUT_taujet    = 10.0;
    103103
     104  // Isolation
    104105  ISOL_PT         = 2.0;      //minimal pt of tracks for isolation criteria
    105106  ISOL_Cone       = 0.5;      //Cone  for isolation criteria
    106 
     107  ISOL_Calo_ET    = 1E99;     //minimal tower energy for isolation criteria. Default off = 1E99
     108  ISOL_Calo_Cone  = 0.5;      //Cone for calorimetric isolation
     109  ISOL_Calo_Grid  = 3;        //Grid size (N x N) for calorimetric isolation
    107110
    108111  // General jet variable
    … …  
    228231  PTCUT_taujet    = DET.PTCUT_taujet;
    229232
    230   ISOL_PT         = DET.ISOL_PT;      //minimal pt of tracks for isolation criteria
    231   ISOL_Cone       = DET.ISOL_Cone;      //Cone  for isolation criteria
     233  // Isolation
     234  ISOL_PT         = DET.ISOL_PT;     // tracking isolation
     235  ISOL_Cone       = DET.ISOL_Cone; 
     236  ISOL_Calo_ET    = DET.ISOL_Calo_ET;  // calorimeter isolation, defaut off
     237  ISOL_Calo_Cone  = DET.ISOL_Calo_Cone;
     238  ISOL_Calo_Grid  = DET.ISOL_Calo_Grid;
    232239
    233240
    … …  
    347354  PTCUT_taujet    = DET.PTCUT_taujet;
    348355
    349   ISOL_PT         = DET.ISOL_PT;      //minimal pt of tracks for isolation criteria
    350   ISOL_Cone       = DET.ISOL_Cone;      //Cone  for isolation criteria
    351 
     356  // Isolation
     357  ISOL_PT         = DET.ISOL_PT;       // tracking isolation
     358  ISOL_Cone       = DET.ISOL_Cone;   
     359  ISOL_Calo_ET    = DET.ISOL_Calo_ET;  // calorimeter isolation, defaut off
     360  ISOL_Calo_Cone  = DET.ISOL_Calo_Cone;
     361  ISOL_Calo_Grid  = DET.ISOL_Calo_Grid;
    352362
    353363  // General jet variable
    … …  
    497507    else if(strstr(temp_string.c_str(),"PTCUT_taujet"))     {curstring >> varname >> value; PTCUT_taujet      = value;}
    498508
    499     else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT        = value;}
    500     else if(strstr(temp_string.c_str(),"ISOL_Cone"))        {curstring >> varname >> value; ISOL_Cone        = value;}
    501 
     509    else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT           = value;}
     510    else if(strstr(temp_string.c_str(),"ISOL_Cone"))        {curstring >> varname >> value; ISOL_Cone         = value;}
     511    else if(strstr(temp_string.c_str(),"ISOL_Calo_ET"))     {curstring >> varname >> value; ISOL_Calo_ET      = value;}
     512    else if(strstr(temp_string.c_str(),"ISOL_Calo_Cone"))   {curstring >> varname >> value; ISOL_Calo_Cone    = value;}
     513    else if(strstr(temp_string.c_str(),"ISOL_Calo_Grid"))   {curstring >> varname >> ivalue; ISOL_Calo_Grid    = ivalue;}
    502514
    503515    else if(strstr(temp_string.c_str(),"JET_coneradius"))   {curstring >> varname >> value; JET_coneradius    = value;}
    … …  
    787799  f_out << left << setw(40) <<"* Cone for isolation criteria: "<<""
    788800        << left << setw(20) <<ISOL_Cone             <<""<< right << setw(10)<<"*"<<"\n";
     801
     802  if(ISOL_Calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied                                  *"<<"\n";
     803  else {
     804   f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
     805         << left << setw(20) <<ISOL_Calo_ET               <<""<< right << setw(10)<<"*"<<"\n";
     806   f_out << left << setw(40) <<"* Cone for calorimetric isolation: "<<""
     807         << left << setw(20) <<ISOL_Calo_Cone             <<""<< right << setw(10)<<"*"<<"\n";
     808   f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
     809         << left << setw(20) <<ISOL_Calo_Grid             <<""<< right << setw(10)<<"*"<<"\n";
     810  }
     811
     812
    789813  f_out<<"*                                                                    *"<<"\n";
    790814  f_out<<"#***************                                                     *"<<"\n";
    … …  
    10721096//***********************Isolation criteria***********************
    10731097//****************************************************************
    1074 bool RESOLution::Isolation(const float phi, const float eta,const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone)
     1098bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso )
    10751099{
    10761100   bool isolated = false;
    1077    float deltar=5000.; // Initial value; should be high; no further repercussion
    1078    // loop on all final charged particles, with p_t >2, close enough from the electron
    1079    for(unsigned int i=0; i < tracks.size(); i++)
    1080       {
    1081        if(tracks[i].PT < pt_second_track)continue;
    1082        float genDeltaR = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta);
     1101   ptiso = 0; // sum of all track pt in isolation cone
     1102   float deltar=1E99; // Initial value; should be high; no further repercussion
     1103
     1104   // loop on all tracks, with p_t above threshold, close enough from the charged lepton
     1105   for(unsigned int i=0; i < tracks.size(); i++) {
     1106       if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
     1107       float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta);
    10831108         if(
    10841109             (genDeltaR > deltar) ||
    1085              (genDeltaR==0)
     1110             (genDeltaR==0)  // rejets the track of the particle itself
    10861111           ) continue ;
    1087            deltar=genDeltaR;
     1112           deltar=genDeltaR;    // finds the closest track
     1113
     1114         // as long as (genDeltaR==0) is put above, the particle itself is not taken into account
     1115         if( genDeltaR < ISOL_Cone) ptiso += tracks[i].PT; // dR cut on tracks
    10881116      }
    10891117   if(deltar > isolCone) isolated = true;
    10901118   return isolated;
     1119}
     1120
     1121// ******* Calorimetric isolation
     1122float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers) {
     1123 // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy
     1124 // 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.
     1125   if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
     1126   float etrat=0;
     1127   // available parameters: ISOL_Calo_ET , ISOL_Calo_Cone ,
     1128/*   for(unsigned int i=0; i < towers.size(); i++) {
     1129     if(towers[i].E > ISOL_Calo_ET) {
     1130        float genDeltaR = DeltaR(part.Phi(),part.Eta(),towers[i].getPhi(),towers[i].getEta());
     1131        if(genDeltaR < ISOL_Calo_Cone) {
     1132           ptiso += towers[i].getET();
     1133        }           
     1134     }
     1135   } // loop on towers
     1136   ptiso -=
     1137*/
     1138  etrat = 100*etrat/part.Pt();
     1139  if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
     1140  else if(etrat>99) cout << "Error: etrat shoud be in [0;99] in CaloIsolation (" << etrat <<")\n";
     1141  return etrat;
    10911142}
    10921143
Note: See TracChangeset for help on using the changeset viewer.