ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionFuncs.cc
Revision: 1.15
Committed: Tue Jan 15 10:35:23 2013 UTC (12 years, 4 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.14: +0 -4 lines
Log Message: 
clean up

File Contents

# User Rev Content
1 dkralph 1.1 #include "SelectionFuncs.h"
2    
3     extern vector<SimpleLepton> failingLeptons;
4     extern vector<SimpleLepton> passingLeptons;
5    
6     extern vector<unsigned> cutvec;
7     extern vector<vector<unsigned> > zcutvec;
8     extern vector<vector<unsigned> > zzcutvec;
9     extern map<unsigned,float> evtrhoMap;
10     extern bool passes_HLT_MC;
11 dkralph 1.12 extern RunLumiRangeMap rlrm;
12 khahn 1.5
13 dkralph 1.13 //----------------------------------------------------------------------------------------
14     void increment(map<TString, map<TString,int>* > &counterMap, TString cutName, int z1type, int z2type)
15     //----------------------------------------------------------------------------------------
16     // increment counters that keep track of event counts at each cut
17     {
18     assert(counterMap.find(cutName) != counterMap.end());
19     (*(counterMap[cutName]))["all"] ++;
20     if(z1type!=0 && z2type!=0) {
21     assert(counterMap.find(cutName) != counterMap.end());
22     (*(counterMap[cutName]))[getChannel(z1type,z2type)] ++;
23     }
24     }
25 dkralph 1.1 //----------------------------------------------------------------------------
26 dkralph 1.12 bool setPV(ControlFlags &ctrl,
27 dkralph 1.1 const mithep::Array<mithep::Vertex> * vtxArr,
28     const mithep::Vertex* &vtx)
29     //----------------------------------------------------------------------------
30 dkralph 1.6 {
31 dkralph 1.1
32     const mithep::Vertex *bestPV = 0;
33     float best_sumpt=-1;
34    
35     // good PV requirements
36     const UInt_t fMinNTracksFit = 0;
37     const Double_t fMinNdof = 4;
38     const Double_t fMaxAbsZ = 24;
39     const Double_t fMaxRho = 2;
40    
41     for(int i=0; i<vtxArr->GetEntries(); ++i) {
42     const mithep::Vertex *pv = (mithep::Vertex*)(vtxArr->At(i));
43     if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
44    
45     // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
46     if(!pv->IsValid()) continue;
47     if(pv->NTracksFit() < fMinNTracksFit) continue;
48     if(pv->Ndof() < fMinNdof) continue;
49     if(fabs(pv->Z()) > fMaxAbsZ) continue;
50     if(pv->Position().Rho() > fMaxRho) continue;
51    
52     // take the first ...
53     bestPV = pv;
54     break;
55    
56     // this never reached ...
57     float tmp_sumpt=0;
58     for( int t=0; t<pv->NTracks(); t++ )
59     tmp_sumpt += pv->Trk(t)->Pt();
60    
61     if( tmp_sumpt > best_sumpt ) {
62     bestPV = pv;
63     best_sumpt = tmp_sumpt;
64     if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
65     }
66     }
67    
68     // sync
69     if(!bestPV)
70     return false;
71     else {
72     vtx = bestPV;
73     return true;
74     }
75     }
76     //----------------------------------------------------------------------------------------
77     void writeEntries(ControlFlags ctrl, unsigned total_unskimmed)
78     {
79     ofstream ofs;
80     TString ofname(ctrl.outputfile);
81     ofs.open(ofname.ReplaceAll(".root",".nevents"));
82     assert(ofs.is_open());
83     ofs << total_unskimmed << endl;
84     ofs.close();
85     }
86     //----------------------------------------------------------------------------------------
87     void getEATargets(ControlFlags &ctrl, mithep::MuonTools::EMuonEffectiveAreaTarget &eraMu, mithep::ElectronTools::EElectronEffectiveAreaTarget &eraEle)
88     {
89     if( !ctrl.mc && ctrl.era == 2011 ) {
90     eraMu = mithep::MuonTools::kMuEAData2011;
91     eraEle = mithep::ElectronTools::kEleEAData2011;
92     } else if( !ctrl.mc && ctrl.era == 2012 ) {
93     eraMu = mithep::MuonTools::kMuEAData2012;
94     eraEle = mithep::ElectronTools::kEleEAData2012;
95     } else if( ctrl.mc && ctrl.era == 2011 ) {
96     eraMu = mithep::MuonTools::kMuEAFall11MC;
97     eraEle = mithep::ElectronTools::kEleEAFall11MC;
98     } else if( ctrl.mc && ctrl.era == 2012 ) {
99     eraMu = mithep::MuonTools::kMuEAData2012;
100     eraEle = mithep::ElectronTools::kEleEAData2012;
101     } else {
102     cerr << "unknown era for effective areas ... quitting." << endl;
103     exit(1);
104     }
105     }
106     //----------------------------------------------------------------------------
107 anlevin 1.10 unsigned makePFnoPUArray(const mithep::Array<PFCandidate> * fPFCandidates,
108 dkralph 1.1 vector<bool> &pfNoPileUpflag,
109     const mithep::Array<mithep::Vertex> * vtxArr )
110     //----------------------------------------------------------------------------
111 dkralph 1.13 //
112     // Make a vector of ints, one for each pf candidate, where we only include it in
113     // the isolation sum if it's 1
114     //
115 dkralph 1.1 {
116 dkralph 1.13 assert(pfNoPileUpflag.size() == 0); // make sure we remembered to clear it after the last event
117 dkralph 1.1 for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
118     const PFCandidate *pf = fPFCandidates->At(i);
119     assert(pf);
120    
121     if(pf->PFType() == PFCandidate::eHadron) {
122     if(pf->HasTrackerTrk() &&
123     vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
124     vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
125    
126     pfNoPileUpflag.push_back(1);
127    
128     } else {
129    
130     Bool_t vertexFound = kFALSE;
131     const Vertex *closestVtx = 0;
132     Double_t dzmin = 10000;
133    
134     // loop over vertices
135     for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
136     const Vertex *vtx = vtxArr->At(j);
137     assert(vtx);
138    
139     if(pf->HasTrackerTrk() &&
140     vtx->HasTrack(pf->TrackerTrk()) &&
141     vtx->TrackWeight(pf->TrackerTrk()) > 0) {
142     vertexFound = kTRUE;
143     closestVtx = vtx;
144     break;
145     }
146     Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
147     if(dz < dzmin) {
148     closestVtx = vtx;
149     dzmin = dz;
150     }
151     }
152    
153     // if(fCheckClosestZVertex) {
154     if(1) {
155     // Fallback: if track is not associated with any vertex,
156     // associate it with the vertex closest in z
157     if(vertexFound || closestVtx != vtxArr->At(0)) {
158     // pfPileUp->Add(pf);
159 dkralph 1.12 pfNoPileUpflag.push_back(0); // i.e., don't use it in the isolation sum
160 dkralph 1.1 } else {
161 dkralph 1.12 pfNoPileUpflag.push_back(1); // use it
162 dkralph 1.1 }
163     } else {
164     if(vertexFound && closestVtx != vtxArr->At(0)) {
165     // pfPileUp->Add(pf);
166     pfNoPileUpflag.push_back(0);
167     } else {
168     // PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
169     pfNoPileUpflag.push_back(1);
170     }
171     }
172     } //hadron & trk stuff
173     } else { // hadron
174     // PFCandidate * ptr = pfNoPileUp->AddNew();
175     pfNoPileUpflag.push_back(1);
176     }
177     } // Loop over PF candidates
178    
179     return pfNoPileUpflag.size();
180     }
181 khahn 1.4 //--------------------------------------------------------------------------------
182 dkralph 1.12 void setEfficiencyWeights(ControlFlags &ctrl,
183     unsigned era,
184 khahn 1.8 EventData & evtdat,
185     std::bitset<1024> triggerBits,
186     mithep::TriggerTable *hltTable,
187     mithep::Array<mithep::TriggerObject> *hltObjArr,
188     mithep::TriggerObjectsTable *fTrigObjs,
189     WeightStruct & weights )
190 khahn 1.4 //--------------------------------------------------------------------------------
191 dkralph 1.1 {
192 dkralph 1.12 //
193     // offline weights
194     //
195 dkralph 1.1 vector<SimpleLepton> lepvec = evtdat.Z1leptons;
196     lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
197 khahn 1.8
198 dkralph 1.12 // get eff for each leg, for each lepton
199     vector< pair <double,double> > wlegs; // pair here is eff & err
200     for( int k=0; k<lepvec.size(); k++ ) {
201     double absEta = fabs(lepvec[k].vec.Eta());
202     double pt = lepvec[k].vec.Pt();
203     if( abs(lepvec[k].type) == 13 )
204     wlegs.push_back( muonPerLegOfflineEfficiencyWeight( ctrl, era, absEta, pt ) );
205     else
206     wlegs.push_back( elePerLegOfflineEfficiencyWeight( ctrl, era, absEta, pt, lepvec[k].isTight ) );
207     }
208     // combine weights and propagate uncertainties
209     pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
210     weights.woff = offpair.first;
211     weights.werroff = offpair.second;
212    
213     //
214     // online weights
215     //
216 khahn 1.8 vector<SimpleLepton> muvec, elvec;
217     if(abs(evtdat.Z1leptons[0].type) == 11 ) {
218     elvec.push_back(evtdat.Z1leptons[0]);
219     elvec.push_back(evtdat.Z1leptons[1]);
220     } else {
221     muvec.push_back(evtdat.Z1leptons[0]);
222     muvec.push_back(evtdat.Z1leptons[1]);
223     }
224     if(abs(evtdat.Z2leptons[0].type) == 11 ) {
225     elvec.push_back(evtdat.Z2leptons[0]);
226     elvec.push_back(evtdat.Z2leptons[1]);
227     } else {
228     muvec.push_back(evtdat.Z2leptons[0]);
229     muvec.push_back(evtdat.Z2leptons[1]);
230     }
231 dkralph 1.12 pair<double,double> onpair = getOnlineEfficiencyWeight( ctrl, triggerBits, hltTable, hltObjArr, fTrigObjs, muvec, elvec );
232 dkralph 1.1 weights.won = onpair.first;
233     weights.werron = onpair.second;
234     }
235     //----------------------------------------------------------------------------
236     void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
237     //----------------------------------------------------------------------------
238     {
239     unsigned evt;
240     float rho;
241    
242     cout << "initialzing EvtRhoMap ";
243     //FILE * f = fopen("evtrho.dat", "r");
244     // FILE * f = fopen("allrho.cali.uniq.dat", "r");
245     // FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
246     assert(fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat","r")); FILE * f = fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
247     int lcount=0;
248     while( fscanf( f, "%u:%f", &evt, &rho ) ) {
249     if( feof(f) ) break;
250     evtrhoMap[evt] = rho;
251     lcount++;
252     if( !(lcount%1000) ) cout << "."; flush(cout);
253     }
254     cout << " done" << endl;
255     }
256    
257 khahn 1.5 //----------------------------------------------------------------------------------------
258 dkralph 1.1 void setEra(string fname, ControlFlags &ctrl)
259     {
260 khahn 1.5 //----------------------------------------------------------------------------------------
261 dkralph 1.1 if( ctrl.mc && (fname.find("f11-") != string::npos) ) ctrl.era=2011 ;
262 dkralph 1.9 else if( ctrl.mc && (fname.find("s12-") != string::npos) ) ctrl.era=2012 ;
263     else if( !ctrl.mc && (fname.find("r11") != string::npos) ) ctrl.era=2011 ;
264     else if( !ctrl.mc && (fname.find("r12") != string::npos) ) ctrl.era=2012 ;
265     else assert(0);
266 dkralph 1.1 }
267 khahn 1.5
268 dkralph 1.1 //----------------------------------------------------------------------------------------
269 dkralph 1.6 unsigned getNPU(mithep::Array<mithep::PileupInfo> *puArr, int bx)
270 khahn 1.5 //----------------------------------------------------------------------------------------
271 dkralph 1.1 {
272    
273     // Identify bunch crossing 0 in PU info
274     int ibx0=-1, ibxm=-1, ibxp=-1;
275     assert(puArr);
276     for(unsigned i=0; i<puArr->GetEntries(); ++i) {
277     if(puArr->At(i)->GetBunchCrossing() == 0) ibx0=i;
278     if(puArr->At(i)->GetBunchCrossing() ==-1) ibxm=i;
279     if(puArr->At(i)->GetBunchCrossing() == 1) ibxp=i;
280     }
281    
282     if(bx==0) return puArr->At(ibx0)->GetPU_NumInteractions();
283     else if(bx==-1) return puArr->At(ibxm)->GetPU_NumInteractions();
284     else if(bx== 1) return puArr->At(ibxp)->GetPU_NumInteractions();
285     else assert(0);
286    
287 dkralph 1.6 return 0;
288 dkralph 1.1 }
289 khahn 1.5
290    
291     //----------------------------------------------------------------------------
292     void initPUWeights()
293     //----------------------------------------------------------------------------
294     {
295 dkralph 1.12 // TFile * puf;
296 khahn 1.5
297 dkralph 1.12 // puf = new TFile("data/pileup/PUWeights_F11To2011.root");
298     // hpu_2011 = (TH1D*)(puf->Get("puWeights"));
299     // hpu_2011->SetDirectory(0);
300     // puf->Close();
301    
302     // puf = new TFile("data/pileup/PUWeights_S12To2012_190456-199011.root");
303     // hpu_2012 = (TH1D*)(puf->Get("puWeights"));
304     // hpu_2012->SetDirectory(0);
305     // puf->Close();
306 khahn 1.5 }
307    
308     //----------------------------------------------------------------------------
309 dkralph 1.12 double getPUWeight(int era, TString fname, unsigned npu)
310 khahn 1.5 //----------------------------------------------------------------------------
311     {
312 dkralph 1.12 if(npu > 59) npu = 59; // hardcoded crap in the pu functions
313    
314     if(era==2011)
315     return weightTrue2011(npu);
316     else assert(era==2012);
317    
318     // assume that if it's in /029/ that we want to reweigh to 53x data
319     TString vers,dataVers,mcVers;
320     if(fname.Contains("/028/")) dataVers = mcVers = vers = "52x";
321     else if(fname.Contains("/029/")) dataVers = mcVers = vers = "53x";
322     else {
323     // if(ctrl.debug)
324     // cout << "WARNING in getPUWeight(): couldn't figure out if this is 52x or 53x, so just using 53x" << endl;
325     dataVers = mcVers = vers = "53x";
326 khahn 1.7 }
327 dkralph 1.12
328     if(dataVers=="52x" && mcVers=="52x")
329     return weightTruePileupV07toIchep52X(npu);
330     else if(dataVers=="53x" && mcVers=="53x")
331     return weightTruePileupV10toHcp53X(npu);
332     // else if(dataVers=="53x" && mcVers=="52x") // umm... can't really do this a.t.m..., but this would be the proper function
333     // return weightTruePileupV07toHcp53X(npu);
334     else { cout << "Couldn't work out which npu to use!" << endl; assert(0); };
335    
336     return 1;
337 khahn 1.5 }
338 dkralph 1.12 //----------------------------------------------------------------------------------------
339     void initJets(ControlFlags ctrl,
340     string cmsswpath,
341     vector<TString> &fJetCorrParsv,
342     vector<JetCorrectorParameters> &correctionParameters,
343     FactorizedJetCorrector *&fJetCorrector,
344     JetIDMVA *&fJetIDMVA)
345     //----------------------------------------------------------------------------------------
346     {
347     //
348     // jet energy corrections
349     //
350     vector<TString> corrTypes;
351     corrTypes.push_back("L1FastJet");
352     corrTypes.push_back("L2Relative");
353     corrTypes.push_back("L3Absolute");
354     if(!ctrl.mc)
355     corrTypes.push_back("L2L3Residual");
356    
357     TString tag;
358     if(ctrl.era==2011) {
359     cout << "\n\nNOTE: dont have the right JEC tags yet!\n\n" << endl;
360     // // from h4l hcp twiki:
361     if(ctrl.mc) tag = "START42_V17"; //START44_V13
362     else tag = "GR_R_42_V23"; //GR_R_44_V15
363     } else if(ctrl.era==2012) {
364     cout << "\n\nNOTE: dont have the right JEC tags yet!\n\n" << endl;
365     // // from h4l hcp twiki:
366     if(ctrl.mc) tag = "START53_V7F"; //START53_V10
367     else tag = "GR_P_V41_AN2"; //GR_P_V41_AN1
368     assert(TString(cmsswpath).Contains("CMSSW_5_3"));
369     } else assert(0);
370    
371     TString dataPath = TString(cmsswpath) + "/MitPhysics/data/";
372     for(unsigned icorr=0; icorr<corrTypes.size(); icorr++) {
373     if(ctrl.era==2011) corrTypes[icorr] = "AK5PF_" + corrTypes[icorr];
374     else if(ctrl.era==2012) corrTypes[icorr] = corrTypes[icorr] + "_AK5PF";
375     else assert(0);
376     TString fname(dataPath + tag + "_" + corrTypes[icorr] + ".txt");
377     if( !fopen(fname, "r") ) { cout << " no such file! " << fname << endl; assert(0); }
378     addJetCorr(fname, fJetCorrParsv);
379     }
380 khahn 1.5
381 dkralph 1.12 for(unsigned icorr=0; icorr<fJetCorrParsv.size(); icorr++)
382     correctionParameters.push_back(JetCorrectorParameters(fJetCorrParsv[icorr].Data()));
383     fJetCorrector = new FactorizedJetCorrector(correctionParameters);
384    
385     //
386     // jet ID
387     //
388     fJetIDMVA = new JetIDMVA();
389     JetIDMVA::MVAType mvaType;
390     TString loPtXml(dataPath+"mva_JetID_lowpt.weights.xml");
391     TString cfi(cmsswpath+"/MitPhysics/Utils/python/JetIdParams_cfi.py");
392     TString Xml;
393     if(ctrl.era==2011) {
394     mvaType = JetIDMVA::k42;
395     Xml = dataPath+"TMVAClassification_PuJetIdOptMVA.weights.xml";
396     } else if(ctrl.era==2012) {
397     mvaType = JetIDMVA::k52;
398     Xml = dataPath+"TMVAClassification_5x_BDT_fullPlusRMS.weights.xml";
399     } else assert(0);
400    
401     fJetIDMVA->Initialize(JetIDMVA::kLoose, loPtXml, Xml, mvaType, cfi);
402     }
403     //----------------------------------------------------------------------------------------
404     int getGenSampleType(TString ofname)
405     {
406     if(ofname.Contains(TRegexp("-zz4[em]-"))) return ESampleType::kVV; // powheg ZZ
407     else if(ofname.Contains(TRegexp("-zz2e2m-"))) return ESampleType::kVV; // ggZZ
408     else if(ofname.Contains(TRegexp("-zz4l-gf"))) return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV; // ggZZ
409     else if(ofname.Contains(TRegexp("-zz2l2l-gf"))) return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV; // ggZZ
410     else if(ofname.Contains(TRegexp("-ggzz4l-"))) return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV; // ggZZ
411     else if(ofname.Contains(TRegexp("-ggzz2l2l-"))) return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV; // ggZZ
412 andrey 1.14 else if(ofname.Contains(TRegexp("-x12[0-9]zz4l-"))) return ESampleType::kHZZJHU; // spin samples
413 dkralph 1.12 else if(ofname.Contains(TRegexp("-h[0-9][0-9]*zz4l-"))) return ESampleType::kHZZ; // powheg HZZ
414 dkralph 1.13 else if(ofname.Contains(TRegexp("-h[0-9][0-9]*-"))) return ESampleType::kHZZ; // powheg HZZ where I fucked up the name
415     else if(ofname.Contains(TRegexp("-zllm50-"))) return ESampleType::kZ; // drell-yan
416     else if(ofname.Contains(TRegexp("-zllm1050-"))) return ESampleType::kZ; // drell-yan
417     else if(ofname.Contains(TRegexp("-zjets-"))) return ESampleType::kZ; // drell-yan
418 dkralph 1.12 else assert(0);
419     }
420     //----------------------------------------------------------------------------------------
421 anlevin 1.11 TString getChannel(int z1type, int z2type)
422     {
423     if( abs(z1type) == 11 && abs(z2type) == 11 ) return "4e";//0;
424     else if( abs(z1type) == 13 && abs(z2type) == 13 ) return "4m";//1;
425     else return "2e2m";//2;
426     }
427 dkralph 1.12 //----------------------------------------------------------------------------------------
428     // START: functions from Mangano
429     //----------------------------------------------------------------------------------------
430     float weightTruePileupV07toIchep52X(float input){
431     float w[60] = {
432     2.59113e-07,
433     1.77107e-05,
434     0.0263017,
435     1.66396,
436     0.16421,
437     0.754166,
438     2.84622,
439     5.57103,
440     8.66558,
441     11.5716,
442     11.8712,
443     12.8102,
444     10.3421,
445     8.91019,
446     7.614,
447     6.10397,
448     4.52745,
449     3.26031,
450     2.39558,
451     1.83297,
452     1.47821,
453     1.26728,
454     1.14716,
455     1.06707,
456     0.98066,
457     0.860877,
458     0.708281,
459     0.539789,
460     0.37652,
461     0.237298,
462     0.1338,
463     0.0671236,
464     0.0299236,
465     0.0118723,
466     0.00420968,
467     0.00134235,
468     0.000389563,
469     0.000104892,
470     2.69214e-05,
471     6.79674e-06,
472     1.73307e-06,
473     4.52553e-07,
474     1.21124e-07,
475     3.29924e-08,
476     9.10616e-09,
477     2.53998e-09,
478     7.16146e-10,
479     2.03786e-10,
480     5.84308e-11,
481     1.68192e-11,
482     4.8434e-12,
483     1.38959e-12,
484     3.96112e-13,
485     1.11358e-13,
486     3.17245e-14,
487     5.34916e-15,
488     0,
489     0,
490     0,
491     0};
492    
493     //return w[(int)floor(input+0.5)];
494     return w[(int)floor(input)];
495    
496     /*
497     TH1F h("boh","boh",60,0.,60.);
498     for(int k=0;k<60;k++){
499     h.SetBinContent(k+1,w[k]);
500     }
501     return h.GetBinContent(h.FindBin(input));
502     */
503     }
504    
505    
506     float weightTruePileupV07toHcp53X(float input){
507     float w[60] = {
508     0.0447136,
509     0.11785,
510     0.23825,
511     1.08447,
512     0.102575,
513     0.454605,
514     1.79761,
515     4.00271,
516     6.83281,
517     9.83701,
518     10.7966,
519     12.2356,
520     10.0247,
521     8.49395,
522     7.1125,
523     5.69527,
524     4.31256,
525     3.19305,
526     2.42035,
527     1.91666,
528     1.58485,
529     1.36297,
530     1.21166,
531     1.09466,
532     0.978941,
533     0.84653,
534     0.699235,
535     0.548996,
536     0.408673,
537     0.288194,
538     0.193367,
539     0.124653,
540     0.0781124,
541     0.0479268,
542     0.0287763,
543     0.0167744,
544     0.00941834,
545     0.00507877,
546     0.00264364,
547     0.00134612,
548     0.000682678,
549     0.000351412,
550     0.0001864,
551     0.00010259,
552     5.87818e-05,
553     3.5033e-05,
554     2.17116e-05,
555     1.39777e-05,
556     9.36123e-06,
557     6.53328e-06,
558     4.76598e-06,
559     3.64139e-06,
560     2.92018e-06,
561     2.4602e-06,
562     2.17291e-06,
563     2.01107e-06,
564     1.94392e-06,
565     1.9598e-06,
566     2.0583e-06,
567     2.24895e-06};
568    
569     return w[(int)floor(input)];
570     }
571    
572    
573     float weightTruePileupV10toIchep53X(float input){
574     float w[60] = {
575     2.35693e-06,
576     7.51928e-05,
577     0.0263529,
578     0.609947,
579     0.737917,
580     1.29365,
581     0.994503,
582     0.85454,
583     1.01559,
584     1.33243,
585     1.72454,
586     2.01264,
587     2.00573,
588     1.80333,
589     1.56328,
590     1.37452,
591     1.24753,
592     1.16481,
593     1.11738,
594     1.09701,
595     1.08843,
596     1.08796,
597     1.09768,
598     1.10763,
599     1.09328,
600     1.0339,
601     0.92408,
602     0.771537,
603     0.59283,
604     0.41266,
605     0.256892,
606     0.14188,
607     0.0692543,
608     0.029902,
609     0.0114564,
610     0.00391383,
611     0.00120625,
612     0.000341485,
613     9.09127e-05,
614     2.34008e-05,
615     5.95438e-06,
616     1.5122e-06,
617     3.82094e-07,
618     9.51794e-08,
619     2.32205e-08,
620     5.51698e-09,
621     1.27267e-09,
622     2.84346e-10,
623     6.12799e-11,
624     1.26731e-11,
625     2.50309e-12,
626     4.69797e-13,
627     8.35153e-14,
628     1.39452e-14,
629     2.24718e-15,
630     2.03841e-16,
631     0,
632     0,
633     0,
634     0};
635    
636     return w[(int)floor(input)];
637     }
638    
639    
640    
641     float weightTruePileupV10toHcp53X(float input){
642     float w[60] = {
643     0.409409,
644     0.527276,
645     0.39328,
646     0.507892,
647     0.48029,
648     0.787701,
649     0.632356,
650     0.618033,
651     0.806089,
652     1.14018,
653     1.5788,
654     1.93507,
655     1.957,
656     1.73004,
657     1.46737,
658     1.28278,
659     1.18189,
660     1.13388,
661     1.12578,
662     1.14415,
663     1.16048,
664     1.1618,
665     1.15318,
666     1.13405,
667     1.09239,
668     1.01915,
669     0.914837,
670     0.786744,
671     0.644879,
672     0.502039,
673     0.371688,
674     0.263586,
675     0.18067,
676     0.120472,
677     0.0780184,
678     0.0486113,
679     0.0289039,
680     0.0163367,
681     0.00879674,
682     0.00456046,
683     0.0023098,
684     0.00115977,
685     0.000583207,
686     0.000294815,
687     0.000149865,
688     7.62892e-05,
689     3.87537e-05,
690     1.96105e-05,
691     9.87744e-06,
692     4.95418e-06,
693     2.47913e-06,
694     1.23919e-06,
695     6.19751e-07,
696     3.10125e-07,
697     1.54934e-07,
698     7.71425e-08,
699     3.8182e-08,
700     1.87455e-08,
701     9.10765e-09,
702     9.19802e-09};
703     return w[(int)floor(input)];
704     }
705    
706    
707    
708    
709    
710    
711    
712    
713    
714     float weightTrue2011(float input){
715     if(input>50)
716     return 1;
717    
718    
719     float w[50];
720    
721    
722     w[0]= 0.212929;
723     w[1]= 0.0208114;
724     w[2]= 0.0584048;
725     w[3]= 0.538898;
726     w[4]= 1.357;
727     w[5]= 1.49913;
728     w[6]= 1.42247;
729     w[7]= 1.35904;
730     w[8]= 1.29946;
731     w[9]= 1.27925;
732     w[10]= 1.37845;
733     w[11]= 1.71246;
734     w[12]= 1.5291;
735     w[13]= 1.35234;
736     w[14]= 1.22215;
737     w[15]= 1.0155;
738     w[16]= 1.01137;
739     w[17]= 0.395465;
740     w[18]= 0.230984;
741     w[19]= 0.109883;
742     w[20]= 0.0433739;
743     w[21]= 0.0111497;
744     w[22]= 0.00408801;
745     w[23]= 0.00115678;
746     w[24]= 0.000365505;
747     w[25]= 0.000112391;
748     w[26]= 3.83894e-05;
749     w[27]= 1.60651e-05;
750     w[28]= 4.81412e-06;
751     w[29]= 1.39717e-06;
752     w[30]= 1.92368e-06;
753     w[31]= 4.10748e-06;
754     w[32]= 2.33157e-05;
755     w[33]= 4.0181e-05;
756     w[34]= 4.87786e-05;
757     w[35]= 0.00194128;
758     w[36]= 8.97414e-05;
759     w[37]= 1;
760     w[38]= 1;
761     w[39]= 0.000162709;
762     w[40]= 1;
763     w[41]= 1;
764     w[42]= 1;
765     w[43]= 1;
766     w[44]= 1;
767     w[45]= 1;
768     w[46]= 1;
769     w[47]= 1;
770     w[48]= 1;
771     w[49]= 1;
772    
773    
774     TH1F h("boh","boh",50,0.,50.);
775    
776     for(int k=0;k<50;k++){
777     h.SetBinContent(k+1,w[k]);
778     }
779    
780     return h.GetBinContent(h.FindBin(input));
781    
782     }
783     //----------------------------------------------------------------------------------------
784     float weightTrue2011to2012(float input){
785     if(input>50)
786     return 1;
787    
788     float w[50];
789    
790     w[0]= 0.000443112;
791     w[1]= 0.000248044;
792     w[2]= 0.000273111;
793     w[3]= 0.00109511;
794     w[4]= 0.00195699;
795     w[5]= 0.00480746;
796     w[6]= 0.027013;
797     w[7]= 0.074795;
798     w[8]= 0.166231;
799     w[9]= 0.309545;
800     w[10]= 0.577657;
801     w[11]= 1.12488;
802     w[12]= 1.36899;
803     w[13]= 1.56925;
804     w[14]= 1.89846;
805     w[15]= 2.20828;
806     w[16]= 3.14112;
807     w[17]= 1.87712;
808     w[18]= 1.97062;
809     w[19]= 2.07067;
810     w[20]= 2.17791;
811     w[21]= 1.7176;
812     w[22]= 2.10953;
813     w[23]= 2.0805;
814     w[24]= 2.29498;
815     w[25]= 2.42189;
816     w[26]= 2.80303;
817     w[27]= 3.94091;
818     w[28]= 3.67917;
819     w[29]= 2.26081;
820     w[30]= 2.99726;
821     w[31]= 3.76553;
822     w[32]= 11.285;
823     w[33]= 10.2781;
824     w[34]= 6.73407;
825     w[35]= 148.182;
826     w[36]= 3.88144;
827     w[37]= 1;
828     w[38]= 1;
829     w[39]= 1.48128;
830     w[40]= 1;
831     w[41]= 1;
832     w[42]= 1;
833     w[43]= 1;
834     w[44]= 1;
835     w[45]= 1;
836     w[46]= 1;
837     w[47]= 1;
838     w[48]= 1;
839     w[49]= 1;
840    
841    
842     TH1F h("boh","boh",50,0.,50.);
843    
844     for(int k=0;k<50;k++){
845     h.SetBinContent(k+1,w[k]);
846     }
847    
848     return h.GetBinContent(h.FindBin(input));
849    
850     }
851     //----------------------------------------------------------------------------------------
852     // END: functions from Mangano
853     //----------------------------------------------------------------------------------------