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