// Recalculate the flavour of the jet // Pass the jet eta, phi and ORIGINAL jet flavour // Returns the new jet flavour, 5 = b, 4 = c, 15 = tau, 0 = light int GetJetFlavour(float jet_eta, float jet_phi, int jet_flav) { float dRmin = 0.4; // Use a cut of 5 GeV for LHCP approval float ptmin = 5000; // float ptmin = 0; if(GetDeltaRNearest(jet_eta, jet_phi, 5, ptmin) < dRmin) return 5; else if(GetDeltaRNearest(jet_eta, jet_phi, 4, ptmin) < dRmin) return 4; else if(jet_flav == 15) return 15; else return 0; } // Check the closest hadron to a jet // Pass 4 to check for c-hardons and 5 to check for b-hadrons float GetDeltaRNearest(float eta, float phi, int flav, float minpt ) { float dR(99.); for (int i = 0; i < mchfpart_n; ++i) { if(flav == 4) { if(!IsCHadron(mchfpart_type->at(i))) continue; } if(flav == 5) { if(!IsBHadron(mchfpart_type->at(i))) continue; } if (mchfpart_pt->at(i) < minpt) continue; float etaTruth = mchfpart_eta->at(i); float phiTruth = mchfpart_phi->at(i); float deta = eta-etaTruth; float dphi = acos(cos(phi-phiTruth)); float tmpdR = sqrt( deta*deta+dphi*dphi); if(tmpdR<dR) { dR=tmpdR; } } return dR; } // Function to check whether particle is a c-hadron bool IsCHadron(int pdg) { pdg=abs(pdg); if((pdg>=400&&pdg<500&&pdg !=443)|| (pdg>=4000 && pdg<5000) || (pdg>=10411 && pdg<=10445) || (pdg>=20411 && pdg<=20445) ) return true; else return false; } // Function to check whether particle is a b-hadron bool IsBHadron(int pdg) { pdg=abs(pdg); if((pdg>=511 && pdg<=545) || (pdg>=10511 && pdg<=10545) || (pdg>=20511 && pdg<=20545) || (pdg>=5112 && pdg<=5554)) { return true; } else return false; }When looping over your chosen D3PD jet container recompute the jet flavour by calling the code as shown in the simplified example below: std::vector<int> jet_flavour_hadron_label;
for(int index = 0; index < jet_n; ++index) { // Recompute jet flavour int newJetFlav = GetJetFlavour(jet_eta->at(index),jet_phi->at(index),jet_flavor_truth_label->at(index)); jet_flavour_hadron_label.push_back(newJetFlav); }
*mv1cEval(double w_IP3D, double w_SV1, double w_pu, double w_pc, double w_pb, double jet_pt, double jet_eta)*
Note: The jet_pt and jet_eta we put in the MV1c script is using the re-calibrated jet pt and eta (Moriond2013 JES calibrated jets).
The usage of CalibrationDataInterface (CDI)
The CDI version we should use is CalibrationDataInterface-00-03-06:
atlasoff/PhysicsAnalysis/JetTagging/JetTagPerformanceCalibration/CalibrationDataInterface/tags/CalibrationDataInterface-00-03-06* (New) The CDI input file is:
In order to consider the b-tagging efficiency sample dependence we found in EPS analysis, this CDI file contains the efficiency maps of different samples, e.g. POWHEG ttbar, Sherpa W+jets, and Sherpa Z+jets with VH(bb) basic selections. In such a case, user needs to specify the name of the 3D efficiency maps in both .env file(below) and the place to get the MC efficiency. An example of the .env file (BTagCalibrationMV1c.2012.continuouseff.env) is shown here:/afs/cern.ch/user/g/giacinto/public/forHSG5/Summer2013_continuous_MV1c_customEffMaps_v2.root
File: WeightFiles_btag/Summer2013_continuous_MV1c_customEffMaps_v2.root FileEff: WeightFiles_btag/Summer2013_continuous_MV1c_customEffMaps_v2.root FileSF: WeightFiles_btag/Summer2013_continuous_MV1c_customEffMaps_v2.root MV1c.EfficiencyCalibrationBName: ttbar_powheg;wjets_sherpa;zjets_sherpa MV1c.EfficiencyCalibrationCName: ttbar_powheg;wjets_sherpa;zjets_sherpa MV1c.EfficiencyCalibrationLightName: ttbar_powheg;wjets_sherpa;zjets_sherpa MV1c.EfficiencyCalibrationTName: ttbar_powheg;wjets_sherpa;zjets_sherpa runEigenVectorMethod: True
//m_nTagBin is the number of b-tagging weight bins
// in the continuous case. MeanBWeightForBin(iwb) is the bin center of the bin-iwb.enum{m_nTagBin = 6}; double m_Xbin[m_nTagBin +1];
m_Xbin[0] = 0.0; m_Xbin[1] = 0.4050; m_Xbin[2] = 0.7028; m_Xbin[3] = 0.8353; m_Xbin[4] = 0.9237; m_Xbin[5] = 0.9848; m_Xbin[6] = 1.0;
float ana_MMN::MeanBWeightForBin( int i ) { //mean b-weight of the bin if( i< m_nXbin ) return ( m_Xbin[i]+m_Xbin[i+1] )/2; else { std::cout<<" HbbEvent::MeanBWeightForBin bin higher than number of bins " << i << " " << m_nXbin<<std::endl; return 0; } }
//** here you need to tell the code which sample you are running.
//** TSample is the index of the sample, i.e. ( ttbar_powheg=1;wjets_sherpa=1;zjets_sherpa=2),
//** The below is just an example, m_HistDirName is a input arguement that you pass to the code.
int tmp_TSample = 0; if (m_HistDirName.Contains("Wl") || m_HistDirName.Contains("Wc") || m_HistDirName.Contains("Wb")) { tmp_TSample = 1; } else if (m_HistDirName.Contains("Zl") || m_HistDirName.Contains("Zc") || m_HistDirName.Contains("Zb")) { tmp_TSample = 2; }
//** In order to generate the truth tagging weight of the jet, we need to store
//** the cumulative tagged efficiency( mcTTeffweight ) and anti-tagged efficiency(mcATTeffweight) ,
std::vector<double> mctteff(m_nTagBin, 0);
for(int iwb=0; iwb<m_nTagBin; iwb++) { float weightVal = MeanBWeightForBin(iwb); ajetTTEff.jetTagWeight=weightVal; mceff = m_calibMV1->getMCEfficiency(ajetTTEff, flav, "continuous", Analysis::None, TSample ); mctteff[iwb]=mceff.first; if(weightVal> m_loose_tag_cut ) mcTTeffweight +=mctteff[iwb]; else mcATTeffweight+=mctteff[iwb]; }2) Randomly generate an randomValTT ( randomValATT) which corresponds to the cumulative pdf of the Tag (Anti-Tag) MV1c PDF. In a word, randomValTT and randomValATT are in an uniform distribution.
//** Here we use the jet index to set the random seed,
//** which may indicate that user should store the original jet index from D3PD.
gRandom->SetSeed( EventNumber +ijet*10000001); //Now assign a MV1c value according to the efficiencies in each bin double randomValTT =gRandom->Uniform(0, mcTTeffweight); double randomValATT=gRandom->Uniform(0, mcATTeffweight);3) Calculate the the truth tagged and ant-tagged weight Scan the b-tagging weight binned PDF from 0 to 1, find the bin that corresponds to the randomValTT and randomValATT and assign the bin center as the truth tagged and ant-tagged weight.
double cumefftriedTT=0; double cumefftriedATT=0; bool doneTT=false; bool doneATT=false; double truTagWt(0), truATagWt(0); for(int i=0; i<m_nTagBin; i++) { float weightVal=MeanBWeightForBin(i); if(weightVal > m_loose_tag_cut ) { cumefftriedTT +=mctteff[i]; if(cumefftriedTT>randomValTT&&!doneTT) { //Here you should set the assigned truth-tagged MV1c value to the jet truTagWt = weightVal; doneTT=true; } } else { cumefftriedATT +=mctteff[i]; if(cumefftriedATT>randomValATT&&!doneATT) { //Here you should set the assigned anti-truth-tagged MV1c value to the jet truATagWt = weightVal; doneATT=true; } } }After we have calculated the b-tagging weight as well as the truth tagging and anti-tagging weight, we need also to get the b-tagging efficiency scale factor to compensate efficiency difference between data and MC.
Analysis::CalibrationDataVariables m_variables; m_variables.jetAuthor = "AntiKt4TopoEMJVF0_5"; m_variables.jetPt = MyJetPt; // in MeV m_variables.jetEta = MyJetEta; m_variables.jetTagWeight = MyJetTagWeightBin; // here you specify the MV1 or MV1c value of the jet std::string OperatingPoint="continuous"; int flav = myBJetFlavor; //here put the truth label of the jet according to the previous section std::pair<double,double> eff_SF = m_btagCalib -> getScaleFactor(m_variables, flav,m_bTaggerCutString , Analysis::None);For the truth tagging case (discussed below), the truth-tagged jets will acquire a randomly generated MV1c value, and you'll therefore ask for the scale factor for these. It is actually even simpler to cache both the SF value if the jet will be tagged or if it will be anti-tagged (tag = 80% efficiency), as explained in the previous section:
ajetTTEff.jetTagWeight = MyJetTagWeightBin; // truth tagged randomly generated b-weight; you need to setup pT and eta as well resnomTTEff = m_calibInterface->getScaleFactor(ajetTTEff, flav, m_bTaggerCutString, Analysis::None); ajetATTEff.jetTagWeight = MyJetTagWeightBin; // truth anti-tagged randomly generated b-weight; you need to setup pT and eta as well resnomATTEff = m_calibInterface->getScaleFactor(ajetATTEff, flav, m_bTaggerCutString, Analysis::None);Here ajetTTEff and ajetATTEff just represent the truth tagged or anti-tagged jets with the generated value of the btagging weight. In all these cases the nominal SF is given by the first value of the pair.
int NumBTagEigenVar = m_calibMV1->getNumVariations( btagjet.jetAuthor, JetFlavor, "continuous", Analysis::SFEigen);The final recommendation is still under investigation. For now we decided to use the largest 10 varations for b and lignt jet, and 15 for c jet. Because the number of variation is ordered from low to high, we have to reverse the order of the variation. The number of variations which is described (here): for b-jet, numVar is from 0 to 9, for c-jet, numVar is from 0-14, and for light jet, numVar is from 0 to 9. Y if (truthLabel== 5 ) NumBTagEigenVar = 30; else if (truthLabel== 4 ) NumBTagEigenVar = 24; else if (truthLabel== 15 ) NumBTagEigenVar = 24; int tmp_numVar = NumBTagEigenVar - (numVar+1); When you consider one of these systematic uncertainties, the jet b-tagging SF is modified by calling:
Analysis::CalibResult res_sys = m_calibMV1->getScaleFactor( btagjet, JetFlavor,"continuous", Analysis::SFEigen, tmp_numVar );Please notice: when considering systematic variations (as in the above line) the first value of the pair is the variation up, the second value is the variation down. This is different from the nominal case.
struct ObjCombination { int TagLevel; std::vector<int> m_iHiggsBJet; std::vector<int> m_iExtraJet; };
std::vector<int> m_iBJet; std::vector<int> m_iNonBJet; ObjCombination PreTag_ObjComb; ObjCombination DirectTag_ObjComb; std::vector<ObjCombination> TruthTag_ObjComb;
void ana_MMN::Select_BJets( TString options = "" ) {
m_iBJet.clear(); m_iNonBJet.clear(); for( int ijet=0; ijet<m_jetslist.size(); ijet++) { // count number of b-tagged jets double btag_Wt = m_jetslist.at(ijet).mv1; if (m_tagger == "MV1c") { btag_Wt = m_jetslist.at(ijet).mv1c; } if ( btag_Wt > m_loose_tag_cut && ( m_consider_3jBtag || ijet<=1) ) { m_iBJet.push_back( ijet ); } else { m_iNonBJet.push_back( ijet ); } } }int ana_MMN::Assign_Obj( TString options = "" ) {
//** Here assign the jets build Higgs and make the tag level //** If having 2 bjets, than the 2 bjets as Higgs jets //** If having 1 bjets, than the bjet and the leading non-bjet as Higgs jets //** If having 0 bjets, than the leading 2 non-bjets as Higgs jets //** default pretag setup DirectTag_ObjComb.TagLevel = -1; DirectTag_ObjComb.m_iHiggsBJet[0] = 0; DirectTag_ObjComb.m_iHiggsBJet[1] = 1; for(int ijet=0; ijet<m_jetslist.size()-2; ijet++) { DirectTag_ObjComb.m_iExtraJet.push_back(ijet+2); } if( m_iBJet.size()==2 ){ //** 2-tag setup DirectTag_ObjComb.m_iHiggsBJet[0] = m_iBJet[0]; DirectTag_ObjComb.m_iHiggsBJet[1] = m_iBJet[1]; DirectTag_ObjComb.m_iExtraJet = m_iNonBJet; DirectTag_ObjComb.TagLevel = 2; } else if( m_iBJet.size()==1 ){ //**!!!! need to fix the order of jets soon !!! //** 1-tag setup: Jet ordered by Pt if( m_jetslist[ m_iBJet[0] ].p.Pt()>m_jetslist[ m_iNonBJet[0] ].p.Pt() ) { DirectTag_ObjComb.m_iHiggsBJet[0] = m_iBJet[0]; DirectTag_ObjComb.m_iHiggsBJet[1] = m_iNonBJet[0]; } else { DirectTag_ObjComb.m_iHiggsBJet[0] = m_iNonBJet[0]; DirectTag_ObjComb.m_iHiggsBJet[1] = m_iBJet[0]; } if(m_jetslist.size()>=3 ) DirectTag_ObjComb.m_iExtraJet.insert( DirectTag_ObjComb.m_iExtraJet.end(), m_iNonBJet.begin()+1, m_iNonBJet.end() ); DirectTag_ObjComb.TagLevel = 1; } else if( m_iBJet.size()==0 ){ //** 0-tag setup DirectTag_ObjComb.m_iHiggsBJet[0] = m_iNonBJet[0]; DirectTag_ObjComb.m_iHiggsBJet[1] = m_iNonBJet[1]; if(m_jetslist.size()>=3 ) DirectTag_ObjComb.m_iExtraJet.insert( DirectTag_ObjComb.m_iExtraJet.end(), m_iNonBJet.begin()+2, m_iNonBJet.end() ); DirectTag_ObjComb.TagLevel = 0; } else if( m_iBJet.size()>=3 ){ //** 3-tag setup : veto or not DirectTag_ObjComb.m_iHiggsBJet[0] = m_iBJet[0]; DirectTag_ObjComb.m_iHiggsBJet[1] = m_iBJet[1]; if(m_jetslist.size()>=3 ) DirectTag_ObjComb.m_iExtraJet.insert( DirectTag_ObjComb.m_iExtraJet.end(), m_iBJet.begin()+2, m_iBJet.end() ); if(m_jetslist.size()>=3 ) DirectTag_ObjComb.m_iExtraJet.insert( DirectTag_ObjComb.m_iExtraJet.end(), m_iNonBJet.begin(), m_iNonBJet.end() ); DirectTag_ObjComb.TagLevel = 2; //*** not veto //DirectTag_ObjComb.TagLevel = -1000; //** If veto } //** The truth tagging combination if( m_jetslist.size()==2 ) { //** 2 jets case should be simple TruthTag_ObjComb.resize(1); TruthTag_ObjComb[0].TagLevel = 2; TruthTag_ObjComb[0].m_iHiggsBJet.resize(2);; TruthTag_ObjComb[0].m_iHiggsBJet[0] = 0; TruthTag_ObjComb[0].m_iHiggsBJet[1] = 1; TruthTag_ObjComb[0].m_iExtraJet.resize(0);; } else if( m_jetslist.size()==3 ) { int nCombFor3Jets(1); if( m_consider_3jBtag ) nCombFor3Jets = 3; //** 3 jets case , I need to several permutations. TruthTag_ObjComb.resize( nCombFor3Jets ); for( int icomb=0; icomb<TruthTag_ObjComb.size(); icomb++ ) { TruthTag_ObjComb[icomb].TagLevel = 2; TruthTag_ObjComb[icomb].m_iHiggsBJet.resize(2);; TruthTag_ObjComb[icomb].m_iExtraJet.resize(1);; } //** store all 3 permutation //**1** TruthTag_ObjComb[0].m_iHiggsBJet[0] = 0; TruthTag_ObjComb[0].m_iHiggsBJet[1] = 1; TruthTag_ObjComb[0].m_iExtraJet[0] = 2; if( m_consider_3jBtag ) { //**2** TruthTag_ObjComb[1].m_iHiggsBJet[0] = 0; TruthTag_ObjComb[1].m_iHiggsBJet[1] = 2; TruthTag_ObjComb[1].m_iExtraJet[0] = 1; //**3** TruthTag_ObjComb[2].m_iHiggsBJet[0] = 1; TruthTag_ObjComb[2].m_iHiggsBJet[1] = 2; TruthTag_ObjComb[2].m_iExtraJet[0] = 0; } } else { cout<<" !!!!!!!!! the number of jets is not correct !!!!!!!"<<endl; exit(EXIT_FAILURE); } return DirectTag_ObjComb.TagLevel; }
counting the number of bjets in each b-tagging cut
void ana_MMN::counting_Bjets( TString options = "" ) { m_nLoosebtagJet = 0; m_nMediumbtagJet = 0; m_nTightbtagJet = 0; m_tru_nLoosebtagJet = 0; m_tru_nMediumbtagJet = 0; m_tru_nTightbtagJet = 0; // SMWANG (08/09/13) double btag_Wt = 0.0; double tru_btag_Wt = 0.0; int NumJet = m_jetslist.size(); // Only consider b-tagging for 1st and 2nd leading jets if (!m_consider_3jBtag) { if (NumJet > 2) {NumJet = 2;} } for (int i = 0; i < NumJet; i++) { // count number of b-tagged jets btag_Wt = m_jetslist.at(i).mv1; if (m_tagger == "MV1c") { btag_Wt = m_jetslist.at(i).mv1c; } // btag wt for truth tagging tru_btag_Wt = m_jetslist.at(i).tru_mvWt; if ( btag_Wt > m_loose_tag_cut ) { m_nLoosebtagJet++; } if ( btag_Wt > m_medium_tag_cut ) { m_nMediumbtagJet++; } if ( btag_Wt > m_tight_tag_cut ) { m_nTightbtagJet++; } if ( tru_btag_Wt > m_loose_tag_cut ) { m_tru_nLoosebtagJet++; } if ( tru_btag_Wt > m_medium_tag_cut ) { m_tru_nMediumbtagJet++; } if ( tru_btag_Wt > m_tight_tag_cut ) { m_tru_nTightbtagJet++; } } // jet loop }
* Assign the LL/MM/TT categories
* Compute the b-tagging weight for the event in both direct tagging case and truth tagging case.int ana_MMN::Assign_LMT_Cateogry( ObjCombination & ObjComb ) {
int tmp_TagLevel = ObjComb.TagLevel; //** assign the tagged category in further: LL/MM/TT if( tmp_TagLevel==2 ) { int tmp_nLoosebtagJet = 0; int tmp_nMediumbtagJet = 0; int tmp_nTightbtagJet = 0; if ( m_truthTag || m_do_all_truth_tag ) { tmp_nLoosebtagJet = m_tru_nLoosebtagJet ; tmp_nMediumbtagJet = m_tru_nMediumbtagJet; tmp_nTightbtagJet = m_tru_nTightbtagJet ; } else { tmp_nLoosebtagJet = m_nLoosebtagJet ; tmp_nMediumbtagJet = m_nMediumbtagJet; tmp_nTightbtagJet = m_nTightbtagJet ; } if ( tmp_nLoosebtagJet == 2 && tmp_nMediumbtagJet < 2 && tmp_nTightbtagJet < 2) { // Loose & not Medium & not Tight tmp_TagLevel=2; } else if ( tmp_nMediumbtagJet == 2 && tmp_nTightbtagJet < 2) { // Medium & not Tight tmp_TagLevel=3; } else if ( tmp_nTightbtagJet == 2 ) { // Tight tmp_TagLevel=4; } } return tmp_TagLevel; }
void ana_MMN::comput_btag_weight( ObjCombination & ObjComb ) { if ( 1 <= m_debug ) { cout<<" In comput_btag_weight: Compute the btag weight-- scale factor and truth effieincy! "<<endl; } int tmp_taglevel = ObjComb.TagLevel; int tmp_iJet1 = ObjComb.m_iHiggsBJet[0]; int tmp_iJet2 = ObjComb.m_iHiggsBJet[1]; int tmp_iJet3 = -1; if( ObjComb.m_iExtraJet.size()!=0 ) { tmp_iJet3 = ObjComb.m_iExtraJet[0]; } m_btagging_weight_SF = 1.; // needed to be initalised at 1 m_btagging_weight_tru_SF = 1.; // needed to be initalised at 1 // m_btagging_weight_SF *= m_jetslist.at(tmp_iJet1).scf; m_btagging_weight_SF *= m_jetslist.at(tmp_iJet2).scf; if( ObjComb.m_iExtraJet.size()!=0 ) { if( m_consider_3jBtag ) m_btagging_weight_SF *= m_jetslist.at(tmp_iJet3).scf; } m_btagging_weight_tru_SF *= m_jetslist.at(tmp_iJet1).trf * m_jetslist.at(tmp_iJet1).tru_scf; m_btagging_weight_tru_SF *= m_jetslist.at(tmp_iJet2).trf * m_jetslist.at(tmp_iJet2).tru_scf; if( ObjComb.m_iExtraJet.size()!=0 ) { if( m_consider_3jBtag ) m_btagging_weight_tru_SF *= ( 1 - m_jetslist.at(tmp_iJet3).trf * m_jetslist.at(tmp_iJet3).tru_inscf ); } m_btagging_weight_TRF_2tag = m_btagging_weight_tru_SF; }