ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc

Comparing UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc (file contents):
Revision 1.2 by bendavid, Thu Jun 2 14:18:07 2011 UTC vs.
Revision 1.7 by fabstoec, Mon Jul 4 13:48:50 2011 UTC

#	Line 1 | Line 1
1	–	// $Id$
2	–
1		#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
2		#include "MitAna/DataTree/interface/PhotonCol.h"
3		#include "MitPhysics/Init/interface/ModNames.h"
4		#include "MitPhysics/Utils/interface/IsolationTools.h"
5		#include "MitPhysics/Utils/interface/PhotonTools.h"
6	+	#include <TNtuple.h>
7	+	#include <TRandom3.h>
8
9		using namespace mithep;
10
#	Line 26 | Line 26 | PhotonCiCMod::PhotonCiCMod(const char *n
26		fPileUpDen(0),
27		fElectrons(0),
28
29	<	fPVName("PrimaryVertexes"),
29	>	fPVName(Names::gkPVBeamSpotBrn),
30		fPV(0),
31	<	fPVFromBranch(kTRUE)
31	>	fPVFromBranch(kTRUE),
32	>
33	>	fConversions(0),
34	>	fConversionName(Names::gkMvfConversionBrn),
35	>
36	>	fBeamspot(0),
37	>
38	>	// May10 ReReco
39	>	//   fDataEnCorr_EB_hR9(-0.0047),
40	>	//   fDataEnCorr_EB_lR9(0.0014),
41	>	//   fDataEnCorr_EE_hR9(0.0076),
42	>	//   fDataEnCorr_EE_lR9(0.0008),
43	>
44	>	// prompt Reload
45	>	fDataEnCorr_EB_hR9(0.0001),
46	>	fDataEnCorr_EB_lR9(0.0052),
47	>	fDataEnCorr_EE_hR9(0.0428),
48	>	fDataEnCorr_EE_lR9(0.0180),
49	>
50	>
51	>	fMCSmear_EB_hR9(0.0089),
52	>	fMCSmear_EB_lR9(0.0199),
53	>	fMCSmear_EE_hR9(0.0409),
54	>	fMCSmear_EE_lR9(0.0246),
55	>
56	>	fIsData(false),
57	>
58	>	rng(new TRandom3()),
59	>
60	>	fMCParticleName(Names::gkMCPartBrn),
61	>	fMCParticles(0),
62	>	fPileUpName         ("PileupInfo"),
63	>	fPileUp             (0)
64
65		{
66		// Constructor.
67		}
68
69	+	PhotonCiCMod::~PhotonCiCMod(){
70	+	if(rng) delete rng;
71	+	}
72	+
73		//--------------------------------------------------------------------------------------------------
74		void PhotonCiCMod::Process()
75		{
#	Line 43 | Line 79 | void PhotonCiCMod::Process()
79
80		PhotonOArr *GoodPhotons = new PhotonOArr;
81		GoodPhotons->SetName(fGoodPhotonsName);
82	<
47	<
48	<	Double_t _tRho = -1.;
49	<	if (fPhotons->GetEntries()>0) {
50	<	LoadEventObject(fTrackBranchName,    fTracks);
51	<	LoadEventObject(fPileUpDenName,      fPileUpDen);
52	<	LoadEventObject(fElectronName,       fElectrons);
53	<	LoadEventObject(fPVName,             fPV);
82	>	GoodPhotons->SetOwner(kTRUE);
83
84	<	if(fPileUpDen->GetEntries() > 0)
85	<	_tRho = (Double_t) fPileUpDen->At(0)->Rho();
86	<
87	<	} else {
88	<	AddObjThisEvt(GoodPhotons);
89	<	return;
84	>	Double_t _tRho = -1.;
85	>	LoadEventObject(fTrackBranchName,    fTracks);
86	>	LoadEventObject(fPileUpDenName,      fPileUpDen);
87	>	LoadEventObject(fElectronName,       fElectrons);
88	>	LoadEventObject(fPVName,             fPV);
89	>	LoadEventObject(fConversionName,     fConversions);
90	>	LoadBranch("BeamSpot");
91	>
92	>	if( !fIsData ) {
93	>	LoadBranch(fMCParticleName);
94	>	LoadBranch(fPileUpName);
95		}
62	–
63	–
96
97	<	std::vector<const Photon*> phCat1;
98	<	std::vector<const Photon*> phCat2;
99	<	std::vector<const Photon*> phCat3;
68	<	std::vector<const Photon*> phCat4;
97	>	Float_t numPU = -1.;
98	>	if( !fIsData )
99	>	numPU = (Float_t) fPileUp->At(0)->GetPU_NumInteractions();
100
101	+	if(fPileUpDen->GetEntries() > 0)
102	+	_tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
103
104	<	const BaseVertex* theVtx = NULL;
105	<	if(fPV->GetEntries() > 0)
106	<	theVtx = fPV->At(0);
107	<	else {
108	<	AddObjThisEvt(GoodPhotons);
109	<	return;
110	<	}
104	>	bool doVtxSelection = true;
105	>	bool doMCSmear      = false;
106	>
107	>	const EventHeader* evtHead = this->GetEventHeader();
108	>
109	>	unsigned int evtNum = evtHead->EvtNum();
110	>	Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
111	>	Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000)  );
112	>
113	>	//double evtNumTest = (int) ( ( (double) _evtNum1 )*10000. + (double) _evtNum2 );
114	>
115	>	Float_t _runNum  = (Float_t) evtHead->RunNum();
116	>	Float_t _lumiSec = (Float_t) evtHead->LumiSec();
117
79	–	float cic4_allcuts_temp_sublead[] = {
80	–	3.77459,     2.18305,     1.76811,     1.30029,
81	–	11.5519,     3.48306,     3.87394,     1.89038,
82	–	3.47103,      2.1822,     2.26931,     1.43769,
83	–	0.0105631,  0.00974116,   0.0282572,   0.0265096,
84	–	0.0834648,   0.0821447,   0.0648449,   0.0476437,
85	–	0.94,    0.355935,    0.94,    0.316358,
86	–	98.9979,     1.94497,     96.2292,     96.2294,
87	–	1.5,         1.5,         1.5,         1.5 };   // THIS IS ????
118
119	+	//unsigned int numVertices = fPV->GetEntries();
120	+
121	+	const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
122	+
123	+	PhotonOArr* preselPh  = new PhotonOArr;
124	+
125	+	// 1. we do the pre-selection; but keep the non-passing photons in a secont container...
126		for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
127	<	const Photon *ph = fPhotons->At(i);
127	>	const Photon *ph = fPhotons->At(i);
128	>
129	>	if(ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
130	>	if(ph->Et() < 20.)     continue;
131	>	if(ph->HadOverEm() > 0.15) continue;
132	>	if(ph->AbsEta() < 1.5) {
133	>	if(ph->CoviEtaiEta() > 0.013) continue;
134	>	//if(ph->EcalRecHitIsoDr03() > 10.) continue;
135	>	} else {
136	>	if(ph->CoviEtaiEta() > 0.03) continue;
137	>	//if(ph->EcalRecHitIsoDr03() > 10.) continue;
138	>	}
139
92	–	if (ph->Pt() <= fPhotonPtMin) continue;
93	–	if (ph->AbsEta() >= fAbsEtaMax) continue;
94	–
95	–	if (  ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)  )
96	–	continue;
97	–
98	–	Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
99	–
140		Bool_t passSpikeRemovalFilter = kTRUE;
141	<
141	>
142		if (ph->SCluster() && ph->SCluster()->Seed()) {
143		if(ph->SCluster()->Seed()->Energy() > 5.0 &&
144	<	ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
145	<	passSpikeRemovalFilter = kFALSE;
146	<	}
107	<
144	>	ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
145	>	passSpikeRemovalFilter = kFALSE;
146	>	}
147		if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
148	+
149	+	preselPh->Add(ph);
150	+	}
151	+
152	+	// sort both by pt... again ;)
153	+	preselPh->Sort();
154	+
155	+	unsigned int numPairs = 0;
156	+	if( preselPh->GetEntries() > 0)
157	+	numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
158	+
159	+	// create all possible pairs of pre-selected photons
160	+
161	+	std::vector<unsigned int> idxFst;
162	+	std::vector<unsigned int> idxSec;
163	+	std::vector<bool> pairPasses;
164	+
165	+	if(numPairs > 0) {
166	+	for(unsigned int iFst = 0; iFst <preselPh->GetEntries() - 1; ++iFst) {
167	+	for(unsigned int iSec = iFst + 1; iSec <preselPh->GetEntries(); ++iSec) {
168	+	idxFst.push_back(iFst);
169	+	idxSec.push_back(iSec);
170	+	pairPasses.push_back(true);
171	+	}
172	+	}
173	+	}
174	+
175	+	// array to store the index of 'chosen Vtx' for each pair
176	+	const Vertex** theVtx = new const Vertex*[numPairs];
177	+
178	+	// arays to store the Vtx 'fixed' photons
179	+	Photon** fixPhFst = new Photon*[numPairs];
180	+	Photon** fixPhSec = new Photon*[numPairs];
181	+
182	+	// sotr pair-indices for pairs passing the selection
183	+	std::vector<unsigned int> passPairs;
184	+	passPairs.resize(0);
185	+
186	+	unsigned int theChosenVtx = 0;
187	+
188	+	float kinPh1[20];
189	+	float kinPh2[20];
190	+
191	+	for(int i=0; i<10;  ++i) {
192	+	kinPh1[i] =-99.;
193	+	kinPh2[i] =-99.;
194	+	}
195	+
196	+	float ptBefore1 = -99.;
197	+	float ptBefore2 = -99.;
198	+
199	+	bool print = false;
200	+	if(evtNum == 17031) {
201	+	std::cout<<" ------------------------------------------- "<<std::endl;
202	+	std::cout<<"   printing info for event #"<<evtNum<<std::endl;
203	+	print = true;
204	+	}
205	+
206	+	for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
207
208	<	// compute all relevant observables first
209	<	double ecalIso = ph->EcalRecHitIsoDr03();
210	<	double hcalIso = ph->HcalTowerSumEtDr03();
211	<	double trackIso1 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
212	<	double trackIso2 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, fPV);
213	<	double trackIso3 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
214	<
215	<	double combIso1 = ecalIso+hcalIso+trackIso1 - 0.17*_tRho;
216	<	double combIso2 = ecalIso+hcalIso+trackIso2 - 0.52*_tRho;
217	<
218	<	double tIso1 = (combIso1)*50./ph->Et();
219	<	double tIso2 = (combIso2)*50./ph->Et();
220	<	double tIso3 = (trackIso3)*50./ph->Et();
221	<
222	<	double covIEtaIEta  =ph->CoviEtaiEta();
223	<	double HoE = ph->HadOverEm();
224	<
225	<	double R9 = ph->R9();
226	<
227	<	double dRTrack = PhotonTools::ElectronVetoCiC(ph, fElectrons);
228	<
229	<	// stupid hard coded, will update later...
230	<	if( tIso1                       < cic4_allcuts_temp_sublead[0+0*4]  &&
231	<	tIso2                       < cic4_allcuts_temp_sublead[0+1*4]  &&
232	<	tIso3                       < cic4_allcuts_temp_sublead[0+2*4]  &&
233	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[0+3*4]  &&
234	<	HoE                         < cic4_allcuts_temp_sublead[0+4*4]  &&
235	<	R9                          > cic4_allcuts_temp_sublead[0+5*4]  &&
236	<	dRTrack*100.                > cic4_allcuts_temp_sublead[0+6*4]  &&
237	<	isbarrel ) phCat1.push_back(ph);
238	<
239	<	if( tIso1                       < cic4_allcuts_temp_sublead[1+0*4]  &&
240	<	tIso2                       < cic4_allcuts_temp_sublead[1+1*4]  &&
241	<	tIso3                       < cic4_allcuts_temp_sublead[1+2*4]  &&
242	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[1+3*4]  &&
243	<	HoE                         < cic4_allcuts_temp_sublead[1+4*4]  &&
244	<	R9                          > cic4_allcuts_temp_sublead[1+5*4]  &&
245	<	dRTrack*100.                > cic4_allcuts_temp_sublead[1+6*4]  &&
246	<	isbarrel ) phCat2.push_back(ph);
247	<
248	<	if( tIso1                       < cic4_allcuts_temp_sublead[2+0*4]  &&
249	<	tIso2                       < cic4_allcuts_temp_sublead[2+1*4]  &&
250	<	tIso3                       < cic4_allcuts_temp_sublead[2+2*4]  &&
251	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[2+3*4]  &&
252	<	HoE                         < cic4_allcuts_temp_sublead[2+4*4]  &&
253	<	R9                          > cic4_allcuts_temp_sublead[2+5*4]  &&
254	<	dRTrack*100.                > cic4_allcuts_temp_sublead[2+6*4]     )  phCat3.push_back(ph);
255	<
256	<	if( tIso1                       < cic4_allcuts_temp_sublead[3+0*4]  &&
257	<	tIso2                       < cic4_allcuts_temp_sublead[3+1*4]  &&
258	<	tIso3                       < cic4_allcuts_temp_sublead[3+2*4]  &&
259	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[3+3*4]  &&
260	<	HoE                         < cic4_allcuts_temp_sublead[3+4*4]  &&
261	<	R9                          > cic4_allcuts_temp_sublead[3+5*4]  &&
262	<	dRTrack*100.                > cic4_allcuts_temp_sublead[3+6*4]     )  phCat4.push_back(ph);
263	<
264	<	}
265	<
266	<	// check which category has two vgalid passes
267	<	bool foundCat=false;
268	<	if(phCat1.size() > 1) {
269	<	// pick the first two guys, assuming they are sorted...
270	<	GoodPhotons->Add(phCat1[0]);
271	<	GoodPhotons->Add(phCat1[1]);
174	<	} else if(phCat2.size() > 1) {
175	<	GoodPhotons->Add(phCat2[0]);
176	<	GoodPhotons->Add(phCat2[1]);
177	<	} else {
178	<	if(phCat3.size() > 1) {
179	<	// one must be a Endcap photons... (stupid, but that's how it is in CiC)
180	<	if(phCat3[0]->SCluster()->AbsEta() > 1.5 || phCat3[1]->SCluster()->AbsEta() > 1.5) {
181	<	GoodPhotons->Add(phCat3[0]);
182	<	GoodPhotons->Add(phCat3[1]);
183	<	foundCat=true;
184	<	}
185	<	if (!foundCat) {
186	<	if(phCat4.size() > 1) {
187	<	if(phCat4[0]->SCluster()->AbsEta() > 1.5 || phCat4[1]->SCluster()->AbsEta() > 1.5) {
188	<	GoodPhotons->Add(phCat4[0]);
189	<	GoodPhotons->Add(phCat4[1]);
190	<	}
191	<	}
208	>	// copy the photons for manipulation
209	>	fixPhFst[iPair] = new Photon(*preselPh->At(idxFst[iPair]));
210	>	fixPhSec[iPair] = new Photon(*preselPh->At(idxSec[iPair]));
211	>
212	>	// if this is Data, scale the energy, if MC smear...
213	>	FourVectorM scMomFst = fixPhFst[iPair]->Mom();
214	>	FourVectorM scMomSec = fixPhSec[iPair]->Mom();
215	>	double scaleFac1 = 1.;
216	>	double scaleFac2 = 1.;
217	>	if (fIsData) {
218	>	if(fixPhFst[iPair]->IsEB())
219	>	if(fixPhFst[iPair]->R9() > 0.94)
220	>	scaleFac1 += fDataEnCorr_EB_hR9;
221	>	else
222	>	scaleFac1 += fDataEnCorr_EB_lR9;
223	>	else
224	>	if(fixPhFst[iPair]->R9() > 0.94)
225	>	scaleFac1 += fDataEnCorr_EE_hR9;
226	>	else
227	>	scaleFac1 += fDataEnCorr_EE_lR9;
228	>	if(fixPhSec[iPair]->IsEB())
229	>	if(fixPhSec[iPair]->R9() > 0.94)
230	>	scaleFac2 += fDataEnCorr_EB_hR9;
231	>	else
232	>	scaleFac2 += fDataEnCorr_EB_lR9;
233	>	else
234	>	if(fixPhSec[iPair]->R9() > 0.94)
235	>	scaleFac2 += fDataEnCorr_EE_hR9;
236	>	else
237	>	scaleFac2 += fDataEnCorr_EE_lR9;
238	>	} else {
239	>	// get the smearing for MC photons..
240	>	UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
241	>	UInt_t seed1    = seedBase + (UInt_t) fixPhFst[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhFst[iPair]->SCluster()->Eta()));
242	>	UInt_t seed2    = seedBase + (UInt_t) fixPhSec[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhSec[iPair]->SCluster()->Eta()));
243	>
244	>	double width1 = 0.;
245	>	double width2 = 0.;
246	>	if(fixPhFst[iPair]->IsEB())
247	>	if(fixPhFst[iPair]->R9() > 0.94)
248	>	width1 = fMCSmear_EB_hR9;
249	>	else
250	>	width1 = fMCSmear_EB_lR9;
251	>	else
252	>	if(fixPhFst[iPair]->R9() > 0.94)
253	>	width1 = fMCSmear_EE_hR9;
254	>	else
255	>	width1 = fMCSmear_EE_lR9;
256	>	if(fixPhSec[iPair]->IsEB())
257	>	if(fixPhSec[iPair]->R9() > 0.94)
258	>	width2 = fMCSmear_EB_hR9;
259	>	else
260	>	width2 = fMCSmear_EB_lR9;
261	>	else
262	>	if(fixPhSec[iPair]->R9() > 0.94)
263	>	width2 = fMCSmear_EE_hR9;
264	>	else
265	>	width2 = fMCSmear_EE_lR9;
266	>
267	>	if(doMCSmear) {
268	>	rng->SetSeed(seed1);
269	>	scaleFac1 = rng->Gaus(1.,width1);
270	>	rng->SetSeed(seed2);
271	>	scaleFac2 = rng->Gaus(1.,width2);
272		}
273		}
274	+
275	+	if(iPair==0) {
276	+	ptBefore1 = fixPhFst[iPair]->Pt();
277	+	ptBefore2 = fixPhSec[iPair]->Pt();
278	+	}
279	+
280	+	if(print && false) {
281	+	std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
282	+	std::cout<<"      Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
283	+	std::cout<<"          py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
284	+	std::cout<<"          pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
285	+	std::cout<<"           E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
286	+	std::cout<<"           M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
287	+	std::cout<<"      Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
288	+	std::cout<<"          py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
289	+	std::cout<<"          pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
290	+	std::cout<<"           E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
291	+	std::cout<<"           M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
292	+	}
293	+
294	+	fixPhFst[iPair]->SetMom(scaleFac1*scMomFst.X(), scaleFac1*scMomFst.Y(), scaleFac1*scMomFst.Z(), scaleFac1*scMomFst.E());
295	+	fixPhSec[iPair]->SetMom(scaleFac2*scMomSec.X(), scaleFac2*scMomSec.Y(), scaleFac2*scMomSec.Z(), scaleFac2*scMomSec.E());
296	+
297	+	if(print && false) {
298	+	std::cout<<"      SF 1 = "<<scaleFac1<<std::endl;
299	+	std::cout<<"      SF 2 = "<<scaleFac2<<std::endl;
300	+
301	+	std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
302	+	std::cout<<"      Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
303	+	std::cout<<"          py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
304	+	std::cout<<"          pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
305	+	std::cout<<"           E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
306	+	std::cout<<"           M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
307	+	std::cout<<"      Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
308	+	std::cout<<"          py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
309	+	std::cout<<"          pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
310	+	std::cout<<"           E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
311	+	std::cout<<"           M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
312	+	}
313	+
314	+	// store the vertex for this pair
315	+	if(doVtxSelection) {
316	+	unsigned int iVtx = findBestVertex(fixPhFst[iPair],fixPhSec[iPair],bsp, print);
317	+	theVtx[iPair] =  fPV->At(iVtx);
318	+	if(iPair == 0) theChosenVtx = iVtx;
319	+	} else
320	+	theVtx[iPair] =  fPV->At(0);
321	+
322	+	// fix the kinematics for both events
323	+	FourVectorM newMomFst = fixPhFst[iPair]->MomVtx(theVtx[iPair]->Position());
324	+	FourVectorM newMomSec = fixPhSec[iPair]->MomVtx(theVtx[iPair]->Position());
325	+	fixPhFst[iPair]->SetMom(newMomFst.X(), newMomFst.Y(), newMomFst.Z(), newMomFst.E());
326	+	fixPhSec[iPair]->SetMom(newMomSec.X(), newMomSec.Y(), newMomSec.Z(), newMomSec.E());
327	+
328	+
329	+	if(print && false) {
330	+	std::cout<<"         Vtx = "<<theChosenVtx<<std::endl;
331	+	std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
332	+	std::cout<<"      Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
333	+	std::cout<<"          py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
334	+	std::cout<<"          pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
335	+	std::cout<<"           E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
336	+	std::cout<<"           M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
337	+	std::cout<<"      Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
338	+	std::cout<<"          py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
339	+	std::cout<<"          pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
340	+	std::cout<<"           E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
341	+	std::cout<<"           M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
342	+	}
343	+
344	+
345	+
346	+	if(iPair != 0) {
347	+	// check if both photons pass the CiC selection
348	+	bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false);
349	+	bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false);
350	+	if( pass1 && pass2 )
351	+	passPairs.push_back(iPair);
352	+	} else {
353	+	bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false, kinPh1);
354	+	bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false, kinPh2);
355	+	if( pass1 && pass2 )
356	+	passPairs.push_back(iPair);
357	+	}
358	+	}
359	+
360	+	// loop over all passing pairs and find the one with the highest sum Et
361	+	Photon* phHard = NULL;
362	+	Photon* phSoft = NULL;
363	+	double maxSumEt = 0.;
364	+	for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
365	+	double sumEt = fixPhFst[passPairs[iPair]]->Et();
366	+	sumEt += fixPhSec[passPairs[iPair]]->Et();
367	+	if( sumEt > maxSumEt ) {
368	+	maxSumEt = sumEt;
369	+	phHard = fixPhFst[passPairs[iPair]];
370	+	phSoft = fixPhSec[passPairs[iPair]];
371	+	}
372	+	}
373	+
374	+	if(phHard && phSoft) {
375	+	GoodPhotons->AddOwned(phHard);
376	+	GoodPhotons->AddOwned(phSoft);
377		}
378
379		// sort according to pt
380		GoodPhotons->Sort();
381	<
381	>
382		// add to event for other modules to use
383	<	AddObjThisEvt(GoodPhotons);
383	>	AddObjThisEvt(GoodPhotons);
384	>
385	>	delete preselPh;
386	>
387	>
388	>	bool doFill = (phHard && phSoft);
389	>	Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
390	>	Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);
391	>
392	>	Float_t _pth = -100.;
393	>	if( !fIsData ) _pth = findHiggsPt();
394	>
395	>	Float_t fillEvent[] = { _tRho,
396	>	_pth,
397	>	numPU,
398	>	_mass,
399	>	_ptgg,
400	>	_evtNum1,
401	>	_evtNum2,
402	>	_runNum,
403	>	_lumiSec,
404	>	(float) theChosenVtx,
405	>	(float) numPairs,
406	>	kinPh1[0],
407	>	kinPh1[1],
408	>	kinPh1[2],
409	>	kinPh1[3],
410	>	kinPh1[4],
411	>	kinPh1[5],
412	>	kinPh1[6],
413	>	kinPh1[7],
414	>	kinPh1[8],
415	>	kinPh1[9],
416	>	kinPh1[10],
417	>	kinPh1[11],
418	>	kinPh1[12],
419	>	kinPh1[13],
420	>	kinPh1[14],
421	>	kinPh1[15],
422	>	kinPh1[16],
423	>	kinPh1[17],
424	>	kinPh1[18],
425	>	kinPh1[19],
426	>	kinPh2[0],
427	>	kinPh2[1],
428	>	kinPh2[2],
429	>	kinPh2[3],
430	>	kinPh2[4],
431	>	kinPh2[5],
432	>	kinPh2[6],
433	>	kinPh2[7],
434	>	kinPh2[8],
435	>	kinPh2[9],
436	>	kinPh2[10],
437	>	kinPh2[11],
438	>	kinPh2[12],
439	>	kinPh2[13],
440	>	kinPh2[14],
441	>	kinPh2[15],
442	>	kinPh2[16],
443	>	kinPh2[17],
444	>	kinPh2[18],
445	>	kinPh2[19],
446	>	ptBefore1,
447	>	ptBefore2
448	>	};
449	>
450	>	hCiCTuple->Fill(fillEvent);
451	>
452	>	return;
453	>
454		}
455
456		//--------------------------------------------------------------------------------------------------
#	Line 211 | Line 464 | void PhotonCiCMod::SlaveBegin()
464		ReqEventObject(fElectronName,       fElectrons, kTRUE);
465		ReqEventObject(fPileUpDenName,      fPileUpDen, kTRUE);
466		ReqEventObject(fPVName,             fPV,        fPVFromBranch);
467	+	ReqEventObject(fConversionName,     fConversions,kTRUE);
468	+	ReqBranch("BeamSpot",fBeamspot);
469	+
470	+	if (!fIsData) {
471	+	ReqBranch(fPileUpName, fPileUp);
472	+	ReqBranch(Names::gkMCPartBrn,fMCParticles);
473	+	}
474	+
475	+	hCiCTuple = new TNtuple("hCiCTuple","hCiCTuple","rho:higgspt:numPU:mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ivtx:npairs:ph1Iso1:ph1Iso2:ph1Iso3:ph1Cov:ph1HoE:ph1R9:ph1DR:ph1Pt:ph1Eta:ph1Phi:ph1Eiso3:ph1Eiso4:ph1Hiso4:ph1TisoA:ph1TisoW:ph1Tiso:ph1Et:ph1E:ph1Pass:ph1Cat:ph2Iso1:ph2Iso2:ph2Iso3:ph2Cov:ph2HoE:ph2R9:ph2DR:ph2Pt:ph2Eta:ph2Phi:ph2Eiso3:ph2Eiso4:ph2Hiso4:ph2TisoA:ph2TisoW:ph2Tiso:ph2Et:ph2E:ph2Pass:ph2Cat:ph1UPt:ph2UPt");
476	+
477	+	AddOutput(hCiCTuple);
478
479		}
480	+
481	+	// return the index of the bext vertex
482	+	unsigned int PhotonCiCMod::findBestVertex(Photon* ph1, Photon* ph2, const BaseVertex *bsp, bool print) {
483	+
484	+	// loop over all vertices and assigne the ranks
485	+	int* ptbal_rank  = new int[fPV->GetEntries()];
486	+	int* ptasym_rank = new int[fPV->GetEntries()];
487	+	int* total_rank  = new int[fPV->GetEntries()];
488	+	double* ptbal = new double[fPV->GetEntries()];
489	+	double* ptasym = new double[fPV->GetEntries()];
490	+
491	+	unsigned int numVertices = fPV->GetEntries();
492	+
493	+	double ptgg = 0.;    // stored for later in the conversion
494	+
495	+	if(print) {
496	+	std::cout<<" --------------------------------- "<<std::endl;
497	+	std::cout<<"   looking for Vtx for photon pair "<<std::endl;
498	+	std::cout<<"                            pt 1 = "<<ph1->Et()<<std::endl;
499	+	std::cout<<"                            pt 2 = "<<ph2->Et()<<std::endl;
500	+	std::cout<<"                    among #"<<numVertices<<" Vertices."<<std::endl;
501	+	}
502	+
503	+
504	+	for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
505	+	if(print)
506	+	std::cout<<std::endl<<"       Vertex #"<<iVtx<<std::endl;
507	+
508	+	const Vertex* tVtx = fPV->At(iVtx);
509	+	ptbal [iVtx] = 0.0;
510	+	ptasym[iVtx] = 0.0;
511	+	ptbal_rank [iVtx] = 1;
512	+	ptasym_rank[iVtx] = 1;
513	+
514	+	// compute the photon momenta with respect to this Vtx
515	+	FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
516	+	FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
517	+
518	+	FourVectorM higgsMom = newMomFst+newMomSec;
519	+
520	+	double ph1Eta = newMomFst.Eta();
521	+	double ph2Eta = newMomSec.Eta();
522	+
523	+	double ph1Phi = newMomFst.Phi();
524	+	double ph2Phi = newMomSec.Phi();
525	+
526	+	if(print && iVtx == 0 && false ) {
527	+	std::cout<<"     new momenta Et1 = "<<newMomFst.Et()<<std::endl;
528	+	std::cout<<"                 Eta = "<<newMomFst.Eta()<<std::endl;
529	+	std::cout<<"                 Phi = "<<newMomFst.Phi()<<std::endl;
530	+	std::cout<<"                  Px = "<<newMomFst.Px()<<std::endl;
531	+	std::cout<<"                  Py = "<<newMomFst.Py()<<std::endl;
532	+	std::cout<<"     new momenta Et2 = "<<newMomSec.Et()<<std::endl;
533	+	std::cout<<"                 Eta = "<<newMomSec.Eta()<<std::endl;
534	+	std::cout<<"                 Phi = "<<newMomSec.Phi()<<std::endl;
535	+	std::cout<<"                  Px = "<<newMomSec.Px()<<std::endl;
536	+	std::cout<<"                  Py = "<<newMomSec.Py()<<std::endl;
537	+	}
538	+
539	+	FourVectorM totTrkMom;
540	+	for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
541	+	const Track* tTrk = tVtx->Trk(iTrk);
542	+	//if(tTrk->Pt()<1.) continue;
543	+	double tEta = tTrk->Eta();
544	+	double tPhi = tTrk->Phi();
545	+	double dEta1 = TMath::Abs(tEta-ph1Eta);
546	+	double dEta2 = TMath::Abs(tEta-ph2Eta);
547	+	double dPhi1 = TMath::Abs(tPhi-ph1Phi);
548	+	double dPhi2 = TMath::Abs(tPhi-ph2Phi);
549	+	if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
550	+	if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
551	+
552	+	double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
553	+	double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
554	+
555	+	if( ( iVtx == 0 || iVtx == 1 ) && print && false) {
556	+	std::cout<<"  Track #"<<iTrk<<std::endl;
557	+	std::cout<<"      pt = "<<tTrk->Pt()<<std::endl;
558	+	std::cout<<"     eta = "<<tTrk->Eta()<<std::endl;
559	+	std::cout<<"     phi = "<<tTrk->Phi()<<std::endl;
560	+	std::cout<<"      px = "<<tTrk->Px()<<std::endl;
561	+	std::cout<<"      py = "<<tTrk->Py()<<std::endl;
562	+	std::cout<<"     dR1 = "<<dR1<<std::endl;
563	+	std::cout<<"     dR2 = "<<dR2<<std::endl;
564	+	}
565	+
566	+	if(dR1 < 0.05 || dR2 < 0.05) continue;
567	+
568	+	if(iTrk == 0) totTrkMom = tTrk->Mom4(0);
569	+	else totTrkMom = totTrkMom + tTrk->Mom4(0);
570	+
571	+	}
572	+
573	+	if(iVtx ==0 && print && false) {
574	+	std::cout<<"   Trk passes cuts: "<<std::endl;
575	+	std::cout<<"            px tot = "<<totTrkMom.Px()<<std::endl;
576	+	std::cout<<"            py tot = "<<totTrkMom.Py()<<std::endl;
577	+	}
578	+
579	+	double ptvtx   = totTrkMom.Pt();
580	+	if(iVtx ==0 && print && false)
581	+	std::cout<<"    Total TkPt = "<<ptvtx<<std::endl;
582	+	double pthiggs = higgsMom.Pt();
583	+	if(iVtx ==0 && print && false)
584	+	std::cout<<"    Total H Pt = "<<pthiggs<<std::endl;
585	+	if(iVtx == 0) ptgg = pthiggs;
586	+	ptbal [iVtx]  = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
587	+	ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
588	+	ptbal [iVtx]  = -ptbal[iVtx]/pthiggs;
589	+	ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
590	+
591	+	if(iVtx ==0 && print && false) {
592	+	std::cout<<"    Results: ptbal = "<<ptbal [iVtx]<<std::endl;
593	+	std::cout<<"            ptasym = "<<ptasym[iVtx]<<std::endl;
594	+	}
595	+
596	+	for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
597	+	if(ptbal [iVtx] > ptbal [cVtx])
598	+	ptbal_rank[cVtx]++;
599	+	else
600	+	ptbal_rank[iVtx]++;
601	+	if(ptasym [iVtx] > ptasym [cVtx])
602	+	ptasym_rank[cVtx]++;
603	+	else
604	+	ptasym_rank[iVtx]++;
605	+	}
606	+	}
607	+
608	+
609	+	// compute the total rank
610	+	for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
611	+	if(print) {
612	+	std::cout<<"     Vertex #"<<iVtx<<"  has rank PTB "<<ptbal_rank[iVtx]<<"   ("<<ptbal[iVtx]<<")"<<std::endl;
613	+	std::cout<<"     Vertex #"<<iVtx<<"  has rank PTSYM "<<ptasym_rank[iVtx]<<"   ("<<ptasym[iVtx]<<")"<<std::endl;
614	+	}
615	+	ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
616	+	total_rank [iVtx] = 0;
617	+	for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
618	+	if(ptasym_rank [iVtx] > ptasym_rank [cVtx])
619	+	total_rank[iVtx]++;
620	+	else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) {
621	+	if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
622	+	total_rank[iVtx]++;
623	+	else
624	+	total_rank[cVtx]++;
625	+	}
626	+	else
627	+	total_rank[cVtx]++;
628	+	}
629	+	}
630	+
631	+	if(print) std::cout<<std::endl;
632	+
633	+	unsigned int bestIdx = 0;
634	+	for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
635	+	if(total_rank[iVtx] == 0) bestIdx = iVtx;
636	+	if(print) {
637	+	std::cout<<"     Vertex #"<<iVtx<<"  has rank "<<total_rank[iVtx]<<std::endl;
638	+	}
639	+	}
640	+
641	+	delete[] ptbal_rank  ;
642	+	delete[] ptasym_rank ;
643	+	delete[] ptbal       ;
644	+	delete[] ptasym      ;
645	+
646	+	//return bestIdx;
647	+
648	+	// check if there's a conversion among the pre-selected photons
649	+	const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, fConversions, 0.1, 0.1, 0.1, print);
650	+	const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, fConversions, 0.1, 0.1, 0.1, print);
651	+
652	+	if(print) {
653	+	if (conv1) {
654	+	std::cout<<" Photon 1 has has conversion with P = "<<conv1->Prob()<<std::endl;
655	+	std::cout<<"                                Rho = "<<conv1->Position().Rho()<<std::endl;
656	+	std::cout<<"                                  Z = "<<conv1->Position().Z()<<std::endl;
657	+	}
658	+	if (conv2) {
659	+	std::cout<<" Photon 2 has has conversion with P = "<<conv2->Prob()<<std::endl;
660	+	std::cout<<"                                Rho = "<<conv2->Position().Rho()<<std::endl;
661	+	std::cout<<"                                  Z = "<<conv2->Position().Z()<<std::endl;
662	+	}
663	+	}
664	+
665	+	if( conv1 && ( conv1->Prob() < 0.0005) ) conv1 = NULL;
666	+	if( conv2 && ( conv2->Prob() < 0.0005) ) conv2 = NULL;
667	+
668	+	double zconv  = 0.;
669	+	double dzconv = 0.;
670	+
671	+	if(conv1 || conv2) {
672	+	if( !conv2 ){
673	+	//const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
674	+	const mithep::ThreeVector caloPos1(ph1->CaloPos());
675	+	zconv  = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
676	+	if( ph1->IsEB() ) {
677	+	double rho = conv1->Position().Rho();
678	+	if     ( rho < 15. ) dzconv = 0.06;
679	+	else if( rho < 60. ) dzconv = 0.67;
680	+	else                 dzconv = 2.04;
681	+	} else {
682	+	double z = conv1->Position().Z();
683	+	if     ( TMath::Abs(z) < 50. )   dzconv = 0.18;
684	+	else if( TMath::Abs(z) < 100.)   dzconv = 0.61;
685	+	else                 dzconv = 0.99;
686	+	}
687	+	} else if( !conv1 ) {
688	+	//const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
689	+	const mithep::ThreeVector caloPos2(ph2->CaloPos());
690	+	zconv  = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
691	+	if( ph2->IsEB() ) {
692	+	double rho = conv2->Position().Rho();
693	+	if     ( rho < 15. ) dzconv = 0.06;
694	+	else if( rho < 60. ) dzconv = 0.67;
695	+	else                 dzconv = 2.04;
696	+	} else {
697	+	double z = conv2->Position().Z();
698	+	if     ( TMath::Abs(z) < 50. )   dzconv = 0.18;
699	+	else if( TMath::Abs(z) < 100.)   dzconv = 0.61;
700	+	else                 dzconv = 0.99;
701	+	}
702	+	} else {
703	+	//const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
704	+	const mithep::ThreeVector caloPos1(ph1->CaloPos());
705	+	double z1  = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
706	+	double dz1 = 0.;
707	+	if( ph1->IsEB() ) {
708	+	double rho = conv1->Position().Rho();
709	+	if     ( rho < 15. ) dz1 = 0.06;
710	+	else if( rho < 60. ) dz1 = 0.67;
711	+	else                 dz1 = 2.04;
712	+	} else {
713	+	double z = conv1->Position().Z();
714	+	if     ( TMath::Abs(z) < 50. )   dz1 = 0.18;
715	+	else if( TMath::Abs(z) < 100.)   dz1 = 0.61;
716	+	else                 dz1 = 0.99;
717	+	}
718	+	//const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
719	+	const mithep::ThreeVector caloPos2(ph2->CaloPos());
720	+	double z2  = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
721	+	double dz2 = 0.;
722	+	if( ph2->IsEB() ) {
723	+	double rho = conv2->Position().Rho();
724	+	if     ( rho < 15. ) dz2 = 0.06;
725	+	else if( rho < 60. ) dz2 = 0.67;
726	+	else                 dz2 = 2.04;
727	+	} else {
728	+	double z = conv2->Position().Z();
729	+	if     ( TMath::Abs(z) < 50. )   dz2 = 0.18;
730	+	else if( TMath::Abs(z) < 100.)   dz2 = 0.61;
731	+	else                 dz2 = 0.99;
732	+	}
733	+
734	+	if(print) {
735	+	std::cout<<"  z1 = "<<z1<<std::endl;
736	+	std::cout<<" dz1 = "<<dz1<<std::endl;
737	+	std::cout<<"  z2 = "<<z2<<std::endl;
738	+	std::cout<<" dz2 = "<<dz2<<std::endl;
739	+	}
740	+
741	+	zconv  = ( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ;  // weighted average
742	+	dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
743	+	}
744	+
745	+	if(print) {
746	+	std::cout<<"  Conversion Z = "<<zconv<<std::endl;
747	+	std::cout<<"            dZ = "<<dzconv<<std::endl;
748	+	}
749	+
750	+
751	+	if(false) {
752	+
753	+	// loop over all ranked Vertices and choose the closest to the Conversion one
754	+	int maxVertices = ( ptgg > 30 ? 3 : 5);
755	+	double minDz = -1.;
756	+
757	+	if(print) std::cout<<std::endl<<"  looping over vertices... "<<ptgg<<"  "<<maxVertices<<std::endl;
758	+
759	+	for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
760	+
761	+	if(print) std::cout<<"     "<<iVtx<<"   has rank  "<<total_rank[iVtx]<<std::endl;
762	+
763	+	if(total_rank[iVtx] < maxVertices) {
764	+	const Vertex* tVtx = fPV->At(iVtx);
765	+	double tDz = TMath::Abs(zconv - tVtx->Z());
766	+	if(print) std::cout<<"     is considered with tDz = "<<tDz<<std::endl;
767	+	if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
768	+	minDz = tDz;
769	+	bestIdx = iVtx;
770	+	if(print) std::cout<<"      and is the best now."<<std::endl;
771	+	}
772	+	}
773	+	}
774	+	} else {
775	+	unsigned int bestIdxTmp = bestIdx;
776	+
777	+	// loop over all ranked Vertices and choose the closest to the Conversion one
778	+	double minDz = -1.;
779	+	int maxVertices = ( ptgg > 30 ? 3 : 5);
780	+
781	+	if(print) std::cout<<std::endl<<"  looping over vertices... "<<ptgg<<"  "<<maxVertices<<std::endl;
782	+
783	+	for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
784	+
785	+	if(print) std::cout<<"     "<<iVtx<<"   has rank  "<<total_rank[iVtx]<<std::endl;
786	+
787	+	const Vertex* tVtx = fPV->At(iVtx);
788	+	double tDz = TMath::Abs(zconv - tVtx->Z());
789	+	if(print) std::cout<<"     is considered with tDz = "<<tDz<<std::endl;
790	+	if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
791	+	minDz = tDz;
792	+	bestIdxTmp = iVtx;
793	+	if(print) std::cout<<"      and is the best now."<<std::endl;
794	+	}
795	+	}
796	+
797	+	// check if best Vtx is among higest ranked ones
798	+	if(total_rank[bestIdxTmp] < maxVertices)
799	+	bestIdx = bestIdxTmp;
800	+	}
801	+	}
802	+
803	+	delete[] total_rank  ;
804	+	return bestIdx;
805	+	}
806	+
807	+	double PhotonCiCMod::findHiggsPt() {
808	+
809	+	// loop over all GEN particles and look for status 1 photons
810	+	for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
811	+	const MCParticle* p = fMCParticles->At(i);
812	+	if( !(p->Is(MCParticle::kH)) ) continue;
813	+	return p->Pt();
814	+	}
815	+
816	+	return -1.0;
817	+	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines