mtw factor

Tool/EffsHeader_MuCS.h

@@ -1,3 +1,4 @@
+  const double ttbar_mtwcorrfactor[8] = {1.00308,1.03938,1.0767,1.0988,1.12942,1.17955,1.25167,1.7567};
   const double ttbar_mus_acc[6] = {0.769261,0.871118,0.900174,0.92814,0.921287,0.945285};
   const double ttbar_mus_recoeff[8][8] = {{0.952769,0.906542,0.939655,0.93617,0.93038,0.914894,0.92,0.903766},{0.962189,0.966667,0.930147,0.963989,0.963768,0.915254,0.952381,0.95539},{0.971917,0.966019,0.97619,0.934911,0.964126,0.945652,0.951456,0.976048},{0.979784,0.980769,0.949239,0.944206,0.940367,0.960526,0.971429,0.975},{0.980216,0.943089,0.992958,0.986486,0.982955,0.962963,0.9625,0.953757},{0.981664,0.965517,0.949495,0.972125,0.929204,0.960452,0.971223,0.951557},{0.96882,0.990566,0.958042,0.985294,0.954545,0.954545,0.934579,0.95053},{0.925926,0.928571,0.92,0.896552,0.913043,0.9,0.918367,0.932773}};
   const double ttbar_mus_isoeff[8][8] = {{0.926544,0.808081,0.77193,0.73743,0.729167,0.752688,0.71831,0.599099},{0.966805,0.884058,0.887597,0.849432,0.833948,0.825,0.762295,0.73622},{0.980271,0.945545,0.936585,0.898773,0.884259,0.8125,0.913462,0.766082},{0.995851,0.993464,0.958115,0.930435,0.925,0.899281,0.913462,0.869792},{0.987544,0.991379,0.962687,0.965686,0.955801,0.935185,0.878378,0.899441},{0.991058,0.982143,0.984127,1,0.990338,0.977528,0.945736,0.930147},{0.997625,1,0.992701,0.995,0.978947,1,0.989796,0.958801},{1,1,1,1,1,1,1,0.996865}};

Tool/LostLepton_MuCS_TTbar.cc

@@ -218,15 +218,15 @@ int main(int argc, char* argv[])
   int nevents_baseline= 0;
   //int nevents_baseline_ref= 0;
 
-  double mtwcorrfactor[8];
-  mtwcorrfactor[0] = 1.01;
-  mtwcorrfactor[1] = 1.04;
-  mtwcorrfactor[2] = 1.07;
-  mtwcorrfactor[3] = 1.20;
-  mtwcorrfactor[4] = 1.15;
-  mtwcorrfactor[5] = 1.20;
-  mtwcorrfactor[6] = 1.28;
-  mtwcorrfactor[7] = 1.83;
+  //double mtwcorrfactor[8];
+  //mtwcorrfactor[0] = 1.01;
+  //mtwcorrfactor[1] = 1.04;
+  //mtwcorrfactor[2] = 1.07;
+  //mtwcorrfactor[3] = 1.20;
+  //mtwcorrfactor[4] = 1.15;
+  //mtwcorrfactor[5] = 1.20;
+  //mtwcorrfactor[6] = 1.28;
+  //mtwcorrfactor[7] = 1.83;
 
   //first loop, to generate Acc, reco and Iso effs and also fill expected histgram
   std::cout<<"First loop begin: "<<std::endl;
@@ -671,6 +671,37 @@ int main(int argc, char* argv[])
         }//loop gen electrons/muons
       }//if no muons
 
+      //mtW correction factor calculation
+      if (nElectrons == 0 && nMuons == 1)
+      {
+        vector<TLorentzVector> muonsLVec = tr.getVec<TLorentzVector>("muonsLVec");
+        vector<double> muonsMiniIso = tr.getVec<double>("muonsMiniIso");
+
+        double reco_mus_pt = 0, reco_mus_eta = 0, reco_mus_phi = 0;
+
+        for(unsigned int im = 0 ; im < muonsLVec.size() ; im++)
+        {
+          if( fabs(muonsLVec[im].Eta()) < 2.4 && muonsMiniIso[im] < 0.2 )
+          {
+            reco_mus_pt  = ( muonsLVec.at(im) ).Pt();
+            reco_mus_eta = ( muonsLVec.at(im) ).Eta();
+            reco_mus_phi = ( muonsLVec.at(im) ).Phi();
+          }
+        }
+
+        double deltaphi_mus = DeltaPhi( reco_mus_phi , metphi );
+        double mtW_mus = std::sqrt( 2.0 * reco_mus_pt * met * ( 1.0 - cos(deltaphi_mus) ) );
+        (myBaseHistgram.h_mtw_mus)->Fill(mtW_mus);
+
+        int ptbin_number_allreco = Set_ptbin_number(reco_mus_pt);
+
+        myAccRecoIsoEffs.mtwall[ptbin_number_allreco]++;        
+        if( mtW_mus < 100 )
+        {
+          myAccRecoIsoEffs.mtw100[ptbin_number_allreco]++;
+        }
+      }
+
       //book the variable we need to check
       int ntopjets = tr.getVar<int>("nTopCandSortedCnt");
       double MT2 = tr.getVar<double>("MT22");
@@ -933,7 +964,7 @@ int main(int argc, char* argv[])
           int acbin_number = Set_acbin_number(activity);
 
 	  //mtwcorrfactor
-	  EventWeight_mus = EventWeight_mus * mtwcorrfactor[ptbin_number];
+	  EventWeight_mus = EventWeight_mus * myAccRecoIsoEffs.mtwcorrfactor[ptbin_number];
 	  //dimuon correction factor
           EventWeight_mus = EventWeight_mus * myAccRecoIsoEffs.corrfactor_di_mus;
 
@@ -987,7 +1018,7 @@ int main(int argc, char* argv[])
           //begin to predict lost electrons from muon CS
           double EventWeight_els = 1.0;
           //mtwcorrfactor
-          EventWeight_els = EventWeight_els * mtwcorrfactor[ptbin_number];
+          EventWeight_els = EventWeight_els * myAccRecoIsoEffs.mtwcorrfactor[ptbin_number];
           //dielectron correction factor
           EventWeight_els = EventWeight_els * myAccRecoIsoEffs.corrfactor_di_els;
 
@@ -1112,6 +1143,12 @@ void AccRecoIsoEffs::NumberstoEffs()
     }
   }
 
+  for(i_cal = 0 ; i_cal < PT_BINS ; i_cal++)
+  {
+    mtwcorrfactor[i_cal] = mtwall[i_cal]/mtw100[i_cal];
+    mtwcorrfactor_err[i_cal] = get_stat_Error(mtw100[i_cal],mtwall[i_cal]);
+  }
+
   return ;
 }
 
@@ -1318,6 +1355,18 @@ void AccRecoIsoEffs::printAccRecoIsoEffs()
   int i_cal = 0;
   int j_cal = 0;
   std::cout.precision(3);
+
+  std::cout << "mtW correction factor: " << std::endl;
+  for( i_cal=0 ; i_cal < PT_BINS ; i_cal++ )
+  {
+    std::cout << mtwcorrfactor[i_cal] << "(" << mtwcorrfactor_err[i_cal] << ")"<< " ";;
+    if( i_cal == PT_BINS-1 )
+    {
+      std::cout << std::endl;
+    }
+  }
+
+
   std::cout << std::endl << "Muon information: " << std::endl;
 
   std::cout << "number of muons from top: " << std::endl;
@@ -1675,6 +1724,15 @@ void AccRecoIsoEffs::printEffsHeader()
   int i_cal = 0;
   int j_cal = 0;
 
+  EffsHeader << "  const double ttbar_mtwcorrfactor[" << PT_BINS << "] = ";
+  for( i_cal = 0 ; i_cal < PT_BINS ; i_cal++ )
+  {
+    if( i_cal == 0 ) { EffsHeader << "{"; }
+    EffsHeader << mtwcorrfactor[i_cal];
+    if( i_cal != PT_BINS-1 ) { EffsHeader << ","; }
+    if( i_cal == PT_BINS-1 ) { EffsHeader << "};" << std::endl; }
+  }
+
   EffsHeader << "  const double ttbar_mus_acc[" << NJETS_BINS << "] = "; 
   for( i_cal = 0 ; i_cal < NJETS_BINS ; i_cal++ )
   {

Tool/LostLepton_MuCS_TTbar.h

@@ -52,6 +52,8 @@ class AccRecoIsoEffs
 
   double mus_isoeff_err_allreco[PT_BINS][AC_BINS] = {{0}}, els_isoeff_err_allreco[PT_BINS][AC_BINS] = {{0}};
 
+  double mtwcorrfactor[PT_BINS] = {0}, mtwcorrfactor_err[PT_BINS] = {0};
+
   //here we define the muon/electron number we need to count in the loop
   double nmus[NJETS_BINS] = {0}, nmus_acc[NJETS_BINS] = {0}, nels[NJETS_BINS] = {0}, nels_acc[NJETS_BINS] = {0};
   double nmus_acc_bin[PT_BINS][AC_BINS] = {{0}}, nels_acc_bin[PT_BINS][AC_BINS] = {{0}};
@@ -61,6 +63,8 @@ class AccRecoIsoEffs
   double nmus_reco_allreco[PT_BINS][AC_BINS] = {{0}}, nels_reco_allreco[PT_BINS][AC_BINS] = {{0}};
   double nmus_iso_allreco[PT_BINS][AC_BINS] = {{0}}, nels_iso_allreco[PT_BINS][AC_BINS] = {{0}};
 
+  double mtwall[PT_BINS] = {0}, mtw100[PT_BINS] = {0};
+
   //here we define the event weight we are going to use in the second loop ( muon/electron CS and prediction plots)
   double mus_EventWeight_iso[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}}, mus_EventWeight_reco[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}}, mus_EventWeight_acc[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}};
   double els_EventWeight_iso[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}}, els_EventWeight_reco[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}}, els_EventWeight_acc[NJETS_BINS][PT_BINS][AC_BINS] = {{{0}}};
@@ -190,6 +194,8 @@ class BaseHistgram
   TH1D *h_b_njets30_4_eta_els, *h_b_njets30_5_eta_els, *h_b_njets30_6_eta_els, *h_b_njets30_7_eta_els, *h_b_njets30_8_eta_els, *h_b_njets30_9_eta_els;
   TH1D *h_b_njets30_4_ht_mus, *h_b_njets30_5_ht_mus, *h_b_njets30_6_ht_mus, *h_b_njets30_7_ht_mus, *h_b_njets30_8_ht_mus, *h_b_njets30_9_ht_mus;
 
+  TH1D *h_mtw_mus;
+
   //closure plots definition
   TH1D *h_pred_mu_iso_met, *h_pred_mu_iso_njets, *h_pred_mu_iso_mt2, *h_pred_mu_iso_topmass, *h_pred_mu_iso_ht, *h_pred_mu_iso_mht, *h_pred_mu_iso_ntopjets;
   TH1D *h_pred_mu_id_met, *h_pred_mu_id_njets, *h_pred_mu_id_mt2, *h_pred_mu_id_topmass, *h_pred_mu_id_ht, *h_pred_mu_id_mht, *h_pred_mu_id_ntopjets;
@@ -304,6 +310,8 @@ void BaseHistgram::BookHistgram(const char *outFileName)
   h_b_njets30_8_eta_els = new TH1D("h_b_njets30_8_eta_els","",60,-3,3);
   h_b_njets30_9_eta_els = new TH1D("h_b_njets30_9_eta_els","",60,-3,3);
 
+  h_mtw_mus = new TH1D("h_mtw_mus","",200,0,200);
+
   h_b_activity_mus = new TH1D("h_b_activity_mus","",1000,0,200);
   h_b_activity_els = new TH1D("h_b_activity_els","",1000,0,200);