ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.61
Committed: Sat Feb 23 14:50:41 2013 UTC (12 years, 2 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.60: +18 -10 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include "MitPhysics/Mods/interface/PhotonPairSelector.h"
2 #include "MitAna/DataTree/interface/PhotonCol.h"
3 #include "MitAna/DataTree/interface/PFCandidateCol.h"
4 #include "MitAna/DataTree/interface/PFJetCol.h"
5 #include "MitAna/DataTree/interface/StableData.h"
6 #include "MitAna/DataTree/interface/StableParticle.h"
7 #include "MitAna/DataTree/interface/DecayParticleFwd.h"
8 #include "MitPhysics/Init/interface/ModNames.h"
9 #include "MitPhysics/Utils/interface/IsolationTools.h"
10 #include "MitPhysics/Utils/interface/PhotonTools.h"
11 #include "MitPhysics/Utils/interface/VertexTools.h"
12 #include "MitPhysics/Utils/interface/PhotonFix.h"
13 #include "MitPhysics/Utils/interface/MVATools.h"
14 #include "TDataMember.h"
15 #include <TNtuple.h>
16 #include <TRandom3.h>
17 #include <TSystem.h>
18 #include <TH1D.h>
19 #include "Math/SMatrix.h"
20 #include "Math/SVector.h"
21
22 using namespace mithep;
23
24 ClassImp(mithep::PhotonPairSelector)
25
26 //--------------------------------------------------------------------------------------------------
27 PhotonPairSelector::PhotonPairSelector(const char *name, const char *title) :
28 // Base Module...
29 BaseMod (name,title),
30 // define all the Branches to load
31 fPhotonBranchName (Names::gkPhotonBrn),
32 fElectronName (Names::gkElectronBrn),
33 fGoodElectronName (Names::gkElectronBrn),
34 fConversionName (Names::gkMvfConversionBrn),
35 fPFConversionName ("PFPhotonConversions"),
36 fTrackBranchName (Names::gkTrackBrn),
37 fPileUpDenName (Names::gkPileupEnergyDensityBrn),
38 fPVName (Names::gkPVBeamSpotBrn),
39 fBeamspotName (Names::gkBeamSpotBrn),
40 fPFCandName (Names::gkPFCandidatesBrn),
41 // MC specific stuff...
42 fMCParticleName (Names::gkMCPartBrn),
43 fPileUpName (Names::gkPileupInfoBrn),
44 fJetsName (Names::gkPFJetBrn),
45 fPFMetName ("PFMet"),
46 fGoodPhotonsName (ModNames::gkGoodPhotonsName),
47 fChosenVtxName ("HggChosenVtx"),
48 // ----------------------------------------
49 // Selection Types
50 fPhotonSelType ("NoSelection"),
51 fVertexSelType ("StdSelection"),
52 fPhSelType (kNoPhSelection),
53 fVtxSelType (kStdVtxSelection),
54 //-----------------------------------------
55 // Id Types fab: shouldn't we have kNone as default ???
56 fIdMVAType ("2011IdMVA"),
57 fIdType (MVATools::k2011IdMVA),
58 //-----------------------------------------
59 // preselection Type
60 fShowerShapeType ("2011ShowerShape"),
61 fSSType (PhotonTools::k2011ShowerShape),
62 // ----------------------------------------
63 fPhotonPtMin (20.0),
64 fPhotonEtaMax (2.5),
65 fLeadingPtMin (100.0/3.0),
66 fTrailingPtMin (100.0/4.0),
67 fIsData (false),
68 fPhotonsFromBranch (true),
69 fPVFromBranch (true),
70 fGoodElectronsFromBranch (kTRUE),
71 fUseSingleLegConversions (kTRUE),
72 // ----------------------------------------
73 // collections....
74 fPhotons (0),
75 fElectrons (0),
76 fConversions (0),
77 fPFConversions (0),
78 fTracks (0),
79 fPileUpDen (0),
80 fPV (0),
81 fBeamspot (0),
82 fPFCands (0),
83 fMCParticles (0),
84 fPileUp (0),
85 fJets (0),
86 fPFMet (0),
87 // ---------------------------------------
88 fDataEnCorr_EBlowEta_hR9central(0.),
89 fDataEnCorr_EBlowEta_hR9gap (0.),
90 fDataEnCorr_EBlowEta_lR9 (0.),
91 fDataEnCorr_EBlowEta_lR9central(0.),
92 fDataEnCorr_EBlowEta_lR9gap (0.),
93 fDataEnCorr_EBhighEta_hR9 (0.),
94 fDataEnCorr_EBhighEta_lR9 (0.),
95 fDataEnCorr_EElowEta_hR9 (0.),
96 fDataEnCorr_EElowEta_lR9 (0.),
97 fDataEnCorr_EEhighEta_hR9 (0.),
98 fDataEnCorr_EEhighEta_lR9 (0.),
99 fRunStart (0),
100 fRunEnd (0),
101
102 fMCSmear_EBlowEta_hR9central (0.),
103 fMCSmear_EBlowEta_hR9gap (0.),
104 fMCSmear_EBlowEta_lR9 (0.),
105 fMCSmear_EBlowEta_lR9central (0.),
106 fMCSmear_EBlowEta_lR9gap (0.),
107 fMCSmear_EBhighEta_hR9 (0.),
108 fMCSmear_EBhighEta_lR9 (0.),
109 fMCSmear_EElowEta_hR9 (0.),
110 fMCSmear_EElowEta_lR9 (0.),
111 fMCSmear_EEhighEta_hR9 (0.),
112 fMCSmear_EEhighEta_lR9 (0.),
113
114 fMCSmearMVA_EBlowEta_hR9central (0.),
115 fMCSmearMVA_EBlowEta_hR9gap (0.),
116 fMCSmearMVA_EBlowEta_lR9 (0.),
117 fMCSmearMVA_EBlowEta_lR9central (0.),
118 fMCSmearMVA_EBlowEta_lR9gap (0.),
119 fMCSmearMVA_EBhighEta_hR9 (0.),
120 fMCSmearMVA_EBhighEta_lR9 (0.),
121 fMCSmearMVA_EElowEta_hR9 (0.),
122 fMCSmearMVA_EElowEta_lR9 (0.),
123 fMCSmearMVA_EEhighEta_hR9 (0.),
124 fMCSmearMVA_EEhighEta_lR9 (0.),
125
126
127 // ---------------------------------------
128 fRng (new TRandom3()),
129 fPhFixString ("4_2"),
130 fEtaCorrections (0),
131 // ---------------------------------------
132 fDoDataEneCorr (true),
133 fDoMCSmear (true),
134 fUseSpecialSmearForDPMVA (false),
135 fDoVtxSelection (true),
136 fApplyEleVeto (true),
137 fInvertElectronVeto (kFALSE),
138 // ------------------------------------------------------------------
139 // this block should eventually be deleted...
140 fVariableType_2011 (10),
141
142 fEndcapWeights_2011 (gSystem->Getenv("CMSSW_BASE")+
143 TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
144 TString("Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
145 TString("weights.xml")),
146 fBarrelWeights_2011 (gSystem->Getenv("CMSSW_BASE")+
147 TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
148 TString("Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
149 TString("weights.xml")),
150 fVariableType_2012_globe (1201),
151 //fEndcapWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
152 //TString("/src/MitPhysics/data/")+
153 // TString("TMVA_EEpf_BDT_globe.")+
154 // TString("weights.xml")),
155 //fBarrelWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
156 // TString("/src/MitPhysics/data/")+
157 // TString("TMVA_EBpf_BDT_globe.")+
158 // TString("weights.xml")),
159 //fEndcapWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
160 // TString("/src/MitPhysics/data/")+
161 // TString("2012ICHEP_PhotonID_Endcap_BDT.")+
162 // TString("weights.xml")),
163 fEndcapWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
164 TString("/src/MitPhysics/data/")+
165 TString("2012ICHEP_PhotonID_Endcap_BDT.")+
166 TString("weights_PSCorr.xml")),
167 fBarrelWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
168 TString("/src/MitPhysics/data/")+
169 TString("2012ICHEP_PhotonID_Barrel_BDT.")+
170 TString("weights.xml")),
171 // --------------------------------------------------------------------------
172
173 fbdtCutBarrel (0.0744), //cuts give same eff (relative to presel) with cic
174 fbdtCutEndcap (0.0959), //cuts give same eff (relative to presel) with cic
175
176 // --------------------------------------------------------------------------
177
178 fDoShowerShapeScaling (kFALSE),
179 //
180 fMCErrScaleEB (1.0),
181 fMCErrScaleEE (1.0),
182 fRelativePtCuts (kFALSE),
183
184 fRhoType (RhoUtilities::CMS_RHO_RHOKT6PFJETS),
185
186 // ---------------------------------------------------------------------------
187 fApplyLeptonTag (kFALSE),
188 fLeptonTagElectronsName ("HggLeptonTagElectrons"),
189 fLeptonTagMuonsName ("HggLeptonTagMuons"),
190 fLeptonTagElectrons (0),
191 fLeptonTagMuons (0),
192
193 // ------------------------------------------------------
194 f2012HCP (kFALSE)
195
196 {
197 // Constructor.
198 }
199
200 PhotonPairSelector::~PhotonPairSelector()
201 {
202 if (fRng)
203 delete fRng;
204 }
205
206 //--------------------------------------------------------------------------------------------------
207 void PhotonPairSelector::Process()
208 {
209 // ------------------------------------------------------------
210 // Process entries of the tree.
211 LoadEventObject(fPhotonBranchName, fPhotons);
212 assert(fPhotons);
213 // lepton tag collections
214 if( fApplyLeptonTag ) {
215 LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
216 LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
217 }
218 //printf("check 0\n");
219 // -----------------------------------------------------------
220 // Output Photon Collection. It will ALWAYS contain either 0 or 2 photons
221 PhotonOArr *GoodPhotons = new PhotonOArr;
222 GoodPhotons->SetName(fGoodPhotonsName);
223 GoodPhotons->SetOwner(kTRUE);
224
225
226 VertexOArr *ChosenVtx = new VertexOArr;
227 ChosenVtx->SetName(fChosenVtxName);
228 ChosenVtx->SetOwner(kTRUE);
229
230 // add to event for other modules to use
231 AddObjThisEvt(GoodPhotons);
232 AddObjThisEvt(ChosenVtx);
233
234 //AddObjThisEvt(ChosenVtx);
235
236 if (fPhotons->GetEntries()<2)
237 return;
238
239 LoadEventObject(fElectronName, fElectrons);
240 LoadEventObject(fGoodElectronName, fGoodElectrons);
241 LoadEventObject(fConversionName, fConversions);
242 if (fUseSingleLegConversions) LoadEventObject(fPFConversionName, fPFConversions);
243 LoadEventObject(fTrackBranchName, fTracks);
244 LoadEventObject(fPileUpDenName, fPileUpDen);
245 LoadEventObject(fPVName, fPV);
246 LoadEventObject(fBeamspotName, fBeamspot);
247 LoadEventObject(fPFCandName, fPFCands);
248 LoadEventObject(fJetsName, fJets);
249 LoadEventObject(fPFMetName, fPFMet);
250
251 if (!fIsData) {
252 LoadBranch(fMCParticleName);
253 }
254
255 // ------------------------------------------------------------
256 // load event based information
257 Float_t rho = -99.;
258 if (fPileUpDen->GetEntries() > 0)
259 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
260 Float_t rho2012 = -99;
261 if (fPileUpDen->At(0)->RhoKt6PFJets()>0.) rho2012 = fPileUpDen->At(0)->RhoKt6PFJets();
262 else rho2012 = fPileUpDen->At(0)->Rho();
263 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
264
265 Float_t theRho = rho;
266 switch (fRhoType) {
267 case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
268 theRho = rho2012;
269 break;
270 case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
271 theRho = ( fPileUpDen->GetEntries() ? fPileUpDen->At(0)->RhoLowEta(): rho );
272 break;
273 case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
274 theRho = ( fPileUpDen->GetEntries() ? fPileUpDen->At(0)->Rho() : rho );
275 break;
276 case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
277 theRho = ( fPileUpDen->GetEntries() ? fPileUpDen->At(0)->RhoRandomLowEta() : rho );
278 break;
279 case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
280 theRho = ( fPileUpDen->GetEntries() ? fPileUpDen->At(0)->RhoRandom() : rho );
281 break;
282 default:
283 theRho = rho;
284 }
285 // ------------------------------------------------------------
286 // Get Event header for Run info etc.
287 const EventHeader* evtHead = this->GetEventHeader();
288 unsigned int evtNum = evtHead->EvtNum();
289 UInt_t runNumber = evtHead->RunNum();
290 Float_t _runNum = (Float_t) runNumber;
291 Float_t _lumiSec = (Float_t) evtHead->LumiSec();
292
293 // ------------------------------------------------------------
294 //if(evtNum==9009 || evtNum==9010){
295 // printf("evtNum:%d\n",evtNum);
296 // }
297 // here we'll store the preselected Photons (and which CiCCategory they are...)
298 PhotonOArr* preselPh = new PhotonOArr;
299 std::vector<PhotonTools::CiCBaseLineCats> preselCat;
300 preselCat.resize(0);
301
302 // 1. we do the pre-selection; but keep the non-passing photons in a second container...
303 for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
304 const Photon *ph = fPhotons->At(i);
305
306 if (ph->SCluster()->AbsEta()>= fPhotonEtaMax ||
307 (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566))
308 continue;
309
310 if (ph->Et() < fPhotonPtMin)
311 continue;
312
313 if (ph->HadOverEm() > 0.15)
314 continue;
315
316 if (ph->IsEB()) {
317 if (ph->CoviEtaiEta() > 0.015)
318 continue;
319 }
320 else {
321 if (ph->CoviEtaiEta() > 0.035)
322 continue;
323 }
324
325 // photon passes the preselection
326 ph->Mark(); // mark for later skimming
327 preselPh->Add(ph);
328 }
329
330 if (preselPh->GetEntries()<2) {
331 this->SkipEvent();
332 delete preselPh;
333 return;
334 }
335
336 //fill conversion collection for vertex selection, adding single leg conversions if needed
337 //note that momentum of single leg conversions needs to be recomputed from the track
338 //as it is not filled properly
339 DecayParticleOArr vtxconversions;
340 if (fUseSingleLegConversions) {
341 vtxconversions.SetOwner(kTRUE);
342 for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
343 DecayParticle *conv = new DecayParticle(*fConversions->At(iconv));
344 vtxconversions.AddOwned(conv);
345 }
346
347 for (UInt_t iconv=0; iconv<fPFConversions->GetEntries(); ++iconv) {
348 const DecayParticle *c = fPFConversions->At(iconv);
349 if (c->NDaughters()!=1) continue;
350
351 //std::cout<<" conv"<<iconv<<" mom = "<<c->Mom().X()<<" "<<c->Mom().Y()<<" "<<c->Mom().Z()<<" "<<std::endl;
352
353 DecayParticle *conv = new DecayParticle(*c);
354 const Track *trk = static_cast<const StableParticle*>(conv->Daughter(0))->Trk();
355 conv->SetMom(trk->Px(), trk->Py(), trk->Pz(), trk->P());
356 vtxconversions.AddOwned(conv);
357 }
358 }
359 else {
360 for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
361 const DecayParticle *c = fConversions->At(iconv);
362 vtxconversions.Add(c);
363 }
364 }
365
366
367 float higgspt = -99.;
368 float higgsz = -99.;
369 float higgsmass = -99.;
370
371 if (!fIsData) FindHiggsPtAndZ(higgspt,higgsz,higgsmass);
372
373 // second loop: sort & assign the right categories..
374 preselPh->Sort();
375 for (unsigned int iPh = 0; iPh <preselPh->GetEntries(); ++iPh) {
376 const Photon* ph = preselPh->At(iPh);
377 preselCat.push_back(PhotonTools::CiCBaseLineCat(ph));
378 }
379
380 // ------------------------------------------------------------
381 // compute how many pairs there are ...
382 unsigned int numPairs = 0;
383 if (preselPh->GetEntries() > 0)
384 numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
385 //printf("ming sync check numPairs:%d\n",numPairs);
386
387 // ... and create all possible pairs of pre-selected photons
388 std::vector<unsigned int> idx1st;
389 std::vector<unsigned int> idx2nd;
390 std::vector<PhotonTools::CiCBaseLineCats> cat1st;
391 std::vector<PhotonTools::CiCBaseLineCats> cat2nd;
392
393 // ... this will be used to store whether a given pair passes the cuts
394 std::vector<bool> pairPasses;
395
396 if (numPairs > 0) {
397 for (unsigned int i1st = 0; i1st <preselPh->GetEntries() - 1; ++i1st) {
398 for (unsigned int i2nd = i1st + 1; i2nd <preselPh->GetEntries(); ++i2nd) {
399 idx1st.push_back(i1st);
400 idx2nd.push_back(i2nd);
401 pairPasses.push_back(true);
402 }
403 }
404 }
405
406 // ------------------------------------------------------------
407 // array to store the index of 'chosen Vtx' for each pair
408 // const Vertex** theVtx = new const Vertex*[numPairs]; // holds the 'chosen' Vtx for each Pair
409 // Photon** fixPh1st = new Photon*[numPairs]; // holds the 1st Photon for each Pair
410 // Photon** fixPh2nd = new Photon*[numPairs]; // holds the 2nd photon for each Pair
411
412 std::vector<const Vertex*> theVtx; // holds the 'chosen' Vtx for each Pair
413 std::vector<Photon*> fixPh1st; // holds the 1st Photon for each Pair
414 std::vector<Photon*> fixPh2nd; // holds the 2nd photon for each Pair
415
416 theVtx.reserve(numPairs);
417 fixPh1st.reserve(numPairs);
418 fixPh2nd.reserve(numPairs);
419
420 // store pair-indices for pairs passing the selection
421 std::vector<unsigned int> passPairs;
422 passPairs.resize(0);
423 //if(evtNum==9009 || evtNum==9010){
424 // printf("check 1\n");
425 //}
426 // ------------------------------------------------------------
427 // Loop over all Pairs and do the 'incredible machine' running....
428 std::vector<int> leptag;
429 leptag.resize(numPairs);
430
431 for (unsigned int iPair = 0; iPair < numPairs; ++iPair) {
432 // first we need a hard copy of the incoming photons
433 fixPh1st.push_back(new Photon(*preselPh->At(idx1st[iPair])));
434 fixPh2nd.push_back(new Photon(*preselPh->At(idx2nd[iPair])));
435 // we also store the category, so we don't have to ask all the time...
436 cat1st.push_back(preselCat[idx1st[iPair]]);
437 cat2nd.push_back(preselCat[idx2nd[iPair]]);
438
439 // do showershape scaling.... (fab: outsourced to PhotonTools)
440 if (fDoShowerShapeScaling && !fIsData) {
441 PhotonTools::ScalePhotonShowerShapes(fixPh1st[iPair],fSSType);
442 PhotonTools::ScalePhotonShowerShapes(fixPh2nd[iPair],fSSType);
443 }
444
445 PhotonTools::eScaleCats escalecat1;
446 PhotonTools::eScaleCats escalecat2;
447
448 if(f2012HCP){
449 escalecat1 = PhotonTools::EScaleCatHCP(fixPh1st[iPair]);
450 escalecat2 = PhotonTools::EScaleCatHCP(fixPh2nd[iPair]);
451 }else{
452 escalecat1 = PhotonTools::EScaleCat(fixPh1st[iPair]);
453 escalecat2 = PhotonTools::EScaleCat(fixPh2nd[iPair]);
454 }
455
456 // now we dicide if we either scale (Data) or Smear (MC) the Photons
457 if (fIsData) {
458 if (fDoDataEneCorr) {
459 // starting with scale = 1.
460 double scaleFac1 = 1.;
461 double scaleFac2 = 1.;
462
463 //eta-dependent corrections
464
465 // checking the run Rangees ...
466 Int_t runRange = FindRunRangeIdx(runNumber);
467
468 if(runRange > -1) {
469 if(f2012HCP){
470 scaleFac1 *= GetDataEnCorrHCP(runRange, escalecat1);
471 scaleFac2 *= GetDataEnCorrHCP(runRange, escalecat2);
472 }else{
473 scaleFac1 *= GetDataEnCorr(runRange, escalecat1);
474 scaleFac2 *= GetDataEnCorr(runRange, escalecat2);
475 }
476 }
477 else {
478 printf("Error: Run Range not found for data energy scale correction\n");
479 assert(0);
480 }
481 PhotonTools::ScalePhoton(fixPh1st[iPair], scaleFac1);
482 PhotonTools::ScalePhoton(fixPh2nd[iPair], scaleFac2);
483 }
484 }
485
486 if (fDoMCSmear) {
487
488 double width1 = 0.;
489 double width2 = 0.;
490
491 if(f2012HCP){
492 width1 = GetMCSmearFacHCP(escalecat1);
493 width2 = GetMCSmearFacHCP(escalecat2);
494 }else {
495 width1 = GetMCSmearFac(escalecat1);
496 width2 = GetMCSmearFac(escalecat2);
497 }
498
499 if (!fIsData) {
500 // get the seed to do deterministic smearing...
501 UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
502 UInt_t seed1 = seedBase + (UInt_t) fixPh1st[iPair]->E() +
503 (UInt_t) (TMath::Abs(10.*fixPh1st[iPair]->SCluster()->Eta()));
504 UInt_t seed2 = seedBase + (UInt_t) fixPh2nd[iPair]->E() +
505 (UInt_t) (TMath::Abs(10.*fixPh2nd[iPair]->SCluster()->Eta()));
506 // get the smearing for MC photons..
507 PhotonTools::SmearPhoton(fixPh1st[iPair], fRng, width1, seed1);
508 PhotonTools::SmearPhoton(fixPh2nd[iPair], fRng, width2, seed2);
509 }
510
511 // in case we want to use specail smearing factors for Error computation, recompute widths
512 if( fUseSpecialSmearForDPMVA ) {
513
514 if(f2012HCP){
515 width1 = GetMCSmearFacHCP(escalecat1, true);
516 width2 = GetMCSmearFacHCP(escalecat2, true);
517 }else {
518 width1 = GetMCSmearFac(escalecat1, true);
519 width2 = GetMCSmearFac(escalecat2, true);
520 }
521
522 }
523
524 PhotonTools::SmearPhotonError(fixPh1st[iPair], width1);
525 PhotonTools::SmearPhotonError(fixPh2nd[iPair], width2);
526 }
527
528 // lepton tag for this pair -- ming
529 leptag[iPair] = -1;
530 if(fApplyLeptonTag){
531 leptag[iPair] = 0;
532 if ( fLeptonTagMuons->GetEntries() > 0 ) {
533 //printf("ming sync check muonpt0:%f\n",fLeptonTagMuons->At(0)->Pt());
534 //printf("ming sync check muonpt1:%f\n",fLeptonTagMuons->At(1)->Pt());
535 //printf("ming sync check ipair:%d deltar1:%f deltar2:%f\n",iPair,MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh1st[iPair]),MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh2nd[iPair]));
536 if( (MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh1st[iPair]) >= 1.0) &&
537 (MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh2nd[iPair]) >= 1.0)
538 ){
539 leptag[iPair] = 2;
540 }
541 }
542
543 if(leptag[iPair] <1){
544 if( fLeptonTagElectrons->GetEntries() > 0 ){
545 if( (MathUtils::DeltaR(*fLeptonTagElectrons->At(0),*fixPh1st[iPair]) >= 1) &&
546 (MathUtils::DeltaR(*fLeptonTagElectrons->At(0),*fixPh2nd[iPair]) >= 1) &&
547 (PhotonTools::ElectronVetoCiC(fixPh1st[iPair],fLeptonTagElectrons) >= 1) &&
548 (PhotonTools::ElectronVetoCiC(fixPh2nd[iPair],fLeptonTagElectrons) >= 1) &&
549 (TMath::Abs( (fixPh1st[iPair]->Mom()+fLeptonTagElectrons->At(0)->Mom()).M()-91.19 ) >= 10) &&
550 (TMath::Abs( (fixPh2nd[iPair]->Mom()+fLeptonTagElectrons->At(0)->Mom()).M()-91.19 ) >= 10)
551 //((fixPh1st[iPair]->Pt()/(fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M())>(45/120)) &&
552 //((fixPh2nd[iPair]->Pt()/(fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M())>(30/120))){
553 ){
554 /*int ph1passeveto=1;
555 int ph2passeveto=1;
556
557 for(UInt_t k=0;k<fElectrons->GetEntries();k++){
558 if(fElectrons->At(k)->BestTrk()->NMissingHits()==0){
559 if((fElectrons->At(k)->SCluster()==fixPh1st[iPair]->SCluster()) && (MathUtils::DeltaR(*fElectrons->At(k)->BestTrk(),*fixPh1st[iPair]) < 1)){
560 ph1passeveto=0;
561 }
562 if((fElectrons->At(k)->SCluster()==fixPh2nd[iPair]->SCluster()) && (MathUtils::DeltaR(*fElectrons->At(k)->BestTrk(),*fixPh2nd[iPair]) < 1)){
563 ph2passeveto=0;
564 }
565 }
566 }
567
568 if(ph1passeveto==1 && ph2passeveto==1){
569 leptag[iPair] = 1;
570 }*/
571
572 /*if(PhotonTools::PassElectronVeto(fixPh1st[iPair], fElectrons) && PhotonTools::PassElectronVeto(fixPh2nd[iPair], fElectrons)){
573 leptag[iPair] = 1;
574 }*/
575
576 /*if(PhotonTools::PassElectronVeto(fixPh1st[iPair], fElectrons) && PhotonTools::PassElectronVeto(fixPh2nd[iPair], fElectrons)){
577 leptag[iPair] = 1;
578 }*/
579
580 if(PhotonTools::ElectronVetoCiC(fixPh1st[iPair], fElectrons)>=1 && PhotonTools::ElectronVetoCiC(fixPh2nd[iPair], fElectrons)>=1){
581 leptag[iPair] = 1;
582 }
583
584 }
585 }
586 }
587 }
588
589 //probability that selected vertex is the correct one
590 Double_t vtxProb = 1.0;
591
592 // store the vertex for this pair
593 switch (fVtxSelType) {
594
595 case kStdVtxSelection:
596 theVtx[iPair] = fPV->At(0);
597 break;
598
599 case kCiCVtxSelection:
600 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
601 &vtxconversions,kFALSE,vtxProb);
602 break;
603
604 case kCiCMVAVtxSelection:
605 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
606 &vtxconversions,kTRUE,vtxProb);
607 break;
608
609 case kMITVtxSelection:
610 // need PFCandidate Collection
611 theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp,
612 mithep::FourVector((fixPh1st[iPair]->Mom()+
613 fixPh2nd[iPair]->Mom())));
614 break;
615
616 case kMetSigVtxSelection: {
617 // need PFCandidate Collection, otherwise use 0
618 if( !fJets || !fPFCands ) {
619 theVtx[iPair] = fPV->At(0);
620 break;
621 }
622
623 PFJetOArr pfjets;
624 for (UInt_t ijet=0; ijet<fJets->GetEntries(); ++ijet) {
625 const PFJet *pfjet = dynamic_cast<const PFJet*>(fJets->At(ijet));
626 if (pfjet && MathUtils::DeltaR(*pfjet,*fixPh1st[iPair])>0.3 && MathUtils::DeltaR(*pfjet,*fixPh2nd[iPair])>0.3) pfjets.Add(pfjet);
627 }
628
629 PFCandidateOArr pfcands;
630 for (UInt_t icand=0; icand<fPFCands->GetEntries(); ++icand) {
631 const PFCandidate *pfcand = fPFCands->At(icand);
632 if (MathUtils::DeltaR(*pfcand,*fixPh1st[iPair])>0.05 && MathUtils::DeltaR(*pfcand,*fixPh2nd[iPair])>0.05) pfcands.Add(pfcand);
633 }
634
635 double minsig = 1e6;
636 for (UInt_t ivtx=0; ivtx<fPV->GetEntries(); ++ivtx) {
637 const Vertex *v = fPV->At(ivtx);
638
639 // Met mmet = fMVAMet.GetMet( false,
640 // fixPh1st[iPair]->Pt(),fixPh1st[iPair]->Phi(),fixPh1st[iPair]->Eta(),
641 // fixPh2nd[iPair]->Pt(),fixPh2nd[iPair]->Phi(),fixPh2nd[iPair]->Eta(),
642 // fPFMet->At(0),
643 // &pfcands,fPV->At(ivtx),fPV, fPileUpDen->At(0)->Rho(),
644 // &pfjets,
645 // int(fPV->GetEntries()),
646 // kTRUE);
647
648 Met mmet = fMVAMet.GetMet( false,
649 0.,0.,0.,
650 0.,0.,0.,
651 fPFMet->At(0),
652 &pfcands,fPV->At(ivtx),fPV, fPileUpDen->At(0)->Rho(),
653 &pfjets,
654 int(fPV->GetEntries()),
655 kFALSE);
656
657 ThreeVector fullmet(mmet.Px() - fixPh1st[iPair]->Px() - fixPh2nd[iPair]->Px(),
658 mmet.Py() - fixPh1st[iPair]->Py() - fixPh2nd[iPair]->Py(),
659 0.);
660
661 // ThreeVector fullmet(mmet.Px(),
662 // mmet.Py(),
663 // 0.);
664
665 TMatrixD *metcov = fMVAMet.GetMetCovariance();
666
667 //Double_t metsigma = sqrt(fullmet.X()*fullmet.X()*(*metcov)(0,0) + 2.*fullmet.X()*fullmet.Y()*(*metcov)(1,0) + fullmet.Y()*fullmet.Y()*(*metcov)(1,1))/fullmet.Rho();
668
669 //Double_t metsig = fullmet.Rho()/metsigma;
670
671 ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
672 mcov(0,0) = (*metcov)(0,0);
673 mcov(0,1) = (*metcov)(0,1);
674 mcov(1,0) = (*metcov)(1,0);
675 mcov(1,1) = (*metcov)(1,1);
676
677 ROOT::Math::SVector<double,2> vmet;
678 vmet(0) = fullmet.X();
679 vmet(1) = fullmet.Y();
680
681 mcov.Invert();
682
683 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
684
685
686 Double_t sumptsq = 0.;
687 for (UInt_t ipfc = 0; ipfc<fPFCands->GetEntries(); ++ipfc) {
688 const PFCandidate *pfc = fPFCands->At(ipfc);
689 if (pfc->PFType()!=PFCandidate::eHadron || !pfc->HasTrackerTrk()) continue;
690 if (TMath::Abs( pfc->TrackerTrk()->DzCorrected(*v) ) > 0.2) continue;
691 if (TMath::Abs( pfc->TrackerTrk()->D0Corrected(*v) ) > 0.1) continue;
692 sumptsq += pfc->Pt()*pfc->Pt();
693 }
694
695 if ( sumptsq < 10.0 ) metsig = -99;
696
697 if (metsig<minsig && metsig>0.) {
698 minsig = metsig;
699 theVtx[iPair] = fPV->At(ivtx);
700 }
701
702 printf("ivtx = %i, sumptsq = %5f, met = %5f, metsig = %5f, dzgen = %5f\n",ivtx,sumptsq, fullmet.Rho(), metsig,fPV->At(ivtx)->Z()-higgsz);
703
704 }
705
706 double testpfmetx = 0.;
707 double testpfmety = 0.;
708 for (UInt_t icand=0; icand<pfcands.GetEntries(); ++icand) {
709 const PFCandidate *pfcand = pfcands.At(icand);
710 testpfmetx -= pfcand->Px();
711 testpfmety -= pfcand->Py();
712 //pfcands.Add(pfcand);
713 }
714
715 testpfmetx -= fixPh1st[iPair]->Px();
716 testpfmetx -= fixPh2nd[iPair]->Px();
717
718 testpfmety -= fixPh1st[iPair]->Py();
719 testpfmety -= fixPh2nd[iPair]->Py();
720
721 double testpfmet = sqrt(testpfmetx*testpfmetx + testpfmety*testpfmety);
722
723 printf("bestvtx: metsig = %5f, dzgen = %5f, diphopt = %5f, pfmet = %5f, testpfmet = %5f\n", minsig, theVtx[iPair]->Z()-higgsz, higgspt,fPFMet->At(0)->Pt(),testpfmet);
724
725 }
726 break;
727 default:
728 theVtx[iPair] = fPV->At(0);
729 }
730
731 /*
732 if(leptag == 1){
733 Double_t distVtx = 999.0;
734 Int_t closestVtx = 0;
735 for(UInt_t nv=0; nv<fPV->GetEntries(); nv++){
736 double dz = TMath::Abs(fLeptonTagMuons->At(0)->BestTrk()->DzCorrected(*fPV->At(nv)));
737 if(dz < distVtx) {
738 distVtx = dz;
739 closestVtx = nv;
740 }
741 }
742 theVtx[iPair] = fPV->At(closestVtx);
743 vtxProb = 1;
744 }
745 if(leptag == 2){
746 Double_t distVtx = 999.0;
747 Int_t closestVtx = 0;
748 for(UInt_t nv=0; nv<fPV->GetEntries(); nv++){
749 double dz = TMath::Abs(fLeptonTagElectrons->At(0)->GsfTrk()->DzCorrected(*fPV->At(nv)));
750 if(dz < distVtx) {
751 distVtx = dz;
752 closestVtx = nv;
753 }
754 }
755 theVtx[iPair] = fPV->At(closestVtx);
756 vtxProb = 1;
757 }
758 */
759
760 //set PV ref in photons
761 fixPh1st[iPair]->SetPV(theVtx[iPair]);
762 fixPh2nd[iPair]->SetPV(theVtx[iPair]);
763 fixPh1st[iPair]->SetVtxProb(vtxProb);
764 fixPh2nd[iPair]->SetVtxProb(vtxProb);
765
766 ThreeVector vtxpos = theVtx[iPair]->Position();
767
768 // fix the kinematics for both events
769 FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(vtxpos);
770 FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(vtxpos);
771 fixPh1st[iPair]->SetMom(newMom1st.X(), newMom1st.Y(), newMom1st.Z(), newMom1st.E());
772 fixPh2nd[iPair]->SetMom(newMom2nd.X(), newMom2nd.Y(), newMom2nd.Z(), newMom2nd.E());
773
774 double pairmass = (fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M();
775
776 double leadptcut = fLeadingPtMin;
777 double trailptcut = fTrailingPtMin;
778
779 if (fixPh2nd[iPair]->Pt() > fixPh1st[iPair]->Pt()) {
780 leadptcut = fTrailingPtMin;
781 trailptcut = fLeadingPtMin;
782 }
783
784
785 if (fRelativePtCuts) {
786 leadptcut = leadptcut*pairmass;
787 trailptcut = trailptcut*pairmass;
788 }
789
790
791 //compute id bdt values
792 Double_t bdt1 = -99.;
793 Double_t bdt2 = -99.;
794 // ---------------------------------------------------------------------------------------------------------------
795 // using new interface letting the MVATools handle the type... (fab)
796 if( fIdType != MVATools::kNone ) { // not strictly needed, but cold spped up things slightly if no MVA is needed...
797 bdt1 = fTool.GetMVAbdtValue(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,theRho,fPFCands,fElectrons,fApplyEleVeto);
798 bdt2 = fTool.GetMVAbdtValue(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,theRho,fPFCands,fElectrons,fApplyEleVeto);
799 }
800 fixPh1st[iPair]->SetIdMva(bdt1);
801 fixPh2nd[iPair]->SetIdMva(bdt2);
802 // ---------------------------------------------------------------------------------------------------------------
803 // superseedded by above code... (fab)
804 // switch (fIdType) {
805 // case MVATools::k2011IdMVA:
806 // bdt1 = fTool.GetMVAbdtValue_2011(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho,fElectrons,fApplyEleVeto);
807 // bdt2 = fTool.GetMVAbdtValue_2011(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho,fElectrons,fApplyEleVeto);
808 // fixPh1st[iPair]->SetIdMva(bdt1);
809 // fixPh2nd[iPair]->SetIdMva(bdt2);
810 // break;
811
812 // case MVATools::k2012IdMVA_globe:
813 // bdt1 = fTool.GetMVAbdtValue_2012_globe(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho2012,fPFCands,fElectrons,fApplyEleVeto);
814 // bdt2 = fTool.GetMVAbdtValue_2012_globe(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho2012,fPFCands,fElectrons,fApplyEleVeto);
815 // fixPh1st[iPair]->SetIdMva(bdt1);
816 // fixPh2nd[iPair]->SetIdMva(bdt2);
817 // break;
818
819 // default:
820 // fixPh1st[iPair]->SetIdMva(-99.);
821 // fixPh2nd[iPair]->SetIdMva(-99.);
822 // }
823 // ---------------------------------------------------------------------------------------------------------------
824
825 //printf("applying id\n");
826
827 // check if both photons pass the CiC selection
828 // FIX-ME: Add other possibilities....
829 bool pass1 = false;
830 bool pass2 = false;
831
832 switch (fPhSelType) {
833
834 // --------------------------------------------------------------------
835 // trivial (no) selection with only pt-cuts
836 case kNoPhSelection:
837 pass1 = (fixPh1st[iPair]->Pt() > leadptcut );
838 pass2 = (fixPh2nd[iPair]->Pt() > trailptcut);
839 break;
840
841 // --------------------------------------------------------------------
842 // CiC4 Selection as used for the 2011 7TeV Baseline analysis
843 case kCiCPhSelection:
844 pass1 = PhotonTools::PassCiCSelection(fixPh1st[iPair], theVtx[iPair], fTracks,
845 fElectrons, fPV, rho, leadptcut, fApplyEleVeto);
846 if (pass1)
847 pass2 = PhotonTools::PassCiCSelection(fixPh2nd[iPair], theVtx[iPair], fTracks,
848 fElectrons, fPV, rho, trailptcut, fApplyEleVeto);
849 break;
850
851 // --------------------------------------------------------------------
852 // PF-CiC4 Selection, as used in the 2012 8TeV Baseline analysis
853 case kCiCPFPhSelection:
854
855 // loose preselection for mva
856 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
857 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh1st[iPair],fElectrons,fConversions,bsp,
858 //PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
859 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
860 fInvertElectronVeto);
861 if (pass1)
862 pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
863 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh2nd[iPair],fElectrons,fConversions,bsp,
864 //PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
865 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
866 fInvertElectronVeto);
867
868 pass1 = pass1 && PhotonTools::PassCiCPFIsoSelection(fixPh1st[iPair], theVtx[iPair], fPFCands,
869 fPV, rho2012, leadptcut);
870 if (pass1)
871 pass2 = pass2 && PhotonTools::PassCiCPFIsoSelection(fixPh2nd[iPair], theVtx[iPair], fPFCands,
872 fPV, rho2012, trailptcut);
873 break;
874
875 // --------------------------------------------------------------------
876 // MVA selection
877 case kMVAPhSelection://MVA
878
879 pass1 = fixPh1st[iPair]->Pt()>leadptcut &&
880 PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
881 fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
882 ( ( fixPh1st[iPair]->SCluster()->AbsEta() < 1.5 && fixPh1st[iPair]->IdMva() > fbdtCutBarrel )
883 || ( fixPh1st[iPair]->IdMva() > fbdtCutEndcap ) );
884
885 // we've already compyted the MVA varible above, not needed to redo...
886 // fTool.PassMVASelection(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho,
887 // fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
888
889
890 if (pass1)
891 pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
892 PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
893 fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
894 ( ( fixPh2nd[iPair]->SCluster()->AbsEta() < 1.5 && fixPh2nd[iPair]->IdMva() > fbdtCutBarrel )
895 || ( fixPh2nd[iPair]->IdMva() > fbdtCutEndcap ) );
896
897 // fTool.PassMVASelection(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho,
898 // fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
899
900 break;
901
902 // --------------------------------------------------------------------
903 // MIT proposed pre-selection as used in the 2011 7TeV MVA analyses
904 case kMITPhSelection:
905 // loose preselection for mva
906 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
907 PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
908 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
909 fInvertElectronVeto);
910 if (pass1)
911 pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
912 PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
913 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
914 fInvertElectronVeto);
915
916 break;
917
918 // --------------------------------------------------------------------
919 // updated (PF absed) pre-selection as used in the 2012 8TeV MVA/Baseline analyses
920 case kMITPFPhSelection:
921 // loose preselection for mva
922 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
923 PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
924 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
925 fInvertElectronVeto);
926
927 if (pass1)
928 pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
929 PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
930 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
931 fInvertElectronVeto);
932
933 break;
934
935 // --------------------------------------------------------------------
936 // updated (PF absed) pre-selection as used in the 2012 8TeV MVA/Baseline analyses
937 case kMITPFPhSelectionNoEcal:
938 // loose preselection for mva
939 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
940 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh1st[iPair],fElectrons,fConversions,bsp,
941 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
942 fInvertElectronVeto);
943
944 if (pass1)
945 pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
946 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh2nd[iPair],fElectrons,fConversions,bsp,
947 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
948 fInvertElectronVeto);
949
950 break;
951
952 default:
953 pass1 = true;
954 pass2 = true;
955 }
956
957 //match to good electrons if requested
958 // if (fInvertElectronVeto) {
959 // pass1 &= !PhotonTools::PassElectronVeto(fixPh1st[iPair],fGoodElectrons);
960 // pass2 &= !PhotonTools::PassElectronVeto(fixPh2nd[iPair],fGoodElectrons);
961 // }
962 // finally, if both Photons pass the selections, add the pair to the passing pairs
963 if (pass1 && pass2)
964 passPairs.push_back(iPair);
965 }
966 //if(evtNum==9009 || evtNum==9010){
967 // printf("check 2\n");
968 // }
969 // ---------------------------------------------------------------
970 // ... we're almost done, stay focused...
971 // loop over all passing pairs and find the one with the highest sum Et
972 const Vertex* vtx = NULL;
973 Photon* phHard = NULL;
974 Photon* phSoft = NULL;
975
976 // not used at all....
977 //PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
978 //PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
979
980 double maxSumEt = 0.;
981 int maxleptag = 0;
982 for (unsigned int iPair=0; iPair<passPairs.size(); ++iPair) {
983 double sumEt = fixPh1st[passPairs[iPair]]->Et()+fixPh2nd[passPairs[iPair]]->Et();
984 //if (((sumEt > maxSumEt) && leptag[iPair] == maxleptag) || (leptag[iPair] > maxleptag)) {
985 if (((sumEt > maxSumEt) && leptag[passPairs[iPair]] == maxleptag) || (leptag[passPairs[iPair]] > maxleptag)) {
986 maxSumEt = sumEt;
987 maxleptag = leptag[passPairs[iPair]];
988 phHard = fixPh1st[passPairs[iPair]];
989 phSoft = fixPh2nd[passPairs[iPair]];
990 //catPh1 = cat1st [passPairs[iPair]];
991 //catPh2 = cat2nd [passPairs[iPair]];
992 vtx = theVtx[iPair];
993 }
994 }
995 //if(evtNum==9009 || evtNum==9010){
996 // printf("check 3\n");
997 // }
998 for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
999 if (fixPh1st[iPair]!=phHard) delete fixPh1st[iPair];
1000 if (fixPh2nd[iPair]!=phSoft) delete fixPh2nd[iPair];
1001 }
1002
1003 // ---------------------------------------------------------------
1004 // we have the Photons (*PARTY*)... compute some useful qunatities
1005
1006 if(phHard && phSoft) {
1007 GoodPhotons->AddOwned(phHard);
1008 GoodPhotons->AddOwned(phSoft);
1009 }
1010
1011 // we also store the chosen Vtx, so later modules can use it
1012 //if(evtNum==9009 || evtNum==9010){
1013 // printf("check 4\n");
1014 //}
1015 Vertex* chosenVtx = NULL;
1016 if ( vtx ) {
1017 chosenVtx = new Vertex( *vtx );
1018 ChosenVtx->AddOwned( chosenVtx );
1019 }
1020 //printf("ChosenVtx->GetEntries():%d\n",ChosenVtx->GetEntries()ChosenVtx->GetEntries());
1021 //if(evtNum==9009 || evtNum==9010){
1022 // printf("check 5\n");
1023 //}
1024 // sort according to pt
1025 GoodPhotons->Sort();
1026
1027 // delete auxiliary photon collection...
1028 delete preselPh;
1029
1030 //delete[] theVtx;
1031 //if(evtNum==9009 || evtNum==9010){
1032 // printf(" GoodPhotons->GetEntries():%d\n",GoodPhotons->GetEntries());
1033 // printf("check 6\n");
1034 //}
1035 return;
1036 }
1037
1038 //---------------------------------------------------------------------------------------------------
1039 void PhotonPairSelector::SlaveBegin()
1040 {
1041 // Run startup code on the computer (slave) doing the actual analysis. Here, we just request the
1042 // photon collection branch.
1043
1044 if( fApplyLeptonTag ) {
1045 ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
1046 ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
1047 }
1048
1049 // load all branches
1050 ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
1051 ReqEventObject(fTrackBranchName, fTracks, true);
1052 ReqEventObject(fElectronName, fElectrons, true);
1053 ReqEventObject(fGoodElectronName, fGoodElectrons,fGoodElectronsFromBranch);
1054 ReqEventObject(fPileUpDenName, fPileUpDen, true);
1055 ReqEventObject(fPVName, fPV, fPVFromBranch);
1056 ReqEventObject(fConversionName, fConversions, true);
1057 if (fUseSingleLegConversions) ReqEventObject(fPFConversionName, fPFConversions, true);
1058 ReqEventObject(fBeamspotName, fBeamspot, true);
1059 ReqEventObject(fPFCandName, fPFCands, true);
1060 ReqEventObject(fJetsName, fJets, false);
1061 ReqEventObject(fPFMetName, fPFMet, true);
1062 if (!fIsData) {
1063 //ReqBranch(fPileUpName, fPileUp);
1064 ReqBranch(fMCParticleName, fMCParticles);
1065 }
1066
1067 // determine photon selection type
1068 if (fPhotonSelType.CompareTo("CiCSelection") == 0)
1069 fPhSelType = kCiCPhSelection;
1070 else if (fPhotonSelType.CompareTo("CiCPFSelection") == 0)
1071 fPhSelType = kCiCPFPhSelection;
1072 else if (fPhotonSelType.CompareTo("MVASelection") == 0) //MVA
1073 fPhSelType = kMVAPhSelection;
1074 else if (fPhotonSelType.CompareTo("MITSelection") == 0)
1075 fPhSelType = kMITPhSelection;
1076 else if (fPhotonSelType.CompareTo("MITPFSelection") == 0)
1077 fPhSelType = kMITPFPhSelection;
1078 else if (fPhotonSelType.CompareTo("MITPFSelectionNoEcal") == 0)
1079 fPhSelType = kMITPFPhSelectionNoEcal;
1080 else if (fPhotonSelType.CompareTo("NoSelection") == 0)
1081 fPhSelType = kNoPhSelection;
1082 else {
1083 std::cerr<<" Photon Seclection "<<fPhotonSelType<<" not implemented."<<std::endl;
1084 return;
1085 }
1086
1087 if (fVertexSelType.CompareTo("CiCSelection") == 0) {
1088 fVtxSelType = kCiCVtxSelection;
1089 fVtxTools.InitP(1);
1090 }
1091 else if (fVertexSelType.CompareTo("MITSelection") == 0)
1092 fVtxSelType = kMITVtxSelection;
1093 else if (fVertexSelType.CompareTo("CiCMVASelection") == 0) {
1094 fVtxSelType = kCiCMVAVtxSelection;
1095 fVtxTools.InitP(1);
1096 }
1097 else if (fVertexSelType.CompareTo("CiCMVA2012Selection") == 0) {
1098 fVtxSelType = kCiCMVAVtxSelection;
1099 fVtxTools.InitP(2);
1100 }
1101 else if (fVertexSelType.CompareTo("MetSigSelection") == 0)
1102 fVtxSelType = kMetSigVtxSelection;
1103 else if (fVertexSelType.CompareTo("ZeroVtxSelection") == 0)
1104 fVtxSelType = kStdVtxSelection;
1105 else {
1106 std::cerr<<" Vertex Seclection "<<fVertexSelType<<" not implemented."<<std::endl;
1107 return;
1108 }
1109
1110 if (fIdMVAType.CompareTo("2011IdMVA") == 0)
1111 fIdType = MVATools::k2011IdMVA;
1112 else if (fIdMVAType.CompareTo("2012IdMVA_globe") == 0)
1113 fIdType = MVATools::k2012IdMVA_globe;
1114 else if (fIdMVAType.CompareTo("2012IdMVA") == 0)
1115 fIdType = MVATools::k2012IdMVA;
1116 else if (fIdMVAType.CompareTo("2011IdMVA_HZg") == 0)
1117 fIdType = MVATools::k2011IdMVA_HZg;
1118 else if (fIdMVAType.CompareTo("None") == 0)
1119 fIdType = MVATools::kNone;
1120 else {
1121 std::cerr<<" Id MVA "<<fIdMVAType<<" not implemented."<<std::endl;
1122 return;
1123 }
1124
1125
1126 if (fShowerShapeType.CompareTo("None") == 0)
1127 fSSType = PhotonTools::kNoShowerShapeScaling;
1128 else if (fShowerShapeType.CompareTo("2011ShowerShape") == 0)
1129 fSSType = PhotonTools::k2011ShowerShape;
1130 else if (fShowerShapeType.CompareTo("2012ShowerShape") == 0)
1131 fSSType = PhotonTools::k2012ShowerShape;
1132 else {
1133 std::cerr<<"shower shape scale "<<fShowerShapeType<<" not implemented."<<std::endl;
1134 return;
1135 }
1136
1137 if (fIsData)
1138 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
1139 else
1140 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
1141
1142 printf("initialize photon pair selector\n");
1143
1144 // ----------------------------------------------------------------
1145 // anopther replace block (fab):
1146 // delegate the choice of the weight files to the MVATools...
1147 fTool.InitializeMVA(fIdType);
1148
1149 // switch (fIdType) {
1150 // case k2011IdMVA:
1151 // fTool.InitializeMVA(fVariableType_2011,fEndcapWeights_2011,fBarrelWeights_2011);
1152 // break;
1153
1154 // case k2012IdMVA_globe:
1155 // fTool.InitializeMVA(fVariableType_2012_globe,fEndcapWeights_2012_globe,fBarrelWeights_2012_globe);
1156 // break;
1157 // }
1158 // ----------------------------------------------------------------
1159
1160 if (fVtxSelType==kMetSigVtxSelection) {
1161 //fMVAMet.Initialize();
1162 // fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1163 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1164 // TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
1165 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_42.root"))),
1166 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_42.root"))),
1167 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1_42.root"))),
1168 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2_42.root")))
1169 // );
1170
1171 fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1172 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1173 TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
1174 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_52.root"))),
1175 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_52.root"))),
1176 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1cov_52.root"))),
1177 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2cov_52.root")))
1178 );
1179
1180 }
1181
1182 }
1183
1184 // ----------------------------------------------------------------------------------------
1185 // some helpfer functions....
1186 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
1187 {
1188 pt = -999.;
1189 decayZ = -999.;
1190 mass = -999.;
1191
1192 // loop over all GEN particles and look for status 1 photons
1193 for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
1194 const MCParticle* p = fMCParticles->At(i);
1195 if( p->Is(MCParticle::kH) || (!fApplyEleVeto && p->AbsPdgId()==23) ) {
1196 pt=p->Pt();
1197 decayZ = p->DecayVertex().Z();
1198 mass = p->Mass();
1199 break;
1200 }
1201 }
1202
1203 return;
1204 }
1205
1206 //---------------------------------------------------------------------------------------------------
1207 Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run)
1208 {
1209 // this routine looks for the idx of the run-range
1210 Int_t runRange=-1;
1211 for (UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
1212 if (run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
1213 runRange = (Int_t) iRun;
1214 return runRange;
1215 }
1216 }
1217 return runRange;
1218 }
1219
1220 //---------------------------------------------------------------------------------------------------
1221 Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::eScaleCats cat)
1222 {
1223 switch (cat) {
1224 case PhotonTools::kEBhighEtaGold:
1225 return fDataEnCorr_EBhighEta_hR9[runRange];
1226 case PhotonTools::kEBhighEtaBad:
1227 return fDataEnCorr_EBhighEta_lR9[runRange];
1228 case PhotonTools::kEBlowEtaGoldCenter:
1229 return fDataEnCorr_EBlowEta_hR9central[runRange];
1230 case PhotonTools::kEBlowEtaGoldGap:
1231 return fDataEnCorr_EBlowEta_hR9gap[runRange];
1232 case PhotonTools::kEBlowEtaBad:
1233 return fDataEnCorr_EBlowEta_lR9[runRange];
1234 case PhotonTools::kEEhighEtaGold:
1235 return fDataEnCorr_EEhighEta_hR9[runRange];
1236 case PhotonTools::kEEhighEtaBad:
1237 return fDataEnCorr_EEhighEta_lR9[runRange];
1238 case PhotonTools::kEElowEtaGold:
1239 return fDataEnCorr_EElowEta_hR9[runRange];
1240 case PhotonTools::kEElowEtaBad:
1241 return fDataEnCorr_EElowEta_lR9[runRange];
1242 default:
1243 return 1.;
1244 }
1245 }
1246
1247 //---------------------------------------------------------------------------------------------------
1248 Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::eScaleCats cat, bool useSpecialSmear)
1249 {
1250
1251 if(!useSpecialSmear) {
1252 switch (cat) {
1253 case PhotonTools::kEBhighEtaGold:
1254 return fMCSmear_EBhighEta_hR9;
1255 case PhotonTools::kEBhighEtaBad:
1256 return fMCSmear_EBhighEta_lR9;
1257 case PhotonTools::kEBlowEtaGoldCenter:
1258 return fMCSmear_EBlowEta_hR9central;
1259 case PhotonTools::kEBlowEtaGoldGap:
1260 return fMCSmear_EBlowEta_hR9gap;
1261 case PhotonTools::kEBlowEtaBad:
1262 return fMCSmear_EBlowEta_lR9;
1263 case PhotonTools::kEEhighEtaGold:
1264 return fMCSmear_EEhighEta_hR9;
1265 case PhotonTools::kEEhighEtaBad:
1266 return fMCSmear_EEhighEta_lR9;
1267 case PhotonTools::kEElowEtaGold:
1268 return fMCSmear_EElowEta_hR9;
1269 case PhotonTools::kEElowEtaBad:
1270 return fMCSmear_EElowEta_lR9;
1271 default:
1272 return 1.;
1273 }
1274 } else {
1275 switch (cat) {
1276 case PhotonTools::kEBhighEtaGold:
1277 return fMCSmearMVA_EBhighEta_hR9;
1278 case PhotonTools::kEBhighEtaBad:
1279 return fMCSmearMVA_EBhighEta_lR9;
1280 case PhotonTools::kEBlowEtaGoldCenter:
1281 return fMCSmearMVA_EBlowEta_hR9central;
1282 case PhotonTools::kEBlowEtaGoldGap:
1283 return fMCSmearMVA_EBlowEta_hR9gap;
1284 case PhotonTools::kEBlowEtaBad:
1285 return fMCSmearMVA_EBlowEta_lR9;
1286 case PhotonTools::kEEhighEtaGold:
1287 return fMCSmearMVA_EEhighEta_hR9;
1288 case PhotonTools::kEEhighEtaBad:
1289 return fMCSmearMVA_EEhighEta_lR9;
1290 case PhotonTools::kEElowEtaGold:
1291 return fMCSmearMVA_EElowEta_hR9;
1292 case PhotonTools::kEElowEtaBad:
1293 return fMCSmearMVA_EElowEta_lR9;
1294 default:
1295 return 1.;
1296 }
1297 }
1298
1299 }
1300
1301 Double_t PhotonPairSelector::GetDataEnCorrHCP(Int_t runRange, PhotonTools::eScaleCats cat)
1302 {
1303 switch (cat) {
1304 case PhotonTools::kEBhighEtaGold:
1305 return fDataEnCorr_EBhighEta_hR9[runRange];
1306 case PhotonTools::kEBhighEtaBad:
1307 return fDataEnCorr_EBhighEta_lR9[runRange];
1308 case PhotonTools::kEBlowEtaGoldCenter:
1309 return fDataEnCorr_EBlowEta_hR9central[runRange];
1310 case PhotonTools::kEBlowEtaGoldGap:
1311 return fDataEnCorr_EBlowEta_hR9gap[runRange];
1312 case PhotonTools::kEBlowEtaBadCenter:
1313 return fDataEnCorr_EBlowEta_lR9central[runRange];
1314 case PhotonTools::kEBlowEtaBadGap:
1315 return fDataEnCorr_EBlowEta_lR9gap[runRange];
1316 case PhotonTools::kEEhighEtaGold:
1317 return fDataEnCorr_EEhighEta_hR9[runRange];
1318 case PhotonTools::kEEhighEtaBad:
1319 return fDataEnCorr_EEhighEta_lR9[runRange];
1320 case PhotonTools::kEElowEtaGold:
1321 return fDataEnCorr_EElowEta_hR9[runRange];
1322 case PhotonTools::kEElowEtaBad:
1323 return fDataEnCorr_EElowEta_lR9[runRange];
1324 default:
1325 return 1.;
1326 }
1327 }
1328
1329 //---------------------------------------------------------------------------------------------------
1330 Double_t PhotonPairSelector::GetMCSmearFacHCP(PhotonTools::eScaleCats cat, bool useSpecialSmear)
1331 {
1332
1333 if(!useSpecialSmear) {
1334
1335 switch (cat) {
1336 case PhotonTools::kEBhighEtaGold:
1337 return fMCSmear_EBhighEta_hR9;
1338 case PhotonTools::kEBhighEtaBad:
1339 return fMCSmear_EBhighEta_lR9;
1340 case PhotonTools::kEBlowEtaGoldCenter:
1341 return fMCSmear_EBlowEta_hR9central;
1342 case PhotonTools::kEBlowEtaGoldGap:
1343 return fMCSmear_EBlowEta_hR9gap;
1344 case PhotonTools::kEBlowEtaBadCenter:
1345 return fMCSmear_EBlowEta_lR9central;
1346 case PhotonTools::kEBlowEtaBadGap:
1347 return fMCSmear_EBlowEta_lR9gap;
1348 case PhotonTools::kEEhighEtaGold:
1349 return fMCSmear_EEhighEta_hR9;
1350 case PhotonTools::kEEhighEtaBad:
1351 return fMCSmear_EEhighEta_lR9;
1352 case PhotonTools::kEElowEtaGold:
1353 return fMCSmear_EElowEta_hR9;
1354 case PhotonTools::kEElowEtaBad:
1355 return fMCSmear_EElowEta_lR9;
1356 default:
1357 return 1.;
1358 }
1359
1360 } else {
1361
1362 switch (cat) {
1363 case PhotonTools::kEBhighEtaGold:
1364 return fMCSmearMVA_EBhighEta_hR9;
1365 case PhotonTools::kEBhighEtaBad:
1366 return fMCSmearMVA_EBhighEta_lR9;
1367 case PhotonTools::kEBlowEtaGoldCenter:
1368 return fMCSmearMVA_EBlowEta_hR9central;
1369 case PhotonTools::kEBlowEtaGoldGap:
1370 return fMCSmearMVA_EBlowEta_hR9gap;
1371 case PhotonTools::kEBlowEtaBadCenter:
1372 return fMCSmearMVA_EBlowEta_lR9central;
1373 case PhotonTools::kEBlowEtaBadGap:
1374 return fMCSmearMVA_EBlowEta_lR9gap;
1375 case PhotonTools::kEEhighEtaGold:
1376 return fMCSmearMVA_EEhighEta_hR9;
1377 case PhotonTools::kEEhighEtaBad:
1378 return fMCSmearMVA_EEhighEta_lR9;
1379 case PhotonTools::kEElowEtaGold:
1380 return fMCSmearMVA_EElowEta_hR9;
1381 case PhotonTools::kEElowEtaBad:
1382 return fMCSmearMVA_EElowEta_lR9;
1383 default:
1384 return 1.;
1385 }
1386
1387 }
1388 }
1389
1390 //---------------------------------------------------------------------------------------------------
1391 Float_t PhotonPairSelector::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
1392 PhotonTools::CiCBaseLineCats cat2)
1393 {
1394 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
1395 bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
1396
1397 bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
1398 bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
1399
1400 if (ph1IsEB && ph2IsEB)
1401 return (ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
1402
1403 return (ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
1404 }
1405