20 #include "TTreeReader.h"
21 #include "TTreeReaderValue.h"
23 #include "TDirectory.h"
27 #ifndef FAKEBKGTOOLS_ATLAS_ENVIRONMENT
28 #define declareProperty(n, p, h) ExtraPropertyManager<LhoodMM_tools>::declareProperty(n, &LhoodMM_tools::p, h)
36 using Clock=std::chrono::high_resolution_clock;
48 declareProperty(
"FixHistogramNormalization",
m_fixNormalization,
"Boolean to determine whether or not histograms are scaled such that their normalization is equal to the fake yield computed for the entire sample (true = yes, do the scaleing)");
95 for (
int ilep = 0; ilep <
s_nLepMax; ilep++) {
98 for (
int ientry = 0; ientry < (0x1 << (ilep+1)); ientry++) {
99 m_coeffs[ilep][ientry].resize((0
x1 << (ilep+1)) );
100 for (
int jentry = 0; jentry < (0x1 << (ilep+1)); jentry++) {
122 for (
int ilep = 0; ilep <
s_nLepMax; ilep++) {
134 for (
int ilep = 0; ilep <
s_nLepMax; ilep++) {
135 m_nrf_mat_vec[ilep] = std::shared_ptr<TMatrixT<double>>(std::make_shared<TMatrixT<double>>(0x1 << (ilep+1),1) );
136 m_MMmatrix_vec[ilep] = std::shared_ptr<TMatrixT<double>>(std::make_shared<TMatrixT<double>> ((0x1 << (ilep+1)),(0x1 << (ilep+1))));
137 m_ntlpred_vec[ilep] = std::shared_ptr<TMatrixT<double>>(std::make_shared< TMatrixT<double>>(0x1 << (ilep+1),1));
149 return StatusCode::SUCCESS;
150 }
else if (
ft ==
"MM") {
152 return StatusCode::SUCCESS;
154 ATH_MSG_ERROR(
"Error in LhoodMM_tools::setFitType: please specify \"MM\" for matrix method or \"FF\" for fake factor method");
155 return StatusCode::FAILURE;
163 if(
sc != StatusCode::SUCCESS)
return sc;
170 int ncells =
h1->GetNcells();
172 fitinfovec->resize(ncells);
181 ATH_MSG_ERROR(
"Error in LhoodMM_tools::register1DHistogram: histogram has already been registered");
182 return StatusCode::FAILURE;
184 return StatusCode::SUCCESS;
191 if(
sc != StatusCode::SUCCESS)
return sc;
197 int ncells = h2->GetNcells();
199 fitinfovec->resize(ncells);
204 ATH_MSG_INFO(
"Registered a 2D histogram "<<h2->GetName());
206 ATH_MSG_ERROR(
"Error in LhoodMM_tools::register2DHistogram: histogram has already been registered");
207 return StatusCode::FAILURE;
209 return StatusCode::SUCCESS;
216 if(
sc != StatusCode::SUCCESS)
return sc;
222 int ncells = h3->GetNcells();
224 fitinfovec->resize(ncells);
229 ATH_MSG_INFO(
"Registered a 3D histogram "<<h3->GetName());
231 ATH_MSG_ERROR(
"Error in LhoodMM_tools::register3DHistogram: histogram has already been registered");
232 return StatusCode::FAILURE;
234 return StatusCode::SUCCESS;
238 const std::vector<Efficiency>& realEff_vals,
239 const std::vector<Efficiency>& fakeEff_vals,
240 const std::vector<int>& charges,
245 int nlep = isTight_vals.size();
253 LhoodMMEvent mmevt(nlep, realEff_vals, fakeEff_vals, isTight_vals, charges, extraweight);
258 return StatusCode::SUCCESS;
266 ATH_MSG_WARNING(
"Attempt to add an event with 0 leptons. This event will be ignored.");
267 return StatusCode::SUCCESS;
269 std::vector<bool> isTight_vals;
270 std::vector<Efficiency> realEff_vals;
271 std::vector<Efficiency> fakeEff_vals;
272 std::vector<int> charges;
273 std::vector<FakeBkgTools::ParticleData>::const_iterator particles_it;
276 isTight_vals.push_back(
p.tight);
277 realEff_vals.push_back(
p.real_efficiency);
278 fakeEff_vals.push_back(
p.fake_efficiency);
279 double r_eff =
p.real_efficiency.value(
this);
280 double f_eff =
p.fake_efficiency.value(
this);
284 for(
const std::pair<short unsigned int, FakeBkgTools::Uncertainty> kv :
p.real_efficiency.uncertainties)
286 ATH_MSG_DEBUG(
"real eff uncertainties for first lepton are " << kv.second.up <<
" " << kv.second.down);
288 for(
const std::pair<short unsigned int, FakeBkgTools::Uncertainty> kv :
p.fake_efficiency.uncertainties)
290 ATH_MSG_DEBUG(
"fake eff uncertainties for first lepton are " << kv.second.up <<
" " << kv.second.down);
293 charges.push_back(
p.charge);
294 if ( r_eff < 0.01 && f_eff< 0.01) {
311 Double_t poserr, negerr;
314 yield =
nfakes(&poserr,&negerr);
316 ATH_MSG_DEBUG(
"Leaving getTotalYield with yield = " << yield);
319 statErrDown = -negerr;
326 int nlep = mmevt.
nlep();
338 const float*
val = map1_iter->second;
339 TH1*
h = map1_iter->first;
346 fitInfo = &histoMap->second[
icell];
352 fitInfo = &histoMap->second[
icell];
360 return StatusCode::FAILURE;
364 std::map<TH2*, std::pair<const float*, const float*> >
::iterator map2_iter;
366 std::pair<const float*, const float*>
val = map2_iter->second;
367 TH2*
h = map2_iter->first;
375 fitInfo = &histoMap->second[
icell];
380 fitInfo = &histoMap->second[
icell];
384 return StatusCode::FAILURE;
389 std::map<TH3*, std::tuple<const float*, const float*, const float*> >
::iterator map3_iter;
391 std::tuple<const float*, const float*, const float*>
val = map3_iter->second;
392 TH3*
h = map3_iter->first;
399 fitInfo = &histoMap->second[
icell];
403 icell =
h->FindBin(*(std::get<0>(
val)), *(std::get<1>(
val)), *(std::get<2>(
val)) );
404 fitInfo = &histoMap->second[
icell];
408 return StatusCode::FAILURE;
412 return StatusCode::SUCCESS;
418 int nlep = mmevt.
nlep();
420 int rank = 0x1 << nlep;
430 unsigned int catIndex = 0;
431 for (
int jlep = 0; jlep < nlep; jlep++) {
432 catIndex += (!mmevt.
isTight(jlep)) << jlep;
449 for (
int icomb = 0; icomb < (0x1 << nlep); icomb++) {
451 std::bitset<s_nLepMax+1> tights(icomb);
452 int ntight = tights.count();
454 for (
int jlep = 0; jlep < nlep; jlep++) {
456 totcharge += mmevt.
charge(jlep);
460 ATH_MSG_VERBOSE(
"Setting OSfrac_denom[" << lepidx <<
"][" << ntight <<
"]");
465 ATH_MSG_VERBOSE(
"Setting OSfrac_num[" << lepidx <<
"][" << ntight <<
"]");
466 if ((TMath::Abs(totcharge) < ntight) || ntight == 1) fitInfo.
OSfrac_num[lepidx][ntight]+=
weight;
470 std::vector<std::vector<FakeBkgTools::Efficiency>>
vals(2, std::vector<FakeBkgTools::Efficiency>(nlep));
471 for (
int ilep = 0; ilep < nlep; ilep++) {
480 vals[0][ilep].setToConst(1.0);
490 for (
int irf = 0; irf < rank; irf++) {
491 for (
int itl = 0; itl < rank; itl++) {
493 for (
int ilep = 0; ilep < nlep; ilep++) {
494 if (itl & (0
x1 << ilep) ) {
495 if (irf & (0
x1 << ilep)) {
503 if (irf & (0
x1 << ilep) ) {
510 ATH_MSG_VERBOSE(
"about to set m_coeffs_num[" << lepidx <<
"][" << itl <<
"][" << irf <<
"]");
512 fitInfo.
normterms[lepidx][(itl<<nlep) + irf].
add(curr_coeff_num);
513 fitInfo.
coeffs_num[lepidx][itl][irf].add(curr_coeff_num);
518 return StatusCode::SUCCESS;
521 return StatusCode::FAILURE;
531 for (
int jlep = 0; jlep < ilep; jlep++) {
540 return StatusCode::SUCCESS;
546 double f = -(obs*TMath::Log(pred)-pred);
547 if (obs > 0)
f += obs*TMath::Log(obs)-obs;
551 #ifdef FAKEBKGTOOLS_ATLAS_ENVIRONMENT
552 #define ASG_MSG_VERBOSE(x) do { if ( m_current_lhoodMM_tool->msgLvl(MSG::VERBOSE)) m_current_lhoodMM_tool->msg(MSG::VERBOSE) << x << endmsg; } while(0)
554 #define ASG_MSG_VERBOSE(x) do { if(m_current_lhoodMM_tool->msgLvl(MSG::VERBOSE)) std::cout << x << std::endl; } while (0)
562 const bool verbose = (
l->m_printLevel > 0);
566 for (
int ipar = 0; ipar < npar; ipar++) {
570 int minnlep =
l->m_minnlep;
571 int maxnlep_loose =
l->m_maxnlep_loose;
577 Double_t sinterm_tot = 1.;
579 int theta_nlep_index = 1 + maxnlep_loose -minnlep;
581 for (
int ilep = minnlep; ilep <= maxnlep_loose; ilep++) {
582 theta_nlep_index +=
l->m_real_indices[ilep-1].size();
583 ASG_MSG_VERBOSE(
"theta_nlep_index for ilep = " << ilep <<
" = " << theta_nlep_index);
591 for (
int ilep = minnlep; ilep <= maxnlep_loose; ilep++) {
592 if (
l->m_current_fitInfo->eventCount[ilep-1] == 0) {
593 ASG_MSG_VERBOSE(
"m_real_indices[" << ilep-1 <<
"].size() = " <<
l->m_real_indices[ilep-1].size());
594 real_index +=
l->m_real_indices[ilep-1].size();
595 for (
int ipar = 0; ipar < (
int)
l->m_fake_indices[ilep-1].size()-1; ipar++) {
600 l->m_curr_nlep = ilep;
601 Int_t npar_thisnlep = 0x1 << ilep;
603 int theta_tot_index =
l->m_theta_tot_start_index+ilep-minnlep;
606 if (
l->m_maxnlep_loose - minnlep > 0) {
609 costerm = TMath::Cos(
par[theta_tot_index]);
614 pars_thisnlep[0] =
rsq*sinterm_tot*costerm*sinterm_tot*costerm;
615 sinterm_tot *= TMath::Sin(
par[theta_tot_index]);
617 pars_thisnlep[0] =
rsq;
621 ASG_MSG_VERBOSE(
"m_real_indices[ilep-1].size() = " <<
l->m_real_indices[ilep-1].size());
622 for (
unsigned ipar = 0; ipar <
l->m_real_indices[ilep-1].size(); ipar++) {
623 pars_thisnlep[par_index] =
par[real_index+ipar];
628 for (
int ipar = 0; ipar < (
int)
l->m_fake_indices[ilep-1].size()-1; ipar++) {
630 pars_thisnlep[par_index+ipar] =
par[theta_nlep_index];
631 if(
verbose)
ASG_MSG_VERBOSE(
"f pars_thisnlep[" << par_index+ipar <<
"] = " << pars_thisnlep[par_index+ipar]);
634 fcn_nlep(npar_thisnlep, gin,
f, pars_thisnlep, iflag);
637 real_index +=
l->m_real_indices[ilep-1].size();
648 const bool verbose = (
l->m_printLevel > 0);
650 int nlep =
l->m_curr_nlep;
652 int rank = 0x1 << nlep;
656 for (
int ipar = 0; ipar < npar; ipar++) {
662 std::string
txt =
"testing variable transform: angle pars in fcn: ";
663 for (
int i = 0;
i < npar;
i++) {
671 std::shared_ptr<TMatrixT<double>> nrf, MMmatrix, ntlpred;
672 nrf =
l->m_nrf_mat_vec[lepidx];
677 MMmatrix =
l->m_MMmatrix_vec[lepidx];
678 ntlpred =
l->m_ntlpred_vec[lepidx];
680 (*nrf)(0,0) =
par[0];
681 double rsq = TMath::Abs(
par[0]);
683 int rsize =
l->m_real_indices[lepidx].size();
684 int fsize =
l->m_fake_indices[lepidx].size();
686 for (
int ipar = 0; ipar < rsize; ipar++) {
687 ASG_MSG_VERBOSE(
"In fcn, setting real par " <<
l->m_real_indices[lepidx][ipar] <<
" to " <<
par[ipar+1]);
688 (*nrf)(
l->m_real_indices[lepidx][ipar], 0) =
par[ipar+1];
694 for (
int ipar = 0; ipar < fsize; ipar++) {
696 if (ipar < fsize-1 ) {
697 double costerm = TMath::Cos(
par[rsize+ipar+1]);
698 if(
verbose)
ASG_MSG_VERBOSE(
"for setting fake parameter, sinterm = " << sinterm <<
" par index = " <<
l->m_real_indices[lepidx].size()+ipar+1);
699 (*nrf)(
l->m_fake_indices[lepidx][ipar],0) =
rsq*sinterm*costerm*sinterm*costerm;
701 (*nrf)(
l->m_fake_indices[lepidx][ipar],0) =
rsq*sinterm*sinterm;
705 ASG_MSG_VERBOSE(
"In fcn, setting fake par " <<
l->m_fake_indices[lepidx][ipar] <<
" to " << (*nrf)(
l->m_fake_indices[lepidx][ipar],0));
707 sinterm *= TMath::Sin(
par[rsize+ipar+1]);
715 *ntlpred = (*MMmatrix)*(*nrf);
724 if (
l->m_doFakeFactor) {
730 for (
int ipar = ipar_start; ipar < rank; ipar++) {
731 if(
verbose)
ASG_MSG_VERBOSE(
"Comparing parameter " << ipar <<
": " <<
l->m_current_fitInfo->event_cat.at(lepidx).at(ipar) <<
" to " << (*ntlpred)(ipar,0));
733 int nobs =
l->m_current_fitInfo->event_cat.at(lepidx).at(ipar);
736 s=
l->m_current_fitInfo->event_sumw2.at(lepidx).at(ipar)/nobs;
762 for (
int ilep = 0; ilep <
s_nLepMax; ilep++) {
766 if (
error.size() > 0) {
785 for (
int icomb = 0; icomb < (0x1 << (ilep+1)); icomb++) {
787 ATH_MSG_VERBOSE(
"ilep " << ilep <<
" (0x1 << ilep) " << std::hex << (0
x1 << ilep) <<
" icomb " << std::hex << icomb << std::dec);
791 if (
m_fsvec[ilep]->accept_selection(tights, charges)) {
793 ATH_MSG_VERBOSE(
"tights = " << std::hex << tights << std::dec <<
" nlep = " << nlep);
807 int maxNlep_proc = 0;
809 for(
unsigned c=0;
c<(1
u<<
n);++
c)
814 for (
int ibit = 0; ibit <
n; ibit++) {
815 reals.set(ibit, ~fakes[ibit]);
818 if (
m_fsvec[
n-1]->accept_process(
n, reals, tights) ) {
819 if(
n < minNlep_proc) minNlep_proc =
n;
820 if (
n > maxNlep_proc) {
832 if(
setup() != StatusCode::SUCCESS)
return 0.;
847 double nfake_fit, nfake_fitErr;
856 std::unique_ptr<TMinuit_LHMM> lhoodFit(
new TMinuit_LHMM(npar));
861 Double_t arglist[10];
865 lhoodFit->mnexcm(
"SET ERR", arglist ,1,ierflg);;
873 init_pars.resize(npar);
875 vector< vector <double> > loc_init_pars;
879 loc_init_pars[ilep-1].resize(0
x1 << ilep);
881 for (
int ipar = 0; ipar < (0x1 << ilep); ipar++) {
882 init_pars[
index+ipar] = loc_init_pars[ilep-1][ipar];
887 Double_t
step = TMath::Max(loc_init_pars[
m_minnlep-1][0]/1000,0.001);
893 vector<TString> parameterName;
896 parameterName.push_back(
"nfake_tot");
897 TString sreal =
"nreal";
900 sprintf(tmpchar,
"_%i", ilep);
901 TString tmpstr = sreal;
902 tmpstr.Append(tmpchar);
903 for (
unsigned isublep = 0; isublep <
m_real_indices[ilep-1].size(); isublep++) {
904 TString locsreal = tmpstr;
906 sprintf(tmpchar2,
"_%u", isublep);
907 locsreal.Append(tmpchar2);
908 parameterName.push_back(locsreal);
910 nreal_start_indices[ilep-1] = (parameterName.size());
911 ATH_MSG_VERBOSE(
"nreal_start_indices[" << ilep-1 <<
"] = " << nreal_start_indices[ilep-1]);
916 TString stheta_tot =
"theta_tot";
920 TString tmpstr = stheta_tot;
921 tmpstr.Append(tmpchar);
922 parameterName.push_back(tmpstr);
926 TString stheta =
"theta_";
928 sprintf(tmpchar,
"%i", ilep);
929 TString tmpstr = stheta;
930 tmpstr.Append(tmpchar);
934 TString locstheta = tmpstr;
936 sprintf(tmpchar2,
"%i", jlep);
937 locstheta.Append(tmpchar2);
938 parameterName.push_back(locstheta);
952 ATH_MSG_VERBOSE(
"nfakes for nlep = " << ilep <<
" used to find theta_tot = " << loc_init_pars[ilep-1][0]);
953 theta_tot[theta_index] = TMath::ACos(TMath::Sqrt(TMath::Max(loc_init_pars[ilep-1][0],0.)/(
m_nfakes_std))/sinterm);
954 if (TMath::IsNaN( theta_tot[theta_index] ) ) {
957 sinterm *= TMath::Sin(theta_tot[theta_index]);
970 upper_limits.resize(npar);
976 for (
unsigned isublep = 0; isublep <
m_real_indices[ilep-1].size(); isublep++) {
982 for (
int ipar = real_index; ipar < npar; ipar++) {
988 init_par_values.resize(npar);
994 for (
unsigned isublep = 0; isublep <
m_real_indices[ilep-1].size(); isublep++) {
995 ATH_MSG_VERBOSE(
"Setting parameter " << glob_index <<
" to " << init_pars[init_index+isublep]);
996 init_par_values[glob_index] = init_pars[init_index+isublep];
999 init_index+=
pow(2,ilep);
1010 theta_start_indices[
i] = theta_start_indices[
i-1] +
m_fake_indices[
i-1].size() - 1;
1028 index+= 0x1 << ilep;
1031 for (ipar = 0; ipar < npar; ipar++) {
1032 lhoodFit->mnparm(ipar, parameterName[ipar], init_par_values[ipar],
step, 0., upper_limits[ipar], ierflg);
1043 int nGoodLeptonMult = 0;
1049 for (
unsigned ipar = nreal_start_indices[ilep-1]; ipar < nreal_start_indices[ilep-1] +
m_real_indices[ilep-1].size(); ipar++) {
1052 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1053 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1058 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1059 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1062 for (
unsigned ipar = theta_start_indices[ilep-1]+1; ipar < theta_start_indices[ilep-1] +
m_fake_indices[ilep-1].size() ; ipar++) {
1065 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1066 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1071 index += (0x1 << ilep) - 2;
1074 if (nGoodLeptonMult == 0) {
1075 ATH_MSG_VERBOSE(
"No possible fake contribution for any lepton multiplicity");
1083 lhoodFit->mnexcm(
"MIGRAD", arglist ,2,ierflg);
1087 Double_t amin, edm, errdef;
1088 Int_t nvpar, nparx, icstat;
1089 lhoodFit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
1091 lhoodFit->GetParameter(0, nfake_fit, nfake_fitErr);
1095 if (poserr && negerr) {
1098 lhoodFit->mnexcm(
"MINOS", arglist ,2,ierflg);
1099 lhoodFit->mnerrs(0, *poserr, *negerr, nfake_fitErr, gcc);
1102 lhoodFit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
1103 lhoodFit->GetParameter(0, nfake_fit, nfake_fitErr);
1107 if ( *poserr < 1.
e-5) {
1108 *poserr =
fixPosErr(nfake_fit, lhoodFit.get());
1110 if (*negerr > -1.
e-5 ) {
1111 *negerr =
fixNegErr(nfake_fit, lhoodFit.get());
1114 if (*negerr > -1.
e-5) {
1115 *negerr = -nfake_fitErr;
1117 if (nfake_fit + *negerr < 0.) {
1118 *negerr = -nfake_fit;
1122 if (*poserr < 1.
e-5) {
1123 *poserr = nfake_fitErr;
1137 auto & fitInfo_vec = map1_iter->second;
1143 auto& fitInfo_vec = map2_iter->second;
1144 TH2*
histogram = (TH2*)map2_iter->first;
1150 auto& fitInfo_vec = map3_iter->second;
1173 int nbins =
h->GetNcells();
1176 Double_t
nf, poserr, negerr, shift = 0.;
1178 for (
int ibin = 0; ibin <
nbins; ibin++) {
1182 if (totEvents > 0) {
1183 ATH_MSG_VERBOSE(
"Filling bin " << ibin <<
" with " << totEvents <<
" events");
1185 h->SetBinContent(ibin,
nf+shift);
1186 if (TMath::IsNaN(
h->GetBinContent(ibin))) {
1187 h->SetBinContent(ibin,0.);
1188 h->SetBinError(ibin, 0.);
1190 h->SetBinError(ibin,TMath::Max(poserr,-negerr));
1193 h->SetBinContent(ibin,0.);
1194 h->SetBinError(ibin, 0.);
1196 ATH_MSG_VERBOSE(
"Result is " <<
h->GetBinContent(ibin) <<
" +/- " <<
h->GetBinError(ibin));
1201 double poserr, negerr;
1205 totFakes =
nfakes(&poserr, &negerr);
1208 double totHistContent = 0.;
1209 for (
int ibin = 0; ibin <
nbins; ibin++) {
1210 totHistContent +=
h->GetBinContent(ibin);
1217 if (totHistContent > 0.) {
1218 scaleFactor = totFakes/totHistContent;
1222 for (
int ibin = 1; ibin <=
nbins; ibin++) {
1223 h->SetBinContent(ibin, scaleFactor*
h->GetBinContent(ibin));
1224 h->SetBinError(ibin, scaleFactor*
h->GetBinError(ibin));
1227 return StatusCode::SUCCESS;
1247 int lepidx = nlep-1;
1250 nrf.resize(0
x1 <<nlep);
1255 std::string
txt =
"Testing variable transform: Initial nrf: ";
1256 for (
auto i = nrf.begin();
i != nrf.end(); ++
i)
1262 double nfakes_std_thisnlep = 0;
1265 for (
int ipar = 0; ipar < 0x1 <<nlep; ipar++) {
1269 for (
int ibit = 0; ibit < nlep; ibit++) {
1270 reals.set(ibit, ~fakes[ibit]);
1273 bool countsAsFake =
false;
1274 for (
int jpar = 0; jpar < 0x1 <<nlep; jpar++) {
1276 for (
int kpar = 0; kpar < 0x1 <<nlep; kpar++) {
1278 if (!countsAsFake &&
1279 m_fsvec[lepidx]->accept_process(nlep, reals, tights) &&
1280 m_fsvec[lepidx]->accept_selection(tights, charges) ) {
1282 nfakes_std_thisnlep += nrf[ipar];
1284 countsAsFake =
true;
1289 if (!countsAsFake) {
1296 init_pars[0] = nfakes_std_thisnlep;
1297 for (
unsigned ipar = 1; ipar <=
m_real_indices[lepidx].size(); ipar++) {
1301 if (nfakes_std_thisnlep > 0.5) {
1302 double sinterm = 1.;
1305 init_pars[ipar] = TMath::ACos(TMath::Sqrt(TMath::Max(nrf[
m_fake_indices[lepidx][ifake] ], 0.)/(nfakes_std_thisnlep))/sinterm);
1306 sinterm *= TMath::Sin(init_pars[ipar]);
1316 txt =
"testing variable transform: Initial pars: ";
1317 for (
int i = 0;
i < (0x1 << nlep);
i++) {
1324 for (
int ipar = 2; ipar < (0x1 << nlep); ipar++) {
1325 if (TMath::IsNaN(init_pars[ipar])) {
1326 init_pars[ipar] = 0.;
1328 ATH_MSG_VERBOSE(
"Setting angle parameter " << ipar <<
" to " << init_pars[ipar]);
1341 int lepidx = nlep-1;
1345 const int rank = 0x1 << nlep;
1347 std::vector<FakeBkgTools::Efficiency> coeff_denom(rank);
1352 for (
int irf = 0; irf < rank; irf++) {
1356 coeff_denom[irf].setToConst(0.);
1358 float chargefactor ;
1361 for (
int ibit = 0; ibit < nlep; ibit++) {
1362 reals.set(ibit, ~fakes[ibit]);
1364 for (
int itl = 0; itl < rank; itl++) {
1368 for (
int ibit = 0; ibit < nlep; ibit++) {
1369 tights.set(ibit, ~antitights[ibit]);
1374 if (
m_fsvec[lepidx]->accept_selection(tights, charges)
1375 &&
m_fsvec[lepidx]->accept_process(nlep, reals, tights) ) {
1380 if (nlep > 2 && tights.count() == 2) {
1386 if (nlep == 2 && tights.count() == 2) {
1392 chargefactor =
m_OSfrac[lepidx][tights.count()];
1394 ATH_MSG_VERBOSE(
"chargefactor = " << chargefactor <<
" for nlep = " << nlep);
1400 coeff_denom[irf].add(tmpEff);
1409 if (coeff_denom[irf].
nominal == 0.) {
1413 for (
int itl = 0; itl < rank; itl++) {
1420 std::shared_ptr<TMatrixT<double>> MMmatrix;
1423 for (
int i = 0;
i < rank;
i++) {
1424 for (
int j = 0; j < rank; j++) {
1432 if(verbose) MMmatrix->Print();
1434 TMatrixT<double> MMmatrix_inv(rank,rank);
1435 MMmatrix_inv = *MMmatrix;
1436 MMmatrix_inv.Invert();
1438 TMatrixT<double> MMmatrix_sqr = MMmatrix_inv;
1439 for (
int i = 0;
i < rank;
i++) {
1440 for (
int j = 0; j < rank; j++) {
1441 MMmatrix_sqr(
i,j) *= MMmatrix_sqr[
i][j];
1445 TMatrixT<double> nevents_mat(rank,1), nfake_mat(rank,1), nfake_err_mat(rank,1);
1446 for (
int i = 0;
i < rank;
i++) {
1451 if(
verbose) nevents_mat.Print();
1453 nfake_mat = MMmatrix_inv*nevents_mat;
1456 if(
verbose) MMmatrix->Print();
1458 if(
verbose) nevents_mat.Print();
1460 if(
verbose) nfake_mat.Print();
1463 nfake_err_mat = MMmatrix_sqr*nevents_mat;
1466 for (
int ipar = 0; ipar < (0x1 <<nlep) ; ipar++) {
1467 nrf[ipar] = nfake_mat(ipar, 0);
1473 for (
int ibit = 0; ibit < nlep; ibit++) {
1474 tights.set(ibit, 1);
1475 reals.set(ibit, ~fakes[ibit]);
1477 if (
m_fsvec[lepidx]->accept_process(nlep, reals, tights) &&
m_fsvec[lepidx]->accept_selection(tights, charges)) {
1478 ATH_MSG_VERBOSE(
"Adding " << nfake_mat(ipar,0) <<
" to m_nfakes_std");
1485 ATH_MSG_VERBOSE(
"Accepted " << n_proc_acc <<
" processes for nlep = " << nlep);
1495 Double_t f_from_fit, junk;
1499 double n_fake_guess_hi = TMath::Max(n_fake_fit*5,1.);
1500 double n_fake_guess_lo = n_fake_fit;
1501 double n_fake_guess = n_fake_guess_hi;
1502 double f_with_guess;
1505 bool stopSearch = 0;
1507 double convergeCriteria = 0.01;
1509 int nfake_tot_index = 1;
1511 Double_t arglist[10];
1515 arglist[0] = nfake_tot_index;
1516 arglist[1] = n_fake_fit;
1517 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1518 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1522 lhoodFit->mnexcm(
"MIGRAD", arglist ,2,ierflg);
1525 lhoodFit->mnstat(f_from_fit, junk, junk, ijunk, ijunk, ijunk);
1528 while (!stopSearch) {
1530 arglist[0] = nfake_tot_index;
1531 arglist[1] = n_fake_guess;
1532 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1533 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1537 lhoodFit->mnexcm(
"MIGRAD", arglist ,2,ierflg);
1540 lhoodFit->mnstat(f_with_guess, junk, junk, ijunk, ijunk, ijunk);
1542 ATH_MSG_VERBOSE(
"nlep, n_fake_guess, n_fake_guess_lo, n_fake_guesss_hi, f_from_fit, f_with_guess: " <<
"?" <<
" " << n_fake_guess <<
" " << n_fake_guess_lo <<
" " << n_fake_guess_hi <<
" " << f_from_fit <<
" " << f_with_guess);
1544 if (TMath::IsNaN(f_with_guess)) {
1545 f_with_guess = f_from_fit + 1.;
1547 if ((f_with_guess - f_from_fit) > 0.5) {
1548 n_fake_guess_hi = n_fake_guess;
1550 n_fake_guess_lo = n_fake_guess;
1553 n_fake_guess = 0.5*(n_fake_guess_lo+n_fake_guess_hi);
1555 ATH_MSG_VERBOSE(
"n_fake_guess_lo, hi = " << n_fake_guess_hi <<
" " << n_fake_guess_lo);
1556 if ((n_fake_guess_hi - n_fake_guess_lo)/n_fake_guess_hi < convergeCriteria) {
1558 ATH_MSG_VERBOSE(
"(n_fake_guess_lo, n_fake_fit = " << n_fake_guess_lo <<
" " << n_fake_fit);
1559 if (n_fake_guess_lo > n_fake_fit) {
1568 arglist[0] = nfake_tot_index;
1569 arglist[1] = n_fake_fit;
1570 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1573 return n_fake_guess - n_fake_fit;
1586 Double_t f_from_fit, junk;
1589 lhoodFit->mnstat(f_from_fit, junk, junk, ijunk, ijunk, ijunk);
1592 double n_fake_guess_hi = n_fake_fit;
1593 double n_fake_guess_lo = 0.;
1594 double n_fake_guess = n_fake_guess_lo;
1595 double f_with_guess;
1598 bool stopSearch = 0;
1599 double convergeCriteria = 0.01;
1600 double min_n_fake_guess = 0.05;
1602 int nfake_tot_index = 1;
1604 Double_t arglist[10];
1608 arglist[0] = nfake_tot_index;
1609 arglist[1] = n_fake_fit;
1610 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1611 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1615 lhoodFit->mnexcm(
"MIGRAD", arglist ,2,ierflg);
1618 lhoodFit->mnstat(f_from_fit, junk, junk, ijunk, ijunk, ijunk);
1620 while (!stopSearch) {
1624 arglist[0] = nfake_tot_index;
1625 arglist[1] = n_fake_guess;
1626 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1627 lhoodFit->mnexcm(
"FIX PAR", arglist, 1, ierflg);
1631 lhoodFit->mnexcm(
"MIGRAD", arglist ,2,ierflg);
1633 lhoodFit->mnstat(f_with_guess, junk, junk, ijunk, ijunk, ijunk);
1635 ATH_MSG_VERBOSE(
"nlep, n_fake_guess, n_fake_guess_lo, n_fake_guesss_hi, f_from_fit, f_with_guess: " <<
"?" <<
" " << n_fake_guess <<
" " << n_fake_guess_lo <<
" " << n_fake_guess_hi <<
" " << f_from_fit <<
" " << f_with_guess);
1637 if ((f_with_guess - f_from_fit) > 0.5) {
1638 n_fake_guess_lo = n_fake_guess;
1640 n_fake_guess_hi = n_fake_guess;
1642 n_fake_guess = 0.5*(n_fake_guess_lo+n_fake_guess_hi);
1644 if (((n_fake_guess_hi - n_fake_guess_lo)/n_fake_guess_hi < convergeCriteria) || (n_fake_guess_hi < min_n_fake_guess) ) {
1646 if (n_fake_guess_hi < n_fake_fit) {
1654 arglist[0] = nfake_tot_index;
1655 arglist[1] = n_fake_fit;
1656 lhoodFit->mnexcm(
"SET PAR", arglist, 2, ierflg);
1659 return n_fake_guess - n_fake_fit;
1672 ATH_MSG_ERROR(
"Multiple calls to saveProgress are not supported");
1673 return StatusCode::FAILURE;
1678 std::unique_ptr<TTree>
t(
new TTree(
"LhoodMM_progress",
"Stores current info from LhoodMM_toos"));
1680 std::unique_ptr<TTree> t_nlep(
new TTree(
"LhoodMM_progress_nlep",
"Stores minimum and maximum lepton multiplicities"));
1683 ATH_MSG_VERBOSE(
"Branch split level is " << fitInfoBranch->GetSplitLevel() );
1697 return StatusCode::SUCCESS;
1706 std::unique_ptr<TFile>
fin(
new TFile(
filename.c_str()));
1707 if (
fin ==
nullptr) {
1709 return StatusCode::FAILURE;
1714 return StatusCode::FAILURE;
1720 TTree *t_nlep = (TTree*)
fin->Get((
prefix +
"LhoodMM_progress_nlep").c_str());
1721 if (t_nlep ==
nullptr) {
1723 return StatusCode::FAILURE;
1726 int merged_maxnlep, merged_maxnlep_prev = -1;
1727 t_nlep->SetBranchAddress(
"maxnlep", &merged_maxnlep);
1729 TTree *
t = (TTree*)
fin->Get((
prefix +
"LhoodMM_progress").c_str());
1732 return StatusCode::FAILURE;
1737 Int_t
nentries = (Int_t)t_nlep->GetEntries();
1738 for (Int_t ievt = 0; ievt <
nentries; ievt++) {
1739 t_nlep->GetEntry(ievt);
1740 if (ievt > 0 && (merged_maxnlep != merged_maxnlep_prev)) {
1741 ATH_MSG_ERROR(
"Attempting to merge files with different lepton multiplicities. This is not supported.");
1742 return StatusCode::FAILURE;
1744 merged_maxnlep_prev = merged_maxnlep;
1755 t->SetBranchAddress(
"glb_fitInfo", &fitInfoPtr );
1761 std::unique_ptr<std::map<TH1*, std::vector< LhoodMMFitInfo > > > tmp_fitInfo_1dhisto_map(
new std::map<TH1*, std::vector< LhoodMMFitInfo > >);
1762 auto *tmp_fitInfo_1dhisto_map_ptr = tmp_fitInfo_1dhisto_map.get();
1763 t->SetBranchAddress(
"fitInfo_1dhisto_map", &tmp_fitInfo_1dhisto_map_ptr);
1764 std::unique_ptr<std::map<TH2*, std::vector< LhoodMMFitInfo > > > tmp_fitInfo_2dhisto_map(
new std::map<TH2*, std::vector< LhoodMMFitInfo > >);
1765 auto *tmp_fitInfo_2dhisto_map_ptr = tmp_fitInfo_2dhisto_map.get();
1766 t->SetBranchAddress(
"fitInfo_2dhisto_map", &tmp_fitInfo_2dhisto_map_ptr);
1767 std::map<TH3*, std::vector< LhoodMMFitInfo > > *tmp_fitInfo_3dhisto_map =
new std::map<TH3*, std::vector< LhoodMMFitInfo > >;
1768 t->SetBranchAddress(
"fitInfo_3dhisto_map", &tmp_fitInfo_3dhisto_map);
1771 for (Int_t ievt = 0; ievt <
nentries; ievt++) {
1778 for (
auto&
im: *tmp_fitInfo_1dhisto_map) {
1781 std::string iname = ihistogram->GetName();
1782 if (
hname == iname) {
1794 for (
auto&
im: *tmp_fitInfo_2dhisto_map) {
1797 std::string iname = ihistogram->GetName();
1798 if (
hname == iname) {
1810 for (
auto&
im: *tmp_fitInfo_3dhisto_map){
1813 std::string iname = ihistogram->GetName();
1814 if (
hname == iname) {
1834 return StatusCode::SUCCESS;