165 {
166
167
168 if (
config.verbose) cout <<
"maxnbaseline = " <<
config.maxnbaseline << endl;
169
171
172 int nevents_thiscase;
173
174
175 TTree *ntuple = new TTree("FakeBkg", "Variables from toy MC");
176 ntuple->Branch("nevents", &nevents_thiscase, "nevents/I");
177 Int_t nevents_sel;
178 ntuple->Branch("nevents_sel", &nevents_sel, "nevents_sel/I");
179 Float_t fake_lep_frac;
180 ntuple->Branch("fake_lep_frac", &fake_lep_frac, "fake_lep_frac/F");
181 Float_t true_fakes;
182 ntuple->Branch("true_fakes", &true_fakes,"true_fakes/F");
183 Float_t asm_fakes;
184 ntuple->Branch("asm_fakes", &asm_fakes,"asm_fakes/F");
185 Float_t asm_err;
186 ntuple->Branch("asm_err", &asm_err,"asm_err/F");
187 Float_t fkf_fakes;
188 ntuple->Branch("fkf_fakes", &fkf_fakes,"fkf_fakes/F");
189 Float_t fkf_err;
190 ntuple->Branch("fkf_err", &fkf_err,"fkf_err/F");
191 Float_t lhoodMM_fakes;
192 ntuple->Branch("lhoodMM_fakes", &lhoodMM_fakes,"lhoodMM_fakes/F");
193 Float_t lhoodMM_poserr;
194 ntuple->Branch("lhoodMM_poserr", &lhoodMM_poserr,"lhoodMM_poserr/F");
195 Float_t lhoodMM_negerr;
196 ntuple->Branch("lhoodMM_negerr", &lhoodMM_negerr,"lhoodMM_negerr/F");
197 Float_t lhoodFF_fakes;
198 ntuple->Branch("lhoodFF_fakes", &lhoodFF_fakes,"lhoodFF_fakes/F");
199 Float_t lhoodFF_poserr;
200 ntuple->Branch("lhoodFF_poserr", &lhoodFF_poserr,"lhoodFF_poserr/F");
201 Float_t lhoodFF_negerr;
202 ntuple->Branch("lhoodFF_negerr", &lhoodFF_negerr,"lhoodFF_negerr/F");
203
204
205 map<CP::SystematicVariation,float> lhoodMM_weight_map, lhoodMM_poserr_map, lhoodMM_negerr_map;
206 map<CP::SystematicVariation,float> lhoodFF_weight_map, lhoodFF_poserr_map, lhoodFF_negerr_map;
207 map<CP::SystematicVariation,float> asm_weight_map, asm_err_map;
208 map<CP::SystematicVariation,float> fkf_weight_map, fkf_err_map;
209
210 TRandom3
rand(242868252);
211
212 Double_t realeff_spread =
config.eff_spread;
213 Double_t fakeeff_spread =
config.eff_spread;
214
215 std::vector<bool> isTight;
216 std::vector<std::string>
input = {
config.outputdirname+
"efficiencies_full" };
217
218 if(false)
219 {
220
221
222
223 }
224 else
225 {
226 input.back().append(
".root");
230 }
231
232
233 for (
unsigned icase(0); icase <
config.ncases; icase++) {
234 nevents_sel = 0;
235
236 TH1F* h_lep_pt_lhoodMM = 0;
237 TH1F* h_lep_eta_lhoodMM = 0;
238 TH2F* h_lep_pt_eta_lhoodMM = 0;
239 TH1F* h_lep_pt_lhoodFF = 0;
240 TH1F* h_lep_eta_lhoodFF = 0;
241 TH2F* h_lep_pt_eta_lhoodFF = 0;
242 TH1F* h_lep_pt_asm = 0;
243 TH1F* h_lep_eta_asm = 0;
244 TH2F* h_lep_pt_eta_asm = 0;
245 TH1F* h_lep_pt_fkf = 0;
246 TH1F* h_lep_eta_fkf = 0;
247 TH2F* h_lep_pt_eta_fkf = 0;
248
251 h_lep_pt_lhoodMM =
new TH1F(
hname.c_str(),
hname.c_str(), 10, 0, 100);
253 h_lep_eta_lhoodMM =
new TH1F(
hname.c_str(),
hname.c_str(), 10, -5, 5);
255 h_lep_pt_eta_lhoodMM =
new TH2F(
hname.c_str(),
hname.c_str(), 10, 0, 100, 10, -5, 5);
257 h_lep_pt_lhoodFF =
new TH1F(
hname.c_str(),
hname.c_str(), 10, 0, 100);
259 h_lep_eta_lhoodFF =
new TH1F(
hname.c_str(),
hname.c_str(), 10, -5, 5);
261 h_lep_pt_eta_lhoodFF =
new TH2F(
hname.c_str(),
hname.c_str(), 10, 0, 100, 10, -5, 5);
263 h_lep_pt_asm =
new TH1F(
hname.c_str(),
hname.c_str(), 10, 0, 100);
265 h_lep_eta_asm =
new TH1F(
hname.c_str(),
hname.c_str(), 10, -5, 5);
267 h_lep_pt_eta_asm =
new TH2F(
hname.c_str(),
hname.c_str(), 10, 0, 100, 10, -5, 5);
269 h_lep_pt_fkf =
new TH1F(
hname.c_str(),
hname.c_str(), 10, 0, 100);
271 h_lep_eta_fkf =
new TH1F(
hname.c_str(),
hname.c_str(), 10, -5, 5);
273 h_lep_pt_eta_fkf =
new TH2F(
hname.c_str(),
hname.c_str(), 10, 0, 100, 10, -5, 5);
274 }
275
276 float lep_pt = 0, lep_eta = 0;
277
279
280 cout << "Starting case " << icase << endl;
282
283 std::vector<CP::LhoodMM_tools*> lhmTool_sav;
285 std::vector<CP::LhoodMM_tools*> lhmTool_FF_sav;
287 std::vector<CP::AsymptMatrixTool*> asmTool_sav;
289 std::vector<CP::ApplyFakeFactor*> fkfTool_sav;
290
291 float asmYield(0);
292 float asmErr(0);
293 float fkfYield (0);
294 float fkfErr(0);
296 ANA_CHECK( lhmTool.register1DHistogram(h_lep_pt_lhoodMM, &lep_pt) );
297 ANA_CHECK( lhmTool.register1DHistogram(h_lep_eta_lhoodMM, &lep_eta) );
298 ANA_CHECK( lhmTool.register2DHistogram(h_lep_pt_eta_lhoodMM, &lep_pt, &lep_eta) );
299 ANA_CHECK( lhmTool_FF.register1DHistogram(h_lep_pt_lhoodFF, &lep_pt) );
300 ANA_CHECK( lhmTool_FF.register1DHistogram(h_lep_eta_lhoodFF, &lep_eta) );
301 ANA_CHECK( lhmTool_FF.register2DHistogram(h_lep_pt_eta_lhoodFF, &lep_pt, &lep_eta) );
302 ANA_CHECK( fkfTool.register1DHistogram(h_lep_pt_fkf, &lep_pt) );
303 ANA_CHECK( fkfTool.register1DHistogram(h_lep_eta_fkf, &lep_eta) );
304 ANA_CHECK( fkfTool.register2DHistogram(h_lep_pt_eta_fkf, &lep_pt, &lep_eta) );
305 ANA_CHECK( asmTool.register1DHistogram(h_lep_pt_asm, &lep_pt) );
306 ANA_CHECK( asmTool.register1DHistogram(h_lep_eta_asm, &lep_eta) );
307 ANA_CHECK( asmTool.register2DHistogram(h_lep_pt_eta_asm, &lep_pt, &lep_eta) );
308 }
309
313 ANA_CHECK( lhmTool_FF.setProperty(
"DoFakeFactorFit",
true) );
314
316
318 for (
int iSave = 0; iSave <=
nSave; iSave++) {
319 std:: string toolName = "LhoodMM_tools_save_" + std::to_string(icase) + "_" + std::to_string(iSave);
321 lhmTool_sav.push_back(lhmTool_is);
327 }
328 toolName = "LhoodMM_tools_FF_save_" + std::to_string(icase) + "_" + std::to_string(iSave);
330 lhmTool_FF_sav.push_back(lhmTool_FF_is);
337 }
338
339 toolName = "AsymptMatrixTool_save_" + std::to_string(icase) + "_" + std::to_string(iSave);
344 }
345 asmTool_sav.push_back(asmTool_is);
347
348 toolName = "ApplyFakeFactor_save_" + std::to_string(icase) + "_" + std::to_string(iSave);
350 fkfTool_sav.push_back(fkfTool_is);
356 }
357 }
358 }
359
360
361 Double_t realeff_mean_thiscase =
rand.Gaus(
config.realeff_mean, realeff_spread);
362 if (realeff_mean_thiscase > 0.99) realeff_mean_thiscase = 0.99;
363 Double_t fakeeff_mean_thiscase =
rand.Gaus(
config.fakeeff_mean, fakeeff_spread);
364 if (fakeeff_mean_thiscase < 0.01) fakeeff_mean_thiscase = 0.01;
365 if (fakeeff_mean_thiscase > realeff_mean_thiscase) fakeeff_mean_thiscase = realeff_mean_thiscase-0.01;
366
367 fake_lep_frac =
rand.Uniform();
368
374
375 std::vector<LhoodMMEvent> mmevts;
376
377 true_fakes = 0;
378
379 std::vector<double> realeff, fakeeff;
380
381 vector<FinalState*>
fs;
383 for (
unsigned ilep(0); ilep <=
config.maxnbaseline; ilep++) {
386 fs.push_back(this_fs);
387 }
388
389 if (
config.poisson_fluctuations) {
390 nevents_thiscase =
rand.Poisson(
config.nevents);
391 } else {
392 nevents_thiscase =
config.nevents;
393 }
394
396
397
398
399 for (int pass(0); pass<2; pass++) {
400 int savIndex(0);
401
402
403 int nEventsMultFactor = 1;
404 if (pass > 0) nEventsMultFactor = 10;
405 for (Long64_t ievt(0); ievt<nevents_thiscase*nEventsMultFactor; ievt++) {
406
407 float nlep_frac = 1./(
config.maxnbaseline-
config.minnbaseline+1);
408 Int_t nlep_select =
config.minnbaseline+
rand.Uniform()/nlep_frac;
409
410 float extraweight(1.0);
411
412 realeff.clear();
413 fakeeff.clear();
414
416 vector<FakeBkgTools::ParticleData> leptons_data;
417
418 vector<bool> lep_real;
419
420
421 vector<float> lep_pts;
422 for (int ilep(0); ilep < nlep_select; ilep++) {
423 lep_pts.push_back(100*
rand.Uniform());
424 }
425
426 sort(lep_pts.begin(), lep_pts.end(), std::greater<float>());
427
428 for (int ilep(0); ilep < nlep_select; ilep++) {
430
431 bool isReal(false);
432 if (
rand.Uniform() > fake_lep_frac) {
433 isReal = true;
434 } else {
435 isReal = false;
436 }
437
438 lep_real.push_back(isReal);
441
442 lep_pt = lep_pts[ilep];
443 lep_eta = -5. + 10*
rand.Uniform();
445
447 ANA_CHECK( lookupEfficiencies(*lepton, lepton_data) );
448
449 if (isReal) {
451 TightAcc(*lepton) = true;
452 } else {
453 TightAcc(*lepton) = false;
454 }
455 } else {
457 TightAcc(*lepton) = true;
458 } else {
459 TightAcc(*lepton) = false;
460 }
461 }
463
464 leptons_data.push_back(lepton_data);
465 }
466
467 if (pass == 0) {
468 ANA_CHECK( asmTool.addEvent(leptons, extraweight) );
469 ANA_CHECK( fkfTool.addEvent(leptons, extraweight) );
470 ANA_CHECK( lhmTool.addEvent(leptons, extraweight) );
471 ANA_CHECK( lhmTool_FF.addEvent(leptons, extraweight) );
473 ANA_CHECK( lhmTool_sav[savIndex]->addEvent(leptons, extraweight) );
474 ANA_CHECK( lhmTool_FF_sav[savIndex]->addEvent(leptons, extraweight) );
475 ANA_CHECK( asmTool_sav[savIndex]->addEvent(leptons, extraweight) );
476 ANA_CHECK( fkfTool_sav[savIndex]->addEvent(leptons, extraweight) );
477 }
478
479 }
480
481 bool all_real(true);
482 if (pass == 1) {
483 for (int ilep(0); ilep < nlep_select; ilep++) {
484 if (!lep_real[ilep]) {
485 all_real = false;
486 }
487 }
488 }
489
490
491 if (pass == 1 && all_real) {
492 ANA_CHECK( fkfTool.addEvent(leptons, -extraweight/nEventsMultFactor) );
493
494
495 ANA_CHECK( lhmTool_FF.addEvent(leptons, -extraweight/nEventsMultFactor) );
496
498 ANA_CHECK( fkfTool_sav[savIndex]->addEvent(leptons, -extraweight/nEventsMultFactor) );
499 ANA_CHECK( lhmTool_FF_sav[savIndex]->addEvent(leptons, -extraweight/nEventsMultFactor) );
500 }
501 }
502
503
504
505 vector<int> lep_indices;
506 std::bitset<64> tights, reals, charges;
507 for (int ilep(0); ilep < nlep_select; ilep++) {
508 lep_indices.push_back(ilep);
509 if (lep_real[ilep]) reals.set(ilep);
510 }
511
512 if (pass == 0) {
513 for (long comb(0); comb < (1<<nlep_select); ++comb) {
514 tights = std::bitset<64>(comb);
515 if (
fs[nlep_select]->accept_selection(tights, charges)) {
516 if (
fs[nlep_select]->accept_process(nlep_select, reals, tights)) {
517 true_fakes += extraweight*
comboProb(leptons_data, tights, reals);
518 }
519 }
520 }
521
522 tights.reset();
523 for (int ilep = 0; ilep < nlep_select; ilep++) {
524 if (TightAcc(*leptons[ilep])) tights.set(ilep);
525 }
526 if (
fs[nlep_select]->accept_selection(tights,charges) ) {
527 nevents_sel += extraweight;
528 }
529
530 float wgt(1.);
532 asmYield += wgt;
533 asmErr += wgt*wgt;
534
536 fkfYield += wgt;
537 fkfErr += wgt*wgt;
538 } else {
539
540 if (all_real) {
541 float wgt(1.0);
543 fkfYield -= wgt/nEventsMultFactor;
544 fkfErr += wgt*wgt/(nEventsMultFactor);
545 }
546 }
547
549 if (pass == 0) {
550 if (ievt > 0 && (((
nSave*ievt)%nevents_thiscase ==0) || ievt == nevents_thiscase -1) ) {
553 std::unique_ptr<TFile> saveFile(TFile::Open(
saveFileName.c_str(),
"RECREATE"));
554 ANA_CHECK( lhmTool_sav[savIndex]->saveProgress(saveFile->mkdir(
"fakes")) );
555 saveFile->Close();
556
558 saveFile = std::make_unique<TFile>(
saveFileName.c_str(),
"RECREATE");
559 ANA_CHECK( asmTool_sav[savIndex]->saveProgress(saveFile.get()) );
560 saveFile->Close();
561
562 savIndex++;
563 }
564 } else {
565 if (ievt > 0 && (((
nSave*ievt)%(nevents_thiscase*nEventsMultFactor) ==0) || ievt == nevents_thiscase*nEventsMultFactor -1) ) {
567
569 std::unique_ptr<TFile> saveFile(TFile::Open(
saveFileName.c_str(),
"RECREATE"));
570 ANA_CHECK( lhmTool_FF_sav[savIndex]->saveProgress(saveFile->mkdir(
"fakes")) );
571 saveFile->Close();
572
574 saveFile = std::make_unique<TFile>(
saveFileName.c_str(),
"RECREATE");
575 ANA_CHECK( fkfTool_sav[savIndex]->saveProgress(saveFile.get()) );
576 saveFile->Close();
577 savIndex++;
578
579 float nfakes_tmp(0), err_tmp(0);
580 ANA_CHECK( fkfTool.getTotalYield(nfakes_tmp, err_tmp, err_tmp) );;
581 }
582 }
583 }
584
586 if (*it !=
nullptr)
delete *
it;
587 }
588 leptons.clear();
589 }
590 }
591 asm_err = sqrt(asmErr);
592 asm_fakes = asmYield;
593 ANA_CHECK( asmTool.getTotalYield(asmYield, asmErr, asmErr) );
594
595
596 ANA_CHECK( lhmTool.getTotalYield(lhoodMM_fakes, lhoodMM_poserr, lhoodMM_negerr) );
597
598 ANA_CHECK( lhmTool_FF.getTotalYield(lhoodFF_fakes, lhoodFF_poserr, lhoodFF_negerr) );
599
600 ANA_CHECK( fkfTool.getTotalYield(fkfYield, fkfErr, fkfErr) );
601 fkf_err = fkfErr;
602 fkf_fakes = fkfYield;
603
604 if (
config.test_systematics) {
605 auto sysvars = asmTool.affectingSystematics();
606 std::string systBrName, systBrNameF;
607
608 for(auto& sysvar : sysvars){
609 if (
config.verbose) asmTool.getSystDescriptor().printUncertaintyDescription(sysvar);
610
611 float asm_syst_weight, asm_syst_err;
612 float fkf_syst_weight, fkf_syst_err;
613 float lhoodMM_syst_weight, lhoodMM_syst_poserr, lhoodMM_syst_negerr;
614 float lhoodFF_syst_weight, lhoodFF_syst_poserr, lhoodFF_syst_negerr;
615
616 if (icase == 0) {
619 ANA_CHECK(
setupSystBranchesAsym(
"lhoodMM", sysvar, lhoodMM_syst_weight, lhoodMM_syst_poserr, lhoodMM_syst_negerr, lhoodMM_weight_map, lhoodMM_poserr_map, lhoodMM_negerr_map, ntuple) );
620 ANA_CHECK(
setupSystBranchesAsym(
"lhoodFF", sysvar, lhoodFF_syst_weight, lhoodFF_syst_poserr, lhoodFF_syst_negerr, lhoodFF_weight_map, lhoodFF_poserr_map, lhoodFF_negerr_map, ntuple) );
621 }
622 ANA_CHECK( asmTool.applySystematicVariation({sysvar}) );
623 ANA_CHECK( asmTool.getTotalYield(asm_weight_map.find(sysvar)->second,
624 asm_err_map.find(sysvar)->second,
625 asm_err_map.find(sysvar)->second));
626
627 ANA_CHECK( fkfTool.applySystematicVariation({sysvar}) );;
628 ANA_CHECK( fkfTool.getTotalYield(fkf_weight_map.find(sysvar)->second,
629 fkf_err_map.find(sysvar)->second,
630 fkf_err_map.find(sysvar)->second) );
631 ANA_CHECK( lhmTool.applySystematicVariation({sysvar}) );;
632 ANA_CHECK( lhmTool.getTotalYield(lhoodMM_weight_map.find(sysvar)->second,
633 lhoodMM_poserr_map.find(sysvar)->second,
634 lhoodMM_negerr_map.find(sysvar)->second) );
635 ANA_CHECK( lhmTool_FF.applySystematicVariation({sysvar}) );;
636 ANA_CHECK( lhmTool_FF.getTotalYield(lhoodFF_weight_map.find(sysvar)->second,
637 lhoodFF_poserr_map.find(sysvar)->second,
638 lhoodFF_negerr_map.find(sysvar)->second) );
639 }
640 }
641
642
643 cout << "OUTPUT true_fakes = " << true_fakes << endl;
644 cout << "OUTPUT nfakes for lhoodMM = " << lhoodMM_fakes << " +" << lhoodMM_poserr << " -" << lhoodMM_negerr << endl;
645 cout << "OUTPUT nfakes for lhoodFF = " << lhoodFF_fakes << " +" << lhoodFF_poserr << " -" << lhoodFF_negerr << endl;
646 cout << "OUTPUT nfakes for asm = " << asm_fakes << " +- " << asm_err << endl;
647 cout << "OUTPUT nfakes for fkf = " << fkf_fakes << " +- " << fkf_err << endl;
648 cout << "OUTPUT true_fakes = " << true_fakes << endl;
649 } else {
650
651 ANA_CHECK( doMerge(input,
"lhm",
config, h_lep_pt_lhoodMM, lep_pt, h_lep_eta_lhoodMM, lep_eta, h_lep_pt_eta_lhoodMM, lhoodMM_fakes, lhoodMM_poserr, lhoodMM_negerr, icase) );
652 ANA_CHECK( doMerge(input,
"lhf",
config, h_lep_pt_lhoodFF, lep_pt, h_lep_eta_lhoodFF, lep_eta, h_lep_pt_eta_lhoodFF, lhoodFF_fakes, lhoodFF_poserr, lhoodFF_negerr, icase) );
653 ANA_CHECK( doMerge(input,
"asm",
config, h_lep_pt_asm, lep_pt, h_lep_eta_asm, lep_eta, h_lep_pt_eta_asm, asm_fakes, asm_err, asm_err, icase) );
654 ANA_CHECK( doMerge(input,
"fkf",
config, h_lep_pt_fkf, lep_pt, h_lep_eta_fkf, lep_eta, h_lep_pt_eta_fkf, fkf_fakes, fkf_err, fkf_err, icase) );
655 }
656
657 ntuple->Fill();
658
659 f_out->cd();
661 h_lep_pt_lhoodMM->Write();
662 h_lep_eta_lhoodMM->Write();
663 h_lep_pt_eta_lhoodMM->Write();
664 h_lep_pt_lhoodFF->Write();
665 h_lep_eta_lhoodFF->Write();
666 h_lep_pt_eta_lhoodFF->Write();
667 h_lep_pt_asm->Write();
668 h_lep_eta_asm->Write();
669 h_lep_pt_eta_asm->Write();
670 h_lep_pt_fkf->Write();
671 h_lep_eta_fkf->Write();
672 h_lep_pt_eta_fkf->Write();
673
674 delete h_lep_pt_lhoodMM;
675 delete h_lep_eta_lhoodMM;
676 delete h_lep_pt_eta_lhoodMM;
677 delete h_lep_pt_lhoodFF;
678 delete h_lep_eta_lhoodFF;
679 delete h_lep_pt_eta_lhoodFF;
680 delete h_lep_pt_asm;
681 delete h_lep_eta_asm;
682 delete h_lep_pt_eta_asm;
683 delete h_lep_pt_fkf;
684 delete h_lep_eta_fkf;
685 delete h_lep_pt_eta_fkf;
686
687
688 }
689 }
690
691 f_out->Write();
692 f_out->Close();
693
694 return StatusCode::SUCCESS;
695}
void setProperty(columnar::PythonToolHandle &self, const std::string &key, nb::object value)
DataModel_detail::iterator< DataVector > iterator
Helper class to provide type-safe access to aux data.
void makePrivateStore()
Create a new (empty) private store for this object.
void setPtEtaPhi(float pt, float eta, float phi)
set the 4-vec
void setCharge(float charge)
set the charge of the object
Class providing the definition of the 4-vector interface.
double comboProb(const vector< FakeBkgTools::ParticleData > &leptons_data, const std::bitset< 64 > &tights, const std::bitset< 64 > &reals)
std::unique_ptr< TFile > openRootFile(fbtTestToyMC_config &config)
StatusCode setupSystBranches(const char *baseName, CP::SystematicVariation sysvar, float &weight, float &weight_err, std::map< CP::SystematicVariation, float > &syst_map, std::map< CP::SystematicVariation, float > &syst_err_map, TTree *ntuple)
StatusCode setupSystBranchesAsym(const char *baseName, CP::SystematicVariation sysvar, float &weight, float &weight_poserr, float &weight_negerr, std::map< CP::SystematicVariation, float > &syst_map, std::map< CP::SystematicVariation, float > &syst_poserr_map, std::map< CP::SystematicVariation, float > &syst_negerr_map, TTree *ntuple)
StatusCode writeROOT(const string &name, int type, float realeff_mean, float fakeeff_mean, float eff_spread, float eff_delta_with_pt)
Double_t n_fake_lhood_tot_negerr
Double_t n_fake_lhood_tot_poserr
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Electron_v1 Electron
Definition of the current "egamma version".
DataVector< IParticle > IParticleContainer
Simple convenience declaration of IParticleContainer.
Efficiency real_efficiency
Efficiency fake_efficiency