438 if (mPOLE<=0) mPOLE=125.6;
440 char TREE_NAME[] =
"TestTree";
443 if (mSample>0) cinput = Form(
"HZZ4lTree_jvbfMELA_H%.0f_%iTeV.root", mSample, erg_tev);
444 else cinput = Form(
"HZZ4lTree_jvbfMELA_HAll_%iTeV.root", erg_tev);
447 if (mSample>0) foutput =
new TFile(Form(
"jvbfMELA_Fits_H%.0f_%iTeV.root", mSample, erg_tev),
"recreate");
448 else foutput =
new TFile(Form(
"jvbfMELA_fits_wide_%iTeV.root", erg_tev),
"recreate");
450 TLorentzVector nullFourVector(0, 0, 0, 0);
452 float MC_weight_noxsec;
455 float pjvbf_VAJHU_first;
456 float pjvbf_VAJHU_second;
458 float jet1Pt, jet2Pt;
459 float jet1Eta, jet2Eta;
460 float jet1Phi, jet2Phi;
462 float jet1Pt_Fake, jet2Pt_Fake;
463 float jet1Eta_Fake, jet2Eta_Fake;
464 float jet1Phi_Fake, jet2Phi_Fake;
465 float jet1E_Fake, jet2E_Fake;
466 float jet1px, jet1py, jet1pz;
467 float jet2px, jet2py, jet2pz;
468 float ZZPx, ZZPy, ZZPz, ZZE, dR;
473 TChain* tree =
new TChain(TREE_NAME);
475 tree->SetBranchAddress(
"MC_weight_noxsec", &MC_weight_noxsec);
476 tree->SetBranchAddress(
"ZZMass", &ZZMass);
477 tree->SetBranchAddress(
"pjvbf_VAJHU", &pjvbf_VAJHU);
478 tree->SetBranchAddress(
"pjvbf_VAJHU_first", &pjvbf_VAJHU_first);
479 tree->SetBranchAddress(
"pjvbf_VAJHU_second", &pjvbf_VAJHU_second);
480 tree->SetBranchAddress(
"NJets30", &NJets30);
481 tree->SetBranchAddress(
"jet1Pt", &jet1Pt);
482 tree->SetBranchAddress(
"jet1Eta", &jet1Eta);
483 tree->SetBranchAddress(
"jet1Phi", &jet1Phi);
484 tree->SetBranchAddress(
"jet1E", &jet1E);
485 tree->SetBranchAddress(
"jet2Pt", &jet2Pt);
486 tree->SetBranchAddress(
"jet2Eta", &jet2Eta);
487 tree->SetBranchAddress(
"jet2Phi", &jet2Phi);
488 tree->SetBranchAddress(
"jet2E", &jet2E);
489 tree->SetBranchAddress(
"jet1_Fake_Pt", &jet1Pt_Fake);
490 tree->SetBranchAddress(
"jet1_Fake_Eta", &jet1Eta_Fake);
491 tree->SetBranchAddress(
"jet1_Fake_Phi", &jet1Phi_Fake);
492 tree->SetBranchAddress(
"jet1_Fake_E", &jet1E_Fake);
493 tree->SetBranchAddress(
"jet2_Fake_Pt", &jet2Pt_Fake);
494 tree->SetBranchAddress(
"jet2_Fake_Eta", &jet2Eta_Fake);
495 tree->SetBranchAddress(
"jet2_Fake_Phi", &jet2Phi_Fake);
496 tree->SetBranchAddress(
"jet2_Fake_E", &jet2E_Fake);
500 double bins_mzz[nbinsx+1]={70, 140, 190, 240, 340, 520, 20000};
501 double sum_w[nbinsx]={ 0 };
502 double sum_mzz[nbinsx][2]={ { 0 } };
503 const int nbins_eta = 8;
504 double bins_eta[nbins_eta+1]={ 0, 0.75, 1.5, 2.25, 3, 4, 5.5, 10.5, 12 };
506 TProfile* hProb[nbinsx];
508 TProfile* hEta[nbinsx+1];
509 TH1D* hDist[nbinsx+1];
510 for (
int bin=0; bin<nbinsx+1; bin++){
511 if (bin<nbinsx) hProb[bin] =
new TProfile(Form(
"njets1_pjvbf_bin_%i", bin+1), Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
512 else hProb_all =
new TH1D(
"njets1_pjvbf_all", Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
513 if (bin<nbinsx) hEta[bin] =
new TProfile(Form(
"njets1_eta_bin_%i", bin+1), Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
514 else hEta[bin] =
new TProfile(
"njets1_eta_all", Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
515 if (bin<nbinsx) hDist[bin] =
new TH1D(Form(
"njets1_dist_bin_%i", bin+1), Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
516 else hDist[bin] =
new TH1D(
"njets1_dist_all", Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
517 if (bin<nbinsx) hProb[bin]->Sumw2();
518 else hProb_all->Sumw2();
523 int nEntries = tree->GetEntries();
524 cout << nEntries << endl;
525 for (
int ev = 0;
ev < nEntries;
ev++){
529 for (
int bin=0; bin<nbinsx; bin++){
530 if (ZZMass>=bins_mzz[bin] && ZZMass<bins_mzz[bin+1]){
535 hProb[massbin]->Fill(fabs(jet2Eta_Fake), pjvbf_VAJHU_second, MC_weight_noxsec);
536 hEta[massbin]->Fill(fabs(jet2Eta_Fake), fabs(jet2Eta_Fake), MC_weight_noxsec);
537 hDist[massbin]->Fill(fabs(jet2Eta_Fake), MC_weight_noxsec);
539 hEta[nbinsx]->Fill(fabs(jet2Eta_Fake), fabs(jet2Eta_Fake), MC_weight_noxsec);
540 hDist[nbinsx]->Fill(fabs(jet2Eta_Fake), MC_weight_noxsec);
544 for (
int bin=0; bin<nbinsx; bin++){
545 hProb[bin]->Scale(1./hProb[bin]->Integral(
"width"));
547 for (
int bin=0; bin<nbinsx+1; bin++){
548 for (
int bineta=1; bineta<=nbins_eta; bineta++){
549 double bincontent = hDist[bin]->GetBinContent(bineta);
550 double binerror = hDist[bin]->GetBinError(bineta);
551 double binwidth = hDist[bin]->GetXaxis()->GetBinWidth(bineta);
553 bincontent /= binwidth;
554 binerror /= binwidth;
556 hDist[bin]->SetBinContent(bineta, bincontent);
557 hDist[bin]->SetBinError(bineta, binerror);
559 hDist[bin]->Scale(1./hDist[bin]->Integral(
"width"));
561 for (
int bineta=1; bineta<=nbins_eta; bineta++){
563 double sum_prob_invwsq=0;
564 for (
int bin=0; bin<nbinsx; bin++){
565 double bincontent = hProb[bin]->GetBinContent(bineta);
566 double binerror = hProb[bin]->GetBinError(bineta);
568 double invwsq = 0, prob_invwsq = 0;
570 invwsq = 1./pow(binerror, 2);
571 prob_invwsq = bincontent*invwsq;
573 sum_invwsq += invwsq;
574 sum_prob_invwsq += prob_invwsq;
578 sum_prob_invwsq /= sum_invwsq;
579 sum_invwsq = 1./sqrt(sum_invwsq);
581 hProb_all->SetBinContent(bineta, sum_prob_invwsq);
582 hProb_all->SetBinError(bineta, sum_invwsq);
584 hProb_all->Scale(1./hProb_all->Integral(
"width"));
587 double eta_center_all[nbins_eta];
588 double prob_center_all[nbins_eta];
589 double eta_center_all_err[nbins_eta];
590 double prob_center_all_err[nbins_eta];
591 int nNonEmpty_all = 0;
592 for (
int bineta=1; bineta<=nbins_eta; bineta++){
593 if (hProb_all->GetBinError(bineta)>0){
594 eta_center_all[nNonEmpty_all] = hEta[nbinsx]->GetBinContent(bineta);
595 eta_center_all_err[nNonEmpty_all] = hEta[nbinsx]->GetBinError(bineta);
596 prob_center_all[nNonEmpty_all] = hProb_all->GetBinContent(bineta);
597 prob_center_all_err[nNonEmpty_all] = hProb_all->GetBinError(bineta);
598 cout << nNonEmpty_all <<
'\t' << prob_center_all[nNonEmpty_all] <<
'\t' << prob_center_all_err[nNonEmpty_all] << endl;
603 TGraphErrors* tg_all =
new TGraphErrors(nNonEmpty_all, eta_center_all, prob_center_all, eta_center_all_err, prob_center_all_err);
604 tg_all->SetNameTitle(
"tg_njets1_pjvbf_all", Form(
"m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]));
605 tg_all->GetXaxis()->SetTitle(
"#eta");
606 tg_all->GetYaxis()->SetTitle(
"<P_{VBF}>");
608 TF1* fit_tg_all =
new TF1(
"fit_tg_njets1_pjvbf_all",
"gaus*(1+gaus(3))", bins_eta[0], bins_eta[nbins_eta]);
610 double mypars[npars]={ 0.155, 5.193, 1.46, 0.362, 3.524, 0.715 };
611 fit_tg_all->SetParameters(mypars);
612 fit_tg_all->SetParLimits(0, 0, 1);
613 fit_tg_all->SetParLimits(3, 0, 1);
614 fit_tg_all->SetParLimits(1, 0, 20);
615 fit_tg_all->SetParLimits(4, 0, 20);
616 fit_tg_all->SetParLimits(2, 0, 20);
617 fit_tg_all->SetParLimits(5, 0, 20);
618 tg_all->Fit(fit_tg_all);
620 double* par_postfit = fit_tg_all->GetParameters();
621 double* parerr_postfit = fit_tg_all->GetParErrors();
622 double parIndex[npars];
623 double parIndexErr[npars]={ 0 };
624 for (
int ip=0; ip<npars; ip++) parIndex[ip] = (
double)ip;
626 double integral = 0.5*par_postfit[0]*(1.+TMath::Erf(par_postfit[1]/par_postfit[2]/TMath::Sqrt(2)));
627 double integral2 = 1./(2.*TMath::Sqrt(2.*TMath::Pi()*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))));
628 integral2 *= par_postfit[0]*par_postfit[3]*TMath::Exp(-0.5*pow(par_postfit[1]-par_postfit[4], 2)/(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2)))*
629 (1.+TMath::Erf((TMath::Power(par_postfit[2], 2)*par_postfit[4]+TMath::Power(par_postfit[5], 2)*par_postfit[1]) / (par_postfit[2]*par_postfit[5]*TMath::Sqrt(2.*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))))));
630 integral += integral2;
632 cout << integral << endl;
633 par_postfit[0] /= integral;
634 parerr_postfit[0] /= integral;
636 TGraphErrors* tg_pars =
new TGraphErrors(npars, parIndex, par_postfit, parIndexErr, parerr_postfit);
637 tg_pars->SetNameTitle(
"njets1_pjvbf_pars",
"gaus*(1+gaus(3))");
638 tg_pars->GetXaxis()->SetTitle(
"Parameter index");
639 tg_pars->GetYaxis()->SetTitle(
"Parameter");
640 foutput->WriteTObject(tg_pars);
643 foutput->WriteTObject(tg_all);
646 for (
int bin=0; bin<nbinsx+1; bin++){
647 if (bin<nbinsx) foutput->WriteTObject(hProb[bin]);
648 else foutput->WriteTObject(hProb_all);
649 foutput->WriteTObject(hEta[bin]);
650 foutput->WriteTObject(hDist[bin]);
653 if (bin<nbinsx)
delete hProb[bin];
654 else delete hProb_all;