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

# User Rev Content
1 fabstoec 1.1 #include "MitPhysics/Mods/interface/PhotonPairSelector.h"
2     #include "MitAna/DataTree/interface/PhotonCol.h"
3     #include "MitAna/DataTree/interface/PFCandidateCol.h"
4 bendavid 1.31 #include "MitAna/DataTree/interface/PFJetCol.h"
5 bendavid 1.8 #include "MitAna/DataTree/interface/StableData.h"
6     #include "MitAna/DataTree/interface/StableParticle.h"
7 bendavid 1.36 #include "MitAna/DataTree/interface/DecayParticleFwd.h"
8 fabstoec 1.1 #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 bendavid 1.6 #include "MitPhysics/Utils/interface/PhotonFix.h"
13 mingyang 1.15 #include "MitPhysics/Utils/interface/MVATools.h"
14 bendavid 1.5 #include "TDataMember.h"
15 fabstoec 1.1 #include <TNtuple.h>
16     #include <TRandom3.h>
17 bendavid 1.6 #include <TSystem.h>
18 bendavid 1.9 #include <TH1D.h>
19 bendavid 1.31 #include "Math/SMatrix.h"
20     #include "Math/SVector.h"
21 fabstoec 1.1
22     using namespace mithep;
23    
24     ClassImp(mithep::PhotonPairSelector)
25 fabstoec 1.47
26 fabstoec 1.1 //--------------------------------------------------------------------------------------------------
27 bendavid 1.48 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 fabstoec 1.49 // Id Types fab: shouldn't we have kNone as default ???
56 bendavid 1.48 fIdMVAType ("2011IdMVA"),
57 fabstoec 1.49 fIdType (MVATools::k2011IdMVA),
58 bendavid 1.48 //-----------------------------------------
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 mingyang 1.53 fUseSingleLegConversions (kTRUE),
72 bendavid 1.48 // ----------------------------------------
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 mingyang 1.52 fDataEnCorr_EBlowEta_lR9central(0.),
92     fDataEnCorr_EBlowEta_lR9gap (0.),
93 bendavid 1.48 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 fabstoec 1.56
102 bendavid 1.48 fMCSmear_EBlowEta_hR9central (0.),
103     fMCSmear_EBlowEta_hR9gap (0.),
104     fMCSmear_EBlowEta_lR9 (0.),
105 mingyang 1.52 fMCSmear_EBlowEta_lR9central (0.),
106     fMCSmear_EBlowEta_lR9gap (0.),
107 bendavid 1.48 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 fabstoec 1.56
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 bendavid 1.48 // ---------------------------------------
128     fRng (new TRandom3()),
129     fPhFixString ("4_2"),
130     fEtaCorrections (0),
131     // ---------------------------------------
132     fDoDataEneCorr (true),
133     fDoMCSmear (true),
134 fabstoec 1.56 fUseSpecialSmearForDPMVA (false),
135 bendavid 1.48 fDoVtxSelection (true),
136     fApplyEleVeto (true),
137     fInvertElectronVeto (kFALSE),
138 fabstoec 1.49 // ------------------------------------------------------------------
139     // this block should eventually be deleted...
140 bendavid 1.48 fVariableType_2011 (10),
141 mingyang 1.52
142 bendavid 1.48 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 fabstoec 1.49 // --------------------------------------------------------------------------
172 mingyang 1.52
173 bendavid 1.48 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 mingyang 1.52
176 fabstoec 1.49 // --------------------------------------------------------------------------
177 mingyang 1.52
178 bendavid 1.48 fDoShowerShapeScaling (kFALSE),
179     //
180     fMCErrScaleEB (1.0),
181     fMCErrScaleEE (1.0),
182 fabstoec 1.50 fRelativePtCuts (kFALSE),
183 mingyang 1.52
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 fabstoec 1.1 {
197     // Constructor.
198     }
199    
200 paus 1.27 PhotonPairSelector::~PhotonPairSelector()
201     {
202     if (fRng)
203     delete fRng;
204 fabstoec 1.1 }
205    
206     //--------------------------------------------------------------------------------------------------
207     void PhotonPairSelector::Process()
208 mingyang 1.57 {
209 paus 1.27 // ------------------------------------------------------------
210     // Process entries of the tree.
211 fabstoec 1.1 LoadEventObject(fPhotonBranchName, fPhotons);
212 mingyang 1.57 assert(fPhotons);
213 mingyang 1.52 // lepton tag collections
214     if( fApplyLeptonTag ) {
215     LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
216     LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
217     }
218 mingyang 1.57 //printf("check 0\n");
219 bendavid 1.8 // -----------------------------------------------------------
220 paus 1.27 // Output Photon Collection. It will ALWAYS contain either 0 or 2 photons
221     PhotonOArr *GoodPhotons = new PhotonOArr;
222 bendavid 1.8 GoodPhotons->SetName(fGoodPhotonsName);
223     GoodPhotons->SetOwner(kTRUE);
224 fabstoec 1.29
225    
226     VertexOArr *ChosenVtx = new VertexOArr;
227     ChosenVtx->SetName(fChosenVtxName);
228     ChosenVtx->SetOwner(kTRUE);
229 fabstoec 1.30
230 bendavid 1.8 // add to event for other modules to use
231 paus 1.27 AddObjThisEvt(GoodPhotons);
232 fabstoec 1.29 AddObjThisEvt(ChosenVtx);
233 paus 1.27
234 fabstoec 1.30 //AddObjThisEvt(ChosenVtx);
235    
236 paus 1.27 if (fPhotons->GetEntries()<2)
237     return;
238    
239 fabstoec 1.1 LoadEventObject(fElectronName, fElectrons);
240 paus 1.27 LoadEventObject(fGoodElectronName, fGoodElectrons);
241 fabstoec 1.1 LoadEventObject(fConversionName, fConversions);
242 bendavid 1.36 if (fUseSingleLegConversions) LoadEventObject(fPFConversionName, fPFConversions);
243 fabstoec 1.1 LoadEventObject(fTrackBranchName, fTracks);
244     LoadEventObject(fPileUpDenName, fPileUpDen);
245 paus 1.27 LoadEventObject(fPVName, fPV);
246 fabstoec 1.1 LoadEventObject(fBeamspotName, fBeamspot);
247     LoadEventObject(fPFCandName, fPFCands);
248 bendavid 1.31 LoadEventObject(fJetsName, fJets);
249     LoadEventObject(fPFMetName, fPFMet);
250 fabstoec 1.1
251 bendavid 1.31 if (!fIsData) {
252     LoadBranch(fMCParticleName);
253     }
254 mingyang 1.57
255 paus 1.27 // ------------------------------------------------------------
256 fabstoec 1.1 // load event based information
257 paus 1.27 Float_t rho = -99.;
258     if (fPileUpDen->GetEntries() > 0)
259     rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
260 bendavid 1.33 Float_t rho2012 = -99;
261     if (fPileUpDen->At(0)->RhoKt6PFJets()>0.) rho2012 = fPileUpDen->At(0)->RhoKt6PFJets();
262     else rho2012 = fPileUpDen->At(0)->Rho();
263 fabstoec 1.1 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
264    
265 fabstoec 1.50 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 paus 1.27 // ------------------------------------------------------------
286 fabstoec 1.1 // Get Event header for Run info etc.
287 paus 1.27 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 mingyang 1.57
293 fabstoec 1.1 // ------------------------------------------------------------
294 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
295     // printf("evtNum:%d\n",evtNum);
296     // }
297 fabstoec 1.1 // 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 paus 1.27
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 fabstoec 1.1 const Photon *ph = fPhotons->At(i);
305    
306 paus 1.27 if (ph->SCluster()->AbsEta()>= fPhotonEtaMax ||
307     (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566))
308     continue;
309 fabstoec 1.49
310 paus 1.27 if (ph->Et() < fPhotonPtMin)
311     continue;
312 fabstoec 1.49
313 paus 1.27 if (ph->HadOverEm() > 0.15)
314     continue;
315    
316     if (ph->IsEB()) {
317 fabstoec 1.49 if (ph->CoviEtaiEta() > 0.015)
318 paus 1.27 continue;
319     }
320     else {
321 fabstoec 1.49 if (ph->CoviEtaiEta() > 0.035)
322 paus 1.27 continue;
323     }
324 fabstoec 1.49
325 paus 1.27 // photon passes the preselection
326     ph->Mark(); // mark for later skimming
327 fabstoec 1.1 preselPh->Add(ph);
328 fabstoec 1.2 }
329 fabstoec 1.49
330 paus 1.27 if (preselPh->GetEntries()<2) {
331     this->SkipEvent();
332     delete preselPh;
333     return;
334     }
335 fabstoec 1.49
336 bendavid 1.36 //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 fabstoec 1.54 //std::cout<<" conv"<<iconv<<" mom = "<<c->Mom().X()<<" "<<c->Mom().Y()<<" "<<c->Mom().Z()<<" "<<std::endl;
352    
353 bendavid 1.37 DecayParticle *conv = new DecayParticle(*c);
354 bendavid 1.36 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 fabstoec 1.49 }
364 bendavid 1.36 }
365    
366 bendavid 1.31
367     float higgspt = -99.;
368     float higgsz = -99.;
369     float higgsmass = -99.;
370    
371     if (!fIsData) FindHiggsPtAndZ(higgspt,higgsz,higgsmass);
372    
373 paus 1.27 // second loop: sort & assign the right categories..
374 fabstoec 1.2 preselPh->Sort();
375 paus 1.27 for (unsigned int iPh = 0; iPh <preselPh->GetEntries(); ++iPh) {
376 fabstoec 1.2 const Photon* ph = preselPh->At(iPh);
377 fabstoec 1.1 preselCat.push_back(PhotonTools::CiCBaseLineCat(ph));
378     }
379 paus 1.27
380 fabstoec 1.1 // ------------------------------------------------------------
381     // compute how many pairs there are ...
382     unsigned int numPairs = 0;
383 paus 1.27 if (preselPh->GetEntries() > 0)
384     numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
385 mingyang 1.61 //printf("ming sync check numPairs:%d\n",numPairs);
386 fabstoec 1.54
387 fabstoec 1.1 // ... and create all possible pairs of pre-selected photons
388 paus 1.27 std::vector<unsigned int> idx1st;
389     std::vector<unsigned int> idx2nd;
390 fabstoec 1.1 std::vector<PhotonTools::CiCBaseLineCats> cat1st;
391     std::vector<PhotonTools::CiCBaseLineCats> cat2nd;
392 fabstoec 1.54
393 paus 1.27 // ... this will be used to store whether a given pair passes the cuts
394 fabstoec 1.1 std::vector<bool> pairPasses;
395 fabstoec 1.54
396 paus 1.27 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 fabstoec 1.1 }
403     }
404     }
405    
406 paus 1.27 // ------------------------------------------------------------
407 fabstoec 1.1 // array to store the index of 'chosen Vtx' for each pair
408 bendavid 1.20 // const Vertex** theVtx = new const Vertex*[numPairs]; // holds the 'chosen' Vtx for each Pair
409 paus 1.27 // Photon** fixPh1st = new Photon*[numPairs]; // holds the 1st Photon for each Pair
410 bendavid 1.20 // Photon** fixPh2nd = new Photon*[numPairs]; // holds the 2nd photon for each Pair
411 paus 1.27
412 bendavid 1.20 std::vector<const Vertex*> theVtx; // holds the 'chosen' Vtx for each Pair
413 paus 1.27 std::vector<Photon*> fixPh1st; // holds the 1st Photon for each Pair
414 bendavid 1.20 std::vector<Photon*> fixPh2nd; // holds the 2nd photon for each Pair
415 fabstoec 1.54
416 bendavid 1.20 theVtx.reserve(numPairs);
417     fixPh1st.reserve(numPairs);
418     fixPh2nd.reserve(numPairs);
419 fabstoec 1.1
420     // store pair-indices for pairs passing the selection
421     std::vector<unsigned int> passPairs;
422     passPairs.resize(0);
423 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
424     // printf("check 1\n");
425     //}
426 paus 1.27 // ------------------------------------------------------------
427     // Loop over all Pairs and do the 'incredible machine' running....
428 mingyang 1.58 std::vector<int> leptag;
429     leptag.resize(numPairs);
430    
431 paus 1.27 for (unsigned int iPair = 0; iPair < numPairs; ++iPair) {
432 fabstoec 1.1 // first we need a hard copy of the incoming photons
433 bendavid 1.20 fixPh1st.push_back(new Photon(*preselPh->At(idx1st[iPair])));
434     fixPh2nd.push_back(new Photon(*preselPh->At(idx2nd[iPair])));
435 fabstoec 1.1 // we also store the category, so we don't have to ask all the time...
436     cat1st.push_back(preselCat[idx1st[iPair]]);
437 paus 1.27 cat2nd.push_back(preselCat[idx2nd[iPair]]);
438    
439 fabstoec 1.47 // do showershape scaling.... (fab: outsourced to PhotonTools)
440 mingyang 1.41 if (fDoShowerShapeScaling && !fIsData) {
441 fabstoec 1.47 PhotonTools::ScalePhotonShowerShapes(fixPh1st[iPair],fSSType);
442     PhotonTools::ScalePhotonShowerShapes(fixPh2nd[iPair],fSSType);
443 bendavid 1.19 }
444 mingyang 1.41
445 mingyang 1.52 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 fabstoec 1.47
456 fabstoec 1.1 // now we dicide if we either scale (Data) or Smear (MC) the Photons
457     if (fIsData) {
458 paus 1.27 if (fDoDataEneCorr) {
459     // starting with scale = 1.
460     double scaleFac1 = 1.;
461     double scaleFac2 = 1.;
462 mingyang 1.52
463 bendavid 1.9 //eta-dependent corrections
464 fabstoec 1.47
465 paus 1.27 // checking the run Rangees ...
466     Int_t runRange = FindRunRangeIdx(runNumber);
467 mingyang 1.52
468 paus 1.27 if(runRange > -1) {
469 mingyang 1.52 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 paus 1.27 }
477 bendavid 1.48 else {
478     printf("Error: Run Range not found for data energy scale correction\n");
479     assert(0);
480     }
481 paus 1.27 PhotonTools::ScalePhoton(fixPh1st[iPair], scaleFac1);
482     PhotonTools::ScalePhoton(fixPh2nd[iPair], scaleFac2);
483 fabstoec 1.1 }
484 paus 1.27 }
485 mingyang 1.57
486 fabstoec 1.56 if (fDoMCSmear) {
487 mingyang 1.52
488 fabstoec 1.56 double width1 = 0.;
489     double width2 = 0.;
490    
491 mingyang 1.52 if(f2012HCP){
492     width1 = GetMCSmearFacHCP(escalecat1);
493     width2 = GetMCSmearFacHCP(escalecat2);
494     }else {
495     width1 = GetMCSmearFac(escalecat1);
496     width2 = GetMCSmearFac(escalecat2);
497     }
498 fabstoec 1.56
499 bendavid 1.9 if (!fIsData) {
500     // get the seed to do deterministic smearing...
501     UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
502 paus 1.27 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 bendavid 1.9 // get the smearing for MC photons..
507 paus 1.27 PhotonTools::SmearPhoton(fixPh1st[iPair], fRng, width1, seed1);
508     PhotonTools::SmearPhoton(fixPh2nd[iPair], fRng, width2, seed2);
509 fabstoec 1.1 }
510 fabstoec 1.56
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 mingyang 1.61
522     }
523 fabstoec 1.56
524 bendavid 1.9 PhotonTools::SmearPhotonError(fixPh1st[iPair], width1);
525     PhotonTools::SmearPhotonError(fixPh2nd[iPair], width2);
526 fabstoec 1.1 }
527 mingyang 1.52
528     // lepton tag for this pair -- ming
529 mingyang 1.55 leptag[iPair] = -1;
530 mingyang 1.52 if(fApplyLeptonTag){
531 mingyang 1.55 leptag[iPair] = 0;
532 mingyang 1.52 if ( fLeptonTagMuons->GetEntries() > 0 ) {
533 mingyang 1.61 //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 mingyang 1.52 if( (MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh1st[iPair]) >= 1.0) &&
537 mingyang 1.61 (MathUtils::DeltaR(*fLeptonTagMuons->At(0),*fixPh2nd[iPair]) >= 1.0)
538     ){
539 mingyang 1.58 leptag[iPair] = 2;
540 mingyang 1.52 }
541 mingyang 1.58 }
542 mingyang 1.61
543 mingyang 1.59 if(leptag[iPair] <1){
544 mingyang 1.58 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 mingyang 1.61 (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 mingyang 1.60 /*int ph1passeveto=1;
555 mingyang 1.59 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 mingyang 1.60 }*/
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 mingyang 1.59 }
583 mingyang 1.60
584 mingyang 1.58 }
585 mingyang 1.52 }
586     }
587     }
588 mingyang 1.58
589 bendavid 1.19 //probability that selected vertex is the correct one
590     Double_t vtxProb = 1.0;
591 fabstoec 1.49
592 fabstoec 1.1 // store the vertex for this pair
593 paus 1.27 switch (fVtxSelType) {
594 fabstoec 1.49
595 fabstoec 1.1 case kStdVtxSelection:
596     theVtx[iPair] = fPV->At(0);
597     break;
598 fabstoec 1.49
599 fabstoec 1.1 case kCiCVtxSelection:
600 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
601 bendavid 1.36 &vtxconversions,kFALSE,vtxProb);
602 fabstoec 1.1 break;
603 fabstoec 1.49
604 bendavid 1.19 case kCiCMVAVtxSelection:
605 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
606 bendavid 1.36 &vtxconversions,kTRUE,vtxProb);
607 paus 1.27 break;
608 fabstoec 1.49
609 fabstoec 1.1 case kMITVtxSelection:
610     // need PFCandidate Collection
611 paus 1.27 theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp,
612     mithep::FourVector((fixPh1st[iPair]->Mom()+
613     fixPh2nd[iPair]->Mom())));
614 fabstoec 1.1 break;
615 fabstoec 1.49
616 bendavid 1.31 case kMetSigVtxSelection: {
617 fabstoec 1.40 // need PFCandidate Collection, otherwise use 0
618     if( !fJets || !fPFCands ) {
619     theVtx[iPair] = fPV->At(0);
620     break;
621     }
622 bendavid 1.31
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 bendavid 1.38 const Vertex *v = fPV->At(ivtx);
638    
639 fabstoec 1.40 // 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 bendavid 1.38
648 fabstoec 1.40 Met mmet = fMVAMet.GetMet( false,
649     0.,0.,0.,
650 ceballos 1.46 0.,0.,0.,
651 fabstoec 1.40 fPFMet->At(0),
652     &pfcands,fPV->At(ivtx),fPV, fPileUpDen->At(0)->Rho(),
653     &pfjets,
654     int(fPV->GetEntries()),
655     kFALSE);
656 bendavid 1.38
657 bendavid 1.31 ThreeVector fullmet(mmet.Px() - fixPh1st[iPair]->Px() - fixPh2nd[iPair]->Px(),
658     mmet.Py() - fixPh1st[iPair]->Py() - fixPh2nd[iPair]->Py(),
659     0.);
660    
661 fabstoec 1.40 // ThreeVector fullmet(mmet.Px(),
662     // mmet.Py(),
663     // 0.);
664 bendavid 1.31
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 fabstoec 1.40
669 bendavid 1.31 //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 fabstoec 1.40
683 bendavid 1.31 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
684    
685 bendavid 1.38
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 fabstoec 1.49 if ( sumptsq < 10.0 ) metsig = -99;
696 bendavid 1.38
697 bendavid 1.31 if (metsig<minsig && metsig>0.) {
698     minsig = metsig;
699     theVtx[iPair] = fPV->At(ivtx);
700     }
701    
702 bendavid 1.38 printf("ivtx = %i, sumptsq = %5f, met = %5f, metsig = %5f, dzgen = %5f\n",ivtx,sumptsq, fullmet.Rho(), metsig,fPV->At(ivtx)->Z()-higgsz);
703 bendavid 1.31
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 fabstoec 1.40
718 bendavid 1.31 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 fabstoec 1.40
725 bendavid 1.31 }
726     break;
727 fabstoec 1.1 default:
728     theVtx[iPair] = fPV->At(0);
729 paus 1.27 }
730 mingyang 1.57
731 mingyang 1.53 /*
732 mingyang 1.52 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 mingyang 1.53 */
759    
760 bendavid 1.8 //set PV ref in photons
761     fixPh1st[iPair]->SetPV(theVtx[iPair]);
762     fixPh2nd[iPair]->SetPV(theVtx[iPair]);
763 bendavid 1.19 fixPh1st[iPair]->SetVtxProb(vtxProb);
764     fixPh2nd[iPair]->SetVtxProb(vtxProb);
765 bendavid 1.48
766     ThreeVector vtxpos = theVtx[iPair]->Position();
767    
768 fabstoec 1.1 // fix the kinematics for both events
769 bendavid 1.48 FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(vtxpos);
770     FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(vtxpos);
771 fabstoec 1.1 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 bendavid 1.19 double pairmass = (fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M();
775 paus 1.27
776 bendavid 1.20 double leadptcut = fLeadingPtMin;
777     double trailptcut = fTrailingPtMin;
778 paus 1.27
779 bendavid 1.20 if (fixPh2nd[iPair]->Pt() > fixPh1st[iPair]->Pt()) {
780     leadptcut = fTrailingPtMin;
781     trailptcut = fLeadingPtMin;
782     }
783 paus 1.27
784 fabstoec 1.49
785 bendavid 1.20 if (fRelativePtCuts) {
786     leadptcut = leadptcut*pairmass;
787     trailptcut = trailptcut*pairmass;
788     }
789 fabstoec 1.49
790 mingyang 1.39
791 bendavid 1.21 //compute id bdt values
792 fabstoec 1.49 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 fabstoec 1.50 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 fabstoec 1.49 }
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 paus 1.27
825 bendavid 1.20 //printf("applying id\n");
826 fabstoec 1.56
827 fabstoec 1.1 // check if both photons pass the CiC selection
828     // FIX-ME: Add other possibilities....
829     bool pass1 = false;
830     bool pass2 = false;
831 fabstoec 1.40
832     switch (fPhSelType) {
833 paus 1.27
834 fabstoec 1.40 // --------------------------------------------------------------------
835     // trivial (no) selection with only pt-cuts
836 fabstoec 1.1 case kNoPhSelection:
837 paus 1.27 pass1 = (fixPh1st[iPair]->Pt() > leadptcut );
838     pass2 = (fixPh2nd[iPair]->Pt() > trailptcut);
839 fabstoec 1.1 break;
840 fabstoec 1.40
841     // --------------------------------------------------------------------
842     // CiC4 Selection as used for the 2011 7TeV Baseline analysis
843 fabstoec 1.1 case kCiCPhSelection:
844 paus 1.27 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 fabstoec 1.1 break;
850 fabstoec 1.40
851     // --------------------------------------------------------------------
852     // PF-CiC4 Selection, as used in the 2012 8TeV Baseline analysis
853 bendavid 1.33 case kCiCPFPhSelection:
854    
855     // loose preselection for mva
856     pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
857 mingyang 1.60 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh1st[iPair],fElectrons,fConversions,bsp,
858 mingyang 1.61 //PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
859 bendavid 1.33 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
860     fInvertElectronVeto);
861     if (pass1)
862     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
863 mingyang 1.60 PhotonTools::PassSinglePhotonPreselPFISONoEcal(fixPh2nd[iPair],fElectrons,fConversions,bsp,
864 mingyang 1.61 //PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
865 bendavid 1.33 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 fabstoec 1.40
875     // --------------------------------------------------------------------
876     // MVA selection
877 mingyang 1.15 case kMVAPhSelection://MVA
878 fabstoec 1.49
879 paus 1.27 pass1 = fixPh1st[iPair]->Pt()>leadptcut &&
880     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
881     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
882 fabstoec 1.49 ( ( 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 paus 1.27 if (pass1)
891     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
892     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
893     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
894 fabstoec 1.49 ( ( 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 paus 1.27
900 mingyang 1.15 break;
901 fabstoec 1.40
902     // --------------------------------------------------------------------
903     // MIT proposed pre-selection as used in the 2011 7TeV MVA analyses
904 fabstoec 1.1 case kMITPhSelection:
905 bendavid 1.20 // loose preselection for mva
906 paus 1.27 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
907     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
908 bendavid 1.33 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
909 paus 1.27 fInvertElectronVeto);
910     if (pass1)
911     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
912     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
913 bendavid 1.33 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
914 paus 1.27 fInvertElectronVeto);
915    
916 fabstoec 1.1 break;
917 fabstoec 1.40
918     // --------------------------------------------------------------------
919     // updated (PF absed) pre-selection as used in the 2012 8TeV MVA/Baseline analyses
920 bendavid 1.33 case kMITPFPhSelection:
921     // loose preselection for mva
922     pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
923     PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
924 fabstoec 1.49 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
925     fInvertElectronVeto);
926    
927 bendavid 1.33 if (pass1)
928     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
929     PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
930 fabstoec 1.49 fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
931     fInvertElectronVeto);
932    
933 mingyang 1.60 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 bendavid 1.33 break;
951 mingyang 1.60
952 fabstoec 1.1 default:
953     pass1 = true;
954     pass2 = true;
955     }
956 fabstoec 1.49
957 bendavid 1.8 //match to good electrons if requested
958 fabstoec 1.49 // if (fInvertElectronVeto) {
959     // pass1 &= !PhotonTools::PassElectronVeto(fixPh1st[iPair],fGoodElectrons);
960     // pass2 &= !PhotonTools::PassElectronVeto(fixPh2nd[iPair],fGoodElectrons);
961     // }
962 paus 1.27 // finally, if both Photons pass the selections, add the pair to the passing pairs
963     if (pass1 && pass2)
964     passPairs.push_back(iPair);
965 fabstoec 1.1 }
966 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
967     // printf("check 2\n");
968     // }
969 fabstoec 1.1 // ---------------------------------------------------------------
970 paus 1.18 // ... we're almost done, stay focused...
971 fabstoec 1.1 // loop over all passing pairs and find the one with the highest sum Et
972 paus 1.27 const Vertex* vtx = NULL;
973 fabstoec 1.1 Photon* phHard = NULL;
974     Photon* phSoft = NULL;
975    
976 fabstoec 1.29 // not used at all....
977     //PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
978     //PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
979 paus 1.27
980 fabstoec 1.1 double maxSumEt = 0.;
981 mingyang 1.55 int maxleptag = 0;
982 paus 1.27 for (unsigned int iPair=0; iPair<passPairs.size(); ++iPair) {
983 mingyang 1.58 double sumEt = fixPh1st[passPairs[iPair]]->Et()+fixPh2nd[passPairs[iPair]]->Et();
984 mingyang 1.61 //if (((sumEt > maxSumEt) && leptag[iPair] == maxleptag) || (leptag[iPair] > maxleptag)) {
985     if (((sumEt > maxSumEt) && leptag[passPairs[iPair]] == maxleptag) || (leptag[passPairs[iPair]] > maxleptag)) {
986 fabstoec 1.1 maxSumEt = sumEt;
987 mingyang 1.61 maxleptag = leptag[passPairs[iPair]];
988 paus 1.27 phHard = fixPh1st[passPairs[iPair]];
989     phSoft = fixPh2nd[passPairs[iPair]];
990 fabstoec 1.29 //catPh1 = cat1st [passPairs[iPair]];
991     //catPh2 = cat2nd [passPairs[iPair]];
992 paus 1.27 vtx = theVtx[iPair];
993 fabstoec 1.1 }
994     }
995 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
996     // printf("check 3\n");
997     // }
998 paus 1.27 for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
999 bendavid 1.20 if (fixPh1st[iPair]!=phHard) delete fixPh1st[iPair];
1000     if (fixPh2nd[iPair]!=phSoft) delete fixPh2nd[iPair];
1001     }
1002 paus 1.27
1003 fabstoec 1.1 // ---------------------------------------------------------------
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 fabstoec 1.29 // we also store the chosen Vtx, so later modules can use it
1012 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
1013     // printf("check 4\n");
1014     //}
1015 fabstoec 1.30 Vertex* chosenVtx = NULL;
1016     if ( vtx ) {
1017     chosenVtx = new Vertex( *vtx );
1018     ChosenVtx->AddOwned( chosenVtx );
1019     }
1020 mingyang 1.57 //printf("ChosenVtx->GetEntries():%d\n",ChosenVtx->GetEntries()ChosenVtx->GetEntries());
1021     //if(evtNum==9009 || evtNum==9010){
1022     // printf("check 5\n");
1023     //}
1024 fabstoec 1.1 // sort according to pt
1025     GoodPhotons->Sort();
1026 mingyang 1.57
1027 fabstoec 1.1 // delete auxiliary photon collection...
1028     delete preselPh;
1029 mingyang 1.57
1030 bendavid 1.20 //delete[] theVtx;
1031 mingyang 1.57 //if(evtNum==9009 || evtNum==9010){
1032     // printf(" GoodPhotons->GetEntries():%d\n",GoodPhotons->GetEntries());
1033     // printf("check 6\n");
1034     //}
1035 fabstoec 1.1 return;
1036     }
1037    
1038 paus 1.27 //---------------------------------------------------------------------------------------------------
1039 fabstoec 1.1 void PhotonPairSelector::SlaveBegin()
1040     {
1041 paus 1.27 // Run startup code on the computer (slave) doing the actual analysis. Here, we just request the
1042     // photon collection branch.
1043 fabstoec 1.1
1044 mingyang 1.52 if( fApplyLeptonTag ) {
1045     ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
1046     ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
1047     }
1048    
1049 paus 1.27 // 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 bendavid 1.36 if (fUseSingleLegConversions) ReqEventObject(fPFConversionName, fPFConversions, true);
1058 paus 1.27 ReqEventObject(fBeamspotName, fBeamspot, true);
1059     ReqEventObject(fPFCandName, fPFCands, true);
1060 bendavid 1.31 ReqEventObject(fJetsName, fJets, false);
1061     ReqEventObject(fPFMetName, fPFMet, true);
1062 fabstoec 1.1 if (!fIsData) {
1063 bendavid 1.31 //ReqBranch(fPileUpName, fPileUp);
1064 fabstoec 1.1 ReqBranch(fMCParticleName, fMCParticles);
1065     }
1066 fabstoec 1.49
1067 paus 1.27 // determine photon selection type
1068     if (fPhotonSelType.CompareTo("CiCSelection") == 0)
1069 fabstoec 1.1 fPhSelType = kCiCPhSelection;
1070 bendavid 1.33 else if (fPhotonSelType.CompareTo("CiCPFSelection") == 0)
1071     fPhSelType = kCiCPFPhSelection;
1072 mingyang 1.15 else if (fPhotonSelType.CompareTo("MVASelection") == 0) //MVA
1073     fPhSelType = kMVAPhSelection;
1074 paus 1.27 else if (fPhotonSelType.CompareTo("MITSelection") == 0)
1075 fabstoec 1.1 fPhSelType = kMITPhSelection;
1076 bendavid 1.33 else if (fPhotonSelType.CompareTo("MITPFSelection") == 0)
1077     fPhSelType = kMITPFPhSelection;
1078 mingyang 1.60 else if (fPhotonSelType.CompareTo("MITPFSelectionNoEcal") == 0)
1079     fPhSelType = kMITPFPhSelectionNoEcal;
1080 bendavid 1.33 else if (fPhotonSelType.CompareTo("NoSelection") == 0)
1081 fabstoec 1.1 fPhSelType = kNoPhSelection;
1082 bendavid 1.33 else {
1083     std::cerr<<" Photon Seclection "<<fPhotonSelType<<" not implemented."<<std::endl;
1084     return;
1085     }
1086 fabstoec 1.1
1087 bendavid 1.42 if (fVertexSelType.CompareTo("CiCSelection") == 0) {
1088 fabstoec 1.1 fVtxSelType = kCiCVtxSelection;
1089 bendavid 1.42 fVtxTools.InitP(1);
1090     }
1091 paus 1.27 else if (fVertexSelType.CompareTo("MITSelection") == 0)
1092 fabstoec 1.1 fVtxSelType = kMITVtxSelection;
1093 bendavid 1.42 else if (fVertexSelType.CompareTo("CiCMVASelection") == 0) {
1094     fVtxSelType = kCiCMVAVtxSelection;
1095     fVtxTools.InitP(1);
1096     }
1097     else if (fVertexSelType.CompareTo("CiCMVA2012Selection") == 0) {
1098 paus 1.27 fVtxSelType = kCiCMVAVtxSelection;
1099 bendavid 1.42 fVtxTools.InitP(2);
1100     }
1101 bendavid 1.31 else if (fVertexSelType.CompareTo("MetSigSelection") == 0)
1102     fVtxSelType = kMetSigVtxSelection;
1103 fabstoec 1.28 else if (fVertexSelType.CompareTo("ZeroVtxSelection") == 0)
1104 paus 1.27 fVtxSelType = kStdVtxSelection;
1105 fabstoec 1.28 else {
1106     std::cerr<<" Vertex Seclection "<<fVertexSelType<<" not implemented."<<std::endl;
1107     return;
1108     }
1109 mingyang 1.34
1110     if (fIdMVAType.CompareTo("2011IdMVA") == 0)
1111 fabstoec 1.49 fIdType = MVATools::k2011IdMVA;
1112 mingyang 1.34 else if (fIdMVAType.CompareTo("2012IdMVA_globe") == 0)
1113 fabstoec 1.49 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 mingyang 1.34 else {
1121     std::cerr<<" Id MVA "<<fIdMVAType<<" not implemented."<<std::endl;
1122     return;
1123     }
1124    
1125 fabstoec 1.47
1126     if (fShowerShapeType.CompareTo("None") == 0)
1127     fSSType = PhotonTools::kNoShowerShapeScaling;
1128     else if (fShowerShapeType.CompareTo("2011ShowerShape") == 0)
1129     fSSType = PhotonTools::k2011ShowerShape;
1130 mingyang 1.41 else if (fShowerShapeType.CompareTo("2012ShowerShape") == 0)
1131 fabstoec 1.47 fSSType = PhotonTools::k2012ShowerShape;
1132 mingyang 1.41 else {
1133     std::cerr<<"shower shape scale "<<fShowerShapeType<<" not implemented."<<std::endl;
1134     return;
1135     }
1136    
1137 paus 1.27 if (fIsData)
1138 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
1139 paus 1.27 else
1140 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
1141 mingyang 1.15
1142 paus 1.27 printf("initialize photon pair selector\n");
1143 mingyang 1.15
1144 fabstoec 1.49 // ----------------------------------------------------------------
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 bendavid 1.31
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 mingyang 1.15
1182 fabstoec 1.1 }
1183    
1184     // ----------------------------------------------------------------------------------------
1185     // some helpfer functions....
1186 paus 1.27 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
1187     {
1188     pt = -999.;
1189 fabstoec 1.1 decayZ = -999.;
1190 paus 1.27 mass = -999.;
1191 bendavid 1.31
1192 fabstoec 1.1 // 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 bendavid 1.8 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 fabstoec 1.1 }
1202 paus 1.27
1203 fabstoec 1.1 return;
1204 paus 1.27 }
1205 fabstoec 1.1
1206 paus 1.27 //---------------------------------------------------------------------------------------------------
1207     Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run)
1208     {
1209     // this routine looks for the idx of the run-range
1210 fabstoec 1.1 Int_t runRange=-1;
1211 paus 1.27 for (UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
1212     if (run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
1213 fabstoec 1.1 runRange = (Int_t) iRun;
1214     return runRange;
1215     }
1216     }
1217     return runRange;
1218     }
1219    
1220 paus 1.27 //---------------------------------------------------------------------------------------------------
1221     Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::eScaleCats cat)
1222     {
1223     switch (cat) {
1224 bendavid 1.19 case PhotonTools::kEBhighEtaGold:
1225     return fDataEnCorr_EBhighEta_hR9[runRange];
1226     case PhotonTools::kEBhighEtaBad:
1227     return fDataEnCorr_EBhighEta_lR9[runRange];
1228 bendavid 1.20 case PhotonTools::kEBlowEtaGoldCenter:
1229     return fDataEnCorr_EBlowEta_hR9central[runRange];
1230     case PhotonTools::kEBlowEtaGoldGap:
1231 paus 1.27 return fDataEnCorr_EBlowEta_hR9gap[runRange];
1232 bendavid 1.19 case PhotonTools::kEBlowEtaBad:
1233 paus 1.27 return fDataEnCorr_EBlowEta_lR9[runRange];
1234 bendavid 1.19 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 paus 1.27 return fDataEnCorr_EElowEta_lR9[runRange];
1242 fabstoec 1.1 default:
1243     return 1.;
1244     }
1245     }
1246    
1247 paus 1.27 //---------------------------------------------------------------------------------------------------
1248 fabstoec 1.56 Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::eScaleCats cat, bool useSpecialSmear)
1249 paus 1.27 {
1250 fabstoec 1.56
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 fabstoec 1.1 }
1298 fabstoec 1.56
1299 fabstoec 1.1 }
1300    
1301 mingyang 1.52 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 fabstoec 1.56 Double_t PhotonPairSelector::GetMCSmearFacHCP(PhotonTools::eScaleCats cat, bool useSpecialSmear)
1331 mingyang 1.52 {
1332 fabstoec 1.56
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 mingyang 1.52 }
1388     }
1389    
1390 paus 1.27 //---------------------------------------------------------------------------------------------------
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 fabstoec 1.1
1397     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
1398     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
1399 paus 1.27
1400     if (ph1IsEB && ph2IsEB)
1401     return (ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
1402    
1403     return (ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
1404 fabstoec 1.1 }
1405 bendavid 1.4