API tools faq
paste
Advertisement
Guest User

C++ Source HEP Analysis2

a guest
Aug 29th, 2014
252
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 28.05 KB | None | 0 0
raw download clone embed print report
  1. #include "NtupleTools3.h"
  2. #include <fstream>
  3. #include <iostream>
  4. #include <vector>
  5. #include <TLorentzVector.h>
  6. #include <TCanvas.h>
  7. #include <TH1F.h>
  8. #include <TFile.h>
  9.  
  10. using namespace std;
  11.  
  12. void JackreaderSummerStd(TString list, TString outname,bool useW=true){
  13.  
  14.   TObjArray* arr = list.Tokenize(" ");
  15.   int size=arr->GetEntries();
  16.   if(size%2!=0) {
  17.     cout<<"unbalance file/weight list: "<<list<<endl;
  18.     exit(0);
  19.   }
  20.       vector<TString> files;
  21.       vector<double> weights;
  22.         for(int i=0;i<size;i+=2){
  23.       files.push_back( arr->At(i)->GetName() );
  24.       weights.push_back( atof( arr->At(i+1)->GetName() ) );
  25.     }
  26.  
  27.  
  28.         // here is the place for histograms --- learn how to determine ranges
  29.         TH1F * hHT             = new TH1F("hHT","H_{T} [GeV]",100,0,5000);
  30.     TH1F * hNel            = new TH1F("hNel","Number of Electrons",5,0.5,5.5);
  31.     TH1F * hNmu            = new TH1F("hNmu","Number of Muons",5,0.5,5.5);
  32.     TH1F * hNlep           = new TH1F("hNlep","Lepton Multiplicity",5,0.5,5.5);
  33.         TH1F * hNjet           = new TH1F("hNjet","Jet Multiplicity",14,0.5,14.5);
  34.     TH1F * hMET            = new TH1F("hMET", "E_{T}^{miss} [GeV]", 40,0,1000);
  35.     TH1F * hDphi           = new TH1F("hDphi", "#Delta#phi", 40,0,4);
  36.     TH1F * hLep1Pt         = new TH1F("hLep1Pt","Leading Lepton P_{T} [GeV]",100,0,600);
  37.     TH1F * hLep2Pt         = new TH1F("hLep2Pt","2^{nd} Leading Lepton P_{T} [GeV]",50,0,600);
  38.     TH1F * hLep3Pt         = new TH1F("hLep3Pt","3^{nd} Leading Lepton P_{T} [GeV]",50,0,600);
  39.     TH1F * hJet1Pt         = new TH1F("hJet1Pt","1^{st} Leading Jet P_{T} [GeV]",50,0,1000);
  40.     TH1F * hJet2Pt         = new TH1F("hJet2Pt","2^{nd} Leading Jet P_{T} [GeV]",50,0,1000);
  41.     TH1F * hJet3Pt         = new TH1F("hJet3Pt","3^{rd} Leading Jet P_{T} [GeV]",50,0,1000);
  42.         TH1F * hDiLMass        = new TH1F("hDiLMass","Lepton Mass [GeV]",25,0,250);
  43.     //        TH1F * hDiLMass2       = new TH1F("hDiLMass2","OS DiLepton mass",30,0,450);
  44.         TH1F * hMultiLMassSF12 = new TH1F("hMultiLMassSF12","Lepton mass [GeV]",30,0,450);
  45.     TH1F * hMETMeff        = new TH1F("hMETMeff","E^{miss}_{T}/E^{miss}_{T}+H_{T}", 30,0,1);
  46.     TH1F * hMETLep         = new TH1F("hMETLep", "E^{miss}_{T}/Lepton P_{T}", 100, 0, 1);
  47.         // TH1F * hMassDifMT      = new TH1F("hMassDifMT","MW-MT/MT",100,0,1);
  48.         // TH1F * hMassDifMET     = new TH1F("hMassDifMET","E^{miss}_{T}-M_{T}/M_{T}",100,0,1);
  49.     TH1F * hBelowZ         = new TH1F("hBelowZ", "Lepton mass Below-Z region [GeV]",25,10,80);
  50.     TH1F * hAboveZ         = new TH1F("hAboveZ", "Lepton Mass Above-Z region [GeV]",35,100,400);
  51.     //  TH1F * hTopMass        = new TH1F("hTopmass", "Top Mass",100,400,3000);
  52.     //  TH1F * hWMass          = new TH1F("hWMass", "W mass", 100,0,2000);
  53.     //  TH1F * hMT3            = new TH1F("hMT3", "M_{T} [GeV]", 100,0,200);
  54.     TH1F * hMeff           = new TH1F("hMeff", "M_{eff} [GeV]", 50,0,5000);
  55.     TH1F * h1LepEta        = new TH1F("h1LepEta", "Leading Lepton #eta", 50,-5,5);
  56.     TH1F * h2LepEta        = new TH1F("h2LepEta", "2^{nd} Leading Lepton #eta", 50,-5,5);
  57.     TH1F * h3LepEta        = new TH1F("h3LepEta", "3^{rd} Leading Lepton #eta", 50,-5,5);
  58.     TH1F * h1JetEta        = new TH1F("h1JetEta", "Leading Jet #eta", 50,-5,5);
  59.     TH1F * h2JetEta        = new TH1F("h2JetEta", "2^{nd} Leading Jet #eta", 50,-5,5);
  60.     TH1F * h3JetEta        = new TH1F("h3JetEta", "3^{rd} Leading Jet #eta", 50,-5,5);
  61.     // TH1F * hMuIso          = new TH1F("hMuIso", "Isolation", 20,0,1);
  62.     // TH1F * hMuIso          = new TH1F("hMuIso", "Isolation", 20,0,1);
  63.     TH1F * hST             = new TH1F("hST","S_{T} [GeV]",100,0,5000);
  64.  
  65.     //  TH1F * hMTpre          = new TH1F("hMTpre","M_{T} pre [GeV]",100,0,200);
  66.     ////////////////////// first check plots
  67.     TH1F * hHT1             = new TH1F("hHT1","H_{T} first check [GeV]",100,0,5000);
  68.     TH1F * hNel1            = new TH1F("hNel1","Number of Electrons first check",5,0.5,5.5);
  69.     TH1F * hNmu1            = new TH1F("hNmu1","Number of Muons first check",5,0.5,5.5);
  70.     TH1F * hNlep1           = new TH1F("hNlep1","Number of leptons first check",5,0.5,5.5);
  71.     TH1F * hNjet1           = new TH1F("hNjet1","Number of Jet first check",14,0.5,14.5);
  72.     TH1F * hMET1            = new TH1F("hMET1", "E^{miss}_{T} first check", 40,0,2000);
  73.     TH1F * hDphi1           = new TH1F("hDphi1", "#Delta#phi first check", 40,0,4);
  74.     TH1F * hLep1Pt1         = new TH1F("hLep1Pt1","Leading lepton P_{T} first check",50,0,600);
  75.     TH1F * hMultiLMassSF1   = new TH1F("hMultiLMassSF1","Lepton mass first check",30,0,450);
  76.     TH1F * hLep2Pt1         = new TH1F("hLep2Pt1","2^{nd} Leading Lepton P_{T} first check",50,0,600);
  77.     TH1F * hLep3Pt1         = new TH1F("hLep3Pt1","3^{nd} Leading Lepton P_{T} first check",50,0,600); 
  78.     TH1F * hJet1Pt1         = new TH1F("hJet1Pt1","1^{st} Leading Jet P_{T} first check",50,0,1000);
  79.     TH1F * hJet2Pt1         = new TH1F("hJet2Pt1","2^{nd} Leading Jet P_{T} first check",50,0,1000);
  80.     TH1F * hJet3Pt1         = new TH1F("hJet3Pt1","3^{rd} Leading Jet P_{T} first check",50,0,1000);
  81.     TH1F * hDiLMass1        = new TH1F("hDiLMass1","OS DiLepton mass first check",25,0,250);
  82.     TH1F * hMETMeff1        = new TH1F("hMETMeff1","E^{miss}_{T}/E^{miss}_{T}+H_{T} first check", 30,0,1);
  83.     TH1F * hMeff1           = new TH1F("hMeff1", "M_{eff} first check", 50,0,5000);
  84.     TH1F * h1LepEta1        = new TH1F("h1LepEta1", "#eta_{1lep} first check", 50,-5,5);
  85.     TH1F * h2LepEta1        = new TH1F("h2LepEta1", "#eta_{2lep} first check", 50,-5,5);
  86.     TH1F * h3LepEta1        = new TH1F("h3LepEta1", "#eta_{3lep} first check", 50,-5,5);
  87.     TH1F * h1JetEta1        = new TH1F("h1JetEta1", "#eta_{1jet} first check", 50,-5,5);
  88.     TH1F * h2JetEta1        = new TH1F("h2JetEta1", "#eta_{2jet} first check", 50,-5,5);
  89.     TH1F * h3JetEta1        = new TH1F("h3JetEta1", "#eta_{3jet} first check", 50,-5,5);
  90.     TH1F * hMuIso1          = new TH1F("hMuIso1", "IsoMu  first check", 20,0,1);
  91.     TH1F * hST1             = new TH1F("hST1","S_{T} first check",100,0,5000); 
  92.     TH1F * hBjet            = new TH1F("hBjet","Bjet multiplicity",10,0.5,10.5);
  93.     TH1F * hMuiso           = new TH1F("hMuiso","isolation mu",20,-1.5,1.5);
  94.     TH1F * hEliso           = new TH1F("hEliso","isolation el",20,-1.5,1.5);
  95.     //  TH1F * hMT31           = new TH1F("hMT31", "M_{T} first check", 100,0,200);
  96.     /////////////////// end of first check histos
  97.     ////////////////// second check histos
  98.     TH1F * hHT2             = new TH1F("hHT2","H_{T} second check [GeV]",100,0,5000);
  99.     TH1F * hNel2            = new TH1F("hNel2","Number of Electrons second check",5,0.5,5.5);
  100.     TH1F * hNmu2            = new TH1F("hNmu2","Number of muons second check",5,0.5,5.5);
  101.     TH1F * hNlep2           = new TH1F("hNlep2","Number of leptons second check",5,0.5,5.5);
  102.     TH1F * hNjet2           = new TH1F("hNjet2","Number of Jet second check",14,0.5,14.5);
  103.     TH1F * hMET2            = new TH1F("hMET2", "E^{miss}_{T} second check [GeV]", 40,0,2000);
  104.     TH1F * hDphi2           = new TH1F("hDphi2", "#Delta#phi second check ", 40,0,4);
  105.     TH1F * hLep1Pt2         = new TH1F("hLep1Pt2","Leading lepton P_{T} second check",50,0,600);
  106.     TH1F * hMultiLMassSF2   = new TH1F("hMultiLMassSF2","Lepton mass  second check [GeV]",25,0,250);
  107.     TH1F * hLep2Pt2         = new TH1F("hLep2Pt2","2^{nd} Leading Lepton P_{T} second check",50,0,600);
  108.     TH1F * hLep3Pt2         = new TH1F("hLep3Pt2","3^{nd} Leading Lepton P_{T} second check",50,0,600);
  109.     TH1F * hJet1Pt2         = new TH1F("hJet1Pt2","1^{st} Leading Jet P_{T} second check",50,0,1000);
  110.     TH1F * hJet2Pt2         = new TH1F("hJet2Pt2","2^{nd} Leading Jet P_{T} second check",50,0,1000);
  111.     TH1F * hJet3Pt2         = new TH1F("hJet3Pt2","3^{rd} Leading Jet P_{T} second check",50,0,1000);
  112.     TH1F * hDiLMass2        = new TH1F("hDiLMass2","Lepton mass second check [GeV]",30,0,250);
  113.     TH1F * hMETMeff2        = new TH1F("hMETMeff2","E^{miss}_{T}/E^{miss}_{T}+H_{T} second check", 100,0,1);
  114.     TH1F * hMeff2           = new TH1F("hMeff2", "M_{eff} second check [GeV]", 100,0,5000);
  115.     TH1F * h1LepEta2        = new TH1F("h1LepEta2", "#eta_{1lep} second check", 50,-5,5);
  116.     TH1F * h2LepEta2        = new TH1F("h2LepEta2", "#eta_{2lep}second check", 50,-5,5);
  117.     TH1F * h3LepEta2        = new TH1F("h3LepEta2", "#eta_{3lep}second check", 50,-5,5);
  118.     TH1F * h1JetEta2        = new TH1F("h1JetEta2", "#eta_{1jet}second check", 50,-5,5);
  119.     TH1F * h2JetEta2        = new TH1F("h2JetEta2", "#eta_{2jet}second check", 50,-5,5);
  120.     TH1F * h3JetEta2        = new TH1F("h3JetEta2", "#eta_{3jet}second check", 50,-5,5);
  121.     //  TH1F * hMuIso2          = new TH1F("hMuIso2", "IsoMu second check", 20,0,1);
  122.     TH1F * hST2             = new TH1F("hST2","S_{T} second check",100,0,5000);
  123.     //  TH1F * hMT32           = new TH1F("hMT32", "M_{T} second check", 100,0,200);
  124.     ///////////////// end of second check histos
  125.     EasyChain* tree    = new EasyChain("delphTree");
  126.    
  127.      for(unsigned i=0;i<files.size();i++){
  128.        tree->AddSmartW(files[i],weights[i]);
  129.        cout<<"add: "<<files[i]<<" "<<weights[i]<<endl;
  130.      }
  131.  
  132.      int Nevents=tree->GetEntries();
  133.     cout<<">>>>>>>>>>>>>>>>>>>>>>> total number of events:\t" << Nevents <<endl;
  134.    
  135.     // Text counter
  136.     const int CutNumb = 8;
  137.     const char * CutList[CutNumb] = {"Preselection", "Trigger", "N#mu=3", "SMZcut"/*, "JetMultiplicity>4"*/, "BJetMultip>4", "dPhi>0.3", "MET>250", "MET>250HT>1000"};
  138.    
  139.         // Counter start here:
  140.         double CFCounter[CutNumb];
  141.     int   iCFCounter[CutNumb];
  142.    
  143.     for (int i=0;i < CutNumb; i++){
  144.       CFCounter[i] = 0;
  145.       iCFCounter[i] = 0;
  146.     }
  147.    
  148.     // Histogram counter
  149.     const int CutNumb1 = 7;
  150.     const char * CutList1[CutNumb1] = {"Preselect.","Trigger", "N#mu = 3", "SM Z Req."/*, "Jet Cleaning & Multip."*/, "Numb. BJets > 4", "#phi > 0.3", "H_{T}>750 E^{miss}_{T}>250"};
  151.  
  152.     double CFCounter1[CutNumb1];
  153.         int   iCFCounter1[CutNumb1];
  154.  
  155.     for (int i=0;i < CutNumb1; i++){
  156.       CFCounter1[i] = 0;
  157.       iCFCounter1[i] = 0;
  158.         }
  159.  
  160.     TH1D * CutFlow     = new TH1D("CutFlow","",CutNumb1,0.5,CutNumb1+0.5); 
  161.     TH1D * CutFlow1    = new TH1D("CutFlow1","",CutNumb,0.5,CutNumb+0.5);  
  162.    
  163.     // label bins
  164.         for(int cj = 0; cj < CutNumb1; cj++)
  165.       CutFlow->GetXaxis()->SetBinLabel(cj+1,CutList1[cj]);
  166.         for(int cj = 0; cj < CutNumb; cj++)
  167.       CutFlow1->GetXaxis()->SetBinLabel(cj+1,CutList[cj]);
  168.    
  169.     for(int entry=0; entry < Nevents; entry+=1)
  170.       {
  171.        
  172.         progressT();
  173.         double fw = tree->GetEntryW(entry); // the applied with AddSmartW for the current file/dir
  174.        
  175.        
  176.         double EvWeight = 1;
  177.         if(useW) EvWeight = tree->Get(EvWeight,"EventWeight");
  178.         EvWeight *= fw * 1000; // 1000 for 1/pb->1/fb
  179.        
  180.        
  181.         // 0. CF presel
  182.         CFCounter[0]+= EvWeight;
  183.         iCFCounter[0]++;
  184.         CFCounter1[0]+= EvWeight;
  185.         iCFCounter1[0]++;
  186.        
  187.         // vectors::::
  188.        
  189.         std::vector<TLorentzVector> Test;// = new vector<TLorentzVector>;
  190.         std::vector<TLorentzVector> &Electrons = tree->Get(&Electrons,"Electrons");
  191.         std::vector<int> &ElCh = tree->Get(&ElCh,"ElectronCh");
  192.         std::vector<TLorentzVector> &Muons = tree->Get(&Muons,"Muons");
  193.         std::vector<int> &MuCh = tree->Get(&MuCh,"MuonCh");
  194.         std::vector<double> &ElIso = tree->Get(&ElIso,"ElectronIso");
  195.         std::vector<double> &MuIso = tree->Get(&MuIso,"MuonIso");
  196.         std::vector<TLorentzVector> &Jets = tree->Get(&Jets,"Jets");
  197.         std::vector<int> &JetB = tree->Get(&JetB, "JetB");
  198.         std::vector<TLorentzVector> Lepton;// = new vector<TLorentzVector>;
  199.         std::vector<double> &JetMETdPhi = tree->Get(&JetMETdPhi, "JetMETdPhi");
  200.         std::vector<TLorentzVector> newJets;         
  201.          
  202.         vector<int> LepFlavour;
  203.        
  204.         int Nel = tree->Get( Nel, "Nel");
  205.         int Nmu = tree->Get( Nmu, "Nmu");
  206.        
  207.         int Nlept = Nel + Nmu;
  208.        
  209.         bool  vetoE[100];for(int ie=0;ie<Nel;ie++)vetoE[ie]=0;
  210.         bool  vetoM[100];for(int im=0;im<Nmu;im++)vetoM[im]=0;
  211.        
  212.         // EWKino leptons
  213.        
  214.         int nlept=0;
  215.         for(int ie=0;ie<Nel;ie++) {
  216.           vetoE[ie]|=(Electrons[ie].Pt()<10)||(fabs(Electrons[ie].Eta())>2.4);
  217.           nlept++;
  218.         }
  219.         for(int im=0;im<Nmu;im++) {
  220.           vetoM[im]|=(Muons[im].Pt()<10)||(fabs(Muons[im].Eta())>2.4);
  221.           nlept++;
  222.         }
  223.         vector<bool> veto;
  224.         for(int elid = 0; elid < Nel; elid++)            
  225.           { veto.push_back(vetoE[elid]);
  226.             Lepton.push_back(Electrons[elid]);
  227.             LepFlavour.push_back(11*ElCh[elid]);
  228.           }
  229.        
  230.         for(int muid = 0; muid < Nmu; muid++)
  231.           { veto.push_back(vetoM[muid]);
  232.             Lepton.push_back(Muons[muid]);
  233.             LepFlavour.push_back(13*MuCh[muid]);
  234.           }
  235.        
  236.         double LeadPt = 0;
  237.         TLorentzVector lv_sum;
  238.         int hasSF = 1; // has same flavour Multilepton (muon)
  239.         int hasZ = 0; // has Z like pair
  240.         double m=0;
  241.        
  242.        
  243.         for(int lid = 0; lid < Nlept; lid++)
  244.           {
  245.             if(veto[lid]) continue;
  246.             for(int sid = 0; (sid < Nlept) && (lid != sid); sid++)
  247.               {
  248.             if(veto[sid]) continue;
  249.             {
  250.               int sumflvr = LepFlavour[lid] + LepFlavour[sid]; // summ of lepton flavour, possible comb: SS: 22,24,26 OS: 0,2 *********
  251.               if(abs(sumflvr) == 26) {
  252.                 hasSF++;
  253.                 m=lv_sum.M();
  254.                
  255.                 if(m>0)hMultiLMassSF12->Fill(m,EvWeight);
  256.               }
  257.              
  258.               if(sumflvr == 0) // may be Z decay ****************
  259.                 {
  260.                   lv_sum = Lepton[lid] + Lepton[sid];
  261.                   if(abs(lv_sum.M()-90.) < 15) hasZ++;
  262.                 }
  263.             }    
  264.               }
  265.            
  266.             if(Lepton[lid].Pt() > LeadPt)
  267.               LeadPt = Lepton[lid].Pt();
  268.           }
  269.        
  270.        
  271.         /////////////////// lepton requirement
  272.         bool trigger = true;
  273.         for (int i = 0; i < Nlept; i++)
  274.           if((Lepton[i].Pt() < 10.) || (fabs(Lepton[i].Eta()) > 2.4) /*&& (Lepton[0].Pt() < 20.)*/) trigger = false;
  275.         if(!trigger) continue;
  276.         for(int i=0; i<Nel;i++)
  277.           {
  278.             if(Electrons[i].Pt()>10)hEliso->Fill(ElIso[0],EvWeight);       
  279.           }
  280.         for(int i=0; i<Nmu;i++)
  281.           {
  282.             if(Muons[i].Pt()>10)hMuiso->Fill(MuIso[0],EvWeight);
  283.           }
  284.         //////////////////// Isolation
  285.         bool isomu = false, isoel = false;
  286.         for(int i=0; i<Nmu; i++)
  287.           {
  288.             if(MuIso[i]/*/Muons[i].Pt()*/ > 0.15) isomu = true;
  289.           }
  290.         for(int i=0; i<Nel; i++)
  291.           {
  292.             if(ElIso[i]/*/Electrons[i].Pt()*/ > 0.15) isoel = true;
  293.           }
  294.         if(isoel || isomu) continue;
  295.        
  296.         int trigger_btag[50];  
  297.         int bind = 0;
  298.         ////////////////////// jet requirement
  299.         std::vector<TLorentzVector> triggerJets;
  300.         int Njet_loose = tree-> Get(Njet_loose,"Njet");
  301.         for(int i=0; i<Njet_loose;i++)
  302.           {            
  303.             if((Jets[i].Pt()>40.) && (fabs(Jets[i].Eta())<2.4))
  304.               { triggerJets.push_back(Jets[i]); trigger_btag[bind]=JetB.at(i); bind++;}
  305.            
  306.           }
  307.        
  308.         bind = 0;
  309.         int clean_btag[50];
  310.         //////////////////// jet cleaning
  311.         bool clean = true;
  312.         for ( int s=0; s<triggerJets.size(); s++)
  313.           {
  314.             clean = true;
  315.             for ( int i=0;i<Lepton.size();++i)
  316.               {
  317.             if (triggerJets[s].DeltaR(Lepton[i]) <= 0.4)
  318.               {
  319.                 clean = false;
  320.                 continue;
  321.               }
  322.               }
  323.             if(clean == true) {newJets.push_back(triggerJets[s]); clean_btag[bind]=trigger_btag[s];bind++;}
  324.            
  325.           }
  326.        
  327.  
  328.  
  329.  
  330.        
  331.  
  332.         ///////////////////// HT
  333.         ///////////////////// NOTE: HTnew <= HT40 <= HTold, HTold40 = HT40
  334.         double HT    =  tree->Get(HT,"HT40");
  335.         double HTnew = 0;
  336.         double HTold = 0;
  337.         double HTold40 = 0;
  338.         for(int i=0; i<Jets.size(); i++)
  339.           {
  340.             HTold += Jets[i].Pt();
  341.             if ( Jets[i].Pt() > 40 ) HTold40+= Jets[i].Pt();
  342.           }
  343.         for (int i=0; i<newJets.size(); i++)
  344.           {
  345.             HTnew+=newJets[i].Pt();
  346.           }
  347.         //      cout << "HT40:\t " << HT << "\t newHT:\t " << HTnew << "\t HTold40: \t" << HTold40 << endl;
  348.  
  349.         double MET   =  tree->Get(MET,"DelphMET");
  350.         double Meff  =  MET + HTnew;
  351.         double LepPt = 0;
  352.         for(int i=0; i<Nlept; i++)
  353.           LepPt +=Lepton[i].Pt();
  354.         double ST    =  LepPt + MET;
  355.        
  356.         int Nbjet = 0;       
  357.         vector<int> BJetInd;
  358.         for(int i = 0; i < newJets.size(); ++i)
  359.           {
  360.             if((clean_btag[i] & (1 << 1)) && (newJets[i].Pt() > 50.) && (fabs(newJets[i].Eta()) < 1.8))
  361.               {
  362.             Nbjet++;
  363.             BJetInd.push_back(i);
  364.               }
  365.           }
  366.         hBjet->Fill(Nbjet,EvWeight);
  367.         CFCounter[1]+= EvWeight;
  368.         iCFCounter[1]++;
  369.         CFCounter1[1]+= EvWeight;
  370.         iCFCounter1[1]++;
  371.        
  372.  
  373.         // TH1F * hNbjet = new TH1F("hNbjet","Bjet multiplicity",15,0.5,15.5);
  374.         // TH1F * hHTnew = new TH1F("hHTnew","HT with cleaned jets",15,0.5,15.5);
  375.         // TH1F * hHT40 = new TH1F("hHT40","normal HT",15,0.5,15.5);
  376.         // hNbjet->Fill(Nbjet,EvWeight);
  377.         // hHTnew->Fill(HTnew,EvWeight);
  378.         // hHT40->Fill(HT,EvWeight);
  379.         // hNbjet->Write();
  380.         // hHTnew->Write();
  381.         // hHT40->Write();
  382.        
  383.        
  384.        
  385.        
  386.         //////////////////////
  387.         // Start of CutFlow //
  388.         //////////////////////
  389.            
  390.            // 1. nLeptons
  391.            if (Nmu != 3) continue;
  392.            if (Nel != 0) continue;
  393.            //if(Nlept != 3) continue;
  394.            if ((Lepton[0].Pt()<25.) || (Lepton[1].Pt()<15.) || (Lepton[2].Pt()<10.)) continue;
  395.            CFCounter[2]+= EvWeight;
  396.            iCFCounter[2]++;
  397.            CFCounter1[2]+= EvWeight;
  398.                    iCFCounter1[2]++;
  399.  
  400.            ///////////////////////   MT Definition
  401.  
  402.            // double MET_Phi = tree->Get(MET_Phi,"DelphMET_Phi");
  403.            // double MT = sqrt(2*MET*Lepton[2].Pt()*(1-cos(fabs(MET_Phi - Lepton[2].Phi()))));         
  404.            // hMTpre -> Fill(MT,EvWeight);
  405.  
  406.            //////////////////////   MT Definition
  407.  
  408.            // hasZ cut
  409.            if(lv_sum.M() < 15) continue;
  410.                    if (hasZ > 0) continue;
  411.            CFCounter[3]+= EvWeight;
  412.            iCFCounter[3]++;
  413.            CFCounter1[3]+= EvWeight;
  414.                    iCFCounter1[3]++;       
  415.  
  416.            // jet multiplicity
  417.                    //if(newJets.size() < 4) continue;
  418.            // CFCounter[4]+=EvWeight;
  419.            // iCFCounter[4]++;
  420.                    // CFCounter1[4]+=EvWeight;
  421.                    // iCFCounter1[4]++;
  422.  
  423.          
  424.            // b selection
  425.            if (Nbjet < 4) continue;
  426.            //if (BJetInd.size() < 4) continue;
  427.            // bool bjet = true;
  428.            // //for(int i=0; i< newJets.size(); i++)
  429.            //   if (JetB.size() > 0) bjet = false;
  430.            //   if(!bjet) continue;
  431.            CFCounter[4]+=EvWeight;         
  432.            iCFCounter[4]++;
  433.            CFCounter1[4]+=EvWeight;
  434.                    iCFCounter1[4]++;
  435.  
  436.            /////////// first chek
  437.            double dPhi01 = TMath::Min(JetMETdPhi[0],JetMETdPhi[1]);
  438.            double dPhi23 = TMath::Min(JetMETdPhi[2],JetMETdPhi[3]);
  439.            double dPhi   = TMath::Min(dPhi01,dPhi23);
  440.            hDphi1       ->  Fill(dPhi01, EvWeight);
  441.            if(Nel>0)hNel1     -> Fill(Nel, EvWeight);
  442.            hHT1      ->  Fill(HTnew,EvWeight);
  443.            hMET1     ->  Fill(MET,EvWeight);
  444.            hNjet1    ->  Fill(newJets.size(),EvWeight);
  445.            hNlep1    -> Fill(Nlept, EvWeight);
  446.            if(Nmu>0)hNmu1     -> Fill(Nmu, EvWeight);
  447.            if(m>0)hMultiLMassSF1->Fill(m,EvWeight);
  448.            hLep1Pt1 -> Fill(Lepton[0].Pt(), EvWeight);             
  449.            hLep2Pt1  -> Fill(Lepton[1].Pt(),EvWeight);
  450.            hLep3Pt1  -> Fill(Lepton[2].Pt(),EvWeight);
  451.            if(newJets.size()>2) hJet1Pt1  -> Fill(newJets[0].Pt(),EvWeight);
  452.            if(newJets.size()>2) hJet2Pt1  -> Fill(newJets[1].Pt(),EvWeight);
  453.            if(newJets.size()>2) hJet3Pt1  -> Fill(newJets[2].Pt(),EvWeight);
  454.            hDiLMass1 -> Fill(lv_sum.M(),EvWeight);
  455.            hMETMeff1 -> Fill(MET/Meff,EvWeight);
  456.            hMeff1    -> Fill(Meff,EvWeight);
  457.            if(Nmu>0)h1LepEta1 -> Fill(Muons[0].Eta(),EvWeight);
  458.            if(Nmu>0)h2LepEta1 -> Fill(Muons[1].Eta(),EvWeight);
  459.            if(Nmu>0)h3LepEta1 -> Fill(Muons[2].Eta(),EvWeight);
  460.            if(newJets.size()>2)h1JetEta1 -> Fill(newJets[0].Eta(),EvWeight);
  461.            if(newJets.size()>2)h2JetEta1 -> Fill(newJets[1].Eta(),EvWeight);
  462.            if(newJets.size()>2)h3JetEta1 -> Fill(newJets[2].Eta(),EvWeight);
  463.            //          if(Nmu>0)hMuIso1 -> Fill(MuIso[0], EvWeight);
  464.            hST1->Fill(ST,EvWeight);
  465.            //          hMT31->Fill(MT, EvWeight);  
  466.            ////////// end of the first check
  467.  
  468.  
  469.            //del phi cut > 0.3
  470.  
  471.            // if (dPhi01 < 0.3) continue;
  472.            // if (dPhi23 < 0.3) continue;
  473.            if (dPhi < 0.3) continue;
  474.            CFCounter[5]+=EvWeight;
  475.            iCFCounter[5]++;
  476.            CFCounter1[5]+=EvWeight;
  477.                    iCFCounter1[5]++;
  478.            
  479.  
  480.            // HT cut
  481.  
  482.            // if (HT < 500) continue;
  483.            // CFCounter[8]+=EvWeight;
  484.            // iCFCounter[8]++;
  485.            // CFCounter1[7]+=EvWeight;
  486.                    // iCFCounter1[7]++;
  487.            
  488.            //Medium MET-HT
  489.            if(MET < 250 /*|| HTnew < 750*/) continue;
  490.               CFCounter1[6]+=EvWeight;
  491.               iCFCounter1[6]++;
  492.               CFCounter[6]+=EvWeight;
  493.               iCFCounter[6]++;
  494.              
  495.  
  496.  
  497.            // MET-HT
  498.            if(MET>250 && HTnew>750) {
  499.              CFCounter[7]+=EvWeight;
  500.              iCFCounter[7]++;
  501.            }
  502.            //          CFCounter1[7]+=EvWeight;
  503.            //                   iCFCounter1[7]++;
  504.  
  505.            
  506.            //          double HT_cut = 500.;
  507.            //          double MET_cut = 250.;
  508.            //          for(HT_cut = 500.; HT_cut<1201 ; HT_cut+=50)
  509.            //            {
  510.            //              for(MET_cut = 250.; MET_cut<751 ; MET_cut+=50)
  511.            //            {
  512.            //              if(MET>MET_cut && HT>HT_cut)
  513.            //{
  514.            //      
  515.            //        }
  516.            //    }
  517.            //  }
  518.            
  519.            // MET cut          
  520.  
  521.            // if (MET < 500) continue;         //////////// <<<<<<<<<<<<
  522.            // CFCounter[11]+=EvWeight;         
  523.            // iCFCounter[11]++;
  524.            // CFCounter1[8]+=EvWeight;
  525.                    // iCFCounter1[8]++;
  526.                        
  527.            ////////////////////// second check
  528.            //          hMT32->Fill(MT, EvWeight);          
  529.            hDphi2       ->  Fill(dPhi01, EvWeight);
  530.                    if(Nel>0)hNel2     -> Fill(Nel, EvWeight);
  531.            hHT2      ->  Fill(HTnew,EvWeight);
  532.                    hMET2     ->  Fill(MET,EvWeight);
  533.                    hNjet2    ->  Fill(newJets.size(),EvWeight);
  534.                    hNlep2    -> Fill(Nlept, EvWeight);
  535.                    if(Nmu>0)hNmu2     -> Fill(Nmu, EvWeight);
  536.            if(m>0)hMultiLMassSF2->Fill(m,EvWeight);
  537.            hLep1Pt2 -> Fill(Lepton[0].Pt(), EvWeight);                  
  538.            hLep2Pt2  -> Fill(Lepton[1].Pt(),EvWeight);
  539.                    hLep3Pt2  -> Fill(Lepton[2].Pt(),EvWeight);
  540.                    if(newJets.size()>2)hJet1Pt2  -> Fill(newJets[0].Pt(),EvWeight);
  541.                    if(newJets.size()>2)hJet2Pt2  -> Fill(newJets[1].Pt(),EvWeight);
  542.                    if(newJets.size()>2)hJet3Pt2  -> Fill(newJets[2].Pt(),EvWeight);
  543.                    hDiLMass2 -> Fill(lv_sum.M(),EvWeight);
  544.            hMETMeff2 -> Fill(MET/Meff,EvWeight);
  545.                    hMeff2    -> Fill(Meff,EvWeight);
  546.            if(Nmu>0)h1LepEta2 -> Fill(Muons[0].Eta(),EvWeight);
  547.                    if(Nmu>0)h2LepEta2 -> Fill(Muons[1].Eta(),EvWeight);
  548.            if(Nmu>0)h3LepEta2 -> Fill(Muons[2].Eta(),EvWeight);
  549.                    if(newJets.size()>2)h1JetEta2 -> Fill(newJets[0].Eta(),EvWeight);
  550.                    if(newJets.size()>2)h2JetEta2 -> Fill(newJets[1].Eta(),EvWeight);
  551.                    if(newJets.size()>2)h3JetEta2 -> Fill(newJets[2].Eta(),EvWeight);
  552.            //          if(Nmu>0)hMuIso2 -> Fill(MuIso[0], EvWeight);
  553.            hST2->Fill(ST,EvWeight);
  554.            ///////////////////// end of the second check
  555.  
  556.  
  557.            //missing trans energy / effective mass
  558.            
  559.            // if (MET/Meff < 0.2) continue;
  560.            // CFCounter[10]+= EvWeight;
  561.            // iCFCounter[10]++;
  562.            // CFCounter1[8]+= EvWeight;
  563.                    // iCFCounter1[8]++;
  564.  
  565.  
  566.            // double LeadBPt = 0;
  567.            // for (int iJet = 0; iJet < BJetInd.size(); iJet++)
  568.            //   {
  569.            //     if (newJets[BJetInd[iJet]].Pt() > LeadBPt)
  570.            //    LeadBPt = newJets[BJetInd[iJet]].Pt();
  571.            //   }
  572.  
  573.            // double MWmass     = MET + Lepton[3].Pt();
  574.            // double MTopmass   = MET + Lepton[3].Pt() + LeadBPt;
  575.            // double MassDifMT  = MWmass - MT3;
  576.            // double MassDifMET = MET - MT3;
  577.  
  578.  
  579.  
  580.            //          hMT3        ->Fill(MT, EvWeight);  
  581.            //          hWMass      ->Fill(MWmass,EvWeight);
  582.            //          hTopMass    ->Fill(MTopmass,EvWeight);
  583.            //          hMassDifMT  ->Fill(MassDifMT/MT3,EvWeight);
  584.            hMETMeff    ->Fill(MET/Meff, EvWeight);
  585.            //          hMETLep     ->Fill(MWmass/ST,EvWeight);
  586.            //          hMassDifMET ->Fill(MassDifMET/MT3,EvWeight);
  587.            hNjet       ->Fill(newJets.size(),EvWeight);    
  588.            hBelowZ     ->Fill(lv_sum.M(),EvWeight);
  589.            hAboveZ     ->Fill(lv_sum.M(),EvWeight);
  590.            hDiLMass    ->Fill(lv_sum.M(),EvWeight);    
  591.            hHT         ->Fill(HTnew,EvWeight);
  592.            hDphi       ->Fill(dPhi01,EvWeight);
  593.            hMET        ->Fill(MET,EvWeight);
  594.            if(Nel>0)hNel        ->Fill(Nel,EvWeight);
  595.            if(Nmu>0)hNmu        ->Fill(Nmu,EvWeight);
  596.            if(Nmu>0 && Nel>0)hNlep       ->Fill(Nel + Nmu,EvWeight);
  597.            hLep1Pt     ->Fill(Lepton[0].Pt(),EvWeight);
  598.            hLep2Pt     ->Fill(Lepton[1].Pt(),EvWeight);
  599.            hLep3Pt     ->Fill(Lepton[2].Pt(),EvWeight);
  600.            if(newJets.size()>2)hJet1Pt     ->Fill(newJets[0].Pt(),EvWeight);
  601.            if(newJets.size()>2)hJet2Pt     ->Fill(newJets[1].Pt(),EvWeight);
  602.            if(newJets.size()>2)  hJet3Pt     ->Fill(newJets[2].Pt(),EvWeight);
  603.            hMeff       ->Fill(Meff,EvWeight);
  604.            if(Nmu>0)h1LepEta -> Fill(Muons[0].Eta(),EvWeight);
  605.                    if(Nmu>0)h2LepEta -> Fill(Muons[1].Eta(),EvWeight);
  606.                    if(Nmu>0)h3LepEta -> Fill(Muons[2].Eta(),EvWeight);
  607.                    if(newJets.size()>2)h1JetEta -> Fill(newJets[0].Eta(),EvWeight);
  608.                    if(newJets.size()>2)h2JetEta -> Fill(newJets[1].Eta(),EvWeight);
  609.                    if(newJets.size()>2)h3JetEta -> Fill(newJets[2].Eta(),EvWeight);
  610.            //          if(Nmu>0)hMuIso -> Fill(MuIso[0], EvWeight);
  611.            hST->Fill(ST,EvWeight);
  612.  
  613.            //  End of Selection
  614.            
  615.       }
  616.     // ^loop end^
  617.    
  618.     ofstream tfile;
  619.     tfile.open("sunum_140PU_3000fb/SummerStd_"+outname+".txt");
  620.     tfile << "########################################" << endl;
  621.     tfile << "Cut efficiency numbers:" << endl;
  622.     for(int ci = 0; ci < CutNumb; ci++)
  623.       {
  624.         tfile << "After cut " << CutList[ci] <<  "\t\t\t"
  625.           << CFCounter[ci]  << "\t events left\t"<< iCFCounter[ci] <<" cnt\t"<< endl;
  626.         CutFlow->SetBinContent(1+ci,CFCounter1[ci]);
  627.         CutFlow1->SetBinContent(1+ci,CFCounter[ci]);
  628.       }
  629.    
  630.     TFile * outf = new TFile("sunum_140PU_3000fb/SummerStd_"+outname+"_his.root","RECREATE");
  631.  
  632.  
  633.     CutFlow           ->Write();
  634.     CutFlow1          ->Write();
  635.     hBjet->Write();
  636.     hMuiso->Write();
  637.     hEliso->Write();
  638.     ////////////////// first check
  639.  
  640.     hNel1             ->Write();
  641.     hNmu1             ->Write();
  642.     hNlep1            ->Write();   
  643.     hNjet1            ->Write();
  644.     hMultiLMassSF1    ->Write();
  645.     hHT1              ->Write();
  646.     hLep1Pt1  -> Write();
  647.     hLep2Pt1  -> Write();
  648.     hLep3Pt1  -> Write();
  649.     hJet1Pt1  -> Write();
  650.     hJet2Pt1  -> Write();
  651.     hJet3Pt1  -> Write();
  652.     hDiLMass1 -> Write();
  653.     hMETMeff1 -> Write();
  654.     hMeff1    -> Write();
  655.     hMET1     -> Write();
  656.     hDphi1    -> Write();
  657.     h1LepEta1 -> Write();  
  658.     h2LepEta1 -> Write();
  659.     h3LepEta1 -> Write();
  660.     h1JetEta1 -> Write();
  661.     h2JetEta1 -> Write();
  662.     h3JetEta1 -> Write();  
  663.     //  hMuIso1 -> Write();
  664.     hST1-> Write();
  665.     ///////////////// second check
  666.    
  667.     hDphi2    -> Write();
  668.     hNel2     -> Write();
  669.     hHT2      -> Write();
  670.     hMET2     -> Write();
  671.     hNjet2    -> Write();
  672.     hNlep2    -> Write();
  673.     hNmu2     -> Write();
  674.     hLep1Pt2  -> Write();
  675.     hMultiLMassSF2-> Write();
  676.     hLep2Pt2  -> Write();
  677.     hLep3Pt2  -> Write();
  678.     hJet1Pt2  -> Write();
  679.     hJet2Pt2  -> Write();
  680.     hJet3Pt2  -> Write();
  681.     hDiLMass2 -> Write();
  682.     hMETMeff2 -> Write();
  683.     hMeff2    -> Write();
  684.         hDiLMass2 -> Write();
  685.     h1LepEta2 -> Write();
  686.         h2LepEta2 -> Write();
  687.         h3LepEta2 -> Write();
  688.         h1JetEta2 -> Write();
  689.         h2JetEta2 -> Write();
  690.         h3JetEta2 -> Write();
  691.     //  hMuIso2 -> Write();
  692.     hST2-> Write();
  693.     /////////////// full requirements
  694.  
  695.     hNel              ->Write();
  696.         hNmu              ->Write();
  697.         hNlep             ->Write();
  698.         hNjet             ->Write();
  699.         hLep1Pt   -> Write();
  700.         hLep2Pt   -> Write();
  701.         hLep3Pt   -> Write();
  702.         hJet1Pt   -> Write();
  703.         hJet2Pt   -> Write();
  704.         hJet3Pt   -> Write();
  705.         hMultiLMassSF12   ->Write();
  706.         hDiLMass          ->Write();
  707.         hBelowZ           ->Write();
  708.         hAboveZ           ->Write();
  709.         hHT               ->Write();
  710.     hMET              ->Write();
  711.     hDphi             ->Write();
  712.     //  hMassDifMT        ->Write();
  713.     //  hMassDifMET       ->Write();
  714.     hMETLep           ->Write();
  715.     hMETMeff          ->Write();
  716.     hMeff             ->Write();
  717.     h1LepEta -> Write();
  718.         h2LepEta -> Write();
  719.         h3LepEta -> Write();
  720.         h1JetEta -> Write();
  721.         h2JetEta -> Write();
  722.         h3JetEta -> Write();
  723.     //  hMuIso -> Write();
  724.     hST-> Write();
  725.     //  hBjet->Write();
  726.     //  hTopMass          ->Write();
  727.     //  hWMass            ->Write();
  728.     // hMT3              ->Write();
  729.     // hMT31             ->Write();
  730.     // hMT32 ->Write();
  731.     // hMTpre->Write();
  732.  
  733.  
  734.  
  735. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!