ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.30
Committed: Thu May 3 09:02:59 2012 UTC (13 years ago) by fabstoec
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028, Mit_027, Mit_027a
Changes since 1.29: +10 -3 lines
Log Message:
bug-fix for Vtx selection

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.8 #include "MitAna/DataTree/interface/StableData.h"
5     #include "MitAna/DataTree/interface/StableParticle.h"
6 fabstoec 1.1 #include "MitPhysics/Init/interface/ModNames.h"
7     #include "MitPhysics/Utils/interface/IsolationTools.h"
8     #include "MitPhysics/Utils/interface/PhotonTools.h"
9     #include "MitPhysics/Utils/interface/VertexTools.h"
10 bendavid 1.6 #include "MitPhysics/Utils/interface/PhotonFix.h"
11 mingyang 1.15 #include "MitPhysics/Utils/interface/MVATools.h"
12 bendavid 1.5 #include "TDataMember.h"
13 fabstoec 1.1 #include <TNtuple.h>
14     #include <TRandom3.h>
15 bendavid 1.6 #include <TSystem.h>
16 bendavid 1.9 #include <TH1D.h>
17 fabstoec 1.1
18     using namespace mithep;
19    
20     ClassImp(mithep::PhotonPairSelector)
21    
22     //--------------------------------------------------------------------------------------------------
23 paus 1.27 PhotonPairSelector::PhotonPairSelector(const char *name, const char *title) :
24 fabstoec 1.1 // Base Module...
25 paus 1.27 BaseMod (name,title),
26 fabstoec 1.1 // define all the Branches to load
27 paus 1.27 fPhotonBranchName (Names::gkPhotonBrn),
28     fElectronName (Names::gkElectronBrn),
29     fGoodElectronName (Names::gkElectronBrn),
30     fConversionName (Names::gkMvfConversionBrn),
31     fTrackBranchName (Names::gkTrackBrn),
32     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
33     fPVName (Names::gkPVBeamSpotBrn),
34     fBeamspotName (Names::gkBeamSpotBrn),
35     fPFCandName (Names::gkPFCandidatesBrn),
36 fabstoec 1.1 // MC specific stuff...
37 paus 1.27 fMCParticleName (Names::gkMCPartBrn),
38     fPileUpName (Names::gkPileupInfoBrn),
39     fGoodPhotonsName (ModNames::gkGoodPhotonsName),
40 fabstoec 1.29 fChosenVtxName ("HggChosenVtx"),
41 fabstoec 1.1 // ----------------------------------------
42     // Selection Types
43 paus 1.27 fPhotonSelType ("NoSelection"),
44     fVertexSelType ("StdSelection"),
45     fPhSelType (kNoPhSelection),
46     fVtxSelType (kStdVtxSelection),
47 fabstoec 1.1 // ----------------------------------------
48 paus 1.27 fPhotonPtMin (20.0),
49     fPhotonEtaMax (2.5),
50     fLeadingPtMin (100.0/3.0),
51     fTrailingPtMin (100.0/4.0),
52     fIsData (false),
53     fPhotonsFromBranch (true),
54     fPVFromBranch (true),
55     fGoodElectronsFromBranch (kTRUE),
56 fabstoec 1.1 // ----------------------------------------
57     // collections....
58 paus 1.27 fPhotons (0),
59     fElectrons (0),
60     fConversions (0),
61     fTracks (0),
62     fPileUpDen (0),
63     fPV (0),
64     fBeamspot (0),
65     fPFCands (0),
66     fMCParticles (0),
67     fPileUp (0),
68 fabstoec 1.1 // ---------------------------------------
69 paus 1.27 fDataEnCorr_EBlowEta_hR9central(0.),
70     fDataEnCorr_EBlowEta_hR9gap (0.),
71     fDataEnCorr_EBlowEta_lR9 (0.),
72     fDataEnCorr_EBhighEta_hR9 (0.),
73     fDataEnCorr_EBhighEta_lR9 (0.),
74     fDataEnCorr_EElowEta_hR9 (0.),
75     fDataEnCorr_EElowEta_lR9 (0.),
76     fDataEnCorr_EEhighEta_hR9 (0.),
77     fDataEnCorr_EEhighEta_lR9 (0.),
78     fRunStart (0),
79     fRunEnd (0),
80     fMCSmear_EBlowEta_hR9central (0.),
81     fMCSmear_EBlowEta_hR9gap (0.),
82     fMCSmear_EBlowEta_lR9 (0.),
83     fMCSmear_EBhighEta_hR9 (0.),
84     fMCSmear_EBhighEta_lR9 (0.),
85     fMCSmear_EElowEta_hR9 (0.),
86     fMCSmear_EElowEta_lR9 (0.),
87     fMCSmear_EEhighEta_hR9 (0.),
88     fMCSmear_EEhighEta_lR9 (0.),
89 fabstoec 1.1 // ---------------------------------------
90 paus 1.27 fRng (new TRandom3()),
91     fPhFixString ("4_2"),
92     fEtaCorrections (0),
93 fabstoec 1.1 // ---------------------------------------
94 paus 1.27 fDoDataEneCorr (true),
95     fDoMCSmear (true),
96     fDoVtxSelection (true),
97     fApplyEleVeto (true),
98     fInvertElectronVeto (kFALSE),
99 mingyang 1.15 //MVA
100 paus 1.27 fVariableType (10), //please use 4 which is the correct type
101     fEndcapWeights (gSystem->Getenv("CMSSW_BASE")+
102     TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
103     TString("Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
104     TString("weights.xml")),
105     fBarrelWeights (gSystem->Getenv("CMSSW_BASE")+
106     TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
107     TString("Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
108     TString("weights.xml")),
109     fbdtCutBarrel (0.0744), //cuts give same eff (relative to presel) with cic
110     fbdtCutEndcap (0.0959), //cuts give same eff (relative to presel) with cic
111     fDoMCR9Scaling (kFALSE),
112     fMCR9ScaleEB (1.0),
113     fMCR9ScaleEE (1.0),
114     fDoMCSigIEtaIEtaScaling (kFALSE),
115     fDoMCWidthScaling (kFALSE),
116     fDoMCErrScaling (kFALSE),
117     fMCErrScaleEB (1.0),
118     fMCErrScaleEE (1.0),
119     fRelativePtCuts (kFALSE)
120 fabstoec 1.1 {
121     // Constructor.
122     }
123    
124 paus 1.27 PhotonPairSelector::~PhotonPairSelector()
125     {
126     if (fRng)
127     delete fRng;
128 fabstoec 1.1 }
129    
130     //--------------------------------------------------------------------------------------------------
131     void PhotonPairSelector::Process()
132     {
133 paus 1.27 // ------------------------------------------------------------
134     // Process entries of the tree.
135 fabstoec 1.1 LoadEventObject(fPhotonBranchName, fPhotons);
136 paus 1.27
137 bendavid 1.8 // -----------------------------------------------------------
138 paus 1.27 // Output Photon Collection. It will ALWAYS contain either 0 or 2 photons
139     PhotonOArr *GoodPhotons = new PhotonOArr;
140 bendavid 1.8 GoodPhotons->SetName(fGoodPhotonsName);
141     GoodPhotons->SetOwner(kTRUE);
142 fabstoec 1.29
143    
144     VertexOArr *ChosenVtx = new VertexOArr;
145     ChosenVtx->SetName(fChosenVtxName);
146     ChosenVtx->SetOwner(kTRUE);
147 fabstoec 1.30
148 bendavid 1.8 // add to event for other modules to use
149 paus 1.27 AddObjThisEvt(GoodPhotons);
150 fabstoec 1.29 AddObjThisEvt(ChosenVtx);
151 paus 1.27
152 fabstoec 1.30 //AddObjThisEvt(ChosenVtx);
153    
154 paus 1.27 if (fPhotons->GetEntries()<2)
155     return;
156    
157 fabstoec 1.1 LoadEventObject(fElectronName, fElectrons);
158 paus 1.27 LoadEventObject(fGoodElectronName, fGoodElectrons);
159 fabstoec 1.1 LoadEventObject(fConversionName, fConversions);
160     LoadEventObject(fTrackBranchName, fTracks);
161     LoadEventObject(fPileUpDenName, fPileUpDen);
162 paus 1.27 LoadEventObject(fPVName, fPV);
163 fabstoec 1.1 LoadEventObject(fBeamspotName, fBeamspot);
164     LoadEventObject(fPFCandName, fPFCands);
165    
166 paus 1.27 // ------------------------------------------------------------
167 fabstoec 1.1 // load event based information
168 paus 1.27 Float_t rho = -99.;
169     if (fPileUpDen->GetEntries() > 0)
170     rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
171 fabstoec 1.1 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
172    
173 paus 1.27 // ------------------------------------------------------------
174 fabstoec 1.1 // Get Event header for Run info etc.
175 paus 1.27 const EventHeader* evtHead = this->GetEventHeader();
176     unsigned int evtNum = evtHead->EvtNum();
177     UInt_t runNumber = evtHead->RunNum();
178     Float_t _runNum = (Float_t) runNumber;
179     Float_t _lumiSec = (Float_t) evtHead->LumiSec();
180 ceballos 1.13
181 fabstoec 1.1 // ------------------------------------------------------------
182     // here we'll store the preselected Photons (and which CiCCategory they are...)
183     PhotonOArr* preselPh = new PhotonOArr;
184     std::vector<PhotonTools::CiCBaseLineCats> preselCat;
185     preselCat.resize(0);
186 paus 1.27
187     // 1. we do the pre-selection; but keep the non-passing photons in a second container...
188     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
189 fabstoec 1.1 const Photon *ph = fPhotons->At(i);
190    
191 paus 1.27 if (ph->SCluster()->AbsEta()>= fPhotonEtaMax ||
192     (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566))
193     continue;
194     if (ph->Et() < fPhotonPtMin)
195     continue;
196     if (ph->HadOverEm() > 0.15)
197     continue;
198    
199     if (ph->IsEB()) {
200     if (ph->CoviEtaiEta() > 0.015)
201     continue;
202     }
203     else {
204     if (ph->CoviEtaiEta() > 0.035)
205     continue;
206     }
207     // photon passes the preselection
208     ph->Mark(); // mark for later skimming
209 fabstoec 1.1 preselPh->Add(ph);
210 fabstoec 1.2 }
211    
212 paus 1.27 if (preselPh->GetEntries()<2) {
213     this->SkipEvent();
214     delete preselPh;
215     return;
216     }
217 bendavid 1.8
218 paus 1.27 // second loop: sort & assign the right categories..
219 fabstoec 1.2 preselPh->Sort();
220 paus 1.27 for (unsigned int iPh = 0; iPh <preselPh->GetEntries(); ++iPh) {
221 fabstoec 1.2 const Photon* ph = preselPh->At(iPh);
222 fabstoec 1.1 preselCat.push_back(PhotonTools::CiCBaseLineCat(ph));
223     }
224 paus 1.27
225 fabstoec 1.1 // ------------------------------------------------------------
226     // compute how many pairs there are ...
227     unsigned int numPairs = 0;
228 paus 1.27 if (preselPh->GetEntries() > 0)
229     numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
230    
231 fabstoec 1.1 // ... and create all possible pairs of pre-selected photons
232 paus 1.27 std::vector<unsigned int> idx1st;
233     std::vector<unsigned int> idx2nd;
234 fabstoec 1.1 std::vector<PhotonTools::CiCBaseLineCats> cat1st;
235     std::vector<PhotonTools::CiCBaseLineCats> cat2nd;
236 paus 1.27
237     // ... this will be used to store whether a given pair passes the cuts
238 fabstoec 1.1 std::vector<bool> pairPasses;
239 paus 1.27
240     if (numPairs > 0) {
241     for (unsigned int i1st = 0; i1st <preselPh->GetEntries() - 1; ++i1st) {
242     for (unsigned int i2nd = i1st + 1; i2nd <preselPh->GetEntries(); ++i2nd) {
243     idx1st.push_back(i1st);
244     idx2nd.push_back(i2nd);
245     pairPasses.push_back(true);
246 fabstoec 1.1 }
247     }
248     }
249    
250 paus 1.27 // ------------------------------------------------------------
251 fabstoec 1.1 // array to store the index of 'chosen Vtx' for each pair
252 bendavid 1.20 // const Vertex** theVtx = new const Vertex*[numPairs]; // holds the 'chosen' Vtx for each Pair
253 paus 1.27 // Photon** fixPh1st = new Photon*[numPairs]; // holds the 1st Photon for each Pair
254 bendavid 1.20 // Photon** fixPh2nd = new Photon*[numPairs]; // holds the 2nd photon for each Pair
255 paus 1.27
256 bendavid 1.20 std::vector<const Vertex*> theVtx; // holds the 'chosen' Vtx for each Pair
257 paus 1.27 std::vector<Photon*> fixPh1st; // holds the 1st Photon for each Pair
258 bendavid 1.20 std::vector<Photon*> fixPh2nd; // holds the 2nd photon for each Pair
259    
260     theVtx.reserve(numPairs);
261     fixPh1st.reserve(numPairs);
262     fixPh2nd.reserve(numPairs);
263 fabstoec 1.1
264     // store pair-indices for pairs passing the selection
265     std::vector<unsigned int> passPairs;
266     passPairs.resize(0);
267    
268 paus 1.27 // ------------------------------------------------------------
269     // Loop over all Pairs and do the 'incredible machine' running....
270     for (unsigned int iPair = 0; iPair < numPairs; ++iPair) {
271 fabstoec 1.1 // first we need a hard copy of the incoming photons
272 bendavid 1.20 fixPh1st.push_back(new Photon(*preselPh->At(idx1st[iPair])));
273     fixPh2nd.push_back(new Photon(*preselPh->At(idx2nd[iPair])));
274 fabstoec 1.1 // we also store the category, so we don't have to ask all the time...
275     cat1st.push_back(preselCat[idx1st[iPair]]);
276 paus 1.27 cat2nd.push_back(preselCat[idx2nd[iPair]]);
277    
278 bendavid 1.19 //scale regression sigmaE in MC if activated
279     if (fDoMCErrScaling && !fIsData) {
280 paus 1.27 if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5)
281     PhotonTools::ScalePhotonError(fixPh1st[iPair],fMCErrScaleEB);
282     else
283     PhotonTools::ScalePhotonError(fixPh1st[iPair],fMCErrScaleEE);
284    
285     if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5)
286     PhotonTools::ScalePhotonError(fixPh2nd[iPair],fMCErrScaleEB);
287     else
288     PhotonTools::ScalePhotonError(fixPh2nd[iPair],fMCErrScaleEE);
289     }
290 bendavid 1.19
291     //scale R9 in Monte Carlo if activated
292     if (fDoMCR9Scaling && !fIsData) {
293 paus 1.27 if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5)
294     PhotonTools::ScalePhotonR9(fixPh1st[iPair],fMCR9ScaleEB);
295     else
296     PhotonTools::ScalePhotonR9(fixPh1st[iPair],fMCR9ScaleEE);
297    
298     if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5)
299     PhotonTools::ScalePhotonR9(fixPh2nd[iPair],fMCR9ScaleEB);
300     else
301     PhotonTools::ScalePhotonR9(fixPh2nd[iPair],fMCR9ScaleEE);
302 bendavid 1.19 }
303 paus 1.27
304 bendavid 1.23 if (fDoMCSigIEtaIEtaScaling && !fIsData) {
305     if (fixPh1st[iPair]->SCluster()->AbsEta()<1.5) fixPh1st[iPair]->SetCoviEtaiEta(0.87*fixPh1st[iPair]->CoviEtaiEta() + 0.0011);
306     else fixPh1st[iPair]->SetCoviEtaiEta(0.99*fixPh1st[iPair]->CoviEtaiEta());
307 paus 1.27
308 bendavid 1.23 if (fixPh2nd[iPair]->SCluster()->AbsEta()<1.5) fixPh2nd[iPair]->SetCoviEtaiEta(0.87*fixPh2nd[iPair]->CoviEtaiEta() + 0.0011);
309 paus 1.27 else fixPh2nd[iPair]->SetCoviEtaiEta(0.99*fixPh2nd[iPair]->CoviEtaiEta());
310 bendavid 1.23 }
311 bendavid 1.19
312    
313 bendavid 1.25 if (fDoMCWidthScaling && !fIsData) {
314     fixPh1st[iPair]->SetEtaWidth(0.99*fixPh1st[iPair]->EtaWidth());
315     fixPh1st[iPair]->SetPhiWidth(0.99*fixPh1st[iPair]->PhiWidth());
316 paus 1.27
317 bendavid 1.25 fixPh2nd[iPair]->SetEtaWidth(0.99*fixPh2nd[iPair]->EtaWidth());
318     fixPh2nd[iPair]->SetPhiWidth(0.99*fixPh2nd[iPair]->PhiWidth());
319     }
320    
321 bendavid 1.19 PhotonTools::eScaleCats escalecat1 = PhotonTools::EScaleCat(fixPh1st[iPair]);
322     PhotonTools::eScaleCats escalecat2 = PhotonTools::EScaleCat(fixPh2nd[iPair]);
323 paus 1.27
324 fabstoec 1.1 // now we dicide if we either scale (Data) or Smear (MC) the Photons
325     if (fIsData) {
326 paus 1.27 if (fDoDataEneCorr) {
327     // starting with scale = 1.
328     double scaleFac1 = 1.;
329     double scaleFac2 = 1.;
330    
331 bendavid 1.9 //eta-dependent corrections
332 paus 1.27
333     // checking the run Rangees ...
334     Int_t runRange = FindRunRangeIdx(runNumber);
335     if(runRange > -1) {
336 bendavid 1.20 scaleFac1 *= GetDataEnCorr(runRange, escalecat1);
337     scaleFac2 *= GetDataEnCorr(runRange, escalecat2);
338 paus 1.27 }
339     PhotonTools::ScalePhoton(fixPh1st[iPair], scaleFac1);
340     PhotonTools::ScalePhoton(fixPh2nd[iPair], scaleFac2);
341 fabstoec 1.1 }
342 paus 1.27 }
343    
344     if (fDoMCSmear) {
345 bendavid 1.19 double width1 = GetMCSmearFac(escalecat1);
346     double width2 = GetMCSmearFac(escalecat2);
347 bendavid 1.9 if (!fIsData) {
348     // get the seed to do deterministic smearing...
349     UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
350 paus 1.27 UInt_t seed1 = seedBase + (UInt_t) fixPh1st[iPair]->E() +
351     (UInt_t) (TMath::Abs(10.*fixPh1st[iPair]->SCluster()->Eta()));
352     UInt_t seed2 = seedBase + (UInt_t) fixPh2nd[iPair]->E() +
353     (UInt_t) (TMath::Abs(10.*fixPh2nd[iPair]->SCluster()->Eta()));
354 bendavid 1.9 // get the smearing for MC photons..
355 paus 1.27 PhotonTools::SmearPhoton(fixPh1st[iPair], fRng, width1, seed1);
356     PhotonTools::SmearPhoton(fixPh2nd[iPair], fRng, width2, seed2);
357 fabstoec 1.1 }
358 bendavid 1.9 PhotonTools::SmearPhotonError(fixPh1st[iPair], width1);
359     PhotonTools::SmearPhotonError(fixPh2nd[iPair], width2);
360 fabstoec 1.1 }
361    
362 bendavid 1.19 //probability that selected vertex is the correct one
363     Double_t vtxProb = 1.0;
364    
365 fabstoec 1.1 // store the vertex for this pair
366 paus 1.27 switch (fVtxSelType) {
367 fabstoec 1.28
368 fabstoec 1.1 case kStdVtxSelection:
369     theVtx[iPair] = fPV->At(0);
370     break;
371 fabstoec 1.28
372 fabstoec 1.1 case kCiCVtxSelection:
373 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
374     fConversions,kFALSE,vtxProb);
375 fabstoec 1.1 break;
376 fabstoec 1.28
377 bendavid 1.19 case kCiCMVAVtxSelection:
378 paus 1.27 theVtx[iPair] = fVtxTools.findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV,
379     fConversions,kTRUE,vtxProb);
380     break;
381 fabstoec 1.28
382 fabstoec 1.1 case kMITVtxSelection:
383     // need PFCandidate Collection
384 paus 1.27 theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp,
385     mithep::FourVector((fixPh1st[iPair]->Mom()+
386     fixPh2nd[iPair]->Mom())));
387 fabstoec 1.1 break;
388     default:
389     theVtx[iPair] = fPV->At(0);
390 paus 1.27 }
391 fabstoec 1.2
392 bendavid 1.8 //set PV ref in photons
393     fixPh1st[iPair]->SetPV(theVtx[iPair]);
394     fixPh2nd[iPair]->SetPV(theVtx[iPair]);
395 bendavid 1.19 fixPh1st[iPair]->SetVtxProb(vtxProb);
396     fixPh2nd[iPair]->SetVtxProb(vtxProb);
397 paus 1.27
398 fabstoec 1.1 // fix the kinematics for both events
399     FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(theVtx[iPair]->Position());
400     FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(theVtx[iPair]->Position());
401     fixPh1st[iPair]->SetMom(newMom1st.X(), newMom1st.Y(), newMom1st.Z(), newMom1st.E());
402     fixPh2nd[iPair]->SetMom(newMom2nd.X(), newMom2nd.Y(), newMom2nd.Z(), newMom2nd.E());
403    
404 bendavid 1.19 double pairmass = (fixPh1st[iPair]->Mom() + fixPh2nd[iPair]->Mom()).M();
405 paus 1.27
406 bendavid 1.20 double leadptcut = fLeadingPtMin;
407     double trailptcut = fTrailingPtMin;
408 paus 1.27
409 bendavid 1.20 if (fixPh2nd[iPair]->Pt() > fixPh1st[iPair]->Pt()) {
410     leadptcut = fTrailingPtMin;
411     trailptcut = fLeadingPtMin;
412     }
413 paus 1.27
414    
415 bendavid 1.20 if (fRelativePtCuts) {
416     leadptcut = leadptcut*pairmass;
417     trailptcut = trailptcut*pairmass;
418     }
419 paus 1.27
420 mingyang 1.15
421 bendavid 1.21 //compute id bdt values
422 paus 1.27 Double_t bdt1 = fTool.GetMVAbdtValue(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho, fElectrons, fApplyEleVeto);
423     Double_t bdt2 = fTool.GetMVAbdtValue(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho, fElectrons, fApplyEleVeto);
424 bendavid 1.21
425     fixPh1st[iPair]->SetIdMva(bdt1);
426     fixPh2nd[iPair]->SetIdMva(bdt2);
427 mingyang 1.15
428 paus 1.27
429 bendavid 1.20 //printf("applying id\n");
430 paus 1.27
431 fabstoec 1.1 // check if both photons pass the CiC selection
432     // FIX-ME: Add other possibilities....
433     bool pass1 = false;
434     bool pass2 = false;
435 paus 1.27
436     switch (fPhSelType) {
437 fabstoec 1.1 case kNoPhSelection:
438 paus 1.27 pass1 = (fixPh1st[iPair]->Pt() > leadptcut );
439     pass2 = (fixPh2nd[iPair]->Pt() > trailptcut);
440 fabstoec 1.1 break;
441     case kCiCPhSelection:
442 paus 1.27 pass1 = PhotonTools::PassCiCSelection(fixPh1st[iPair], theVtx[iPair], fTracks,
443     fElectrons, fPV, rho, leadptcut, fApplyEleVeto);
444     if (pass1)
445     pass2 = PhotonTools::PassCiCSelection(fixPh2nd[iPair], theVtx[iPair], fTracks,
446     fElectrons, fPV, rho, trailptcut, fApplyEleVeto);
447 fabstoec 1.1 break;
448 mingyang 1.15 case kMVAPhSelection://MVA
449 paus 1.27 pass1 = fixPh1st[iPair]->Pt()>leadptcut &&
450     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
451     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
452     fTool.PassMVASelection(fixPh1st[iPair],theVtx[iPair],fTracks,fPV,rho,
453     fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
454     if (pass1)
455     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
456     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
457     fTracks,theVtx[iPair],rho,fApplyEleVeto) &&
458     fTool.PassMVASelection(fixPh2nd[iPair],theVtx[iPair],fTracks,fPV,rho,
459     fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyEleVeto);
460    
461 mingyang 1.15 break;
462 fabstoec 1.1 case kMITPhSelection:
463 bendavid 1.20 // loose preselection for mva
464 paus 1.27 pass1 = fixPh1st[iPair]->Pt() > leadptcut &&
465     PhotonTools::PassSinglePhotonPresel(fixPh1st[iPair],fElectrons,fConversions,bsp,
466     fTracks,theVtx[iPair],rho,fApplyEleVeto,
467     fInvertElectronVeto);
468     if (pass1)
469     pass2 = fixPh2nd[iPair]->Pt() > trailptcut &&
470     PhotonTools::PassSinglePhotonPresel(fixPh2nd[iPair],fElectrons,fConversions,bsp,
471     fTracks,theVtx[iPair],rho,fApplyEleVeto,
472     fInvertElectronVeto);
473    
474 fabstoec 1.1 break;
475     default:
476     pass1 = true;
477     pass2 = true;
478     }
479 paus 1.27
480 bendavid 1.8 //match to good electrons if requested
481     if (fInvertElectronVeto) {
482     pass1 &= !PhotonTools::PassElectronVeto(fixPh1st[iPair],fGoodElectrons);
483     pass2 &= !PhotonTools::PassElectronVeto(fixPh2nd[iPair],fGoodElectrons);
484     }
485 paus 1.27 // finally, if both Photons pass the selections, add the pair to the passing pairs
486     if (pass1 && pass2)
487     passPairs.push_back(iPair);
488 fabstoec 1.1 }
489 paus 1.27
490    
491 fabstoec 1.1 // ---------------------------------------------------------------
492 paus 1.18 // ... we're almost done, stay focused...
493 fabstoec 1.1 // loop over all passing pairs and find the one with the highest sum Et
494 paus 1.27 const Vertex* vtx = NULL;
495 fabstoec 1.1 Photon* phHard = NULL;
496     Photon* phSoft = NULL;
497    
498 fabstoec 1.29 // not used at all....
499     //PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
500     //PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
501 paus 1.27
502 fabstoec 1.1 double maxSumEt = 0.;
503 paus 1.27 for (unsigned int iPair=0; iPair<passPairs.size(); ++iPair) {
504     double sumEt = fixPh1st[passPairs[iPair]]->Et()
505     + fixPh2nd[passPairs[iPair]]->Et();
506     if (sumEt > maxSumEt) {
507 fabstoec 1.1 maxSumEt = sumEt;
508 paus 1.27 phHard = fixPh1st[passPairs[iPair]];
509     phSoft = fixPh2nd[passPairs[iPair]];
510 fabstoec 1.29 //catPh1 = cat1st [passPairs[iPair]];
511     //catPh2 = cat2nd [passPairs[iPair]];
512 paus 1.27 vtx = theVtx[iPair];
513 fabstoec 1.1 }
514     }
515 bendavid 1.20
516 paus 1.27 for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
517 bendavid 1.20 if (fixPh1st[iPair]!=phHard) delete fixPh1st[iPair];
518     if (fixPh2nd[iPair]!=phSoft) delete fixPh2nd[iPair];
519     }
520 paus 1.27
521 fabstoec 1.1 // ---------------------------------------------------------------
522     // we have the Photons (*PARTY*)... compute some useful qunatities
523    
524     if(phHard && phSoft) {
525     GoodPhotons->AddOwned(phHard);
526     GoodPhotons->AddOwned(phSoft);
527     }
528    
529 fabstoec 1.29 // we also store the chosen Vtx, so later modules can use it
530 fabstoec 1.30
531     Vertex* chosenVtx = NULL;
532     if ( vtx ) {
533     chosenVtx = new Vertex( *vtx );
534     ChosenVtx->AddOwned( chosenVtx );
535     }
536 fabstoec 1.29
537 fabstoec 1.1 // sort according to pt
538     GoodPhotons->Sort();
539 fabstoec 1.30
540 fabstoec 1.1 // delete auxiliary photon collection...
541     delete preselPh;
542 bendavid 1.20 //delete[] theVtx;
543 paus 1.27
544 fabstoec 1.1 return;
545     }
546    
547 paus 1.27 //---------------------------------------------------------------------------------------------------
548 fabstoec 1.1 void PhotonPairSelector::SlaveBegin()
549     {
550 paus 1.27 // Run startup code on the computer (slave) doing the actual analysis. Here, we just request the
551     // photon collection branch.
552 fabstoec 1.1
553 paus 1.27 // load all branches
554     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
555     ReqEventObject(fTrackBranchName, fTracks, true);
556     ReqEventObject(fElectronName, fElectrons, true);
557     ReqEventObject(fGoodElectronName, fGoodElectrons,fGoodElectronsFromBranch);
558     ReqEventObject(fPileUpDenName, fPileUpDen, true);
559     ReqEventObject(fPVName, fPV, fPVFromBranch);
560     ReqEventObject(fConversionName, fConversions, true);
561     ReqEventObject(fBeamspotName, fBeamspot, true);
562     ReqEventObject(fPFCandName, fPFCands, true);
563 fabstoec 1.1 if (!fIsData) {
564     ReqBranch(fPileUpName, fPileUp);
565     ReqBranch(fMCParticleName, fMCParticles);
566     }
567 paus 1.27 // determine photon selection type
568     if (fPhotonSelType.CompareTo("CiCSelection") == 0)
569 fabstoec 1.1 fPhSelType = kCiCPhSelection;
570 mingyang 1.15 else if (fPhotonSelType.CompareTo("MVASelection") == 0) //MVA
571     fPhSelType = kMVAPhSelection;
572 paus 1.27 else if (fPhotonSelType.CompareTo("MITSelection") == 0)
573 fabstoec 1.1 fPhSelType = kMITPhSelection;
574 paus 1.27 else
575 fabstoec 1.1 fPhSelType = kNoPhSelection;
576    
577 paus 1.27 if (fVertexSelType.CompareTo("CiCSelection") == 0)
578 fabstoec 1.1 fVtxSelType = kCiCVtxSelection;
579 paus 1.27 else if (fVertexSelType.CompareTo("MITSelection") == 0)
580 fabstoec 1.1 fVtxSelType = kMITVtxSelection;
581 paus 1.27 else if (fVertexSelType.CompareTo("CiCMVASelection") == 0)
582     fVtxSelType = kCiCMVAVtxSelection;
583 fabstoec 1.28 else if (fVertexSelType.CompareTo("ZeroVtxSelection") == 0)
584 paus 1.27 fVtxSelType = kStdVtxSelection;
585 fabstoec 1.28 else {
586     std::cerr<<" Vertex Seclection "<<fVertexSelType<<" not implemented."<<std::endl;
587     return;
588     }
589 fabstoec 1.1
590 paus 1.27 if (fIsData)
591 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
592 paus 1.27 else
593 bendavid 1.9 fPhFixFile = gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
594 mingyang 1.15
595 paus 1.27 printf("initialize photon pair selector\n");
596 mingyang 1.15
597     fTool.InitializeMVA(fVariableType,fEndcapWeights,fBarrelWeights);
598 paus 1.27 fVtxTools.InitP();
599 mingyang 1.15
600 fabstoec 1.1 }
601    
602     // ----------------------------------------------------------------------------------------
603     // some helpfer functions....
604 paus 1.27 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
605     {
606     pt = -999.;
607 fabstoec 1.1 decayZ = -999.;
608 paus 1.27 mass = -999.;
609 fabstoec 1.1
610     // loop over all GEN particles and look for status 1 photons
611     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
612     const MCParticle* p = fMCParticles->At(i);
613 bendavid 1.8 if( p->Is(MCParticle::kH) || (!fApplyEleVeto && p->AbsPdgId()==23) ) {
614     pt=p->Pt();
615     decayZ = p->DecayVertex().Z();
616     mass = p->Mass();
617     break;
618     }
619 fabstoec 1.1 }
620 paus 1.27
621 fabstoec 1.1 return;
622 paus 1.27 }
623 fabstoec 1.1
624 paus 1.27 //---------------------------------------------------------------------------------------------------
625     Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run)
626     {
627     // this routine looks for the idx of the run-range
628 fabstoec 1.1 Int_t runRange=-1;
629 paus 1.27 for (UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
630     if (run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
631 fabstoec 1.1 runRange = (Int_t) iRun;
632     return runRange;
633     }
634     }
635     return runRange;
636     }
637    
638 paus 1.27 //---------------------------------------------------------------------------------------------------
639     Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::eScaleCats cat)
640     {
641     switch (cat) {
642 bendavid 1.19 case PhotonTools::kEBhighEtaGold:
643     return fDataEnCorr_EBhighEta_hR9[runRange];
644     case PhotonTools::kEBhighEtaBad:
645     return fDataEnCorr_EBhighEta_lR9[runRange];
646 bendavid 1.20 case PhotonTools::kEBlowEtaGoldCenter:
647     return fDataEnCorr_EBlowEta_hR9central[runRange];
648     case PhotonTools::kEBlowEtaGoldGap:
649 paus 1.27 return fDataEnCorr_EBlowEta_hR9gap[runRange];
650 bendavid 1.19 case PhotonTools::kEBlowEtaBad:
651 paus 1.27 return fDataEnCorr_EBlowEta_lR9[runRange];
652 bendavid 1.19 case PhotonTools::kEEhighEtaGold:
653     return fDataEnCorr_EEhighEta_hR9[runRange];
654     case PhotonTools::kEEhighEtaBad:
655     return fDataEnCorr_EEhighEta_lR9[runRange];
656     case PhotonTools::kEElowEtaGold:
657     return fDataEnCorr_EElowEta_hR9[runRange];
658     case PhotonTools::kEElowEtaBad:
659 paus 1.27 return fDataEnCorr_EElowEta_lR9[runRange];
660 fabstoec 1.1 default:
661     return 1.;
662     }
663     }
664    
665 paus 1.27 //---------------------------------------------------------------------------------------------------
666     Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::eScaleCats cat)
667     {
668     switch (cat) {
669 bendavid 1.19 case PhotonTools::kEBhighEtaGold:
670     return fMCSmear_EBhighEta_hR9;
671     case PhotonTools::kEBhighEtaBad:
672     return fMCSmear_EBhighEta_lR9;
673 bendavid 1.20 case PhotonTools::kEBlowEtaGoldCenter:
674     return fMCSmear_EBlowEta_hR9central;
675     case PhotonTools::kEBlowEtaGoldGap:
676 paus 1.27 return fMCSmear_EBlowEta_hR9gap;
677 bendavid 1.19 case PhotonTools::kEBlowEtaBad:
678 paus 1.27 return fMCSmear_EBlowEta_lR9;
679 bendavid 1.19 case PhotonTools::kEEhighEtaGold:
680     return fMCSmear_EEhighEta_hR9;
681     case PhotonTools::kEEhighEtaBad:
682     return fMCSmear_EEhighEta_lR9;
683     case PhotonTools::kEElowEtaGold:
684     return fMCSmear_EElowEta_hR9;
685     case PhotonTools::kEElowEtaBad:
686 paus 1.27 return fMCSmear_EElowEta_lR9;
687 fabstoec 1.1 default:
688     return 1.;
689     }
690     }
691    
692 paus 1.27 //---------------------------------------------------------------------------------------------------
693     Float_t PhotonPairSelector::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
694     PhotonTools::CiCBaseLineCats cat2)
695     {
696     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
697     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
698 fabstoec 1.1
699     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
700     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
701 paus 1.27
702     if (ph1IsEB && ph2IsEB)
703     return (ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
704    
705     return (ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
706 fabstoec 1.1 }
707 bendavid 1.4