LeptonSelection/src/FSR.cc

#include <vector>
#include "MathUtils.h"
#include "FSR.h"
#include "IsolationSelection.h"
#include "CommonDefs.h"
#include "TLorentzVector.h"

using namespace std;

extern vector<bool> PFnoPUflag;

//--------------------------------------------------------------------------------------------------
// typeI = PF IDed photons.  NB : repurpose PFnoPUflag, flip for recovered photons  
// so that they are skipped in the isolation calculation
//--------------------------------------------------------------------------------------------------
bool recover_typeI_Photon( ControlFlags & ctrl,
                           mithep::Electron * el, 
                           const int electronIndex,
                           vector<SimpleLepton> &lepvec,
                           const mithep::Array<mithep::PFCandidate> * pfArr,
                           const mithep::Array<mithep::Electron> *eleArr,
                           TLorentzVector * Zvec,
                           vector<const mithep::PFCandidate*> &photonsToVeto ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[electronIndex].fsrRecoveryAttempted ) return false;

  vector<int> photonIndices; 
  for( int i=0; i<pfArr->GetEntries(); i++ ) { 
    if( !(PFnoPUflag[i])) continue; // my PF no PU hack
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
        pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {

      if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
      //      float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
      float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                           lepvec[electronIndex].vec.Phi(), lepvec[electronIndex].vec.Eta());
      if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;

      //
      // veto if close to an electron SC 
      //
      bool flagEleSC = false;
      for( int j=0; j<lepvec.size(); j++ ) { 
        if( !(abs(lepvec[j].type) == 11 ) )      continue;
        if( !(lepvec[j].status.looseIDAndPre()) ) continue;
        double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi(); 
        float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
        float dEta = fabs(pf->Eta()-eeta);
        float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
        if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi 
                            << "\tdEta: " << dEta 
                            << "\tetaPH: " << pf->Eta()
                            << "\tetaELH: " << eeta
                            << "\tdR:" << dR << endl;
        if( (dPhi<2.&& dEta<0.05) || dR<0.15 ) { 
            flagEleSC = true;
            break;
        }
        if( flagEleSC ) break;
      }
      if( flagEleSC ) continue;
      if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;


      //
      // check that input electron is the closest lepton to this photon
      //
      bool found_closer_lepton=false;
      for( int j=0; j<lepvec.size(); j++ ) { 
        if( j == electronIndex ) continue;
        if( !(lepvec[j].status.looseIDAndPre()) ) continue;
        float tmp_dR =  mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                                  lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
        if( tmp_dR < dR ) {
          if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
                              <<tmp_dR<<" vs "<<dR<<") skipping..." << endl;   
        found_closer_lepton=true;
        break;
        }
      }
      if( found_closer_lepton ) continue;


      //
      // Z mass OK?
      //
      TLorentzVector pvec;
      pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
      float newMass = (pvec + *Zvec).M(); 
      if( !( newMass > 4.   && 
             newMass < 100. && 
             (fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
             ) ) continue;
      if( ctrl.debug ) std::cerr << "FSR :: improved Zmass  ... " <<
        Zvec->M() << " -> " << newMass << std::endl;


      //
      // "keep all photons close to one of the 4L electrons ..."
      //
      if( dR < 0.07 ) { 
        if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing  ... " << std::endl;
        photonIndices.push_back(i);
      }

      //      
      // "need tighter cuts for other photons ..."
      //
      if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
      //      if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
        if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing  ... " << std::endl;
        photonIndices.push_back(i);
      }
    }
  }    

  float highest_pt  = -1;   int highest_pt_index=-1;
  float smallest_dR = 999.; int smallest_dR_index=-1;
  for( int i=0; i<photonIndices.size(); i++ ) {
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)(pfArr->At(photonIndices[i]));
    float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
    if( pf->Pt() > highest_pt ) { 
      highest_pt_index = photonIndices[i];
      highest_pt       = pf->Pt();
    }
    if( dR < smallest_dR ) { 
      smallest_dR_index = photonIndices[i];
      smallest_dR       = dR;
    }
  }

  const mithep::PFCandidate * thepf;
  if( highest_pt > 4. ) { 
    thepf  = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
    // "... remove it from lepton isolation ..."
    //    PFnoPUflag[highest_pt_index] = 0;
    // TMP, commented flip above for FSR study 
    // gammaMatches[highest_pt_index].push_back(lepvec[electronIndex].index);
    photonsToVeto.push_back(thepf);
  } else if( smallest_dR != 999. ) {
    thepf  = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
    // "... remove it from lepton isolation ..."
    // PFnoPUflag[smallest_dR_index] = 0;
    // TMP, commented flip above for FSR study 
    // gammaMatches[smallest_dR_index].push_back(lepvec[electronIndex].index);
    photonsToVeto.push_back(thepf);
  } else { 
    return false;
  }
  
  if( thepf != NULL ) { 
    // add to the electron
    TLorentzVector elvec,phvec,newelvec;
    elvec.SetPtEtaPhiM( el->Pt(), el->Eta(), el->Phi(), ELECTRON_MASS);
    phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
    newelvec = elvec+phvec;
    // don't update the electron object, just simplelepton
    //     el->SetPtEtaPhi        (newelvec.Pt(),
    //                      newelvec.Eta(),
    //                      newelvec.Phi());
    lepvec[electronIndex].vec += phvec;
    lepvec[electronIndex].fsrRecoveryAttempted = true;
    return true;      
  }
  return false;
}


//--------------------------------------------------------------------------------------------------
// typeI = PF IDed photons.  NB : repurpose PFnoPUflag, flip for recovered photons  
// so that they are skipped in the isolation calculation
//--------------------------------------------------------------------------------------------------
bool recover_typeI_Photon( ControlFlags & ctrl,
                           mithep::Muon * mu, 
                           const int muonIndex,
                           vector<SimpleLepton> &lepvec,
                           const mithep::Array<mithep::PFCandidate> * pfArr,
                           const mithep::Array<mithep::Electron> *eleArr,
                           TLorentzVector * Zvec, 
                           vector<const mithep::PFCandidate*> &photonsToVeto ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;

  vector<int> photonIndices; 
  for( int i=0; i<pfArr->GetEntries(); i++ ) { 
    if( !(PFnoPUflag[i])) continue; // my PF no PU hack
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
        pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {

      if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
      //      float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
      float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                           lepvec[muonIndex].vec.Phi(), lepvec[muonIndex].vec.Eta());
      if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;

      //
      // veto if close to an electron SC 
      //
      bool flagEleSC = false;
      for( int j=0; j<lepvec.size(); j++ ) { 
        if( !(abs(lepvec[j].type) == 11 ) )      continue;
        if( !(lepvec[j].status.looseIDAndPre()) ) continue;
        double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi(); 
        float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
        float dEta = fabs(pf->Eta()-eeta);
        float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
        if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi 
                            << "\tdEta: " << dEta 
                            << "\tetaPH: " << pf->Eta()
                            << "\tetaELH: " << eeta
                            << "\tdR:" << dR << endl;
        if( (dPhi<2.&& dEta<0.05) || dR<0.15 ) { 
            flagEleSC = true;
            break;
        }
        if( flagEleSC ) break;
      }
      if( flagEleSC ) continue;
      if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;


      //
      // check that input muon is the closest lepton to this photon
      //
      bool found_closer_lepton=false;
      for( int j=0; j<lepvec.size(); j++ ) { 
        if( j == muonIndex ) continue;
        if( !(lepvec[j].status.looseIDAndPre()) ) continue;
        float tmp_dR =  mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                                  lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
        if( tmp_dR < dR ) {
          if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
                              <<tmp_dR<<" vs "<<dR<<") skipping..." << endl;   
          found_closer_lepton=true;
          break;
        }
      }
      if( found_closer_lepton ) continue;


      //
      // Z mass OK?
      //
      TLorentzVector pvec;
      pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
      float newMass = (pvec + *Zvec).M(); 
      if( !( newMass > 4.   && 
             newMass < 100. && 
             (fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
             ) ) continue;
      if( ctrl.debug ) std::cerr << "FSR :: improved Zmass  ... " <<
        Zvec->M() << " -> " << newMass << std::endl;

      //
      // "keep all photons close to one of the 4L muons ..."
      //
      if( dR < 0.07 ) { 
        if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing  ... " << std::endl;
        photonIndices.push_back(i);
      }

      //      
      // "need tighter cuts for other photons ..."
      //
      if( ctrl.debug ) std::cerr << "FSR :: pass tighter?, pT: " << pf->Pt() << std::endl;
      if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
      //      if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
        if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing index= " << i  << std::endl;
        photonIndices.push_back(i);
      }
    }
  }    

  float highest_pt  = -1;   int highest_pt_index=-1;
  float smallest_dR = 999.; int smallest_dR_index=-1;
  for( int i=0; i<photonIndices.size(); i++ ) {
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)(pfArr->At(photonIndices[i]));
    float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
    if( pf->Pt() > highest_pt ) { 
      highest_pt_index = photonIndices[i];
      highest_pt       = pf->Pt();
    }
    if( dR < smallest_dR ) { 
      smallest_dR_index = photonIndices[i];
      smallest_dR       = dR;
    }
  }

  const mithep::PFCandidate * thepf;
  if( highest_pt > 4. ) { 
    if(ctrl.debug) std::cerr << "FSR :: taking highest pt gamma, index = " <<  highest_pt_index << endl;
    thepf  = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
    // "... remove it from lepton isolation ..."
    //    PFnoPUflag[highest_pt_index] = 0;
    // TMP, commented flip above for FSR study 
    // gammaMatches[highest_pt_index].push_back(lepvec[muonIndex].index);
    photonsToVeto.push_back(thepf);
  } else if( smallest_dR != 999. ) {
    if(ctrl.debug) std::cerr << "FSR :: taking smallest dR gamma, index = " <<  highest_pt_index << endl;
    thepf  = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
    // "... remove it from lepton isolation ..."
    //    PFnoPUflag[smallest_dR_index] = 0;
    // TMP, commented flip above for FSR study 
    //gammaMatches[smallest_dR_index].push_back(lepvec[muonIndex].index);
    photonsToVeto.push_back(thepf);
  } else { 
    return false;
  }

  TLorentzVector pvec;
  if( thepf != NULL ) { 
    // add to the muon
    if( ctrl.debug ) cerr << "FSR :: before return, oldpT=" << mu->Pt() << endl;
    TLorentzVector muvec,phvec,newmuvec;
    muvec.SetPtEtaPhiM( mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS);
    phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
    pvec = phvec;
    newmuvec = muvec+phvec;
    // don't update the muon object, just simplelepton
    //     mu->SetPtEtaPhi        (newmuvec.Pt(),
    //                      newmuvec.Eta(),
    //                      newmuvec.Phi());
    lepvec[muonIndex].vec += phvec;
    lepvec[muonIndex].fsrRecoveryAttempted = true;
    return true;      
  }
  return false;
}


//--------------------------------------------------------------------------------------------------
// typeII = "PFClusters linked to muons"
//--------------------------------------------------------------------------------------------------
bool recover_typeII_Photon( ControlFlags & ctrl,
                            mithep::Muon * mu, 
                            const int muonIndex,
                            vector<SimpleLepton> &lepvec,
                            const mithep::Array<mithep::PFCandidate> * pfArr ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;

  bool foundPF=false;
  const mithep::PFCandidate * thepf;
  for( int i=0; i<pfArr->GetEntries(); i++ ) { 
    if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
    if( abs(pf->PFType()) == mithep::PFCandidate::eMuon 
        && (mu->TrackerTrk()==pf->TrackerTrk())) 
      {
        if(ctrl.debug) cout << "FSR :: t2, found pf muon, pt " <<  mu->Pt() << endl;
        foundPF = true;
        thepf = pf;
        break;
      }
  }
  
  if( foundPF ) {
    double sintet = thepf->Pt()/thepf->E();
    double phpt = thepf->EECal() * sintet;
    if ( thepf->EECal() >= 2.0 && phpt >= 2.0 ) {
      // don't update the muon object, just simplelepton
      //       mu->SetPtEtaPhi        (mu->Pt()+phpt,
      //                              mu->Eta(),
      //                              mu->Phi());
      TLorentzVector pvec;
      pvec.SetPtEtaPhiM( phpt,mu->Eta(),mu->Phi(),0.);
      lepvec[muonIndex].vec += pvec;
      if(ctrl.debug) cout << "FSR :: t2, new pt " <<  lepvec[muonIndex].vec.Pt() << endl;
      lepvec[muonIndex].fsrRecoveryAttempted = true;  
      return true;
    }
  }
  
  return false;
}
Revision:	1.4
Committed:	Tue Jun 5 19:35:47 2012 UTC (12 years, 11 months ago) by khahn
Content type:	text/plain
Branch:	MAIN
CVS Tags:	synced_FSR
Changes since 1.3:	+74 -60 lines
Log Message:	final sync for FSR
#	User	Rev	Content
1	khahn	1.1	#include <vector>
2			#include "MathUtils.h"
3			#include "FSR.h"
4	khahn	1.3	#include "IsolationSelection.h"
5	khahn	1.1	#include "CommonDefs.h"
6			#include "TLorentzVector.h"
7
8			using namespace std;
9
10			extern vector<bool> PFnoPUflag;
11
12			//--------------------------------------------------------------------------------------------------
13	khahn	1.3	// typeI = PF IDed photons. NB : repurpose PFnoPUflag, flip for recovered photons
14			// so that they are skipped in the isolation calculation
15			//--------------------------------------------------------------------------------------------------
16			bool recover_typeI_Photon( ControlFlags & ctrl,
17			mithep::Electron * el,
18			const int electronIndex,
19			vector<SimpleLepton> &lepvec,
20			const mithep::Array<mithep::PFCandidate> * pfArr,
21			const mithep::Array<mithep::Electron> *eleArr,
22	khahn	1.4	TLorentzVector * Zvec,
23			vector<const mithep::PFCandidate*> &photonsToVeto )
24	khahn	1.2	//--------------------------------------------------------------------------------------------------
25	khahn	1.3	{
26	khahn	1.4	if( lepvec[electronIndex].fsrRecoveryAttempted ) return false;
27
28	khahn	1.3	vector<int> photonIndices;
29			for( int i=0; i<pfArr->GetEntries(); i++ ) {
30			if( !(PFnoPUflag[i])) continue; // my PF no PU hack
31			const mithep::PFCandidate pf = (mithep::PFCandidate)((*pfArr)[i]);
32			if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
33			pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {
34	khahn	1.2
35	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
36	khahn	1.4	// float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
37			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
38			lepvec[electronIndex].vec.Phi(), lepvec[electronIndex].vec.Eta());
39	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;
40	khahn	1.2
41	khahn	1.3	//
42			// veto if close to an electron SC
43			//
44			bool flagEleSC = false;
45			for( int j=0; j<lepvec.size(); j++ ) {
46			if( !(abs(lepvec[j].type) == 11 ) ) continue;
47			if( !(lepvec[j].status.looseIDAndPre()) ) continue;
48			double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi();
49			float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
50			float dEta = fabs(pf->Eta()-eeta);
51			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
52			if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi
53			<< "\tdEta: " << dEta
54			<< "\tetaPH: " << pf->Eta()
55			<< "\tetaELH: " << eeta
56			<< "\tdR:" << dR << endl;
57			if( (dPhi<2.&& dEta<0.05) \|\| dR<0.15 ) {
58			flagEleSC = true;
59			break;
60			}
61			if( flagEleSC ) break;
62			}
63			if( flagEleSC ) continue;
64			if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;
65	khahn	1.2
66	khahn	1.4
67			//
68			// check that input electron is the closest lepton to this photon
69			//
70			bool found_closer_lepton=false;
71			for( int j=0; j<lepvec.size(); j++ ) {
72			if( j == electronIndex ) continue;
73			if( !(lepvec[j].status.looseIDAndPre()) ) continue;
74			float tmp_dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
75			lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
76			if( tmp_dR < dR ) {
77			if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
78			<<tmp_dR<<" vs "<<dR<<") skipping..." << endl;
79			found_closer_lepton=true;
80			break;
81			}
82			}
83			if( found_closer_lepton ) continue;
84
85
86	khahn	1.2	//
87	khahn	1.3	// Z mass OK?
88	khahn	1.2	//
89	khahn	1.3	TLorentzVector pvec;
90			pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
91			float newMass = (pvec + *Zvec).M();
92			if( !( newMass > 4. &&
93			newMass < 100. &&
94			(fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
95			) ) continue;
96			if( ctrl.debug ) std::cerr << "FSR :: improved Zmass ... " <<
97			Zvec->M() << " -> " << newMass << std::endl;
98
99
100	khahn	1.2	//
101	khahn	1.3	// "keep all photons close to one of the 4L electrons ..."
102	khahn	1.2	//
103	khahn	1.3	if( dR < 0.07 ) {
104			if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing ... " << std::endl;
105			photonIndices.push_back(i);
106	khahn	1.2	}
107	khahn	1.3
108			//
109			// "need tighter cuts for other photons ..."
110	khahn	1.2	//
111	khahn	1.4	if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
112			// if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
113	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing ... " << std::endl;
114			photonIndices.push_back(i);
115	khahn	1.2	}
116			}
117	khahn	1.3	}
118
119			float highest_pt = -1; int highest_pt_index=-1;
120			float smallest_dR = 999.; int smallest_dR_index=-1;
121			for( int i=0; i<photonIndices.size(); i++ ) {
122			const mithep::PFCandidate pf = (mithep::PFCandidate)(pfArr->At(photonIndices[i]));
123			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
124			if( pf->Pt() > highest_pt ) {
125			highest_pt_index = photonIndices[i];
126			highest_pt = pf->Pt();
127			}
128			if( dR < smallest_dR ) {
129			smallest_dR_index = photonIndices[i];
130			smallest_dR = dR;
131			}
132	khahn	1.2	}
133
134	khahn	1.3	const mithep::PFCandidate * thepf;
135			if( highest_pt > 4. ) {
136	khahn	1.4	thepf = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
137	khahn	1.3	// "... remove it from lepton isolation ..."
138	khahn	1.4	// PFnoPUflag[highest_pt_index] = 0;
139	khahn	1.3	// TMP, commented flip above for FSR study
140			// gammaMatches[highest_pt_index].push_back(lepvec[electronIndex].index);
141	khahn	1.4	photonsToVeto.push_back(thepf);
142	khahn	1.3	} else if( smallest_dR != 999. ) {
143	khahn	1.4	thepf = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
144	khahn	1.3	// "... remove it from lepton isolation ..."
145	khahn	1.4	// PFnoPUflag[smallest_dR_index] = 0;
146	khahn	1.3	// TMP, commented flip above for FSR study
147			// gammaMatches[smallest_dR_index].push_back(lepvec[electronIndex].index);
148	khahn	1.4	photonsToVeto.push_back(thepf);
149	khahn	1.3	} else {
150			return false;
151			}
152
153			if( thepf != NULL ) {
154			// add to the electron
155			TLorentzVector elvec,phvec,newelvec;
156			elvec.SetPtEtaPhiM( el->Pt(), el->Eta(), el->Phi(), ELECTRON_MASS);
157			phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
158			newelvec = elvec+phvec;
159			// don't update the electron object, just simplelepton
160			// el->SetPtEtaPhi (newelvec.Pt(),
161			// newelvec.Eta(),
162			// newelvec.Phi());
163			lepvec[electronIndex].vec += phvec;
164	khahn	1.4	lepvec[electronIndex].fsrRecoveryAttempted = true;
165	khahn	1.3	return true;
166			}
167			return false;
168	khahn	1.2	}
169
170
171
172			//--------------------------------------------------------------------------------------------------
173	khahn	1.1	// typeI = PF IDed photons. NB : repurpose PFnoPUflag, flip for recovered photons
174			// so that they are skipped in the isolation calculation
175			//--------------------------------------------------------------------------------------------------
176	khahn	1.2	bool recover_typeI_Photon( ControlFlags & ctrl,
177			mithep::Muon * mu,
178	khahn	1.3	const int muonIndex,
179			vector<SimpleLepton> &lepvec,
180			const mithep::Array<mithep::PFCandidate> * pfArr,
181			const mithep::Array<mithep::Electron> *eleArr,
182	khahn	1.4	TLorentzVector * Zvec,
183			vector<const mithep::PFCandidate*> &photonsToVeto )
184	khahn	1.1	//--------------------------------------------------------------------------------------------------
185			{
186	khahn	1.4	if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;
187
188	khahn	1.1	vector<int> photonIndices;
189			for( int i=0; i<pfArr->GetEntries(); i++ ) {
190			if( !(PFnoPUflag[i])) continue; // my PF no PU hack
191			const mithep::PFCandidate pf = (mithep::PFCandidate)((*pfArr)[i]);
192			if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
193			pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {
194
195	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
196	khahn	1.4	// float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
197			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
198			lepvec[muonIndex].vec.Phi(), lepvec[muonIndex].vec.Eta());
199	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;
200
201			//
202			// veto if close to an electron SC
203			//
204			bool flagEleSC = false;
205			for( int j=0; j<lepvec.size(); j++ ) {
206			if( !(abs(lepvec[j].type) == 11 ) ) continue;
207			if( !(lepvec[j].status.looseIDAndPre()) ) continue;
208			double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi();
209			float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
210			float dEta = fabs(pf->Eta()-eeta);
211			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
212			if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi
213			<< "\tdEta: " << dEta
214			<< "\tetaPH: " << pf->Eta()
215			<< "\tetaELH: " << eeta
216			<< "\tdR:" << dR << endl;
217			if( (dPhi<2.&& dEta<0.05) \|\| dR<0.15 ) {
218			flagEleSC = true;
219			break;
220			}
221			if( flagEleSC ) break;
222			}
223			if( flagEleSC ) continue;
224			if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;
225
226	khahn	1.4
227			//
228			// check that input muon is the closest lepton to this photon
229			//
230			bool found_closer_lepton=false;
231			for( int j=0; j<lepvec.size(); j++ ) {
232			if( j == muonIndex ) continue;
233			if( !(lepvec[j].status.looseIDAndPre()) ) continue;
234			float tmp_dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
235			lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
236			if( tmp_dR < dR ) {
237			if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
238			<<tmp_dR<<" vs "<<dR<<") skipping..." << endl;
239			found_closer_lepton=true;
240			break;
241			}
242			}
243			if( found_closer_lepton ) continue;
244
245
246
247	khahn	1.3	//
248			// Z mass OK?
249			//
250			TLorentzVector pvec;
251			pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
252			float newMass = (pvec + *Zvec).M();
253			if( !( newMass > 4. &&
254			newMass < 100. &&
255			(fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
256			) ) continue;
257			if( ctrl.debug ) std::cerr << "FSR :: improved Zmass ... " <<
258			Zvec->M() << " -> " << newMass << std::endl;
259
260			//
261	khahn	1.1	// "keep all photons close to one of the 4L muons ..."
262	khahn	1.3	//
263			if( dR < 0.07 ) {
264			if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing ... " << std::endl;
265			photonIndices.push_back(i);
266			}
267
268			//
269	khahn	1.1	// "need tighter cuts for other photons ..."
270	khahn	1.3	//
271			if( ctrl.debug ) std::cerr << "FSR :: pass tighter?, pT: " << pf->Pt() << std::endl;
272	khahn	1.4	if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
273			// if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
274	khahn	1.3	if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing index= " << i << std::endl;
275	khahn	1.2	photonIndices.push_back(i);
276	khahn	1.3	}
277	khahn	1.1	}
278	khahn	1.3	}
279
280			float highest_pt = -1; int highest_pt_index=-1;
281			float smallest_dR = 999.; int smallest_dR_index=-1;
282			for( int i=0; i<photonIndices.size(); i++ ) {
283			const mithep::PFCandidate pf = (mithep::PFCandidate)(pfArr->At(photonIndices[i]));
284			float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
285			if( pf->Pt() > highest_pt ) {
286			highest_pt_index = photonIndices[i];
287			highest_pt = pf->Pt();
288			}
289			if( dR < smallest_dR ) {
290			smallest_dR_index = photonIndices[i];
291			smallest_dR = dR;
292	khahn	1.1	}
293			}
294	khahn	1.3
295			const mithep::PFCandidate * thepf;
296			if( highest_pt > 4. ) {
297			if(ctrl.debug) std::cerr << "FSR :: taking highest pt gamma, index = " << highest_pt_index << endl;
298	khahn	1.4	thepf = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
299	khahn	1.3	// "... remove it from lepton isolation ..."
300	khahn	1.4	// PFnoPUflag[highest_pt_index] = 0;
301	khahn	1.3	// TMP, commented flip above for FSR study
302			// gammaMatches[highest_pt_index].push_back(lepvec[muonIndex].index);
303	khahn	1.4	photonsToVeto.push_back(thepf);
304	khahn	1.3	} else if( smallest_dR != 999. ) {
305			if(ctrl.debug) std::cerr << "FSR :: taking smallest dR gamma, index = " << highest_pt_index << endl;
306	khahn	1.4	thepf = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
307	khahn	1.3	// "... remove it from lepton isolation ..."
308	khahn	1.4	// PFnoPUflag[smallest_dR_index] = 0;
309	khahn	1.3	// TMP, commented flip above for FSR study
310			//gammaMatches[smallest_dR_index].push_back(lepvec[muonIndex].index);
311	khahn	1.4	photonsToVeto.push_back(thepf);
312	khahn	1.3	} else {
313			return false;
314			}
315
316			TLorentzVector pvec;
317			if( thepf != NULL ) {
318			// add to the muon
319			if( ctrl.debug ) cerr << "FSR :: before return, oldpT=" << mu->Pt() << endl;
320			TLorentzVector muvec,phvec,newmuvec;
321			muvec.SetPtEtaPhiM( mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS);
322			phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
323			pvec = phvec;
324			newmuvec = muvec+phvec;
325			// don't update the muon object, just simplelepton
326			// mu->SetPtEtaPhi (newmuvec.Pt(),
327			// newmuvec.Eta(),
328			// newmuvec.Phi());
329			lepvec[muonIndex].vec += phvec;
330	khahn	1.4	lepvec[muonIndex].fsrRecoveryAttempted = true;
331	khahn	1.3	return true;
332			}
333	khahn	1.1	return false;
334			}
335
336	khahn	1.3
337
338	khahn	1.1	//--------------------------------------------------------------------------------------------------
339			// typeII = "PFClusters linked to muons"
340			//--------------------------------------------------------------------------------------------------
341	khahn	1.2	bool recover_typeII_Photon( ControlFlags & ctrl,
342			mithep::Muon * mu,
343	khahn	1.3	const int muonIndex,
344			vector<SimpleLepton> &lepvec,
345	khahn	1.1	const mithep::Array<mithep::PFCandidate> * pfArr )
346			//--------------------------------------------------------------------------------------------------
347			{
348	khahn	1.3	if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;
349	khahn	1.1
350			bool foundPF=false;
351			const mithep::PFCandidate * thepf;
352			for( int i=0; i<pfArr->GetEntries(); i++ ) {
353			if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
354			const mithep::PFCandidate pf = (mithep::PFCandidate)((*pfArr)[i]);
355	khahn	1.3	if( abs(pf->PFType()) == mithep::PFCandidate::eMuon
356			&& (mu->TrackerTrk()==pf->TrackerTrk()))
357			{
358			if(ctrl.debug) cout << "FSR :: t2, found pf muon, pt " << mu->Pt() << endl;
359			foundPF = true;
360			thepf = pf;
361			break;
362			}
363	khahn	1.1	}
364
365			if( foundPF ) {
366			double sintet = thepf->Pt()/thepf->E();
367			double phpt = thepf->EECal() * sintet;
368			if ( thepf->EECal() >= 2.0 && phpt >= 2.0 ) {
369	khahn	1.3	// don't update the muon object, just simplelepton
370			// mu->SetPtEtaPhi (mu->Pt()+phpt,
371			// mu->Eta(),
372			// mu->Phi());
373			TLorentzVector pvec;
374			pvec.SetPtEtaPhiM( phpt,mu->Eta(),mu->Phi(),0.);
375			lepvec[muonIndex].vec += pvec;
376			if(ctrl.debug) cout << "FSR :: t2, new pt " << lepvec[muonIndex].vec.Pt() << endl;
377	khahn	1.4	lepvec[muonIndex].fsrRecoveryAttempted = true;
378	khahn	1.1	return true;
379			}
380			}
381
382			return false;
383			}