LeptonSelection/src/FSR.cc

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

extern vector<bool> PFnoPUflag;

void addPhoton(EventData &ret, TLorentzVector &pvec)
{
  assert(ret.fsrPhotons.size() < 3);
  SimpleLepton photon;
  photon.vec.SetPtEtaPhiM( pvec.Pt(), pvec.Eta(), pvec.Phi(), 0);

  // don't add the same photon twice...
  if(ret.fsrPhotons.size() > 0)
    if(dr(ret.fsrPhotons[0], photon) < 0.01) return;
  if(ret.fsrPhotons.size() > 1)
    if(dr(ret.fsrPhotons[1], photon) < 0.01) return;

  // cout << "  possibly adding photon: "
  //      << setw(12) << pvec.Pt()
  //      << " existing: "
  //      << setw(12) << ( (ret.fsrPhotons.size() > 0) ? ret.fsrPhotons[0].vec.Pt() : -1 )
  //      << setw(12) << ( (ret.fsrPhotons.size() > 1) ? ret.fsrPhotons[1].vec.Pt() : -1 )
  //      << endl;
    
  if(ret.fsrPhotons.size() == 2) {
    // only want to save two, so if there's already two in the vector, store the highest two out of the three photons
    sort( ret.fsrPhotons.begin(), ret.fsrPhotons.end(), SimpleLepton::lep_pt_sort );
    assert(ret.fsrPhotons[0].vec.Pt() > ret.fsrPhotons[1].vec.Pt());
    if(ret.fsrPhotons[1].vec.Pt() < photon.vec.Pt()) {
      ret.fsrPhotons.erase(ret.fsrPhotons.begin() + 1);
      ret.fsrPhotons.push_back(photon);
    }
  } else
    ret.fsrPhotons.push_back(photon);
  sort( ret.fsrPhotons.begin(), ret.fsrPhotons.end(), SimpleLepton::lep_pt_sort );
  // cout << "  possibly added photon: "
  //      << setw(12) << pvec.Pt()
  //      << " existing: "
  //      << setw(12) << ( (ret.fsrPhotons.size() > 0) ? ret.fsrPhotons[0].vec.Pt() : -1 )
  //      << setw(12) << ( (ret.fsrPhotons.size() > 1) ? ret.fsrPhotons[1].vec.Pt() : -1 )
  //      << endl;
}    
//--------------------------------------------------------------------------------------------------
// 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,
                           EventData &ret,
                           Electron * el, 
                           const int electronIndex,
                           vector<SimpleLepton> &lepvec,
                           const Array<PFCandidate> * pfArr,
                           const Array<Electron> *eleArr,
                           TLorentzVector * Zvec,
                           vector<const PFCandidate*> &photonsToVeto ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[electronIndex].fsrRecoveryAttempted ) return false;

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

      if( ctrl.debug ) cout << "    FSR :: pass preselection ... pt: "<< pf->Pt() << endl;
      //      float dR = MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
      float dR = MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                           lepvec[electronIndex].vec.Phi(), lepvec[electronIndex].vec.Eta());
      if( ctrl.debug ) cout << "    FSR :: dR = " << dR << 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(MathUtils::DeltaPhi(pf->Phi(),ephi));
        float dEta = fabs(pf->Eta()-eeta);
        float dR = 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 ) cout << "    FSR :: not matched to an ele SC ... " << 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 =  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 ) cout << "    FSR :: improved Zmass  ... " <<
        Zvec->M() << " -> " << newMass << endl;


      //
      // "keep all photons close to one of the 4L electrons ..."
      //
      if( dR < 0.07 ) { 
        if( ctrl.debug ) cout << "    FSR :: dR < 0.07, pushing  ... " << 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 ) cout << "    FSR :: tighter cuts, pushing  ... " << 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 PFCandidate *pf = (PFCandidate*)(pfArr->At(photonIndices[i]));
    float dR = 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 PFCandidate * thepf;
  if( highest_pt > 4. ) { 
    thepf  = (const 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 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.);
    addPhoton(ret,phvec);
    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,
                           EventData &ret,
                           Muon * mu, 
                           const int muonIndex,
                           vector<SimpleLepton> &lepvec,
                           const Array<PFCandidate> * pfArr,
                           const Array<Electron> *eleArr,
                           TLorentzVector * Zvec, 
                           vector<const PFCandidate*> &photonsToVeto ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;

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

      if( ctrl.debug ) cout << "    FSR :: pass preselection ... pt: "<< pf->Pt() << endl;
      //      float dR = MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
      float dR = MathUtils::DeltaR(pf->Phi(),pf->Eta(), 
                                           lepvec[muonIndex].vec.Phi(), lepvec[muonIndex].vec.Eta());
      if( ctrl.debug ) cout << "    FSR :: dR = " << dR << 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(MathUtils::DeltaPhi(pf->Phi(),ephi));
        float dEta = fabs(pf->Eta()-eeta);
        float dR = 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 ) cout << "    FSR :: not matched to an ele SC ... " << 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 =  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 ) cout << "    FSR :: improved Zmass  ... " <<
        Zvec->M() << " -> " << newMass << endl;

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

      //      
      // "need tighter cuts for other photons ..."
      //
      if( ctrl.debug ) cout << "    FSR :: pass tighter?, pT: " << pf->Pt() << 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 ) cout << "    FSR :: tighter cuts, pushing index= " << i  << 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 PFCandidate *pf = (PFCandidate*)(pfArr->At(photonIndices[i]));
    float dR = 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 PFCandidate * thepf;
  if( highest_pt > 4. ) { 
    if(ctrl.debug) cout << "    FSR :: taking highest pt gamma, index = " <<  highest_pt_index << endl;
    thepf  = (const 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) cout << "    FSR :: taking smallest dR gamma, index = " <<  highest_pt_index << endl;
    thepf  = (const 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 ) cout << "    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.);
    addPhoton(ret,phvec);
    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,
                            EventData &ret,
                            Muon * mu, 
                            const int muonIndex,
                            vector<SimpleLepton> &lepvec,
                            const Array<PFCandidate> * pfArr ) 
//--------------------------------------------------------------------------------------------------
{
  if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;

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