Mods/src/PhotonCiCMod.cc

// $Id: PhotonCiCMod.cc,v 1.2 2011/06/02 14:18:07 bendavid Exp $

#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
#include "MitAna/DataTree/interface/PhotonCol.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitPhysics/Utils/interface/PhotonTools.h"

using namespace mithep;

ClassImp(mithep::PhotonCiCMod)

//--------------------------------------------------------------------------------------------------
PhotonCiCMod::PhotonCiCMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fPhotonBranchName  (Names::gkPhotonBrn),
  fGoodPhotonsName   (ModNames::gkGoodPhotonsName),
  fTrackBranchName   (Names::gkTrackBrn),
  fPileUpDenName     (Names::gkPileupEnergyDensityBrn),
  fElectronName      ("Electrons"),
  fPhotonPtMin(20.0),
  fApplySpikeRemoval(kFALSE),
  fAbsEtaMax(999.99),
  fPhotons(0),
  fTracks(0),
  fPileUpDen(0),
  fElectrons(0),
  
  fPVName("PrimaryVertexes"),
  fPV(0),
  fPVFromBranch(kTRUE)

{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void PhotonCiCMod::Process()
{
  // Process entries of the tree. 

  LoadEventObject(fPhotonBranchName,   fPhotons);

  PhotonOArr *GoodPhotons = new PhotonOArr;
  GoodPhotons->SetName(fGoodPhotonsName);
  
  Double_t _tRho = -1.;
  if (fPhotons->GetEntries()>0) {
    LoadEventObject(fTrackBranchName,    fTracks);
    LoadEventObject(fPileUpDenName,      fPileUpDen);
    LoadEventObject(fElectronName,       fElectrons);
    LoadEventObject(fPVName,             fPV);    

    if(fPileUpDen->GetEntries() > 0)
      _tRho = (Double_t) fPileUpDen->At(0)->Rho();    
  
  } else {
    AddObjThisEvt(GoodPhotons);
    return;
  }

  int numVertices = fPV->GetEntries();
  if(numVertices == 0) {
    AddObjThisEvt(GoodPhotons);
    return;
  }    
  
  PhotonOArr* preselPh  = new PhotonOArr;
  PhotonOArr* failingPh = new PhotonOArr;
  
  const Vertex** vtxRanked = new const Vertex*[numVertices];
  
  // 1. we do the pre-selection; but keep the non-passing photons in a secont container...
  for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {    
    const Photon *ph = fPhotons->At(i);
    if(ph->AbsEta() > 2.5) {
      failingPh->Add(ph);
      continue;
    }
    if(ph->Et() < 30.) {
      failingPh->Add(ph);
      continue;
    }
    if(ph->HadOverEm() > 0.15) {
      failingPh->Add(ph);
      continue;
    }
    if(ph->AbsEta() < 1.5) {
      if(ph->CoviEtaiEta() > 0.013) {
        failingPh->Add(ph);
        continue;
      }
      if(ph->EcalRecHitIsoDr03() > 10.) {
        failingPh->Add(ph);
        continue;
      }
    } else {
      if(ph->CoviEtaiEta() > 0.03) {
        failingPh->Add(ph);
        continue;
      }
      if(ph->EcalRecHitIsoDr03() > 10.) {
        failingPh->Add(ph);
        continue;
      }
    }
    preselPh->Add(ph);
  }

  // sort both by pt
  preselPh->Sort();
  failingPh->Sort();
  
  // if we don't have two photons in total, nothing to do...
  if( (preselPh->GetEntries() + failingPh->GetEntries() ) < 2) {
    AddObjThisEvt(GoodPhotons);
    delete preselPh;
    delete failingPh;
    return;
  }

  // now we assign the two hardest photons to be the ones to select the Vertex...
  double ptggV[2];
  ptggV[0] = 0.;
  ptggV[1] = 0.;  
  for(unsigned int iPh =0; iPh < preselPh->GetEntries() && iPh < 2 ; ++iPh) {
    ptggV[0] += preselPh->At(iPh)->Px();
    ptggV[1] += preselPh->At(iPh)->Px();
  }
  for(unsigned int iPh =preselPh->GetEntries(); iPh < (preselPh->GetEntries()+failingPh->GetEntries()) && iPh < 2 ; ++iPh) {
    ptggV[0] += failingPh->At(iPh)->Px();
    ptggV[1] += failingPh->At(iPh)->Px();
  }  
  double ptgg = TMath::Sqrt(ptggV[0]*ptggV[0]+ptggV[1]+ptggV[1]);
  
  // loop over all vertices and assigne the ranks
  int* ptbal_rank  = new int[fPV->GetEntries()];
  int* ptasym_rank = new int[fPV->GetEntries()];
  int* total_rank  = new int[fPV->GetEntries()];
  double* ptbal = new double[fPV->GetEntries()];
  double* ptasym = new double[fPV->GetEntries()];

  for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
    const Vertex* _tVtx = fPV->At(iVtx);
    ptbal [iVtx] = 0.0;
    ptasym[iVtx] = 0.0;
    ptbal_rank [iVtx] = 1;
    ptasym_rank[iVtx] = 1;
    double ptvtx = 0.;
    for(unsigned int iTrk = 0; iTrk < _tVtx->NTracks(); ++iTrk) {
      const Track* _tTrk = _tVtx->Trk(iTrk);
      ptvtx += _tTrk->Pt();
      ptbal[iVtx]+= (_tTrk->Px()*ptggV[0])+(_tTrk->Py()*ptggV[1]);
    }
    ptbal [iVtx] = -ptbal[iVtx]/ptgg;
    ptasym[iVtx] = (ptvtx - ptgg)*(ptvtx+ptgg);
    for(int cVtx =0; cVtx < iVtx; ++cVtx) {
      if(ptbal [iVtx] > ptbal [cVtx])
        ptbal_rank [cVtx]++;
      else
        ptbal_rank [iVtx]++;
      if(ptasym [iVtx] > ptasym [cVtx])
        ptasym_rank [cVtx]++;
      else
        ptasym_rank [iVtx]++;
    }
  }

  // compute the total rank
  for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
    ptbal_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);  
    total_rank [iVtx] = 0;
    for(int cVtx =0; cVtx < iVtx; ++cVtx) {
      if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
        total_rank[iVtx]++;
      else
        total_rank[cVtx]++;
    }
  }

  // sort the Vtx by total rank
  for(int iVtx = 0; iVtx < numVertices; ++iVtx)
    vtxRanked[total_rank[iVtx]] = fPV->At(iVtx);
  
  // pick the higest ranked vertex
  const BaseVertex* theVtx = vtxRanked[0];
  
  // delete all utility arrays
  delete ptbal_rank  ;
  delete ptasym_rank ;
  delete total_rank  ;
  delete ptbal       ;
  delete ptasym      ;

  // fix all the photon kinematics using the new Vtx
  PhotonOArr *fixedPhotons = new PhotonOArr;
  for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {    
    const Photon *ph = fPhotons->At(i);
    // first copy this guy
    Photon* fPh = new Photon(*ph);
    const SuperCluster* sc = ph->SCluster();
    // compute the new Photon momentum (assuming zero mass)
    ThreeVectorC tPh_3mom  = sc->Point() - theVtx->Position();
    tPh_3mom = tPh_3mom/tPh_3mom.R();
    fPh->SetMom(sc->Energy()*tPh_3mom.X(), sc->Energy()*tPh_3mom.Y(), sc->Energy()*tPh_3mom.Z(), sc->Energy());
    fixedPhotons->Add(fPh);
  }
  
  // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
  float cic4_allcuts_temp_sublead[] = { 
    3.8,         2.2,         1.77,        1.29,
    11.7,        3.4,         3.9,         1.84,
    3.5,         2.2,         2.3,         1.45,
    0.0106,      0.0097,      0.028,       0.027,
    0.082,       0.062,       0.065,       0.048,
    0.94,        0.36,        0.94,        0.32,
    1.,          0.062,       0.97,        0.97,
    1.5,         1.5,         1.5,         1.5 };  // the last line is POixelmatchVeto and un-used
  
  for (UInt_t i=0; i<fixedPhotons->GetEntries(); ++i) {    
    const Photon *ph = fixedPhotons->At(i);        

    // cut on Et instead of Pt???
    if ( ph->Pt() <= fPhotonPtMin   ) continue;
    if ( ph->AbsEta() >= fAbsEtaMax ) continue;
    
    if (  ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)  ) 
      continue;
    
    Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
    
    Bool_t passSpikeRemovalFilter = kTRUE;
    
    if (ph->SCluster() && ph->SCluster()->Seed()) {
      if(ph->SCluster()->Seed()->Energy() > 5.0 && 
         ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
        passSpikeRemovalFilter = kFALSE;
    }
    
    if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
    
    // compute all relevant observables first
    double ecalIso = ph->EcalRecHitIsoDr03();
    double hcalIso = ph->HcalTowerSumEtDr03();
    double trackIso1 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
    double trackIso2 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, fPV);
    double trackIso3 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
    
    double combIso1 = ecalIso+hcalIso+trackIso1 - 0.17*_tRho;
    double combIso2 = ecalIso+hcalIso+trackIso2 - 0.52*_tRho;
    
    double tIso1 = (combIso1) *50./ph->Et();
    double tIso2 = (combIso2) *50./ph->Et();
    double tIso3 = (trackIso3)*50./ph->Et();    
    
    double covIEtaIEta  =ph->CoviEtaiEta();
    double HoE = ph->HadOverEm();
    
    double R9 = ph->R9();
    
    double dRTrack = PhotonTools::ElectronVetoCiC(ph, fElectrons);

    // check which category it is ...
    int _tCat = 1;
    if ( !isbarrel ) _tCat = 3;
    if ( R9 < 0.94 ) _tCat++;
    
    if( tIso1                       < cic4_allcuts_temp_sublead[_tCat-1+0*4]  && 
        tIso2                       < cic4_allcuts_temp_sublead[_tCat-1+1*4]  && 
        tIso3                       < cic4_allcuts_temp_sublead[_tCat-1+2*4]  && 
        covIEtaIEta                 < cic4_allcuts_temp_sublead[_tCat-1+3*4]  && 
        HoE                         < cic4_allcuts_temp_sublead[_tCat-1+4*4]  && 
        R9                          > cic4_allcuts_temp_sublead[_tCat-1+5*4]  && 
        dRTrack*100.                > cic4_allcuts_temp_sublead[_tCat-1+6*4]  && 
        isbarrel ) GoodPhotons->Add(ph);
  }

  // sort according to pt
  GoodPhotons->Sort();
  
  // add to event for other modules to use
  AddObjThisEvt(GoodPhotons);  
}

//--------------------------------------------------------------------------------------------------
void PhotonCiCMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we just request the photon collection branch.

  ReqEventObject(fPhotonBranchName,   fPhotons,   kTRUE);
  ReqEventObject(fTrackBranchName,    fTracks,    kTRUE);
  ReqEventObject(fElectronName,       fElectrons, kTRUE);  
  ReqEventObject(fPileUpDenName,      fPileUpDen, kTRUE);
  ReqEventObject(fPVName,             fPV,        fPVFromBranch);

}

Revision:	1.3
Committed:	Wed Jun 15 10:46:51 2011 UTC (13 years, 10 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_023, Mit_022a, Mit_022
Changes since 1.2:	+175 -92 lines
Log Message:	updated CiC selection
#	User	Rev	Content
1	fabstoec	1.3	// $Id: PhotonCiCMod.cc,v 1.2 2011/06/02 14:18:07 bendavid Exp $
2	fabstoec	1.1
3			#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
4			#include "MitAna/DataTree/interface/PhotonCol.h"
5			#include "MitPhysics/Init/interface/ModNames.h"
6			#include "MitPhysics/Utils/interface/IsolationTools.h"
7			#include "MitPhysics/Utils/interface/PhotonTools.h"
8
9			using namespace mithep;
10
11			ClassImp(mithep::PhotonCiCMod)
12
13			//--------------------------------------------------------------------------------------------------
14			PhotonCiCMod::PhotonCiCMod(const char name, const char title) :
15			BaseMod(name,title),
16			fPhotonBranchName (Names::gkPhotonBrn),
17			fGoodPhotonsName (ModNames::gkGoodPhotonsName),
18			fTrackBranchName (Names::gkTrackBrn),
19			fPileUpDenName (Names::gkPileupEnergyDensityBrn),
20			fElectronName ("Electrons"),
21	bendavid	1.2	fPhotonPtMin(20.0),
22	fabstoec	1.1	fApplySpikeRemoval(kFALSE),
23			fAbsEtaMax(999.99),
24			fPhotons(0),
25			fTracks(0),
26			fPileUpDen(0),
27			fElectrons(0),
28
29			fPVName("PrimaryVertexes"),
30			fPV(0),
31			fPVFromBranch(kTRUE)
32
33			{
34			// Constructor.
35			}
36
37			//--------------------------------------------------------------------------------------------------
38			void PhotonCiCMod::Process()
39			{
40			// Process entries of the tree.
41
42			LoadEventObject(fPhotonBranchName, fPhotons);
43	bendavid	1.2
44			PhotonOArr *GoodPhotons = new PhotonOArr;
45			GoodPhotons->SetName(fGoodPhotonsName);
46
47	fabstoec	1.1	Double_t _tRho = -1.;
48			if (fPhotons->GetEntries()>0) {
49			LoadEventObject(fTrackBranchName, fTracks);
50			LoadEventObject(fPileUpDenName, fPileUpDen);
51			LoadEventObject(fElectronName, fElectrons);
52			LoadEventObject(fPVName, fPV);
53
54			if(fPileUpDen->GetEntries() > 0)
55			_tRho = (Double_t) fPileUpDen->At(0)->Rho();
56
57	bendavid	1.2	} else {
58			AddObjThisEvt(GoodPhotons);
59			return;
60			}
61	fabstoec	1.3
62			int numVertices = fPV->GetEntries();
63			if(numVertices == 0) {
64			AddObjThisEvt(GoodPhotons);
65			return;
66			}
67
68			PhotonOArr* preselPh = new PhotonOArr;
69			PhotonOArr* failingPh = new PhotonOArr;
70
71			const Vertex** vtxRanked = new const Vertex*[numVertices];
72	fabstoec	1.1
73	fabstoec	1.3	// 1. we do the pre-selection; but keep the non-passing photons in a secont container...
74			for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
75			const Photon *ph = fPhotons->At(i);
76			if(ph->AbsEta() > 2.5) {
77			failingPh->Add(ph);
78			continue;
79			}
80			if(ph->Et() < 30.) {
81			failingPh->Add(ph);
82			continue;
83			}
84			if(ph->HadOverEm() > 0.15) {
85			failingPh->Add(ph);
86			continue;
87			}
88			if(ph->AbsEta() < 1.5) {
89			if(ph->CoviEtaiEta() > 0.013) {
90			failingPh->Add(ph);
91			continue;
92			}
93			if(ph->EcalRecHitIsoDr03() > 10.) {
94			failingPh->Add(ph);
95			continue;
96			}
97			} else {
98			if(ph->CoviEtaiEta() > 0.03) {
99			failingPh->Add(ph);
100			continue;
101			}
102			if(ph->EcalRecHitIsoDr03() > 10.) {
103			failingPh->Add(ph);
104			continue;
105			}
106			}
107			preselPh->Add(ph);
108			}
109	fabstoec	1.1
110	fabstoec	1.3	// sort both by pt
111			preselPh->Sort();
112			failingPh->Sort();
113
114			// if we don't have two photons in total, nothing to do...
115			if( (preselPh->GetEntries() + failingPh->GetEntries() ) < 2) {
116			AddObjThisEvt(GoodPhotons);
117			delete preselPh;
118			delete failingPh;
119			return;
120			}
121	fabstoec	1.1
122	fabstoec	1.3	// now we assign the two hardest photons to be the ones to select the Vertex...
123			double ptggV[2];
124			ptggV[0] = 0.;
125			ptggV[1] = 0.;
126			for(unsigned int iPh =0; iPh < preselPh->GetEntries() && iPh < 2 ; ++iPh) {
127			ptggV[0] += preselPh->At(iPh)->Px();
128			ptggV[1] += preselPh->At(iPh)->Px();
129			}
130			for(unsigned int iPh =preselPh->GetEntries(); iPh < (preselPh->GetEntries()+failingPh->GetEntries()) && iPh < 2 ; ++iPh) {
131			ptggV[0] += failingPh->At(iPh)->Px();
132			ptggV[1] += failingPh->At(iPh)->Px();
133			}
134			double ptgg = TMath::Sqrt(ptggV[0]*ptggV[0]+ptggV[1]+ptggV[1]);
135
136			// loop over all vertices and assigne the ranks
137			int* ptbal_rank = new int[fPV->GetEntries()];
138			int* ptasym_rank = new int[fPV->GetEntries()];
139			int* total_rank = new int[fPV->GetEntries()];
140			double* ptbal = new double[fPV->GetEntries()];
141			double* ptasym = new double[fPV->GetEntries()];
142
143			for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
144			const Vertex* _tVtx = fPV->At(iVtx);
145			ptbal [iVtx] = 0.0;
146			ptasym[iVtx] = 0.0;
147			ptbal_rank [iVtx] = 1;
148			ptasym_rank[iVtx] = 1;
149			double ptvtx = 0.;
150			for(unsigned int iTrk = 0; iTrk < _tVtx->NTracks(); ++iTrk) {
151			const Track* _tTrk = _tVtx->Trk(iTrk);
152			ptvtx += _tTrk->Pt();
153			ptbal[iVtx]+= (_tTrk->Px()ptggV[0])+(_tTrk->Py()ptggV[1]);
154			}
155			ptbal [iVtx] = -ptbal[iVtx]/ptgg;
156			ptasym[iVtx] = (ptvtx - ptgg)*(ptvtx+ptgg);
157			for(int cVtx =0; cVtx < iVtx; ++cVtx) {
158			if(ptbal [iVtx] > ptbal [cVtx])
159			ptbal_rank [cVtx]++;
160			else
161			ptbal_rank [iVtx]++;
162			if(ptasym [iVtx] > ptasym [cVtx])
163			ptasym_rank [cVtx]++;
164			else
165			ptasym_rank [iVtx]++;
166			}
167			}
168	fabstoec	1.1
169	fabstoec	1.3	// compute the total rank
170			for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
171			ptbal_rank [iVtx] = ptbal_rank [iVtx]ptasym_rank [iVtx](iVtx+1);
172			total_rank [iVtx] = 0;
173			for(int cVtx =0; cVtx < iVtx; ++cVtx) {
174			if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
175			total_rank[iVtx]++;
176			else
177			total_rank[cVtx]++;
178			}
179	bendavid	1.2	}
180	fabstoec	1.1
181	fabstoec	1.3	// sort the Vtx by total rank
182			for(int iVtx = 0; iVtx < numVertices; ++iVtx)
183			vtxRanked[total_rank[iVtx]] = fPV->At(iVtx);
184
185			// pick the higest ranked vertex
186			const BaseVertex* theVtx = vtxRanked[0];
187
188			// delete all utility arrays
189			delete ptbal_rank ;
190			delete ptasym_rank ;
191			delete total_rank ;
192			delete ptbal ;
193			delete ptasym ;
194	fabstoec	1.1
195	fabstoec	1.3	// fix all the photon kinematics using the new Vtx
196			PhotonOArr *fixedPhotons = new PhotonOArr;
197	fabstoec	1.1	for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
198	fabstoec	1.3	const Photon *ph = fPhotons->At(i);
199			// first copy this guy
200			Photon* fPh = new Photon(*ph);
201			const SuperCluster* sc = ph->SCluster();
202			// compute the new Photon momentum (assuming zero mass)
203			ThreeVectorC tPh_3mom = sc->Point() - theVtx->Position();
204			tPh_3mom = tPh_3mom/tPh_3mom.R();
205			fPh->SetMom(sc->Energy()tPh_3mom.X(), sc->Energy()tPh_3mom.Y(), sc->Energy()*tPh_3mom.Z(), sc->Energy());
206			fixedPhotons->Add(fPh);
207			}
208
209			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
210			float cic4_allcuts_temp_sublead[] = {
211			3.8, 2.2, 1.77, 1.29,
212			11.7, 3.4, 3.9, 1.84,
213			3.5, 2.2, 2.3, 1.45,
214			0.0106, 0.0097, 0.028, 0.027,
215			0.082, 0.062, 0.065, 0.048,
216			0.94, 0.36, 0.94, 0.32,
217			1., 0.062, 0.97, 0.97,
218			1.5, 1.5, 1.5, 1.5 }; // the last line is POixelmatchVeto and un-used
219
220			for (UInt_t i=0; i<fixedPhotons->GetEntries(); ++i) {
221			const Photon *ph = fixedPhotons->At(i);
222	fabstoec	1.1
223	fabstoec	1.3	// cut on Et instead of Pt???
224			if ( ph->Pt() <= fPhotonPtMin ) continue;
225			if ( ph->AbsEta() >= fAbsEtaMax ) continue;
226	fabstoec	1.1
227	bendavid	1.2	if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566) )
228			continue;
229
230	fabstoec	1.1	Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
231
232			Bool_t passSpikeRemovalFilter = kTRUE;
233
234			if (ph->SCluster() && ph->SCluster()->Seed()) {
235			if(ph->SCluster()->Seed()->Energy() > 5.0 &&
236			ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
237			passSpikeRemovalFilter = kFALSE;
238			}
239
240			if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
241
242			// compute all relevant observables first
243			double ecalIso = ph->EcalRecHitIsoDr03();
244			double hcalIso = ph->HcalTowerSumEtDr03();
245			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
246			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, fPV);
247			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
248	fabstoec	1.3
249	fabstoec	1.1	double combIso1 = ecalIso+hcalIso+trackIso1 - 0.17*_tRho;
250			double combIso2 = ecalIso+hcalIso+trackIso2 - 0.52*_tRho;
251
252	fabstoec	1.3	double tIso1 = (combIso1) *50./ph->Et();
253			double tIso2 = (combIso2) *50./ph->Et();
254	fabstoec	1.1	double tIso3 = (trackIso3)*50./ph->Et();
255
256			double covIEtaIEta =ph->CoviEtaiEta();
257			double HoE = ph->HadOverEm();
258
259			double R9 = ph->R9();
260
261			double dRTrack = PhotonTools::ElectronVetoCiC(ph, fElectrons);
262
263	fabstoec	1.3	// check which category it is ...
264			int _tCat = 1;
265			if ( !isbarrel ) _tCat = 3;
266			if ( R9 < 0.94 ) _tCat++;
267
268			if( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4] &&
269			tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4] &&
270			tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4] &&
271			covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4] &&
272			HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4] &&
273			R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4] &&
274			dRTrack100. > cic4_allcuts_temp_sublead[_tCat-1+64] &&
275			isbarrel ) GoodPhotons->Add(ph);
276	fabstoec	1.1	}
277
278			// sort according to pt
279			GoodPhotons->Sort();
280	fabstoec	1.3
281	fabstoec	1.1	// add to event for other modules to use
282			AddObjThisEvt(GoodPhotons);
283			}
284
285			//--------------------------------------------------------------------------------------------------
286			void PhotonCiCMod::SlaveBegin()
287			{
288			// Run startup code on the computer (slave) doing the actual analysis. Here,
289			// we just request the photon collection branch.
290
291			ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
292			ReqEventObject(fTrackBranchName, fTracks, kTRUE);
293			ReqEventObject(fElectronName, fElectrons, kTRUE);
294			ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
295			ReqEventObject(fPVName, fPV, fPVFromBranch);
296
297			}
298	fabstoec	1.3