ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.26
Committed: Mon Dec 19 23:45:00 2011 UTC (13 years, 4 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025d
Changes since 1.25: +3 -3 lines
Log Message: 
updated id mva

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