Advertisement
Guest User

C++ Source HEP Analysis2

a guest
Aug 29th, 2014
225
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 28.05 KB | None | 0 0
  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