ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionFuncs.cc
Revision: 1.13
Committed: Mon Dec 17 12:34:02 2012 UTC (12 years, 5 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled
Changes since 1.12: +22 -0 lines
Log Message:
*** empty log message ***

File Contents

# Content
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 extern RunLumiRangeMap rlrm;
12
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 //----------------------------------------------------------------------------
26 bool setPV(ControlFlags &ctrl,
27 const mithep::Array<mithep::Vertex> * vtxArr,
28 const mithep::Vertex* &vtx)
29 //----------------------------------------------------------------------------
30 {
31
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 unsigned makePFnoPUArray(const mithep::Array<PFCandidate> * fPFCandidates,
108 vector<bool> &pfNoPileUpflag,
109 const mithep::Array<mithep::Vertex> * vtxArr )
110 //----------------------------------------------------------------------------
111 //
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 {
116 assert(pfNoPileUpflag.size() == 0); // make sure we remembered to clear it after the last event
117 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 pfNoPileUpflag.push_back(0); // i.e., don't use it in the isolation sum
160 } else {
161 pfNoPileUpflag.push_back(1); // use it
162 }
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 //--------------------------------------------------------------------------------
182 void setEfficiencyWeights(ControlFlags &ctrl,
183 unsigned era,
184 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 //--------------------------------------------------------------------------------
191 {
192 //
193 // offline weights
194 //
195 vector<SimpleLepton> lepvec = evtdat.Z1leptons;
196 lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
197
198 // 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 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 pair<double,double> onpair = getOnlineEfficiencyWeight( ctrl, triggerBits, hltTable, hltObjArr, fTrigObjs, muvec, elvec );
232 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 //----------------------------------------------------------------------------------------
258 void setEra(string fname, ControlFlags &ctrl)
259 {
260 //----------------------------------------------------------------------------------------
261 if( ctrl.mc && (fname.find("f11-") != string::npos) ) ctrl.era=2011 ;
262 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 }
267
268 //----------------------------------------------------------------------------------------
269 unsigned getNPU(mithep::Array<mithep::PileupInfo> *puArr, int bx)
270 //----------------------------------------------------------------------------------------
271 {
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 return 0;
288 }
289
290
291 //----------------------------------------------------------------------------
292 void initPUWeights()
293 //----------------------------------------------------------------------------
294 {
295 // TFile * puf;
296
297 // 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 }
307
308 //----------------------------------------------------------------------------
309 double getPUWeight(int era, TString fname, unsigned npu)
310 //----------------------------------------------------------------------------
311 {
312 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 }
327
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 }
338 //----------------------------------------------------------------------------------------
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
381 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 else if(ofname.Contains(TRegexp("-x12[0-9]zz4l-"))) return -1; // spin samples
413 else if(ofname.Contains(TRegexp("-h[0-9][0-9]*zz4l-"))) return ESampleType::kHZZ; // powheg HZZ
414 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 else assert(0);
419 }
420 //----------------------------------------------------------------------------------------
421 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 //----------------------------------------------------------------------------------------
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 //----------------------------------------------------------------------------------------