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

Diff Legend

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