ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.39
Committed: Tue May 29 21:21:46 2012 UTC (12 years, 11 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028a
Changes since 1.38: +8 -3 lines
Log Message:
add new scale

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    
26     //--------------------------------------------------------------------------------------------------
27 paus 1.27 PhotonPairSelector::PhotonPairSelector(const char *name, const char *title) :
28 fabstoec 1.1 // Base Module...
29 paus 1.27 BaseMod (name,title),
30 fabstoec 1.1 // define all the Branches to load
31 paus 1.27 fPhotonBranchName (Names::gkPhotonBrn),
32     fElectronName (Names::gkElectronBrn),
33     fGoodElectronName (Names::gkElectronBrn),
34     fConversionName (Names::gkMvfConversionBrn),
35 bendavid 1.36 fPFConversionName ("PFPhotonConversions"),
36 paus 1.27 fTrackBranchName (Names::gkTrackBrn),
37     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
38     fPVName (Names::gkPVBeamSpotBrn),
39     fBeamspotName (Names::gkBeamSpotBrn),
40     fPFCandName (Names::gkPFCandidatesBrn),
41 fabstoec 1.1 // MC specific stuff...
42 paus 1.27 fMCParticleName (Names::gkMCPartBrn),
43     fPileUpName (Names::gkPileupInfoBrn),
44 bendavid 1.31 fJetsName (Names::gkPFJetBrn),
45     fPFMetName ("PFMet"),
46 paus 1.27 fGoodPhotonsName (ModNames::gkGoodPhotonsName),
47 fabstoec 1.29 fChosenVtxName ("HggChosenVtx"),
48 fabstoec 1.1 // ----------------------------------------
49     // Selection Types
50 paus 1.27 fPhotonSelType ("NoSelection"),
51     fVertexSelType ("StdSelection"),
52     fPhSelType (kNoPhSelection),
53     fVtxSelType (kStdVtxSelection),
54 mingyang 1.34 //-----------------------------------------
55     // Id Types
56     fIdMVAType ("2011IdMVA"),
57     fIdType (k2011IdMVA),
58 fabstoec 1.1 // ----------------------------------------
59 paus 1.27 fPhotonPtMin (20.0),
60     fPhotonEtaMax (2.5),
61     fLeadingPtMin (100.0/3.0),
62     fTrailingPtMin (100.0/4.0),
63     fIsData (false),
64     fPhotonsFromBranch (true),
65     fPVFromBranch (true),
66     fGoodElectronsFromBranch (kTRUE),
67 bendavid 1.36 fUseSingleLegConversions (kFALSE),
68 fabstoec 1.1 // ----------------------------------------
69     // collections....
70 paus 1.27 fPhotons (0),
71     fElectrons (0),
72     fConversions (0),
73 bendavid 1.36 fPFConversions (0),
74 paus 1.27 fTracks (0),
75     fPileUpDen (0),
76     fPV (0),
77     fBeamspot (0),
78     fPFCands (0),
79     fMCParticles (0),
80     fPileUp (0),
81 bendavid 1.31 fJets (0),
82     fPFMet (0),
83 fabstoec 1.1 // ---------------------------------------
84 paus 1.27 fDataEnCorr_EBlowEta_hR9central(0.),
85     fDataEnCorr_EBlowEta_hR9gap (0.),
86     fDataEnCorr_EBlowEta_lR9 (0.),
87     fDataEnCorr_EBhighEta_hR9 (0.),
88     fDataEnCorr_EBhighEta_lR9 (0.),
89     fDataEnCorr_EElowEta_hR9 (0.),
90     fDataEnCorr_EElowEta_lR9 (0.),
91     fDataEnCorr_EEhighEta_hR9 (0.),
92     fDataEnCorr_EEhighEta_lR9 (0.),
93     fRunStart (0),
94     fRunEnd (0),
95     fMCSmear_EBlowEta_hR9central (0.),
96     fMCSmear_EBlowEta_hR9gap (0.),
97     fMCSmear_EBlowEta_lR9 (0.),
98     fMCSmear_EBhighEta_hR9 (0.),
99     fMCSmear_EBhighEta_lR9 (0.),
100     fMCSmear_EElowEta_hR9 (0.),
101     fMCSmear_EElowEta_lR9 (0.),
102     fMCSmear_EEhighEta_hR9 (0.),
103     fMCSmear_EEhighEta_lR9 (0.),
104 fabstoec 1.1 // ---------------------------------------
105 paus 1.27 fRng (new TRandom3()),
106     fPhFixString ("4_2"),
107     fEtaCorrections (0),
108 fabstoec 1.1 // ---------------------------------------
109 paus 1.27 fDoDataEneCorr (true),
110     fDoMCSmear (true),
111     fDoVtxSelection (true),
112     fApplyEleVeto (true),
113     fInvertElectronVeto (kFALSE),
114 mingyang 1.15 //MVA
115 mingyang 1.34 fVariableType_2011 (10),
116     fEndcapWeights_2011 (gSystem->Getenv("CMSSW_BASE")+
117     TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
118     TString("Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
119 mingyang 1.32 TString("weights.xml")),
120 mingyang 1.34 fBarrelWeights_2011 (gSystem->Getenv("CMSSW_BASE")+
121     TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
122     TString("Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
123     TString("weights.xml")),
124     fVariableType_2012_globe (1201),
125     fEndcapWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
126     TString("/src/MitPhysics/data/")+
127     TString("TMVA_EEpf_BDT_globe.")+
128     TString("weights.xml")),
129     fBarrelWeights_2012_globe (gSystem->Getenv("CMSSW_BASE")+
130     TString("/src/MitPhysics/data/")+
131     TString("TMVA_EBpf_BDT_globe.")+
132 mingyang 1.32 TString("weights.xml")),
133 paus 1.27 fbdtCutBarrel (0.0744), //cuts give same eff (relative to presel) with cic
134     fbdtCutEndcap (0.0959), //cuts give same eff (relative to presel) with cic
135     fDoMCR9Scaling (kFALSE),
136     fMCR9ScaleEB (1.0),
137     fMCR9ScaleEE (1.0),
138     fDoMCSigIEtaIEtaScaling (kFALSE),
139     fDoMCWidthScaling (kFALSE),
140 mingyang 1.39 fDoMCNewScaling (kFALSE),
141 mingyang 1.35 fMCSigIEtaIEtaScaleEB (0.87),
142     fMCSigIEtaIEtaShiftEB (0.0011),
143     fMCSigIEtaIEtaScaleEE (0.99),
144     fMCSigIEtaIEtaShiftEE (0),
145     fMCWidthScale (0.99),
146     fMCWidthShift (0),
147     //
148 paus 1.27 fDoMCErrScaling (kFALSE),
149     fMCErrScaleEB (1.0),
150     fMCErrScaleEE (1.0),
151     fRelativePtCuts (kFALSE)
152 fabstoec 1.1 {
153     // Constructor.
154     }
155    
156 paus 1.27 PhotonPairSelector::~PhotonPairSelector()
157     {
158     if (fRng)
159     delete fRng;
160 fabstoec 1.1 }
161    
162     //--------------------------------------------------------------------------------------------------
163     void PhotonPairSelector::Process()
164     {
165 paus 1.27 // ------------------------------------------------------------
166     // Process entries of the tree.
167 fabstoec 1.1 LoadEventObject(fPhotonBranchName, fPhotons);
168 paus 1.27
169 bendavid 1.8 // -----------------------------------------------------------
170 paus 1.27 // Output Photon Collection. It will ALWAYS contain either 0 or 2 photons
171     PhotonOArr *GoodPhotons = new PhotonOArr;
172 bendavid 1.8 GoodPhotons->SetName(fGoodPhotonsName);
173     GoodPhotons->SetOwner(kTRUE);
174 fabstoec 1.29
175    
176     VertexOArr *ChosenVtx = new VertexOArr;
177     ChosenVtx->SetName(fChosenVtxName);
178     ChosenVtx->SetOwner(kTRUE);
179 fabstoec 1.30
180 bendavid 1.8 // add to event for other modules to use
181 paus 1.27 AddObjThisEvt(GoodPhotons);
182 fabstoec 1.29 AddObjThisEvt(ChosenVtx);
183 paus 1.27
184 fabstoec 1.30 //AddObjThisEvt(ChosenVtx);
185    
186 paus 1.27 if (fPhotons->GetEntries()<2)
187     return;
188    
189 fabstoec 1.1 LoadEventObject(fElectronName, fElectrons);
190 paus 1.27 LoadEventObject(fGoodElectronName, fGoodElectrons);
191 fabstoec 1.1 LoadEventObject(fConversionName, fConversions);
192 bendavid 1.36 if (fUseSingleLegConversions) LoadEventObject(fPFConversionName, fPFConversions);
193 fabstoec 1.1 LoadEventObject(fTrackBranchName, fTracks);
194     LoadEventObject(fPileUpDenName, fPileUpDen);
195 paus 1.27 LoadEventObject(fPVName, fPV);
196 fabstoec 1.1 LoadEventObject(fBeamspotName, fBeamspot);
197     LoadEventObject(fPFCandName, fPFCands);
198 bendavid 1.31 LoadEventObject(fJetsName, fJets);
199     LoadEventObject(fPFMetName, fPFMet);
200 fabstoec 1.1
201 bendavid 1.31 if (!fIsData) {
202     LoadBranch(fMCParticleName);
203     }
204    
205    
206 paus 1.27 // ------------------------------------------------------------
207 fabstoec 1.1 // load event based information
208 paus 1.27 Float_t rho = -99.;
209     if (fPileUpDen->GetEntries() > 0)
210     rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
211 bendavid 1.33 Float_t rho2012 = -99;
212     if (fPileUpDen->At(0)->RhoKt6PFJets()>0.) rho2012 = fPileUpDen->At(0)->RhoKt6PFJets();
213     else rho2012 = fPileUpDen->At(0)->Rho();
214 fabstoec 1.1 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
215    
216 paus 1.27 // ------------------------------------------------------------
217 fabstoec 1.1 // Get Event header for Run info etc.
218 paus 1.27 const EventHeader* evtHead = this->GetEventHeader();
219     unsigned int evtNum = evtHead->EvtNum();
220     UInt_t runNumber = evtHead->RunNum();
221     Float_t _runNum = (Float_t) runNumber;
222     Float_t _lumiSec = (Float_t) evtHead->LumiSec();
223 ceballos 1.13
224 fabstoec 1.1 // ------------------------------------------------------------
225     // here we'll store the preselected Photons (and which CiCCategory they are...)
226     PhotonOArr* preselPh = new PhotonOArr;
227     std::vector<PhotonTools::CiCBaseLineCats> preselCat;
228     preselCat.resize(0);
229 paus 1.27
230     // 1. we do the pre-selection; but keep the non-passing photons in a second container...
231     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
232 fabstoec 1.1 const Photon *ph = fPhotons->At(i);
233    
234 paus 1.27 if (ph->SCluster()->AbsEta()>= fPhotonEtaMax ||
235     (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566))
236     continue;
237     if (ph->Et() < fPhotonPtMin)
238     continue;
239     if (ph->HadOverEm() > 0.15)
240     continue;
241    
242     if (ph->IsEB()) {
243     if (ph->CoviEtaiEta() > 0.015)
244     continue;
245     }
246     else {
247     if (ph->CoviEtaiEta() > 0.035)
248     continue;
249     }
250     // photon passes the preselection
251     ph->Mark(); // mark for later skimming
252 fabstoec 1.1 preselPh->Add(ph);
253 fabstoec 1.2 }
254    
255 paus 1.27 if (preselPh->GetEntries()<2) {
256     this->SkipEvent();
257     delete preselPh;
258     return;
259     }
260 bendavid 1.8
261 bendavid 1.36 //fill conversion collection for vertex selection, adding single leg conversions if needed
262     //note that momentum of single leg conversions needs to be recomputed from the track
263     //as it is not filled properly
264     DecayParticleOArr vtxconversions;
265     if (fUseSingleLegConversions) {
266     vtxconversions.SetOwner(kTRUE);
267     for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
268     DecayParticle *conv = new DecayParticle(*fConversions->At(iconv));
269     vtxconversions.AddOwned(conv);
270     }
271    
272     for (UInt_t iconv=0; iconv<fPFConversions->GetEntries(); ++iconv) {
273     const DecayParticle *c = fPFConversions->At(iconv);
274     if (c->NDaughters()!=1) continue;
275    
276 bendavid 1.37 DecayParticle *conv = new DecayParticle(*c);
277 bendavid 1.36 const Track *trk = static_cast<const StableParticle*>(conv->Daughter(0))->Trk();
278     conv->SetMom(trk->Px(), trk->Py(), trk->Pz(), trk->P());
279     vtxconversions.AddOwned(conv);
280     }
281     }
282     else {
283     for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
284     const DecayParticle *c = fConversions->At(iconv);
285     vtxconversions.Add(c);
286     }
287     }
288    
289 bendavid 1.31
290     float higgspt = -99.;
291     float higgsz = -99.;
292     float higgsmass = -99.;
293    
294     if (!fIsData) FindHiggsPtAndZ(higgspt,higgsz,higgsmass);
295    
296 paus 1.27 // second loop: sort & assign the right categories..
297 fabstoec 1.2 preselPh->Sort();
298 paus 1.27 for (unsigned int iPh = 0; iPh <preselPh->GetEntries(); ++iPh) {
299 fabstoec 1.2 const Photon* ph = preselPh->At(iPh);
300 fabstoec 1.1 preselCat.push_back(PhotonTools::CiCBaseLineCat(ph));
301     }
302 paus 1.27
303 fabstoec 1.1 // ------------------------------------------------------------
304     // compute how many pairs there are ...
305     unsigned int numPairs = 0;
306 paus 1.27 if (preselPh->GetEntries() > 0)
307     numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
308    
309 fabstoec 1.1 // ... and create all possible pairs of pre-selected photons
310 paus 1.27 std::vector<unsigned int> idx1st;
311     std::vector<unsigned int> idx2nd;
312 fabstoec 1.1 std::vector<PhotonTools::CiCBaseLineCats> cat1st;
313     std::vector<PhotonTools::CiCBaseLineCats> cat2nd;
314 paus 1.27
315     // ... this will be used to store whether a given pair passes the cuts
316 fabstoec 1.1 std::vector<bool> pairPasses;
317 paus 1.27
318     if (numPairs > 0) {
319     for (unsigned int i1st = 0; i1st <preselPh->GetEntries() - 1; ++i1st) {
320     for (unsigned int i2nd = i1st + 1; i2nd <preselPh->GetEntries(); ++i2nd) {
321     idx1st.push_back(i1st);
322     idx2nd.push_back(i2nd);
323     pairPasses.push_back(true);
324 fabstoec 1.1 }
325     }
326     }
327    
328 paus 1.27 // ------------------------------------------------------------
329 fabstoec 1.1 // array to store the index of 'chosen Vtx' for each pair
330 bendavid 1.20 // const Vertex** theVtx = new const Vertex*[numPairs]; // holds the 'chosen' Vtx for each Pair
331 paus 1.27 // Photon** fixPh1st = new Photon*[numPairs]; // holds the 1st Photon for each Pair
332 bendavid 1.20 // Photon** fixPh2nd = new Photon*[numPairs]; // holds the 2nd photon for each Pair
333 paus 1.27
334 bendavid 1.20 std::vector<const Vertex*> theVtx; // holds the 'chosen' Vtx for each Pair
335 paus 1.27 std::vector<Photon*> fixPh1st; // holds the 1st Photon for each Pair
336 bendavid 1.20 std::vector<Photon*> fixPh2nd; // holds the 2nd photon for each Pair
337    
338     theVtx.reserve(numPairs);
339     fixPh1st.reserve(numPairs);
340     fixPh2nd.reserve(numPairs);
341 fabstoec 1.1
342     // store pair-indices for pairs passing the selection
343     std::vector<unsigned int> passPairs;
344     passPairs.resize(0);
345    
346 paus 1.27 // ------------------------------------------------------------
347     // Loop over all Pairs and do the 'incredible machine' running....
348     for (unsigned int iPair = 0; iPair < numPairs; ++iPair) {
349 fabstoec 1.1 // first we need a hard copy of the incoming photons
350 bendavid 1.20 fixPh1st.push_back(new Photon(*preselPh->At(idx1st[iPair])));
351     fixPh2nd.push_back(new Photon(*preselPh->At(idx2nd[iPair])));
352 fabstoec 1.1 // we also store the category, so we don't have to ask all the time...
353     cat1st.push_back(preselCat[idx1st[iPair]]);
354 paus 1.27 cat2nd.push_back(preselCat[idx2nd[iPair]]);
355    
356 bendavid 1.19 //scale regression sigmaE in MC if activated
357     if (fDoMCErrScaling && !fIsData) {
358 paus 1.27 if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5)
359     PhotonTools::ScalePhotonError(fixPh1st[iPair],fMCErrScaleEB);
360     else
361     PhotonTools::ScalePhotonError(fixPh1st[iPair],fMCErrScaleEE);
362    
363     if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5)
364     PhotonTools::ScalePhotonError(fixPh2nd[iPair],fMCErrScaleEB);
365     else
366     PhotonTools::ScalePhotonError(fixPh2nd[iPair],fMCErrScaleEE);
367     }
368 bendavid 1.19
369     //scale R9 in Monte Carlo if activated
370     if (fDoMCR9Scaling && !fIsData) {
371 paus 1.27 if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5)
372     PhotonTools::ScalePhotonR9(fixPh1st[iPair],fMCR9ScaleEB);
373     else
374     PhotonTools::ScalePhotonR9(fixPh1st[iPair],fMCR9ScaleEE);
375    
376     if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5)
377     PhotonTools::ScalePhotonR9(fixPh2nd[iPair],fMCR9ScaleEB);
378     else
379     PhotonTools::ScalePhotonR9(fixPh2nd[iPair],fMCR9ScaleEE);
380 bendavid 1.19 }
381 paus 1.27
382 bendavid 1.23 if (fDoMCSigIEtaIEtaScaling && !fIsData) {
383 mingyang 1.35 if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5) fixPh1st[iPair]->SetCoviEtaiEta(fMCSigIEtaIEtaScaleEB*fixPh1st[iPair]->CoviEtaiEta()+fMCSigIEtaIEtaShiftEB);
384     else fixPh1st[iPair]->SetCoviEtaiEta(fMCSigIEtaIEtaScaleEE*fixPh1st[iPair]->CoviEtaiEta()+fMCSigIEtaIEtaShiftEE);
385 paus 1.27
386 mingyang 1.35 if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5) fixPh2nd[iPair]->SetCoviEtaiEta(fMCSigIEtaIEtaScaleEB*fixPh2nd[iPair]->CoviEtaiEta()+fMCSigIEtaIEtaShiftEB);
387     else fixPh2nd[iPair]->SetCoviEtaiEta(fMCSigIEtaIEtaScaleEE*fixPh2nd[iPair]->CoviEtaiEta()+fMCSigIEtaIEtaShiftEE);
388 bendavid 1.23 }
389 bendavid 1.19
390    
391 bendavid 1.25 if (fDoMCWidthScaling && !fIsData) {
392 mingyang 1.35 fixPh1st[iPair]->SetEtaWidth(fMCWidthScale*fixPh1st[iPair]->EtaWidth()+fMCWidthShift);
393     fixPh1st[iPair]->SetPhiWidth(fMCWidthScale*fixPh1st[iPair]->PhiWidth()+fMCWidthShift);
394 paus 1.27
395 mingyang 1.35 fixPh2nd[iPair]->SetEtaWidth(fMCWidthScale*fixPh2nd[iPair]->EtaWidth()+fMCWidthShift);
396     fixPh2nd[iPair]->SetPhiWidth(fMCWidthScale*fixPh2nd[iPair]->PhiWidth()+fMCWidthShift);
397 bendavid 1.25 }
398    
399 bendavid 1.19 PhotonTools::eScaleCats escalecat1 = PhotonTools::EScaleCat(fixPh1st[iPair]);
400     PhotonTools::eScaleCats escalecat2 = PhotonTools::EScaleCat(fixPh2nd[iPair]);
401 paus 1.27
402 fabstoec 1.1 // now we dicide if we either scale (Data) or Smear (MC) the Photons
403     if (fIsData) {
404 paus 1.27 if (fDoDataEneCorr) {
405     // starting with scale = 1.
406     double scaleFac1 = 1.;
407     double scaleFac2 = 1.;
408    
409 bendavid 1.9 //eta-dependent corrections
410 paus 1.27
411     // checking the run Rangees ...
412     Int_t runRange = FindRunRangeIdx(runNumber);
413     if(runRange > -1) {
414 bendavid 1.20 scaleFac1 *= GetDataEnCorr(runRange, escalecat1);
415     scaleFac2 *= GetDataEnCorr(runRange, escalecat2);
416 paus 1.27 }
417     PhotonTools::ScalePhoton(fixPh1st[iPair], scaleFac1);
418     PhotonTools::ScalePhoton(fixPh2nd[iPair], scaleFac2);
419 fabstoec 1.1 }
420 paus 1.27 }
421    
422     if (fDoMCSmear) {
423 bendavid 1.19 double width1 = GetMCSmearFac(escalecat1);
424     double width2 = GetMCSmearFac(escalecat2);
425 bendavid 1.9 if (!fIsData) {
426     // get the seed to do deterministic smearing...
427     UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
428 paus 1.27 UInt_t seed1 = seedBase + (UInt_t) fixPh1st[iPair]->E() +
429     (UInt_t) (TMath::Abs(10.*fixPh1st[iPair]->SCluster()->Eta()));
430     UInt_t seed2 = seedBase + (UInt_t) fixPh2nd[iPair]->E() +
431     (UInt_t) (TMath::Abs(10.*fixPh2nd[iPair]->SCluster()->Eta()));
432 bendavid 1.9 // get the smearing for MC photons..
433 paus 1.27 PhotonTools::SmearPhoton(fixPh1st[iPair], fRng, width1, seed1);
434     PhotonTools::SmearPhoton(fixPh2nd[iPair], fRng, width2, seed2);
435 fabstoec 1.1 }
436 bendavid 1.9 PhotonTools::SmearPhotonError(fixPh1st[iPair], width1);
437     PhotonTools::SmearPhotonError(fixPh2nd[iPair], width2);
438 fabstoec 1.1 }
439    
440 bendavid 1.19 //probability that selected vertex is the correct one
441     Double_t vtxProb = 1.0;
442    
443 fabstoec 1.1 // store the vertex for this pair
444 paus 1.27 switch (fVtxSelType) {
445 fabstoec 1.28
446 fabstoec 1.1 case kStdVtxSelection:
447     theVtx[iPair] = fPV->At(0);
448     break;
449 fabstoec 1.28
450 fabstoec 1.1 case kCiCVtxSelection:
451 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
452 bendavid 1.36 &vtxconversions,kFALSE,vtxProb);
453 fabstoec 1.1 break;
454 fabstoec 1.28
455 bendavid 1.19 case kCiCMVAVtxSelection:
456 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
457 bendavid 1.36 &vtxconversions,kTRUE,vtxProb);
458 paus 1.27 break;
459 fabstoec 1.28
460 fabstoec 1.1 case kMITVtxSelection:
461     // need PFCandidate Collection
462 paus 1.27 theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp,
463     mithep::FourVector((fixPh1st[iPair]->Mom()+
464     fixPh2nd[iPair]->Mom())));
465 fabstoec 1.1 break;
466 bendavid 1.31 case kMetSigVtxSelection: {
467     // need PFCandidate Collection
468    
469     PFJetOArr pfjets;
470     for (UInt_t ijet=0; ijet<fJets->GetEntries(); ++ijet) {
471     const PFJet *pfjet = dynamic_cast<const PFJet*>(fJets->At(ijet));
472     if (pfjet && MathUtils::DeltaR(*pfjet,*fixPh1st[iPair])>0.3 && MathUtils::DeltaR(*pfjet,*fixPh2nd[iPair])>0.3) pfjets.Add(pfjet);
473     }
474    
475     PFCandidateOArr pfcands;
476     for (UInt_t icand=0; icand<fPFCands->GetEntries(); ++icand) {
477     const PFCandidate *pfcand = fPFCands->At(icand);
478     if (MathUtils::DeltaR(*pfcand,*fixPh1st[iPair])>0.05 && MathUtils::DeltaR(*pfcand,*fixPh2nd[iPair])>0.05) pfcands.Add(pfcand);
479     }
480    
481     double minsig = 1e6;
482     for (UInt_t ivtx=0; ivtx<fPV->GetEntries(); ++ivtx) {
483 bendavid 1.38 const Vertex *v = fPV->At(ivtx);
484    
485    
486    
487    
488 bendavid 1.31 // Met mmet = fMVAMet.GetMet( false,
489     // fixPh1st[iPair]->Pt(),fixPh1st[iPair]->Phi(),fixPh1st[iPair]->Eta(),
490     // fixPh2nd[iPair]->Pt(),fixPh2nd[iPair]->Phi(),fixPh2nd[iPair]->Eta(),
491     // fPFMet->At(0),
492     // &pfcands,fPV->At(ivtx),fPV, fPileUpDen->At(0)->Rho(),
493     // &pfjets,
494     // int(fPV->GetEntries()),
495     // kTRUE);
496    
497     Met mmet = fMVAMet.GetMet( false,
498     0.,0.,0.,
499     fPFMet->At(0),
500     &pfcands,fPV->At(ivtx),fPV, fPileUpDen->At(0)->Rho(),
501     &pfjets,
502     int(fPV->GetEntries()),
503     kFALSE);
504    
505     ThreeVector fullmet(mmet.Px() - fixPh1st[iPair]->Px() - fixPh2nd[iPair]->Px(),
506     mmet.Py() - fixPh1st[iPair]->Py() - fixPh2nd[iPair]->Py(),
507     0.);
508    
509     // ThreeVector fullmet(mmet.Px(),
510     // mmet.Py(),
511     // 0.);
512    
513     TMatrixD *metcov = fMVAMet.GetMetCovariance();
514    
515     //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();
516    
517     //Double_t metsig = fullmet.Rho()/metsigma;
518    
519     ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
520     mcov(0,0) = (*metcov)(0,0);
521     mcov(0,1) = (*metcov)(0,1);
522     mcov(1,0) = (*metcov)(1,0);
523     mcov(1,1) = (*metcov)(1,1);
524    
525     ROOT::Math::SVector<double,2> vmet;
526     vmet(0) = fullmet.X();
527     vmet(1) = fullmet.Y();
528    
529     mcov.Invert();
530    
531     Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
532    
533 bendavid 1.38
534     Double_t sumptsq = 0.;
535     for (UInt_t ipfc = 0; ipfc<fPFCands->GetEntries(); ++ipfc) {
536     const PFCandidate *pfc = fPFCands->At(ipfc);
537     if (pfc->PFType()!=PFCandidate::eHadron || !pfc->HasTrackerTrk()) continue;
538     if (TMath::Abs( pfc->TrackerTrk()->DzCorrected(*v) ) > 0.2) continue;
539     if (TMath::Abs( pfc->TrackerTrk()->D0Corrected(*v) ) > 0.1) continue;
540     sumptsq += pfc->Pt()*pfc->Pt();
541     }
542    
543     if (sumptsq<10.0) metsig = -99;
544    
545 bendavid 1.31 if (metsig<minsig && metsig>0.) {
546     minsig = metsig;
547     theVtx[iPair] = fPV->At(ivtx);
548     }
549    
550 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);
551 bendavid 1.31
552     }
553    
554     double testpfmetx = 0.;
555     double testpfmety = 0.;
556     for (UInt_t icand=0; icand<pfcands.GetEntries(); ++icand) {
557     const PFCandidate *pfcand = pfcands.At(icand);
558     testpfmetx -= pfcand->Px();
559     testpfmety -= pfcand->Py();
560     //pfcands.Add(pfcand);
561     }
562    
563     testpfmetx -= fixPh1st[iPair]->Px();
564     testpfmetx -= fixPh2nd[iPair]->Px();
565    
566     testpfmety -= fixPh1st[iPair]->Py();
567     testpfmety -= fixPh2nd[iPair]->Py();
568    
569     double testpfmet = sqrt(testpfmetx*testpfmetx + testpfmety*testpfmety);
570    
571     printf("bestvtx: metsig = %5f, dzgen = %5f, diphopt = %5f, pfmet = %5f, testpfmet = %5f\n", minsig, theVtx[iPair]->Z()-higgsz, higgspt,fPFMet->At(0)->Pt(),testpfmet);
572    
573     }
574     break;
575 fabstoec 1.1 default:
576     theVtx[iPair] = fPV->At(0);
577 paus 1.27 }
578 fabstoec 1.2
579 bendavid 1.8 //set PV ref in photons
580     fixPh1st[iPair]->SetPV(theVtx[iPair]);
581     fixPh2nd[iPair]->SetPV(theVtx[iPair]);
582 bendavid 1.19 fixPh1st[iPair]->SetVtxProb(vtxProb);
583     fixPh2nd[iPair]->SetVtxProb(vtxProb);
584 paus 1.27
585 fabstoec 1.1 // fix the kinematics for both events
586     FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(theVtx[iPair]->Position());
587     FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(theVtx[iPair]->Position());
588     fixPh1st[iPair]->SetMom(newMom1st.X(), newMom1st.Y(), newMom1st.Z(), newMom1st.E());
589     fixPh2nd[iPair]->SetMom(newMom2nd.X(), newMom2nd.Y(), newMom2nd.Z(), newMom2nd.E());
590    
591 bendavid 1.19 double pairmass = (fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M();
592 paus 1.27
593 bendavid 1.20 double leadptcut = fLeadingPtMin;
594     double trailptcut = fTrailingPtMin;
595 paus 1.27
596 bendavid 1.20 if (fixPh2nd[iPair]->Pt() > fixPh1st[iPair]->Pt()) {
597     leadptcut = fTrailingPtMin;
598     trailptcut = fLeadingPtMin;
599     }
600 paus 1.27
601    
602 bendavid 1.20 if (fRelativePtCuts) {
603     leadptcut = leadptcut*pairmass;
604     trailptcut = trailptcut*pairmass;
605     }
606 paus 1.27
607 mingyang 1.39
608 bendavid 1.21 //compute id bdt values
609 mingyang 1.34 Double_t bdt1;
610     Double_t bdt2;
611 mingyang 1.39 Bool_t applyNewScale=kFALSE;
612     if ((!fIsData) && fDoMCNewScaling ){
613     applyNewScale=kTRUE;
614     }
615 mingyang 1.15
616 mingyang 1.34 switch (fIdType) {
617     case k2011IdMVA:
618     bdt1 = fTool.GetMVAbdtValue_2011(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho,fElectrons,fApplyEleVeto);
619     bdt2 = fTool.GetMVAbdtValue_2011(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho,fElectrons,fApplyEleVeto);
620     fixPh1st[iPair]->SetIdMva(bdt1);
621     fixPh2nd[iPair]->SetIdMva(bdt2);
622     break;
623    
624     case k2012IdMVA_globe:
625 mingyang 1.39 bdt1 = fTool.GetMVAbdtValue_2012_globe(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho2012,fPFCands,applyNewScale,fElectrons,fApplyEleVeto);
626     bdt2 = fTool.GetMVAbdtValue_2012_globe(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho2012,fPFCands,applyNewScale,fElectrons,fApplyEleVeto);
627 mingyang 1.34 fixPh1st[iPair]->SetIdMva(bdt1);
628     fixPh2nd[iPair]->SetIdMva(bdt2);
629     break;
630     }
631 paus 1.27
632 bendavid 1.20 //printf("applying id\n");
633 paus 1.27
634 fabstoec 1.1 // check if both photons pass the CiC selection
635     // FIX-ME: Add other possibilities....
636     bool pass1 = false;
637     bool pass2 = false;
638 paus 1.27
639     switch (fPhSelType) {
640 fabstoec 1.1 case kNoPhSelection:
641 paus 1.27 pass1 = (fixPh1st[iPair]->Pt() > leadptcut );
642     pass2 = (fixPh2nd[iPair]->Pt() > trailptcut);
643 fabstoec 1.1 break;
644     case kCiCPhSelection:
645 paus 1.27 pass1 = PhotonTools::PassCiCSelection(fixPh1st[iPair], theVtx[iPair], fTracks,
646     fElectrons, fPV, rho, leadptcut, fApplyEleVeto);
647     if (pass1)
648     pass2 = PhotonTools::PassCiCSelection(fixPh2nd[iPair], theVtx[iPair], fTracks,
649     fElectrons, fPV, rho, trailptcut, fApplyEleVeto);
650 fabstoec 1.1 break;
651 bendavid 1.33 case kCiCPFPhSelection:
652    
653     // loose preselection for mva
654     pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
655     PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
656     fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
657     fInvertElectronVeto);
658     if (pass1)
659     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
660     PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
661     fTracks,theVtx[iPair],rho2012,fPFCands,fApplyEleVeto,
662     fInvertElectronVeto);
663    
664     pass1 = pass1 && PhotonTools::PassCiCPFIsoSelection(fixPh1st[iPair], theVtx[iPair], fPFCands,
665     fPV, rho2012, leadptcut);
666     if (pass1)
667     pass2 = pass2 && PhotonTools::PassCiCPFIsoSelection(fixPh2nd[iPair], theVtx[iPair], fPFCands,
668     fPV, rho2012, trailptcut);
669     break;
670 mingyang 1.15 case kMVAPhSelection://MVA
671 paus 1.27 pass1 = fixPh1st[iPair]->Pt()>leadptcut &&
672     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
673     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
674     fTool.PassMVASelection(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho,
675     fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
676     if (pass1)
677     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
678     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
679     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
680     fTool.PassMVASelection(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho,
681     fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
682    
683 mingyang 1.15 break;
684 fabstoec 1.1 case kMITPhSelection:
685 bendavid 1.20 // loose preselection for mva
686 paus 1.27 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
687     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
688 bendavid 1.33 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
689 paus 1.27 fInvertElectronVeto);
690     if (pass1)
691     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
692     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
693 bendavid 1.33 fTracks,theVtx[iPair],rho2012,fApplyEleVeto,
694 paus 1.27 fInvertElectronVeto);
695    
696 fabstoec 1.1 break;
697 bendavid 1.33 case kMITPFPhSelection:
698     // loose preselection for mva
699     pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
700     PhotonTools::PassSinglePhotonPreselPFISO(fixPh1st[iPair],fElectrons,fConversions,bsp,
701     fTracks,theVtx[iPair],rho,fPFCands,fApplyEleVeto,
702     fInvertElectronVeto);
703     if (pass1)
704     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
705     PhotonTools::PassSinglePhotonPreselPFISO(fixPh2nd[iPair],fElectrons,fConversions,bsp,
706     fTracks,theVtx[iPair],rho,fPFCands,fApplyEleVeto,
707     fInvertElectronVeto);
708    
709     break;
710 fabstoec 1.1 default:
711     pass1 = true;
712     pass2 = true;
713     }
714 paus 1.27
715 bendavid 1.8 //match to good electrons if requested
716 bendavid 1.33 // if (fInvertElectronVeto) {
717     // pass1 &= !PhotonTools::PassElectronVeto(fixPh1st[iPair],fGoodElectrons);
718     // pass2 &= !PhotonTools::PassElectronVeto(fixPh2nd[iPair],fGoodElectrons);
719     // }
720 paus 1.27 // finally, if both Photons pass the selections, add the pair to the passing pairs
721     if (pass1 && pass2)
722     passPairs.push_back(iPair);
723 fabstoec 1.1 }
724 paus 1.27
725    
726 fabstoec 1.1 // ---------------------------------------------------------------
727 paus 1.18 // ... we're almost done, stay focused...
728 fabstoec 1.1 // loop over all passing pairs and find the one with the highest sum Et
729 paus 1.27 const Vertex* vtx = NULL;
730 fabstoec 1.1 Photon* phHard = NULL;
731     Photon* phSoft = NULL;
732    
733 fabstoec 1.29 // not used at all....
734     //PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
735     //PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
736 paus 1.27
737 fabstoec 1.1 double maxSumEt = 0.;
738 paus 1.27 for (unsigned int iPair=0; iPair<passPairs.size(); ++iPair) {
739     double sumEt = fixPh1st[passPairs[iPair]]->Et()
740     + fixPh2nd[passPairs[iPair]]->Et();
741     if (sumEt > maxSumEt) {
742 fabstoec 1.1 maxSumEt = sumEt;
743 paus 1.27 phHard = fixPh1st[passPairs[iPair]];
744     phSoft = fixPh2nd[passPairs[iPair]];
745 fabstoec 1.29 //catPh1 = cat1st [passPairs[iPair]];
746     //catPh2 = cat2nd [passPairs[iPair]];
747 paus 1.27 vtx = theVtx[iPair];
748 fabstoec 1.1 }
749     }
750 bendavid 1.20
751 paus 1.27 for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
752 bendavid 1.20 if (fixPh1st[iPair]!=phHard) delete fixPh1st[iPair];
753     if (fixPh2nd[iPair]!=phSoft) delete fixPh2nd[iPair];
754     }
755 paus 1.27
756 fabstoec 1.1 // ---------------------------------------------------------------
757     // we have the Photons (*PARTY*)... compute some useful qunatities
758    
759     if(phHard && phSoft) {
760     GoodPhotons->AddOwned(phHard);
761     GoodPhotons->AddOwned(phSoft);
762     }
763    
764 fabstoec 1.29 // we also store the chosen Vtx, so later modules can use it
765 fabstoec 1.30
766     Vertex* chosenVtx = NULL;
767     if ( vtx ) {
768     chosenVtx = new Vertex( *vtx );
769     ChosenVtx->AddOwned( chosenVtx );
770     }
771 fabstoec 1.29
772 fabstoec 1.1 // sort according to pt
773     GoodPhotons->Sort();
774 fabstoec 1.30
775 fabstoec 1.1 // delete auxiliary photon collection...
776     delete preselPh;
777 bendavid 1.20 //delete[] theVtx;
778 paus 1.27
779 fabstoec 1.1 return;
780     }
781    
782 paus 1.27 //---------------------------------------------------------------------------------------------------
783 fabstoec 1.1 void PhotonPairSelector::SlaveBegin()
784     {
785 paus 1.27 // Run startup code on the computer (slave) doing the actual analysis. Here, we just request the
786     // photon collection branch.
787 fabstoec 1.1
788 paus 1.27 // load all branches
789     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
790     ReqEventObject(fTrackBranchName, fTracks, true);
791     ReqEventObject(fElectronName, fElectrons, true);
792     ReqEventObject(fGoodElectronName, fGoodElectrons,fGoodElectronsFromBranch);
793     ReqEventObject(fPileUpDenName, fPileUpDen, true);
794     ReqEventObject(fPVName, fPV, fPVFromBranch);
795     ReqEventObject(fConversionName, fConversions, true);
796 bendavid 1.36 if (fUseSingleLegConversions) ReqEventObject(fPFConversionName, fPFConversions, true);
797 paus 1.27 ReqEventObject(fBeamspotName, fBeamspot, true);
798     ReqEventObject(fPFCandName, fPFCands, true);
799 bendavid 1.31 ReqEventObject(fJetsName, fJets, false);
800     ReqEventObject(fPFMetName, fPFMet, true);
801 fabstoec 1.1 if (!fIsData) {
802 bendavid 1.31 //ReqBranch(fPileUpName, fPileUp);
803 fabstoec 1.1 ReqBranch(fMCParticleName, fMCParticles);
804     }
805 paus 1.27 // determine photon selection type
806     if (fPhotonSelType.CompareTo("CiCSelection") == 0)
807 fabstoec 1.1 fPhSelType = kCiCPhSelection;
808 bendavid 1.33 else if (fPhotonSelType.CompareTo("CiCPFSelection") == 0)
809     fPhSelType = kCiCPFPhSelection;
810 mingyang 1.15 else if (fPhotonSelType.CompareTo("MVASelection") == 0) //MVA
811     fPhSelType = kMVAPhSelection;
812 paus 1.27 else if (fPhotonSelType.CompareTo("MITSelection") == 0)
813 fabstoec 1.1 fPhSelType = kMITPhSelection;
814 bendavid 1.33 else if (fPhotonSelType.CompareTo("MITPFSelection") == 0)
815     fPhSelType = kMITPFPhSelection;
816     else if (fPhotonSelType.CompareTo("NoSelection") == 0)
817 fabstoec 1.1 fPhSelType = kNoPhSelection;
818 bendavid 1.33 else {
819     std::cerr<<" Photon Seclection "<<fPhotonSelType<<" not implemented."<<std::endl;
820     return;
821     }
822 fabstoec 1.1
823 paus 1.27 if (fVertexSelType.CompareTo("CiCSelection") == 0)
824 fabstoec 1.1 fVtxSelType = kCiCVtxSelection;
825 paus 1.27 else if (fVertexSelType.CompareTo("MITSelection") == 0)
826 fabstoec 1.1 fVtxSelType = kMITVtxSelection;
827 paus 1.27 else if (fVertexSelType.CompareTo("CiCMVASelection") == 0)
828     fVtxSelType = kCiCMVAVtxSelection;
829 bendavid 1.31 else if (fVertexSelType.CompareTo("MetSigSelection") == 0)
830     fVtxSelType = kMetSigVtxSelection;
831 fabstoec 1.28 else if (fVertexSelType.CompareTo("ZeroVtxSelection") == 0)
832 paus 1.27 fVtxSelType = kStdVtxSelection;
833 fabstoec 1.28 else {
834     std::cerr<<" Vertex Seclection "<<fVertexSelType<<" not implemented."<<std::endl;
835     return;
836     }
837 mingyang 1.34
838     if (fIdMVAType.CompareTo("2011IdMVA") == 0)
839     fIdType = k2011IdMVA;
840     else if (fIdMVAType.CompareTo("2012IdMVA_globe") == 0)
841     fIdType = k2012IdMVA_globe;
842     else {
843     std::cerr<<" Id MVA "<<fIdMVAType<<" not implemented."<<std::endl;
844     return;
845     }
846    
847 paus 1.27 if (fIsData)
848 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
849 paus 1.27 else
850 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
851 mingyang 1.15
852 paus 1.27 printf("initialize photon pair selector\n");
853 mingyang 1.15
854 mingyang 1.34 switch (fIdType) {
855     case k2011IdMVA:
856     fTool.InitializeMVA(fVariableType_2011,fEndcapWeights_2011,fBarrelWeights_2011);
857     break;
858    
859     case k2012IdMVA_globe:
860     fTool.InitializeMVA(fVariableType_2012_globe,fEndcapWeights_2012_globe,fBarrelWeights_2012_globe);
861     break;
862     }
863    
864 paus 1.27 fVtxTools.InitP();
865 bendavid 1.31
866     if (fVtxSelType==kMetSigVtxSelection) {
867     //fMVAMet.Initialize();
868     // fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
869     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
870     // TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
871     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_42.root"))),
872     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_42.root"))),
873     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1_42.root"))),
874     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2_42.root")))
875     // );
876    
877     fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
878     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
879     TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
880     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_52.root"))),
881     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_52.root"))),
882     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1cov_52.root"))),
883     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2cov_52.root")))
884     );
885    
886     }
887 mingyang 1.15
888 fabstoec 1.1 }
889    
890     // ----------------------------------------------------------------------------------------
891     // some helpfer functions....
892 paus 1.27 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
893     {
894     pt = -999.;
895 fabstoec 1.1 decayZ = -999.;
896 paus 1.27 mass = -999.;
897 bendavid 1.31
898 fabstoec 1.1 // loop over all GEN particles and look for status 1 photons
899     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
900     const MCParticle* p = fMCParticles->At(i);
901 bendavid 1.8 if( p->Is(MCParticle::kH) || (!fApplyEleVeto && p->AbsPdgId()==23) ) {
902     pt=p->Pt();
903     decayZ = p->DecayVertex().Z();
904     mass = p->Mass();
905     break;
906     }
907 fabstoec 1.1 }
908 paus 1.27
909 fabstoec 1.1 return;
910 paus 1.27 }
911 fabstoec 1.1
912 paus 1.27 //---------------------------------------------------------------------------------------------------
913     Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run)
914     {
915     // this routine looks for the idx of the run-range
916 fabstoec 1.1 Int_t runRange=-1;
917 paus 1.27 for (UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
918     if (run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
919 fabstoec 1.1 runRange = (Int_t) iRun;
920     return runRange;
921     }
922     }
923     return runRange;
924     }
925    
926 paus 1.27 //---------------------------------------------------------------------------------------------------
927     Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::eScaleCats cat)
928     {
929     switch (cat) {
930 bendavid 1.19 case PhotonTools::kEBhighEtaGold:
931     return fDataEnCorr_EBhighEta_hR9[runRange];
932     case PhotonTools::kEBhighEtaBad:
933     return fDataEnCorr_EBhighEta_lR9[runRange];
934 bendavid 1.20 case PhotonTools::kEBlowEtaGoldCenter:
935     return fDataEnCorr_EBlowEta_hR9central[runRange];
936     case PhotonTools::kEBlowEtaGoldGap:
937 paus 1.27 return fDataEnCorr_EBlowEta_hR9gap[runRange];
938 bendavid 1.19 case PhotonTools::kEBlowEtaBad:
939 paus 1.27 return fDataEnCorr_EBlowEta_lR9[runRange];
940 bendavid 1.19 case PhotonTools::kEEhighEtaGold:
941     return fDataEnCorr_EEhighEta_hR9[runRange];
942     case PhotonTools::kEEhighEtaBad:
943     return fDataEnCorr_EEhighEta_lR9[runRange];
944     case PhotonTools::kEElowEtaGold:
945     return fDataEnCorr_EElowEta_hR9[runRange];
946     case PhotonTools::kEElowEtaBad:
947 paus 1.27 return fDataEnCorr_EElowEta_lR9[runRange];
948 fabstoec 1.1 default:
949     return 1.;
950     }
951     }
952    
953 paus 1.27 //---------------------------------------------------------------------------------------------------
954     Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::eScaleCats cat)
955     {
956     switch (cat) {
957 bendavid 1.19 case PhotonTools::kEBhighEtaGold:
958     return fMCSmear_EBhighEta_hR9;
959     case PhotonTools::kEBhighEtaBad:
960     return fMCSmear_EBhighEta_lR9;
961 bendavid 1.20 case PhotonTools::kEBlowEtaGoldCenter:
962     return fMCSmear_EBlowEta_hR9central;
963     case PhotonTools::kEBlowEtaGoldGap:
964 paus 1.27 return fMCSmear_EBlowEta_hR9gap;
965 bendavid 1.19 case PhotonTools::kEBlowEtaBad:
966 paus 1.27 return fMCSmear_EBlowEta_lR9;
967 bendavid 1.19 case PhotonTools::kEEhighEtaGold:
968     return fMCSmear_EEhighEta_hR9;
969     case PhotonTools::kEEhighEtaBad:
970     return fMCSmear_EEhighEta_lR9;
971     case PhotonTools::kEElowEtaGold:
972     return fMCSmear_EElowEta_hR9;
973     case PhotonTools::kEElowEtaBad:
974 paus 1.27 return fMCSmear_EElowEta_lR9;
975 fabstoec 1.1 default:
976     return 1.;
977     }
978     }
979    
980 paus 1.27 //---------------------------------------------------------------------------------------------------
981     Float_t PhotonPairSelector::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
982     PhotonTools::CiCBaseLineCats cat2)
983     {
984     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
985     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
986 fabstoec 1.1
987     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
988     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
989 paus 1.27
990     if (ph1IsEB && ph2IsEB)
991     return (ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
992    
993     return (ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
994 fabstoec 1.1 }
995 bendavid 1.4