Skip to content
Snippets Groups Projects
Commit f1a2684b authored by Dag Gillberg's avatar Dag Gillberg
Browse files

Added very first implmenation of WG1 quark mass uncertaitnies + updated...

Added very first implmenation of WG1 quark mass uncertaitnies + updated pertubative STXS scheme uncertainties
parent 0d434cc4
No related branches found
No related tags found
9 merge requests!58791DataQualityConfigurations: Modify L1Calo config for web display,!46784MuonCondInterface: Enable thread-safety checking.,!46776Updated LArMonitoring config file for WD to match new files produced using MT,!45419WIP: GNNC tool,!45405updated ART test cron job,!42417Draft: DIRE and VINCIA Base Fragments for Pythia 8.3,!28528Revert 63f845ae,!27054Atr20369 210,!26342Monopole: Handle fractionally charged particles
......@@ -201,7 +201,7 @@ namespace xAOD {
updateWeights(hw.nominal,hw.ct14nlo,hw.ct14nlo_0118);
updateWeights(hw.nominal,hw.qcd_wg1_mu,hw.qcd_wg1_res);
updateWeights(hw.nominal,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12);
updateWeights(hw.nominal,hw.qcd_wg1_pTH,hw.qcd_wg1_qm);
updateWeights(hw.nominal,hw.qcd_wg1_pTH,hw.qcd_wg1_qm_b,hw.qcd_wg1_qm_t);
updateWeights(hw.nominal,hw.qcd_nnlops_nnlo,hw.qcd_nnlops_pow);
updateWeights(hw.nominal,hw.qcd_stxs);
}
......@@ -362,13 +362,23 @@ namespace xAOD {
// High pT uncertainty
static double y1_1 = 0.88, y2_1 = 1.16, x2_1 = 150;
static double y1_ge2 = 0.88, y2_ge2 = 1.16, x2_ge2 = 225;
double pTH_unc = 1.0;
if (Njets>=2) pTH_unc = pTH>x2_ge2?y2_ge2:y1_ge2+(y2_ge2-y1_ge2)*pTH/x2_ge2;
else if (Njets==1) pTH_unc = pTH>x2_1?y2_1:y1_1+(y2_1-y1_1)*pTH/x2_1;
hw.qcd_wg1_pTH = pTH_unc*w_nom;
double pTH_uncSF = 1.0;
if (Njets==1) pTH_uncSF = linInter(pTH,0,y1_1,x2_1,y2_1);
else if (Njets>=2) pTH_uncSF = linInter(pTH,0,y1_ge2,x2_ge2,y2_ge2);
//pTH>x2_ge2?y2_ge2:y1_ge2+(y2_ge2-y1_ge2)*pTH/x2_ge2;
//else if (Njets==1) pTH_uncSF = pTH>x2_1?y2_1:y1_1+(y2_1-y1_1)*pTH/x2_1;
hw.qcd_wg1_pTH = pTH_uncSF*w_nom;
double qmSF=1.0;
if (Njets==0) qmSF = linInter(pTH,6,0.92,24,1.05);
else if (Njets==1) qmSF = linInter(pTH,35,0.92,62.5,1.05);
else if (Njets>=2) qmSF = linInter(pTH,40,0.92,120,1.05);
hw.qcd_wg1_qm_b = w_nom * qmSF;
hw.qcd_wg1_qm_t = w_nom * linInter(pTH,160,1.0,500,1.37);
// Quark-mass uncdertainty - TODO
hw.qcd_wg1_qm = w_nom;
//if (pTH>500) qm=1.5; else if (pTH>160) qm=1.0+0.5*(pTH-160)/340;
//hw.qcd_wg1_qm = w_nom;
/*********
* Tackmann STXS uncertainty scheme
......@@ -383,12 +393,13 @@ namespace xAOD {
// This is followed by Dsig60, Dsig120 and Dsig200
// These are extracted from Powheg NNLOPS scale variations (using this tool!)
// As the envelope of hw.qcd_nnlops, which gives (all x-sec below have Njets>=1):
// sig60 = 9.687 +/- 1.566 pb
// sig120 = 2.534 +/- 0.526 pb
// sig200 = 0.562 +/- 0.115 pb
static double sig0_60=8.458, sig60_200=9.125, sig120_200=1.972,
sig0_120=sig0_60+sig60_200-sig120_200, sig200_plus=0.562;
static double Dsig60_200=1.457, Dsig120_200=0.411, Dsig200_plus=0.115;
// sig(60,200) = 9.095 +/- 1.445 pb, BLPTW 10.9%
// sig(120,200) = 1.961 +/- 0.401 pb, BLPTW 13.1%
// sig(200,inf) = 0.582 +/- 0.121 pb, BLPTW 15.1%
static double sig0_60=8.719, sig60_200=9.095, sig120_200=1.961,
sig0_120=sig0_60+sig60_200-sig120_200, sig200_plus=0.582; // 0.121 (-) 0.151*0.582
// static double Dsig60_200=1.457, Dsig120_200=0.411, Dsig200_plus=0.115; // with 40k events, and no BLPTW subraction
static double Dsig60_200=1.055, Dsig120_200=0.206, Dsig200_plus=0.0832; // with 2M evts, and subtraction
double dsig60=0, dsig120=0, dsig200=0;
if (Njets>=1) {
if (pTH<60) dsig60=-Dsig60_200/sig0_60; // -17.2%
......@@ -397,7 +408,7 @@ namespace xAOD {
if (pTH<120) dsig120 = -Dsig120_200/sig0_120; // -2.6%
else if (pTH<200) dsig120 = Dsig120_200/sig120_200; // +20.8%
if (pTH>200) dsig200=Dsig200_plus/sig200_plus; // +20.5%
if (pTH>200) dsig200=Dsig200_plus/sig200_plus; // +14.3%
}
hw.qcd_stxs.push_back((1.0+dsig60)*w_nom);
hw.qcd_stxs.push_back((1.0+dsig120)*w_nom);
......
......@@ -35,7 +35,7 @@ namespace xAOD {
/// WG1 proposed QCD uncertainty scheme
double qcd_wg1_mu, qcd_wg1_res, qcd_wg1_mig01, qcd_wg1_mig12;
double qcd_wg1_pTH, qcd_wg1_qm;
double qcd_wg1_pTH, qcd_wg1_qm_b, qcd_wg1_qm_t;
/// Tackmann proposed QCD uncertainty scheme, TODO
std::vector<double> qcd_stxs;
......@@ -51,6 +51,8 @@ namespace xAOD {
double pTH;
int Njets30, STXS;
std::vector<double> qcd_wg1() { return {qcd_wg1_mu,qcd_wg1_res,qcd_wg1_mig01,qcd_wg1_mig12,qcd_wg1_pTH,qcd_wg1_qm_b,qcd_wg1_qm_t}; }
/// methods to print weights to the screen
char *uncStr(double var, double nom) { return var==0?Form(" N/A"):Form("%s%.1f%%",var>=nom?"+":"",(var-nom)/nom*100); }
void print() {
......@@ -84,8 +86,8 @@ namespace xAOD {
printf("\n WG1 proposed QCD uncertainty scheme\n");
printf(" mu: %s, res: %s, mig01: %s, mig12: %s\n",
uncStr(qcd_wg1_mu,n),uncStr(qcd_wg1_res,n),uncStr(qcd_wg1_mig01,n),uncStr(qcd_wg1_mig12,n));
printf(" pTH: %s, quark-mass: %s\n",
uncStr(qcd_wg1_pTH,n),uncStr(qcd_wg1_qm,n));
printf(" pTH: %s, quark-mass, b: %s, t: %s\n",
uncStr(qcd_wg1_pTH,n),uncStr(qcd_wg1_qm_b,n),uncStr(qcd_wg1_qm_b,n));
}
if (qcd.size()) {
......@@ -151,6 +153,13 @@ namespace xAOD {
private:
/// linear interpolation
double linInter(double x, double x1, double y1, double x2, double y2) {
if (x<x1) return y1;
if (x>x2) return y2;
return y1+(y2-y1)*(x-x1)/(x2-x1);
}
/// Access the HiggsWeights
HiggsWeights getHiggsWeightsInternal(int HTXS_Njets30=-1, double HTXS_pTH=-99.0, int HTXS_cat=-1);
......
......@@ -85,7 +85,7 @@ int main( int argc, char* argv[] ) {
HistV h_pTH_qcd = makeHistos(8,"pTH_qcd",Nbins,min,max,ptTit); // Default Powheg QCD variations (NLO)
HistV h_pTH_nnlops_qcd = makeHistos(26,"pTH_nnlops_qcd",Nbins,min,max,ptTit); // NNLOPS internal QCD vars
HistV h_pTH_nnlops_qcd2 = makeHistos(2,"pTH_nnlo_qcd",Nbins,min,max,ptTit); // NNLO and Powheg vars for NNLOPS
HistV h_pTH_wg1qcd = makeHistos(6,"pTH_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
HistV h_pTH_wg1qcd = makeHistos(7,"pTH_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
Nbins=10; min=-0.5; max=9.5; Str tit=";#it{N}_{jets}";
TH1F *h_Njets = new TH1F("Njets30",tit,Nbins,min,max);
......@@ -96,7 +96,7 @@ int main( int argc, char* argv[] ) {
HistV h_Njets_qcd = makeHistos(8,"Njets30_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_Njets_nnlops_qcd = makeHistos(26,"Njets30_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_Njets_nnlops_qcd2 = makeHistos(2,"Njets30_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
HistV h_Njets_wg1qcd = makeHistos(6,"Njets30_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
HistV h_Njets_wg1qcd = makeHistos(7,"Njets30_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
Nbins=52; min=1; max=53; tit=";STXS fine index";
TH1F *h_STXS = new TH1F("STXS",tit,Nbins,min,max);
......@@ -107,7 +107,7 @@ int main( int argc, char* argv[] ) {
HistV h_STXS_qcd = makeHistos(8,"STXS_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_STXS_nnlops_qcd = makeHistos(26,"STXS_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_STXS_nnlops_qcd2 = makeHistos(2,"STXS_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
HistV h_STXS_wg1qcd = makeHistos(6,"STXS_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
HistV h_STXS_wg1qcd = makeHistos(7,"STXS_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
Nbins=80; min=-4; max=4; tit=";#it{y_{H}}";
TH1F *h_yH = new TH1F("yH",tit,Nbins,min,max);
......@@ -118,7 +118,7 @@ int main( int argc, char* argv[] ) {
HistV h_yH_qcd = makeHistos(8,"yH_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_yH_nnlops_qcd = makeHistos(26,"yH_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_yH_nnlops_qcd2 = makeHistos(2,"yH_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
HistV h_yH_wg1qcd = makeHistos(6,"yH_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
HistV h_yH_wg1qcd = makeHistos(7,"yH_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
HistV h_yH_other = makeHistos(8,"yH_other",Nbins,min,max,tit); // WG1 propsed scheme for ggF
// Initialise the application:
......@@ -258,7 +258,8 @@ int main( int argc, char* argv[] ) {
fillHistos(h_STXS_nnlops_qcd2,HTXS_index,nnlops_2np_qcd);
// 3.d WG1 propsed uncertainty scheme. Recommended.
NumV wg1_qcd={hw.qcd_wg1_mu,hw.qcd_wg1_res,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12,hw.qcd_wg1_pTH,hw.qcd_wg1_qm};
//NumV wg1_qcd={hw.qcd_wg1_mu,hw.qcd_wg1_res,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12,hw.qcd_wg1_pTH,hw.qcd_wg1_qm};
NumV wg1_qcd=hw.qcd_wg1();
fillHistos(h_pTH_wg1qcd,h.Pt(),wg1_qcd);
fillHistos(h_yH_wg1qcd,h.Rapidity(),wg1_qcd);
fillHistos(h_Njets_wg1qcd,HTXS_Njets30,wg1_qcd);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment