Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:28 

0001 #include "PhysicsTools/Heppy/interface/BTagSF.h"
0002 #include "TMath.h"
0003 
0004 namespace heppy {
0005 
0006   BTagSF::BTagSF(int seed) { randm = new TRandom3(seed); }
0007 
0008   BTagSF::~BTagSF() { delete randm; }
0009 
0010   Bool_t BTagSF::isbtagged(Float_t pt,
0011                            Float_t eta,
0012                            Float_t csv,
0013                            Int_t jetflavor,
0014                            Bool_t isdata,
0015                            UInt_t btagsys,
0016                            UInt_t mistagsys,
0017                            Bool_t is2012) {
0018     randm->SetSeed((int)((eta + 5) * 100000));
0019 
0020     Bool_t btagged = kFALSE;
0021 
0022     if (isdata) {
0023       if (csv > 0.679)
0024         btagged = kTRUE;
0025       else
0026         btagged = kFALSE;
0027       return btagged;
0028     }
0029 
0030     Double_t SFb = 0.0;
0031     Double_t eff_b = 0.719;
0032 
0033     SFb = getSFb(pt, btagsys, is2012);
0034 
0035     Double_t promoteProb_btag = 0;  // ~probability to promote to tagged
0036     Double_t demoteProb_btag = 0;   // ~probability to demote from tagged
0037 
0038     if (SFb < 1) {
0039       demoteProb_btag = std::abs(1.0 - SFb);
0040     } else {
0041       promoteProb_btag = std::abs(SFb - 1.0) / ((SFb / eff_b) - 1.0);
0042     }
0043 
0044     if (std::abs(jetflavor) == 5) {  // real b-jet
0045       if (csv > 0.679) {             // if tagged
0046         if (demoteProb_btag > 0 && randm->Uniform() > demoteProb_btag)
0047           btagged = kTRUE;  // leave it tagged
0048         else
0049           btagged = kFALSE;  // demote it to untagged
0050       } else {
0051         if (promoteProb_btag > 0 && randm->Uniform() < promoteProb_btag)
0052           btagged = kTRUE;  // promote it to tagged
0053         else
0054           btagged = kFALSE;  // leave it untagged
0055       }
0056       return btagged;
0057     }
0058 
0059     // not a real b-jet, apply mistag
0060 
0061     Double_t SFl = 0, eff_l = 0;
0062 
0063     // c or light
0064     if (std::abs(jetflavor) == 4) {
0065       // SFc = SFb with twice the quoted uncertainty
0066       SFl = getSFc(pt, btagsys, is2012);
0067       eff_l =
0068           0.192 *
0069           SFl;  // eff_c in MC for CSVM = (-1.5734604211*x*x*x*x +  1.52798999269*x*x*x +  0.866697059943*x*x +  -1.66657942274*x +  0.780639301724), x = 0.679
0070     } else {
0071       SFl = getSFl(pt, eta, mistagsys, is2012);
0072       eff_l = getMistag(pt, eta);
0073     }
0074 
0075     Double_t promoteProb_mistag = 0;  // ~probability to promote to tagged
0076     Double_t demoteProb_mistag = 0;   // ~probability to demote from tagged
0077 
0078     if (SFl > 1) {
0079       promoteProb_mistag = std::abs(SFl - 1.0) / ((SFl / eff_l) - 1.0);
0080     } else {
0081       demoteProb_mistag = SFl;
0082     }
0083 
0084     if (csv > 0.679) {  // if tagged
0085       if (demoteProb_mistag > 0 && randm->Uniform() > demoteProb_mistag)
0086         btagged = kFALSE;  // demote it to untagged
0087       else
0088         btagged = kTRUE;  // leave it tagged
0089     } else {              // not tagged
0090       if (promoteProb_mistag > 0 && randm->Uniform() < promoteProb_mistag)
0091         btagged = kTRUE;  // promote it to tagged
0092       else
0093         btagged = kFALSE;  // leave it untagged
0094     }
0095 
0096     return btagged;
0097   }
0098 
0099   Double_t BTagSF::getSFb(Float_t pt, UInt_t btagsys, Bool_t is2012) {
0100     // pT dependent scale factors
0101     // Tagger: CSVM within 30 < pt < 670 GeV, abs(eta) < 2.4, x = pt
0102     // SFb = 0.6981*((1.+(0.414063*x))/(1.+(0.300155*x))); (2011)
0103     // SFb = (0.938887+(0.00017124*x))+(-2.76366e-07*(x*x)); (2012)
0104     // for pt > 670 (800) GeV: use the SFb value at 670 (800) GeV with twice the quoted uncertainty for 2011 (2012)
0105     // for pt < 30 (20) GeV: use the SFb value at 30 (20) GeV with a +/-0.12 absolute uncertainty (twice the quoted uncertainty) for 2011 (2012)
0106     // i.e SFb(pt<30) = SFb(pt=30) +/- 0.12, so the relative uncertainty is 0.12/SFb(pt=30) for 2011
0107 
0108     Float_t x = pt;
0109     if (!is2012 && pt >= 670.0)
0110       x = 669.9;
0111     if (!is2012 && pt < 30.0)
0112       x = 30.0;
0113     if (is2012 && pt >= 800.0)
0114       x = 799.9;
0115     if (is2012 && pt < 20.0)
0116       x = 20.0;
0117 
0118     Double_t SFb = 1.0;
0119     if (!is2012) {
0120       SFb = 0.6981 * ((1. + (0.414063 * x)) / (1. + (0.300155 * x)));
0121     } else {
0122       SFb = (0.938887 + (0.00017124 * x)) + (-2.76366e-07 * (x * x));
0123     }
0124     if (btagsys == kNo)
0125       return SFb;
0126 
0127     Double_t SFb_error_2011[] = {0.0295675,
0128                                  0.0295095,
0129                                  0.0210867,
0130                                  0.0219349,
0131                                  0.0227033,
0132                                  0.0204062,
0133                                  0.0185857,
0134                                  0.0256242,
0135                                  0.0383341,
0136                                  0.0409675,
0137                                  0.0420284,
0138                                  0.0541299,
0139                                  0.0578761,
0140                                  0.0655432};
0141     Float_t ptmin_2011[] = {30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500};
0142     Float_t ptmax_2011[] = {40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 670};
0143     Double_t SFb_error_2012[] = {0.0415707,
0144                                  0.0204209,
0145                                  0.0223227,
0146                                  0.0206655,
0147                                  0.0199325,
0148                                  0.0174121,
0149                                  0.0202332,
0150                                  0.0182446,
0151                                  0.0159777,
0152                                  0.0218531,
0153                                  0.0204688,
0154                                  0.0265191,
0155                                  0.0313175,
0156                                  0.0415417,
0157                                  0.0740446,
0158                                  0.0596716};
0159     Float_t ptmin_2012[] = {20, 30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 600};
0160     Float_t ptmax_2012[] = {30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 600, 800};
0161 
0162     Double_t SFb_error_x = 0.0;
0163 
0164     UInt_t nbins = is2012 ? 16 : 14;
0165 
0166     for (UInt_t ibin = 0; ibin < nbins; ibin++) {
0167       if (!is2012) {
0168         if (x >= ptmin_2011[ibin] && x < ptmax_2011[ibin])
0169           SFb_error_x = SFb_error_2011[ibin];
0170       } else {
0171         if (x >= ptmin_2012[ibin] && x < ptmax_2012[ibin])
0172           SFb_error_x = SFb_error_2012[ibin];
0173       }
0174     }
0175     if (!is2012) {
0176       if (pt >= 670.0)
0177         SFb_error_x = 2.0 * 0.0655432;
0178       if (pt < 30.0)
0179         SFb_error_x = 0.12;
0180     } else {
0181       if (pt >= 800.0)
0182         SFb_error_x = 2.0 * 0.0717567;
0183       if (pt < 20.0)
0184         SFb_error_x = 2.0 * 0.0554504;
0185     }
0186 
0187     Double_t scalefactor = SFb;
0188 
0189     if (btagsys == kDown)
0190       scalefactor = (SFb - SFb_error_x);
0191     if (btagsys == kUp)
0192       scalefactor = (SFb + SFb_error_x);
0193 
0194     return scalefactor;
0195   }
0196 
0197   Double_t BTagSF::getSFc(Float_t pt, UInt_t btagsys, Bool_t is2012) {
0198     // SFc = SFb with twice the quoted uncertainty
0199 
0200     Float_t x = pt;
0201     if (!is2012 && pt >= 670.0)
0202       x = 669.9;
0203     if (!is2012 && pt < 30.0)
0204       x = 30.0;
0205     if (is2012 && pt >= 800.0)
0206       x = 799.9;
0207     if (is2012 && pt < 20.0)
0208       x = 20.0;
0209 
0210     Double_t SFc = 1.0;
0211     if (!is2012) {
0212       SFc = 0.6981 * ((1. + (0.414063 * x)) / (1. + (0.300155 * x)));
0213     } else {
0214       SFc = (0.938887 + (0.00017124 * x)) + (-2.76366e-07 * (x * x));
0215     }
0216     if (btagsys == kNo)
0217       return SFc;
0218 
0219     Double_t SFb_error_2011[] = {0.0295675,
0220                                  0.0295095,
0221                                  0.0210867,
0222                                  0.0219349,
0223                                  0.0227033,
0224                                  0.0204062,
0225                                  0.0185857,
0226                                  0.0256242,
0227                                  0.0383341,
0228                                  0.0409675,
0229                                  0.0420284,
0230                                  0.0541299,
0231                                  0.0578761,
0232                                  0.0655432};
0233     Float_t ptmin_2011[] = {30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500};
0234     Float_t ptmax_2011[] = {40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 670};
0235     Double_t SFb_error_2012[] = {0.0415707,
0236                                  0.0204209,
0237                                  0.0223227,
0238                                  0.0206655,
0239                                  0.0199325,
0240                                  0.0174121,
0241                                  0.0202332,
0242                                  0.0182446,
0243                                  0.0159777,
0244                                  0.0218531,
0245                                  0.0204688,
0246                                  0.0265191,
0247                                  0.0313175,
0248                                  0.0415417,
0249                                  0.0740446,
0250                                  0.0596716};
0251     Float_t ptmin_2012[] = {20, 30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 600};
0252     Float_t ptmax_2012[] = {30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 600, 800};
0253 
0254     Double_t SFc_error_x = 0.0;
0255 
0256     UInt_t nbins = is2012 ? 16 : 14;
0257 
0258     for (UInt_t ibin = 0; ibin < nbins; ibin++) {
0259       if (!is2012) {
0260         if (x >= ptmin_2011[ibin] && x < ptmax_2011[ibin])
0261           SFc_error_x = 2.0 * SFb_error_2011[ibin];
0262       } else {
0263         if (x >= ptmin_2012[ibin] && x < ptmax_2012[ibin])
0264           SFc_error_x = 2.0 * SFb_error_2012[ibin];
0265       }
0266     }
0267     if (!is2012) {
0268       if (pt >= 670.0)
0269         SFc_error_x = 2.0 * 2.0 * 0.0655432;
0270       if (pt < 30.0)
0271         SFc_error_x = 2.0 * 0.12;
0272     } else {
0273       if (pt >= 800.0)
0274         SFc_error_x = 2.0 * 2.0 * 0.0717567;
0275       if (pt < 20.0)
0276         SFc_error_x = 2.0 * 2.0 * 0.0554504;
0277     }
0278 
0279     Double_t scalefactor = SFc;
0280 
0281     if (btagsys == kDown)
0282       scalefactor = (SFc - SFc_error_x);
0283     if (btagsys == kUp)
0284       scalefactor = (SFc + SFc_error_x);
0285 
0286     return scalefactor;
0287   }
0288 
0289   Double_t BTagSF::getSFl(Float_t pt, Float_t eta, UInt_t mistagsys, Bool_t is2012) {
0290     Float_t x = std::min(double(pt), is2012 ? 670.0 : 800.0);
0291 
0292     Double_t SFl = 0;
0293 
0294     if (!is2012) {
0295       if (std::abs(eta) >= 0.0 && std::abs(eta) < 0.8) {
0296         if (mistagsys == kNo)
0297           SFl = ((1.06182 + (0.000617034 * x)) + (-1.5732e-06 * (x * x))) + (3.02909e-10 * (x * (x * x)));
0298         else if (mistagsys == kDown)
0299           SFl = ((0.972455 + (7.51396e-06 * x)) + (4.91857e-07 * (x * x))) + (-1.47661e-09 * (x * (x * x)));
0300         else if (mistagsys == kUp)
0301           SFl = ((1.15116 + (0.00122657 * x)) + (-3.63826e-06 * (x * x))) + (2.08242e-09 * (x * (x * x)));
0302       } else if (std::abs(eta) >= 0.8 && std::abs(eta) < 1.6) {
0303         if (mistagsys == kNo)
0304           SFl = ((1.111 + (-9.64191e-06 * x)) + (1.80811e-07 * (x * x))) + (-5.44868e-10 * (x * (x * x)));
0305         else if (mistagsys == kDown)
0306           SFl = ((1.02055 + (-0.000378856 * x)) + (1.49029e-06 * (x * x))) + (-1.74966e-09 * (x * (x * x)));
0307         else if (mistagsys == kUp)
0308           SFl = ((1.20146 + (0.000359543 * x)) + (-1.12866e-06 * (x * x))) + (6.59918e-10 * (x * (x * x)));
0309       } else if (std::abs(eta) >= 1.6 && std::abs(eta) < 2.4) {
0310         if (mistagsys == kNo)
0311           SFl = ((1.08498 + (-0.000701422 * x)) + (3.43612e-06 * (x * x))) + (-4.11794e-09 * (x * (x * x)));
0312         else if (mistagsys == kDown)
0313           SFl = ((0.983476 + (-0.000607242 * x)) + (3.17997e-06 * (x * x))) + (-4.01242e-09 * (x * (x * x)));
0314         else if (mistagsys == kUp)
0315           SFl = ((1.18654 + (-0.000795808 * x)) + (3.69226e-06 * (x * x))) + (-4.22347e-09 * (x * (x * x)));
0316       }
0317     } else {
0318       if (std::abs(eta) >= 0.0 && std::abs(eta) < 0.8) {
0319         if (mistagsys == kNo)
0320           SFl = ((1.07541 + (0.00231827 * x)) + (-4.74249e-06 * (x * x))) + (2.70862e-09 * (x * (x * x)));
0321         else if (mistagsys == kDown)
0322           SFl = ((0.964527 + (0.00149055 * x)) + (-2.78338e-06 * (x * x))) + (1.51771e-09 * (x * (x * x)));
0323         else if (mistagsys == kUp)
0324           SFl = ((1.18638 + (0.00314148 * x)) + (-6.68993e-06 * (x * x))) + (3.89288e-09 * (x * (x * x)));
0325       } else if (std::abs(eta) >= 0.8 && std::abs(eta) < 1.6) {
0326         if (mistagsys == kNo)
0327           SFl = ((1.05613 + (0.00114031 * x)) + (-2.56066e-06 * (x * x))) + (1.67792e-09 * (x * (x * x)));
0328         else if (mistagsys == kDown)
0329           SFl = ((0.946051 + (0.000759584 * x)) + (-1.52491e-06 * (x * x))) + (9.65822e-10 * (x * (x * x)));
0330         else if (mistagsys == kUp)
0331           SFl = ((1.16624 + (0.00151884 * x)) + (-3.59041e-06 * (x * x))) + (2.38681e-09 * (x * (x * x)));
0332       } else if (std::abs(eta) >= 1.6 && std::abs(eta) < 2.4) {
0333         if (mistagsys == kNo)
0334           SFl = ((1.05625 + (0.000487231 * x)) + (-2.22792e-06 * (x * x))) + (1.70262e-09 * (x * (x * x)));
0335         else if (mistagsys == kDown)
0336           SFl = ((0.956736 + (0.000280197 * x)) + (-1.42739e-06 * (x * x))) + (1.0085e-09 * (x * (x * x)));
0337         else if (mistagsys == kUp)
0338           SFl = ((1.15575 + (0.000693344 * x)) + (-3.02661e-06 * (x * x))) + (2.39752e-09 * (x * (x * x)));
0339       }
0340     }
0341 
0342     return SFl;
0343   }
0344 
0345   Double_t BTagSF::getMistag(Float_t pt, Float_t eta) {
0346     Float_t x = std::min(double(pt), 670.0);
0347 
0348     Double_t eff_l = 0.0;
0349 
0350     if (std::abs(eta) >= 0.0 && std::abs(eta) < 0.8) {
0351       eff_l = ((0.00967751 + (2.54564e-05 * x)) + (-6.92256e-10 * (x * x)));
0352     } else if (std::abs(eta) >= 0.8 && std::abs(eta) < 1.6) {
0353       eff_l = ((0.00974141 + (5.09503e-05 * x)) + (2.0641e-08 * (x * x)));
0354     } else if (std::abs(eta) >= 1.6 && std::abs(eta) < 2.4) {
0355       eff_l = ((0.013595 + (0.000104538 * x)) + (-1.36087e-08 * (x * x)));
0356     }
0357 
0358     return eff_l;
0359   }
0360 
0361 }  // namespace heppy
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!