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