Changeset 22 in huonglan

Timestamp:

Sep 2, 2011, 3:07:15 PM (13 years ago)

Author:

huonglan

Message:

sPlots_NewRel improved

Location:

SPLOT/NEWREALEASE

Files:

• EFFICIENCY/CatEPS.tar (deleted)
• EFFICIENCY/EfficiencyDependOnCut.C (modified) (17 diffs)
• EFFICIENCY/EfficiencyDependOnCut_C.d (deleted)
• MERGEFILES/MakeFakeData.C (modified) (18 diffs)
• sPlots_NewRel.C (modified) (38 diffs)

SPLOT/NEWREALEASE/EFFICIENCY/EfficiencyDependOnCut.C

@@ -10 +10 @@
 #include <TLegend.h>
 #include <TROOT.h>
+#include <TMinuit.h>
 
 #include <iostream>
@@ -15 +16 @@
 using namespace std;
 
+//Set the cut values
+double Cut1 = 0.7;
+double Cut2 = 0.7;
+double Cut3 = 48;
+double Cut4 = 250;
+double Cut5 = 620;
+
+//Define class Event
+class EventC {
+public:
+  EventC() { m_w = 1; m_fs = 1; m_fb = 0; m_cat = 1; }
+  EventC(double w0, double fs0, double fb0) { m_w = w0; m_fs = fs0; m_fb = fb0; }
+  EventC(double w0, int cat0) { m_w = w0; m_cat = cat0; m_fs = 1; m_fb = 0; }
+  double w() { return m_w; }
+  double fs() { return m_fs; }
+  double fb() { return m_fb; }
+  int cat() { return m_cat; }
+  void fs(double fs0) { m_fs = fs0; }
+  void fb(double fb0) { m_fb = fb0; }
+
+private:
+  double m_w, m_fs, m_fb;
+  int m_cat;
+};
+
+//global variables
+std::vector<EventC> m_EventList;
+double m_F1[3][31][6];
+double m_F0[3][31][6];
+double m_N1exp[3][31];
+double m_N0exp[3][31];
+
+void fcnN0N1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
+  double L = 0;
+  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
+    EventC evt = m_EventList[iev];
+    double smt = par[0]*evt.fb() + par[1]*evt.fs();
+    if (smt > 0) L += evt.w()*log(smt);
+  }
+  f = par[0] + par[1] - L; 
+};
+
+vector<double> vFe(double cos1,double cos2,double dMLep, int i, int j)
+{
+  double f0e = 0;
+  double f1e = 0;
+  
+  //cout << "cos1 = " << cos1 << " cos2 = " << cos2 << " dMLep = " << dMLep << endl;
+  
+  int it = -1;
+  if (cos1<Cut1 && cos2<Cut2){
+    if (dMLep<Cut3) it = 0;
+    if (dMLep>Cut3) it = 1;
+  }
+  else if (cos1>Cut1 && cos2>Cut2){
+    if (dMLep<Cut3) it = 2;
+    if (dMLep>Cut3) it = 3;
+  }
+  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
+    if (dMLep<Cut3) it = 4;
+    if (dMLep>Cut3) it = 5;
+  }
+  
+  f0e = m_F0[i][j][it];
+  f1e = m_F1[i][j][it];
+  
+  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
+  
+  vector<double> vfe;
+  vfe.push_back(f0e);
+  vfe.push_back(f1e);
+  return vfe;
+}
+
+//Make list of events in Data file
+void MakeEventListf0f1(int i, int j){
+  m_EventList.clear();  
+
+  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
+  TTree *tD = (TTree*)fD->Get("T"); 
+  int nD = tD->GetEntries();
+  for (int iev = 0; iev < nD; iev++) {
+    tD->GetEntry(iev);
+    
+    double mt    = tD->GetLeaf("mtophad")->GetValue();
+    double HT    = tD->GetLeaf("HT")->GetValue();
+    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
+    double y     = tD->GetLeaf("cosWqq")->GetValue();
+    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
+
+    int hasBJ = (int)tD->GetLeaf("hasBadJet")->GetValue();
+    int RunN  = (int)tD->GetLeaf("RunNumber_ana")->GetValue();
+    
+    if (RunN >= 180614 && hasBJ==1) continue;
+
+    double cos1 = fabs(x);
+    double cos2 = fabs(y);
+    
+    //sPlots total
+    if (mt > Cut4 && HT > Cut5){
+      vector<double> fsfb =  vFe(cos1,cos2,dMLep,i,j);
+      double f0 = fsfb.at(0);
+      double f1 = fsfb.at(1);
+      //cout << "f0 = " << f0 << " f1 = " << f1 << endl;
+      EventC eventC(1,f1,f0);
+      m_EventList.push_back(eventC);
+    } 
+  }
+}
+
+//Fit the number N0 and N1 using the events from data files
+std::pair<pair<double,double>,double> Fit(double N0_exp, double N1_exp, int i, int j){
+  // Test
+  MakeEventListf0f1(i,j);
+  cout << "Number of events in the LH " << m_EventList.size() << endl;
+
+  int ifl = 0;
+  double step = 0.01;
+  
+  double fitPar0 = N0_exp;
+  double nmax0   = fitPar0 + 4*sqrt(N0_exp);
+  double nmin0   = fitPar0 - 4*sqrt(N0_exp);
+  
+  double fitPar1 = N1_exp;
+  double nmax1   = fitPar1 + 4*sqrt(N1_exp);
+  double nmin1   =         - 4*sqrt(N1_exp);
+  
+  double n0best, err0best;
+  double n1best, err1best;
+  
+  TMinuit* myMinuit = new TMinuit(2);
+  //myMinuit->Command("SET PRI -1");
+  myMinuit->Command("SET NOW");
+  myMinuit->Command("SET ERR 0.5");
+  myMinuit->Command("SET STR 2");   // try to improve mimimum search
+  myMinuit->SetFCN(fcnN0N1);
+  myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl); 
+  myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl); 
+  //myMinuit->FixParameter(0);
+  myMinuit->Command("MIGRAD 500 1");
+  myMinuit->GetParameter(0,n0best,err0best); 
+  myMinuit->GetParameter(1,n1best,err1best); 
+  
+  std::cout << "Expected N0 = " << N0_exp << " fitted = " << n0best << " +- " << err0best << std::endl;
+  std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
+  
+  double N0 = n0best;
+  double N1 = n1best;
+  double sig2 = err1best*err1best;
+
+  pair<double,double> N0N1 = make_pair<double,double> (N0,N1);
+  
+  return make_pair<pair<double,double>,double> (N0N1,sig2);
+};
+
+
+//##############################################################
+//####################### MAIN PROGRAM #########################
+//##############################################################
+
 void EfficiencyDependOnCut(){
-  TTree *t[3];
+  TTree *t[4];
 
   TFile *fB = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_AllBkg.root");
@@ -24 +185 @@
   t[1] = (TTree*)fS->Get("T");
 
+  TFile *fBS = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_FakeData.root");
+  t[2] = (TTree*)fBS->Get("T");
+
   TFile *fD = TFile::Open("/exp/atlas/huonglan/DATA2011/resData.root");
-  t[2] = (TTree*)fD->Get("T");
-
-  double Cut1 = 0.7;
-  double Cut2 = 0.7;
-  double Cut3 = 48;
-
-  double   e1[3][3][31];
-  double   e2[3][3][31];   
-  double   e3[3][3][31];   
-  double   e4[3][3][31];   
-  double   e5[3][3][31];   
-  double   e6[3][3][31];   
-  double iTot[3][3][31];
+  t[3] = (TTree*)fD->Get("T");
+
+  double   e1[3][4][31];
+  double   e2[3][4][31];   
+  double   e3[3][4][31];   
+  double   e4[3][4][31];   
+  double   e5[3][4][31];   
+  double   e6[3][4][31];   
+  double iTot[3][4][31];
 
   for (int k = 0; k < 3; k++){
-    for (int i = 0; i < 3; i++){
+    for (int i = 0; i < 4; i++){
       for (int j = 0; j < 31; j++){
         iTot[k][i][j] =0 ;
@@ -56 +216 @@
   double mjjjErr[31]; for (int i = 0; i < 31; i++) { mjjjErr[i] = 0;}
 
-  for (int i = 0; i < 3; i++){
+  for (int i = 0; i < 4; i++){
     int nevt = t[i]->GetEntries();
     for (int iev = 0; iev < nevt; iev++) {
@@ -64 +224 @@
       double wei = 1;
       double cos1 = 1, cos2 = 1;
+      int hfor = -1;
 
       double mt   = t[i]->GetLeaf("mtophad")->GetValue();
       double ht   = t[i]->GetLeaf("HT")->GetValue();
-      if (i==0||i==1) {
+      int hasBJ   = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      if (i==0||i==1||i==2) {
         wei  = t[i]->GetLeaf("wei")->GetValue();
         cos1 = fabs(t[i]->GetLeaf("cos1")->GetValue());
         cos2 = fabs(t[i]->GetLeaf("cos2")->GetValue());
-      } else if (i==2){
+        hfor = (int)t[i]->GetLeaf("hfor")->GetValue();
+      } else if (i==3){
         cos1 = fabs(t[i]->GetLeaf("cos_tlepCMtpair")->GetValue());
         cos2 = fabs(t[i]->GetLeaf("cosWqq")->GetValue());
       }
       double dM   = t[i]->GetLeaf("MinLep1")->GetValue();
+
+      if (hasBJ!=0 || hfor==4) continue;
 
       int it = -1;
@@ -110 +275 @@
   }
 
-  double    eB[3][31][6],    eS[3][31][6],    eD[3][31][6]; 
-  double eBErr[3][31][6], eSErr[3][31][6], eDErr[3][31][6]; 
+  double    eB[3][31][6],    eS[3][31][6],    eBS[3][31][6],    eD[3][31][6]; 
+  double eBErr[3][31][6], eSErr[3][31][6], eBSErr[3][31][6], eDErr[3][31][6]; 
 
   for (int i = 0; i < 3; i++){
@@ -136 +301 @@
       if (iTot[i][1][j]==0) { for (int k = 0; k < 6; k++) { eS[i][j][k] = 0; } }
       
-      eD[i][j][0] = e1[i][2][j]/iTot[i][2][j]; 
-      eD[i][j][1] = e2[i][2][j]/iTot[i][2][j]; 
-      eD[i][j][2] = e3[i][2][j]/iTot[i][2][j]; 
-      eD[i][j][3] = e4[i][2][j]/iTot[i][2][j]; 
-      eD[i][j][4] = e5[i][2][j]/iTot[i][2][j]; 
-      eD[i][j][5] = e6[i][2][j]/iTot[i][2][j]; 
-
-      if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
+      eBS[i][j][0] = e1[i][2][j]/iTot[i][2][j];
+      eBS[i][j][1] = e2[i][2][j]/iTot[i][2][j];
+      eBS[i][j][2] = e3[i][2][j]/iTot[i][2][j];
+      eBS[i][j][3] = e4[i][2][j]/iTot[i][2][j];
+      eBS[i][j][4] = e5[i][2][j]/iTot[i][2][j];
+      eBS[i][j][5] = e6[i][2][j]/iTot[i][2][j];
+
+      if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eBS[i][j][k] = 0; } }
+      
+      eD[i][j][0] = e1[i][3][j]/iTot[i][3][j]; 
+      eD[i][j][1] = e2[i][3][j]/iTot[i][3][j]; 
+      eD[i][j][2] = e3[i][3][j]/iTot[i][3][j]; 
+      eD[i][j][3] = e4[i][3][j]/iTot[i][3][j]; 
+      eD[i][j][4] = e5[i][3][j]/iTot[i][3][j]; 
+      eD[i][j][5] = e6[i][3][j]/iTot[i][3][j]; 
+
+      if (iTot[i][3][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
       
       //cout << "e1S[" << i << "][" << j << "] = " << e1S[i][j];
@@ -150 +324 @@
   }
   
+  //3 is the number of regions
+  //31 is the number of the cut for mjjj
+  //6 is the number of categories
   for (int i = 0; i < 3; i++){
     for (int j = 0; j < 31; j++){
       for (int k = 0; k < 6; k++){
-        eBErr[i][j][k] = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
-        eSErr[i][j][k] = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
-        eDErr[i][j][k] = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][2][j]);
-      }
+        eBErr[i][j][k]  = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
+        eSErr[i][j][k]  = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
+        eBSErr[i][j][k] = sqrt(eBS[i][j][k]*(1-eBS[i][j][k])/iTot[i][2][j]);
+        eDErr[i][j][k]  = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][3][j]);
+
+        m_F1[i][j][k] = eS[i][j][k];
+        m_F0[i][j][k] = eB[i][j][k];
+        m_N1exp[i][j] = iTot[i][1][j];
+        m_N0exp[i][j] = iTot[i][0][j];
+
+        cout << "region " << i << " cut mjjj " << mjjj[j] << " cate " << k+1 ;
+        cout << " iTotB = " << iTot[i][0][j] ;
+        cout << " iTotS = " << iTot[i][1][j] << endl;
+      }
+    }
+  }
+
+  //3 is the number of regions
+  //31 is the number of the cut for mjjj
+  double N0fit[3][31], N1fit[3][31], N1err[3][31];
+  for (int i = 0; i < 3; i++){
+    for (int j = 0; j < 31; j++){
+      pair<pair<double,double>,double> N0N1sig2 = Fit(m_N0exp[i][j],m_N1exp[i][j],i,j);
+      pair<double,double> N0N1 = N0N1sig2.first;
+      double N0   = N0N1.first; 
+      double N1   = N0N1.second; 
+      double sig2 = N0N1sig2.second;
+      N0fit[i][j] = N0;
+      N1fit[i][j] = N1;
+      N1err[i][j] = sig2;
     }
   }
@@ -175 +378 @@
   TString cat[6] = {"cat1","cat2","cat3","cat4","cat5","cat6"};
 
-  TGraphErrors *heffB[3][6];
-  TGraphErrors *heffS[3][6];
-  TGraphErrors *heffD[3][6];
-
+//   TGraphErrors *heffB[3][6];
+//   TGraphErrors *heffS[3][6];
+//   TGraphErrors *heffBS[3][6];
+//   TGraphErrors *heffD[3][6];
+
+  TGraph *heffB[3][6];
+  TGraph *heffS[3][6];
+  TGraph *heffBS[3][6];
+  TGraph *heffD[3][6];
+
+  //*************Define histograms on N1************
+  TGraphErrors *hN1[3];
+  for (int i = 0; i < 3; i++){
+    double N1fittemp[31], N1errtemp[31];
+    for (int j = 0; j < 31; j++){ 
+      N1fittemp[j] = N1fit[i][j];
+      N1errtemp[j] = N1err[i][j]; 
+      cout << "mjjj = " << mjjj[j] << endl;
+      cout << " m_N1exp[" << i << "][" << j << "]= " << m_N1exp[i][j] ;
+      cout << " m_N0exp[" << i << "][" << j << "]= " << m_N0exp[i][j] ;
+      cout << " N1fit[" << i << "][" << j << "]= " << N1fit[i][j] ;
+      cout << " N0fit[" << i << "][" << j << "]= " << N0fit[i][j] ;
+      cout << " N1err[" << i << "][" << j << "]= " << N1err[i][j] << endl;
+    }
+    hN1[i] = new TGraphErrors(31,mjjj,N1fittemp,mjjjErr,N1errtemp);
+  }
+
+  //*************Define histograms on efficiencies****************
   for (int i = 0; i < 3; i++){
     for (int j = 0; j < 6; j++){
@@ -189 +416 @@
         effErr[ic] = eBErr[i][ic][j];
       }
-      heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
-      heffB[i][j]->SetLineColor(4);
+      //heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      heffB[i][j] = new TGraph(31, mjjj, eff);
+      heffB[i][j]->SetLineColor(8);
+      heffB[i][j]->SetLineWidth(2);
       
       name = "heffS"; name += region[i]; name += cat[j];
@@ -197 +426 @@
         effErr[ic] = eSErr[i][ic][j];
       }
-      heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      //heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      heffS[i][j] = new TGraph(31, mjjj, eff);
       heffS[i][j]->SetLineColor(2);
+      heffS[i][j]->SetLineWidth(2);
+      
+      name = "heffBS"; name += region[i]; name += cat[j];
+      for (int ic = 0; ic < 31; ic++) { 
+        eff[ic]    = eBS[i][ic][j];
+        effErr[ic] = eBSErr[i][ic][j];
+      }
+      //heffBS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      heffBS[i][j] = new TGraph(31, mjjj, eff);
+      heffBS[i][j]->SetLineColor(kMagenta);
+      heffBS[i][j]->SetLineWidth(2);
       
       name = "heffD"; name += region[i]; name += cat[j];
@@ -205 +446 @@
         effErr[ic] = eDErr[i][ic][j];
       }
-      heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      //heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
+      heffD[i][j] = new TGraph(31, mjjj, eff);
+      heffD[i][j]->SetLineWidth(2);
     }
   }
@@ -211 +454 @@
   legend->AddEntry(heffB[0][1],"BkgMC");
   legend->AddEntry(heffS[0][1],"SigMC");
+  legend->AddEntry(heffBS[0][1],"Sig+BkgMC");
   legend->AddEntry(heffD[0][1],"Data");
+
+  //****************** DRAWING *********************
+
+  TCanvas *cN1 = new TCanvas("cN1", "cN1", 1000, 400);
+  cN1->Divide(3);
+  
+  cN1->cd(1);
+  hN1[0]->Draw("ALP");
+  cN1->cd(2);
+  hN1[1]->Draw("ALP");
+  cN1->cd(3);
+  hN1[2]->Draw("ALP");
+
+  cN1->Print("N1fit.eps");
 
   TCanvas *c[6];
@@ -223 +481 @@
     frame1->SetTitle(name);
     frame1->Draw();
-    heffB[0][i]->Draw("LP");
-    heffS[0][i]->Draw();
+    heffS[0][i]->Draw("LP");
     heffD[0][i]->Draw();
+    heffBS[0][i]->Draw();
+    heffB[0][i]->Draw();
     legend->Draw("same");
 
@@ -232 +491 @@
     frame2->SetTitle(name);
     frame2->Draw();
-    heffB[1][i]->Draw("LP");
-    heffS[1][i]->Draw();
+    heffS[1][i]->Draw("LP");
     heffD[1][i]->Draw();
+    heffBS[1][i]->Draw();
+    heffB[1][i]->Draw();
     legend->Draw("same");
 
@@ -241 +501 @@
     frame3->SetTitle(name);
     frame3->Draw();
-    heffB[2][i]->Draw("LP");
-    heffS[2][i]->Draw();
+    heffS[2][i]->Draw("LP");
     heffD[2][i]->Draw();
+    heffBS[2][i]->Draw();
+    heffB[2][i]->Draw();
     legend->Draw("same");
 
@@ -250 +511 @@
   }
 
-  //********************************************
-
-  //TCanvas *cN = new TCanvas("cN", "", 500, 400);
-  //TGraph *hB = new TGraph(31, mjjj, N1);
-  //hB->Draw("ALP");
 }

SPLOT/NEWREALEASE/MERGEFILES/MakeFakeData.C

@@ -39 +39 @@
 
   TString s3 = "v3.root";
-  TString ver = "ToTestEfficiency";
+  TString ver = "v3";
 
   //set the name of 78 files
@@ -144 +144 @@
   //###################################
 
-  TFile *fData    = new TFile("res_Data"+ver+".root","RECREATE");
+  TFile *fData   = new TFile("res_FakeData"+ver+".root","RECREATE");
   TTree *tWeiD   = new TTree("T","T");
-  float b_weiD  = -100;
-  float b_mtophad= -100;    
-  float b_HT= -100;         
-  float b_cos1= -100;       
-  float b_cos2= -100;       
-  float b_MinLep1= -100;    
-  short b_hfor= -100;       
-  short b_hasBadJet= -100;  
-  short b_isEleEvent= -100; 
-  int b_RunNumber= -100;  
-  tWeiD->Branch("wei",        &b_weiD,               "wei/F");
+  float b_wei        = -100;
+  float b_mtophad    = -100;    
+  float b_HT         = -100;         
+  float b_cos1       = -100;       
+  float b_cos2       = -100;       
+  float b_MinLep1    = -100;    
+  short b_hfor       = -100;       
+  short b_hasBadJet  = -100;  
+  short b_isEleEvent = -100; 
+  int b_RunNumber    = -100;  
+  tWeiD->Branch("wei",        &b_wei,                "wei/F");
   tWeiD->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiD->Branch("HT",         &b_HT,                  "HT/F");
@@ -172 +172 @@
     countD++;
     tData->GetEntry(ifi);
-    b_mtophad = tSig->GetLeaf("mtophad")->GetValue();
-    b_HT = tSig->GetLeaf("HT")->GetValue();
+    b_mtophad = tData->GetLeaf("mtophad")->GetValue();
+    b_HT      = tData->GetLeaf("HT")->GetValue();
 
     //if (b_mtophad<250 || b_HT<620) continue; 
-    if (ifi%1000==0) cout << "Reading event " << ifi << " from Signal tree" << endl;
-
-    //TString currentFile = tSig->GetTree()->GetCurrentFile()->GetName();
-
-    b_weiS = Wei[0]*tSig->GetLeaf("WeiEvent")->GetValue();
-    b_cos1 = tSig->GetLeaf("cos_tlepCMtpair")->GetValue();
-    b_cos2 = tSig->GetLeaf("cosWqq")->GetValue();
-    b_MinLep1 = tSig->GetLeaf("MinLep1")->GetValue();           
-    b_hfor = (short)tSig->GetLeaf("hfor")->GetValue();
-    b_isEleEvent = (short)tSig->GetLeaf("isEleEvent")->GetValue();
-    b_hasBadJet = (short)tSig->GetLeaf("hasBadJet")->GetValue();
-    b_RunNumber = (int)tSig->GetLeaf("RunNumber_ana")->GetValue();
-
-    tWeiS->Fill();
-  }
-
-  fSig->cd();
-  tWeiS->Write();
-  fSig->Close();
-
+    if (ifi%1000==0) cout << "Reading event " << ifi << " from Data tree" << endl;
+
+    TString currentFile = tData->GetTree()->GetCurrentFile()->GetName();
+
+    b_hfor       = (short)tData->GetLeaf("hfor")->GetValue();
+    for (int i = 0; i < 78; i++){
+      double kfac = 1, ksE = 1, ksM = 1;
+      if (i>=2 && i<=32)  {kfac = 1.2; ksE = 0.906; ksM = 0.814;}//k factor for W+jets
+      if (i>=40 && i<=69) {kfac = 1.25;}//Zjets
+
+      if (fname[i]==currentFile)
+        b_wei = Wei[i]*(tData->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
+    }
+    b_cos1       = tData->GetLeaf("cos_tlepCMtpair")->GetValue();
+    b_cos2       = tData->GetLeaf("cosWqq")->GetValue();
+    b_MinLep1    = tData->GetLeaf("MinLep1")->GetValue();           
+    b_hfor       = (short)tData->GetLeaf("hfor")->GetValue();
+    b_isEleEvent = (short)tData->GetLeaf("isEleEvent")->GetValue();
+    b_hasBadJet  = (short)tData->GetLeaf("hasBadJet")->GetValue();
+    b_RunNumber  = (int)tData->GetLeaf("RunNumber_ana")->GetValue();
+
+    tWeiD->Fill();
+  }
+
+  fData->cd();
+  tWeiD->Write();
+  fData->Close();
 
   //###################################
@@ -201 +208 @@
   //###################################
 
-  TFile *fSig    = new TFile("res_Signal"+ver+".root","RECREATE");
+  TFile *fSig    = new TFile("res_Signal"+s2[0]+ver+".root","RECREATE");
+
   TTree *tWeiS   = new TTree("T","T");
-  float b_weiS  = -100;
-  float b_mtophad= -100;    
-  float b_HT= -100;         
-  float b_cos1= -100;       
-  float b_cos2= -100;       
-  float b_MinLep1= -100;    
-  short b_hfor= -100;       
-  short b_hasBadJet= -100;  
-  short b_isEleEvent= -100; 
-  int b_RunNumber= -100;  
-  tWeiS->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiS->Branch("wei",        &b_wei,                "wei/F");
   tWeiS->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiS->Branch("HT",         &b_HT,                  "HT/F");
@@ -230 +228 @@
     tSig->GetEntry(ifi);
     b_mtophad = tSig->GetLeaf("mtophad")->GetValue();
-    b_HT = tSig->GetLeaf("HT")->GetValue();
+    b_HT      = tSig->GetLeaf("HT")->GetValue();
 
     //if (b_mtophad<250 || b_HT<620) continue; 
@@ -237 +235 @@
     //TString currentFile = tSig->GetTree()->GetCurrentFile()->GetName();
 
-    b_weiS = Wei[0]*tSig->GetLeaf("WeiEvent")->GetValue();
-    b_cos1 = tSig->GetLeaf("cos_tlepCMtpair")->GetValue();
-    b_cos2 = tSig->GetLeaf("cosWqq")->GetValue();
-    b_MinLep1 = tSig->GetLeaf("MinLep1")->GetValue();           
-    b_hfor = (short)tSig->GetLeaf("hfor")->GetValue();
+    b_wei        = Wei[0]*tSig->GetLeaf("WeiEvent")->GetValue();
+    b_cos1       = tSig->GetLeaf("cos_tlepCMtpair")->GetValue();
+    b_cos2       = tSig->GetLeaf("cosWqq")->GetValue();
+    b_MinLep1    = tSig->GetLeaf("MinLep1")->GetValue();           
+    b_hfor       = (short)tSig->GetLeaf("hfor")->GetValue();
     b_isEleEvent = (short)tSig->GetLeaf("isEleEvent")->GetValue();
-    b_hasBadJet = (short)tSig->GetLeaf("hasBadJet")->GetValue();
-    b_RunNumber = (int)tSig->GetLeaf("RunNumber_ana")->GetValue();
+    b_hasBadJet  = (short)tSig->GetLeaf("hasBadJet")->GetValue();
+    b_RunNumber  = (int)tSig->GetLeaf("RunNumber_ana")->GetValue();
 
     tWeiS->Fill();
@@ -271 +269 @@
 
   TTree *tWeiB   = new TTree("T","T");
-  tWeiB->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiB->Branch("wei",        &b_wei,               "wei/F");
   tWeiB->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiB->Branch("HT",         &b_HT,                  "HT/F");
@@ -303 +301 @@
 
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tAllBkg->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
+        b_wei = Wei[i]*(tAllBkg->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
     }
     b_cos1       = tAllBkg->GetLeaf("cos_tlepCMtpair")->GetValue();
@@ -337 +335 @@
 
   TTree *tWeiT   = new TTree("T","T");
-  tWeiT->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiT->Branch("wei",        &b_wei,               "wei/F");
   tWeiT->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiT->Branch("HT",         &b_HT,                  "HT/F");
@@ -364 +362 @@
     for (int i = 0; i < 78; i++){
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tTop->GetLeaf("WeiEvent")->GetValue());
+        b_wei = Wei[i]*(tTop->GetLeaf("WeiEvent")->GetValue());
     }
     b_cos1       = tTop->GetLeaf("cos_tlepCMtpair")->GetValue();
@@ -386 +384 @@
 
   TTree *tWeiW   = new TTree("T","T");
-  tWeiW->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiW->Branch("wei",        &b_wei,               "wei/F");
   tWeiW->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiW->Branch("HT",         &b_HT,                  "HT/F");
@@ -417 +415 @@
 
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tAllW->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
+        b_wei = Wei[i]*(tAllW->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
     }
     b_cos1       = tAllW->GetLeaf("cos_tlepCMtpair")->GetValue();
@@ -435 +433 @@
   //*********** File Single Top  *****************
   //##############################################
-  TFile *fSingT   = new TFile("res_SingleTop"+ver+".root","RECREATE");
+  TFile *fSingT   = new TFile("res_SingT"+ver+".root","RECREATE");
 
   TTree *tWeiST   = new TTree("T","T");
-  tWeiST->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiST->Branch("wei",        &b_wei,               "wei/F");
   tWeiST->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiST->Branch("HT",         &b_HT,                  "HT/F");
@@ -465 +463 @@
     for (int i = 0; i < 78; i++){
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tSingT->GetLeaf("WeiEvent")->GetValue());
+        b_wei = Wei[i]*(tSingT->GetLeaf("WeiEvent")->GetValue());
     }
     b_hfor       = (short)tSingT->GetLeaf("hfor")->GetValue();
@@ -487 +485 @@
 
   TTree *tWeiZ   = new TTree("T","T");
-  tWeiZ->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiZ->Branch("wei",        &b_wei,               "wei/F");
   tWeiZ->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiZ->Branch("HT",         &b_HT,                  "HT/F");
@@ -518 +516 @@
 
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tAllZ->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
+        b_wei = Wei[i]*(tAllZ->GetLeaf("WeiEvent")->GetValue())*WJet_SF(b_hfor)*kfac*ksE*ksM;
     }
     b_cos1       = tAllZ->GetLeaf("cos_tlepCMtpair")->GetValue();
@@ -539 +537 @@
 
   TTree *tWeiQCD   = new TTree("T","T");
-  tWeiQCD->Branch("wei",        &b_weiS,               "wei/F");
+  tWeiQCD->Branch("wei",        &b_wei,               "wei/F");
   tWeiQCD->Branch("mtophad",    &b_mtophad,        "mtophad/F");
   tWeiQCD->Branch("HT",         &b_HT,                  "HT/F");
@@ -566 +564 @@
     for (int i = 0; i < 78; i++){
       if (fname[i]==currentFile)
-        b_weiS = Wei[i]*(tQCD->GetLeaf("WeiEvent")->GetValue());
+        b_wei = Wei[i]*(tQCD->GetLeaf("WeiEvent")->GetValue());
     }
     b_cos1       = tQCD->GetLeaf("cos_tlepCMtpair")->GetValue();

SPLOT/NEWREALEASE/sPlots_NewRel.C

@@ -12 +12 @@
 #include <vector>
 
+#include <iostream>
+
+using namespace std;
+
+void sPlots_NewRel_N0known();
+
 //Set the cut values
 double Cut1 = 0.7;
@@ -18 +24 @@
 double Cut4 = 250;
 double Cut5 = 620;
-double sigmaMin;
-int icount = 0;
-double iTot[2];
-double iCatS[6], iCatB[6];
-double F1[6];
-double F0[6];
-
-double N1 = 0., N0 = 0., sig2 = 0.;
+
+//global variables 
+int    m_icount = 0;
+double m_iTot[2];              //for (int i = 0; i < 2; i++){ m_iTot[i] = 0;}
+double m_iCatS[6], m_iCatB[6]; //for (int i = 0; i < 6; i++){ m_iCatS[i] = 0; m_iCatB[i] = 0;}
+double m_F1[6];                //for (int i = 0; i < 6; i++){ m_F1[i] = 0; }
+double m_F0[6];                //for (int i = 0; i < 6; i++){ m_F0[i] = 0; }
+double m_N1 = 0., m_N0 = 0., m_sig2 = 0.;
+
+TH1F *h[7];
 
 //Define class Event
@@ -46 +54 @@
 
 std::vector<EventC> m_EventList;
-
-void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
+bool m_ListAlreadyDone = false;
+
+void fcnN1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
   double L = 0;
-  for (int iev = 0; iev < m_EventList.size(); iev++) {
+  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
+    EventC evt = m_EventList[iev];
+    double smt = m_N0*evt.fb() + par[0]*evt.fs();
+    if (smt > 0) L += evt.w()*log(smt);
+  }
+  f = par[0] - L; 
+};
+
+void fcnN0N1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
+  double L = 0;
+  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
     EventC evt = m_EventList[iev];
     double smt = par[0]*evt.fb() + par[1]*evt.fs();
@@ -56 +75 @@
   f = par[0] + par[1] - L; 
 };
-
 
 vector<double> InvMatrix(double a, double b, double c, double d){
@@ -78 +96 @@
   //######## COMPUTE the sWeights
   TString filename = "";
-  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
-  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
+  TString type[2] = {"AllBkgTopMCAtNLO","Signal115420"};
   TFile *file[2];
   TTree *t[2];
   
   for (int i = 0; i < 2; i++){
-    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
+    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
     file[i] = TFile::Open(filename);
     t[i]    = (TTree*)file[i]->Get("T");
   }
 
-  for (int i = 0; i < 2; i++){
-    iTot[i] = 0;
-  }
-  for (int i = 0; i < 6; i++){
-    iCatS[i] = 0;
-    iCatB[i] = 0;
-  }
-  
   for (int i = 0; i < 2; i++){
     int nevt = t[i]->GetEntries();
@@ -111 +120 @@
       //cout << "WeiTot = " << WeiTot << endl;
 
-      int hasBJ = t[i]->GetLeaf("hasBadJet")->GetValue();
-      double hf = t[i]->GetLeaf("hfor")->GetValue();
+      int hasBJ = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int hf    = (int)t[i]->GetLeaf("hfor")->GetValue();
 
       if (hf!=4 && hasBJ==0){
@@ -120 +129 @@
       
         if (mt > Cut4 && HT > Cut5) {
-          icount++;
-          iTot[i] += WeiTot;
+          m_iTot[i] += WeiTot;
           int it = -1;
           if (cos1<Cut1 && cos2<Cut2){
@@ -138 +146 @@
           }
           if (i == 0) 
-            iCatB[it] += WeiTot;
+            m_iCatB[it] += WeiTot;
           else
-            iCatS[it] += WeiTot;
+            m_iCatS[it] += WeiTot;
         }
       }
@@ -146 +154 @@
   }//End of reading 2 files
 
-  cout << "iTot0 = " << iTot[0] << endl;
-  cout << "iTot1 = " << iTot[1] << endl;
-
-  N1 = iTot[1];
-  N0 = iTot[0];
+  m_N0 = m_iTot[0];
+  m_N1 = m_iTot[1];
+
+  cout << "To calculate efficiencies using " << type[1] << endl;
+  cout << "iTot0 = " << m_iTot[0] << endl;
+  cout << "iTot1 = " << m_iTot[1] << endl;
 
   double sigmaInv = 0;
   double sumeffS = 0, sumeffB = 0;
   for (int i = 0; i < 6; i++){
-    F1[i] = iCatS[i]/iTot[1];
-    F0[i] = iCatB[i]/iTot[0];
-    sumeffS += F1[i];
-    sumeffB += F0[i];
-    cout << "F1[" << i << "] = " << F1[i] << endl;
-    cout << "F0[" << i << "] = " << F0[i] << endl;
-
-    double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
-    if (F1[i]==0 && F0[i]==0)  term = 0;
+    m_F1[i] = m_iCatS[i]/m_iTot[1];
+    m_F0[i] = m_iCatB[i]/m_iTot[0];
+    sumeffS += m_F1[i];
+    sumeffB += m_F0[i];
+    cout << "m_F1[" << i << "] = " << m_F1[i] << endl;
+    cout << "m_F0[" << i << "] = " << m_F0[i] << endl;
+
+    double term = m_F1[i]*m_F1[i]/(m_N1*m_F1[i] + m_N0*m_F0[i]);
+    if (m_F1[i]==0 && m_F0[i]==0)  term = 0;
     sigmaInv += term;
   }
-
-  double sigma2 = 1/N1/N1/sigmaInv;
-  sigmaMin = sqrt(sigma2);
-  
-  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
-  cout << "sigmaMin = " << sigmaMin << endl;
 }
 
@@ -195 +198 @@
   }
   
-  f0e = F0[it];
-  f1e = F1[it];
+  f0e = m_F0[it];
+  f1e = m_F1[it];
   
   //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
@@ -207 +210 @@
 
 void MakeEventListf0f1(){
+  m_ListAlreadyDone = true;
   //######## COMPUTE the sWeights
   TString filename = "";
-  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
-  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
+  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
   TFile *file[2];
   TTree *t[2];
   
   for (int i = 0; i < 2; i++){
-    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
+    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
     file[i] = TFile::Open(filename);
     t[i]    = (TTree*)file[i]->Get("T");
@@ -231 +234 @@
       double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
       double WeiTot = t[i]->GetLeaf("wei")->GetValue();
+
+      // Yes some events have 0 weight due to PU reweighting
+      // (just to save time, remove them)
+      if (fabs(WeiTot) < 1e-9) continue; 
+      //cout << "Reading event " << iev << " of type " << i << " w = " << WeiTot << endl; 
       
       //cout << "WeiTot = " << WeiTot << endl;
       
-      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
-      double hf = t[i]->GetLeaf("hfor")->GetValue();
+      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int hf     = (int)t[i]->GetLeaf("hfor")->GetValue();
 
       if (hf!=4 && hasBJ==0){
@@ -255 +263 @@
 }
 
-Fit(){
+void Fit(bool FitN0){
   // Test
-  MakeEventListf0f1();
-
-  int ifl = 0;
-  double step = 0.01;
-  double fitPar0 = iTot[0];
-  double fitPar1 = iTot[1];
-  double nmax0   = fitPar0 + 4*sqrt(iTot[0]);
-  double nmin0   = fitPar0 - 4*sqrt(iTot[0]);
-  double nmax1   = fitPar1 + 4*sqrt(iTot[1]);
-  double nmin1   = - 4*sqrt(iTot[1]);
-  double n0best, err0best;
-  double n1best, err1best;
-  TMinuit* myMinuit = new TMinuit(2);
-  //myMinuit->Command("SET PRI -1");
-  myMinuit->Command("SET NOW");
-  myMinuit->Command("SET ERR 0.5");
-  myMinuit->Command("SET STR 2");   // try to improve mimimum search
-  myMinuit->SetFCN(fcn);
-  myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl); 
-  myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl); 
-  //myMinuit->FixParameter(0);
-  myMinuit->Command("MIGRAD 500 1");
-  myMinuit->GetParameter(0,n0best,err0best); 
-  myMinuit->GetParameter(1,n1best,err1best); 
-  std::cout << "Expected N0 = " << iTot[0] << " fitted = " << n0best << " +- " << err0best << std::endl;
-  std::cout << "Expected N1 = " << iTot[1] << " fitted = " << n1best << " +- " << err1best << std::endl;
-
-  N0 = n0best;
-  N1 = n1best;
-
-  cout << "N0 = " << N0 << " N1 = " << N1 << endl;
-  sig2 = sigmaMin*sigmaMin*N1*N1;
+  if (!m_ListAlreadyDone) MakeEventListf0f1();
+  cout << "Number of events in the LH " << m_EventList.size() << endl;
+
+  if (FitN0){
+
+    double N0_exp = m_iTot[0];
+    double N1_exp = m_iTot[1];
+
+    cout << "FitN0N1 " << " N0_exp " << N0_exp << " N1_exp " << N1_exp << endl;
+
+    int ifl = 0;
+    double step = 0.01;
+
+    double fitPar0 = N0_exp;
+    double nmax0   = fitPar0 + 4*sqrt(N0_exp);
+    double nmin0   = fitPar0 - 4*sqrt(N0_exp);
+
+    double fitPar1 = N1_exp;
+    double nmax1   = fitPar1 + 4*sqrt(N1_exp);
+    double nmin1   =         - 4*sqrt(N1_exp);
+
+    double n0best, err0best;
+    double n1best, err1best;
+
+    TMinuit* myMinuit = new TMinuit(2);
+    //myMinuit->Command("SET PRI -1");
+    myMinuit->Command("SET NOW");
+    myMinuit->Command("SET ERR 0.5");
+    myMinuit->Command("SET STR 2");   // try to improve mimimum search
+    myMinuit->SetFCN(fcnN0N1);
+    myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl); 
+    myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl); 
+    //myMinuit->FixParameter(0);
+    myMinuit->Command("MIGRAD 500 1");
+    myMinuit->GetParameter(0,n0best,err0best); 
+    myMinuit->GetParameter(1,n1best,err1best); 
+
+    std::cout << "Expected N0 = " << N0_exp << " fitted = " << n0best << " +- " << err0best << std::endl;
+    std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
+    
+    m_N0 = n0best;
+    m_N1 = n1best;
+  }
+  else {
+    double N0_exp = m_iTot[0];
+    double N1_exp = m_iTot[1];
+
+    cout << "FitN1 " << " N0_exp " << N0_exp << " N1_exp " << N1_exp << endl;
+
+    int ifl = 0;
+    double step = 0.01;
+    double fitPar1 = N1_exp;
+    double nmax1   = fitPar1 + 4*sqrt(N1_exp);
+    double nmin1   = - 4*sqrt(N1_exp);
+    double n1best, err1best;
+    TMinuit* myMinuitN1 = new TMinuit(1);
+    //myMinuitN1->Command("SET PRI -1");
+    myMinuitN1->Command("SET NOW");
+    myMinuitN1->Command("SET ERR 0.5");
+    myMinuitN1->Command("SET STR 2");   // try to improve mimimum search
+    myMinuitN1->SetFCN(fcnN1);
+    myMinuitN1->mnparm(0,"N1",fitPar1,step,nmin1,nmax1,ifl); 
+    //myMinuitN1->FixParameter(0);
+    myMinuitN1->Command("MIGRAD 500 1");
+    myMinuitN1->GetParameter(0,n1best,err1best); 
+    std::cout << "Expected N0 = " << N0_exp << endl;
+    std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
+   
+    m_N1 = n1best;
+    m_sig2 = err1best*err1best;
+  }
+
+  cout << "bool FitN0 " << FitN0 << " m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
+  //cout << "m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
 }
 
@@ -296 +346 @@
 
 void sPlots_NewRel(){
-  //gROOT->SetStyle("Plain");
-  //Make the table of sP values first
-
   cout << "Cut1 = " << Cut1 << endl;
   cout << "Cut2 = " << Cut2 << endl;
@@ -305 +352 @@
   cout << "Cut5 = " << Cut5 << endl;
 
+  //Compute the efficiencies first
   CalF0F1();
-  cout << "icount = " << icount << endl;
+
+  //Reset m_iTot[0] and m_iTot[1]
+  m_iTot[0] = 0.; m_iTot[1] = 0.;
+  //double F0N[6]
+
+  //Now compute m_N0, m_N1 
+  TString filename = "";
+  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
+  TFile *file[2];
+  TTree *t[2];
+  
+  for (int i = 0; i < 2; i++){
+    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
+    file[i] = TFile::Open(filename);
+    t[i]    = (TTree*)file[i]->Get("T");
+  }
+
+  for (int i = 0; i < 2; i++){
+    int nevt = t[i]->GetEntries();
+    for (int iev = 0; iev < nevt; iev++) {
+      t[i]->GetEntry(iev);
+            
+      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
+      double HT     = t[i]->GetLeaf("HT")->GetValue();
+      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
+
+      int hasBJ = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int    hf = (int)t[i]->GetLeaf("hfor")->GetValue();
+
+      if (hf!=4 && hasBJ==0){
+        if (mt > Cut4 && HT > Cut5) {
+          m_icount++;
+          m_iTot[i] += WeiTot;
+        }
+      }
+    }//Read each file
+  }//End of reading 2 files
+
+  cout << "In main macro" << endl;
+  cout << "iTot0 = " << m_iTot[0] << endl;
+  cout << "iTot1 = " << m_iTot[1] << endl;
+
+  m_N1 = m_iTot[1];
+  m_N0 = m_iTot[0];
+
+  cout << "m_icount = " << m_icount << endl;
 
   //sPlots N0 known
@@ -314 +407 @@
   sPlots_NewRel_N0known();
 
-  //Get right N0 and N1 for standard sPlot
   cout << "***********************************************" << endl;
   cout << "************* sPlot N0 unknown ****************" << endl;
   cout << "***********************************************" << endl;
 
-  Fit();
+  //Get right N0 and N1 for standard sPlot
+  Fit(true);
  
-  TString filename = "";
-  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
-  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
-  TFile *file[2];
-  TTree *t[2];
-  
-  for (int i = 0; i < 2; i++){
-    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
-    file[i] = TFile::Open(filename);
-    t[i]    = (TTree*)file[i]->Get("T");
-  }
-
   double VI00 = 0, VI01 = 0, VI10 = 0, VI11 = 0;
   //start reading 2 files
-  //int icount = 0;
+  //int m_icount = 0;
   for (int i = 0; i < 2; i++){
     int nevt = t[i]->GetEntries();
@@ -351 +432 @@
       double cos2 = fabs(y);
       
-      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
-      double hf = t[i]->GetLeaf("hfor")->GetValue();
+      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
 
       if (hf!=4 && hasBJ==0){
@@ -368 +449 @@
           double vi11;
           
-          vi00 = WeiTot*F0e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
-          vi01 = WeiTot*F0e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
-          vi10 = WeiTot*F1e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
-          vi11 = WeiTot*F1e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
+          vi00 = WeiTot*F0e*F0e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
+          vi01 = WeiTot*F0e*F1e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
+          vi10 = WeiTot*F1e*F0e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
+          vi11 = WeiTot*F1e*F1e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
 
           VI00 += vi00;
@@ -382 +463 @@
   }//End of reading 2 files
 
-  cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
+  cout << "Main macro : m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
 
   cout << "VI00 = " << VI00 << endl;
@@ -402 +483 @@
   cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl; 
 
-  TH1F *hSigN  = new TH1F("hSigN",     "", 35, 250, 600);
-  TH1F *hBkgN  = new TH1F("hBkgN",     "", 35, 250, 600);
-  TH2F *h2DSig = new TH2F("h2DSig",    "", 35, 250, 600, 50, 500, 1500);
-  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 35, 250, 600, 50, 500, 1500);
+  TH1F *hSigN  = new TH1F("hSigN",     "", 150, 0, 1500);
+  TH1F *hBkgN  = new TH1F("hBkgN",     "", 150, 0, 1500);
+  TH2F *h2DSig = new TH2F("h2DSig",    "", 150, 0, 1500, 50, 500, 1500);
+  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 150, 0, 1500, 50, 500, 1500);
 
   //Try standard sPlots
   double SumsP1 = 0.;
   double SumsP0 = 0.;
-  double *sPN = new double[icount];
-  double *mtN = new double[icount];
+  double *sPN = new double[m_icount];
+  double *mtN = new double[m_icount];
   int icc = 0;
   for (int i = 0; i < 2; i++){
@@ -424 +505 @@
       double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
       double WeiTot = t[i]->GetLeaf("wei")->GetValue();
-      
-      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
-      double hf = t[i]->GetLeaf("hfor")->GetValue();
+      if (fabs(WeiTot) < 1e-9) continue; 
+      
+      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
 
       if (hf!=4 && hasBJ==0){
@@ -439 +521 @@
           double F1e = veff.at(1);
           
-          double sP1 = (V11*F1e + V10*F0e)/(N1*F1e + N0*F0e);
-          double sP0 = (V00*F0e + V10*F1e)/(N1*F1e + N0*F0e);
+          double sP1 = (V11*F1e + V10*F0e)/(m_N1*F1e + m_N0*F0e);
+          double sP0 = (V00*F0e + V10*F1e)/(m_N1*F1e + m_N0*F0e);
           
           hSigN->Fill(mt,sP1*WeiTot);
@@ -447 +529 @@
           h2DBkg->Fill(mt,HT,sP0*WeiTot);
           mtN[icc] = mt;
-          if (WeiTot==0) sPN[icc] = 0;
-          else sPN[icc] = sP1*sP1*WeiTot;
+          sPN[icc] = sP1*sP1*WeiTot;
           icc++;
           
@@ -465 +546 @@
   cout << "hBkgN integral  " << hBkgN->Integral(0,1500)  << endl;
   
-  for (int i = 1; i < 36; i++){
-    double vb1 = 250 + (i-1)*10;
-    double vb2 = 250 + i*10;
+  for (int i = 1; i < 151; i++){
+    double vb1 = (i-1)*10;
+    double vb2 = i*10;
     double ErrBin = 0.;
-    for (int j = 0; j < icount; j++){
+    for (int j = 0; j < m_icount; j++){
       double mt = mtN[j];
       double sp = sPN[j];
@@ -476 +557 @@
       }
     }
-    hSigN->SetBinError(i,sqrt(ErrBin));
-  }
-
-  TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
-
-  TCanvas *c = new TCanvas("c", "", 900, 400);
+    if (ErrBin<0) cout << "ErrBin is negative in sPlot standard " << i << " " << ErrBin << endl;
+    hSigN->SetBinError(i,sqrt(fabs(ErrBin)));
+  }
+
+  //TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
+
+  TCanvas *c = new TCanvas("c", "cStandard", 900, 400);
   c->Divide(2);
   c->cd(1); 
   //frame->Draw();
+  double x1 = 250;
+  double x2 = 600 - 0.01;
+  TAxis *ax = hSigN->GetXaxis();
+  int b1 = ax->FindBin(x1);
+  int b2 = ax->FindBin(x2);
+  ax->SetRange(b1,b2);
   hSigN->SetLineColor(2);
   hSigN->SetLineWidth(1);
   hSigN->DrawCopy("E1");
   hSigN->SetLineColor(1);
-  hSigN->SetLineWidth(2.5);
+  hSigN->SetLineWidth(2);
   hSigN->DrawCopy("hist,same");
   c->cd(2);
+  TAxis *axB = hBkgN->GetXaxis();
+  int b1B = axB->FindBin(x1);
+  int b2B = axB->FindBin(x2);
+  axB->SetRange(b1B,b2B);
   hBkgN->SetLineColor(4);
-  hBkgN->SetLineWidth(2.5);
+  hBkgN->SetLineWidth(2);
   hBkgN->Draw();
 
@@ -510 +602 @@
 
 void sPlots_NewRel_N0known(){
-  cout << "in function N0 known " << icount << endl;
-  N0 = iTot[0];
-  N1 = iTot[1];
-  sig2 = sigmaMin*sigmaMin*N1*N1;
-  cout << "sigmaMin = " << sigmaMin << " sig2 = " << sig2 << endl;
-  cout << "N0 = " << N0 << " N1 = " << N1 << endl;
-
-  double k = (sig2 - N1)/(N0 - (sig2 - N1));
+
+  Fit(false);
+
+  cout << "in function N0 known " << m_icount << endl;
+  cout << "m_sig2 = " << m_sig2 << endl;
+  cout << "m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
+
+  double k = (m_sig2 - m_N1)/(m_N0 - (m_sig2 - m_N1));
   
   TH2F *twoDim_mt_HT_Sig     = new TH2F("mt_HT_Sig",     "", 75, 250, 1500, 50, 500, 2000); 
   TH2F *twoDim_mt_HT_Sig_esP = new TH2F("mt_HT_Sig_esP", "", 75, 250, 1500, 50, 500, 2000); 
   
-  TH1F *h[7];
   h[0]  = new TH1F("All1fb",       "",  150,  0,  1500);
   h[1]  = new TH1F("Sig1fb",       "",  150,  0,  1500);
@@ -531 +622 @@
   h[6]  = new TH1F("esPlotSig",    "",  150,  0,  1500);
   
-  double *mtV  = new double[icount];
-  double *err  = new double[icount];
-  double *errN = new double[icount];
+  double *mtV  = new double[m_icount];
+  double *err  = new double[m_icount];
+  double *errN = new double[m_icount];
   int icc = 0;
   
   TString filename = "";
-  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
+  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
   TFile *file[2];
   TTree *t[2];
   for (int i = 0; i < 2; i++){
-    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
+    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
     file[i] = TFile::Open(filename);
     t[i]    = (TTree*)file[i]->Get("T");
@@ -558 +649 @@
       double WeiTot = t[i]->GetLeaf("wei")->GetValue();
 
-      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
-      double hf = t[i]->GetLeaf("hfor")->GetValue();
+      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
+      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
       
       if (hf!=4 && hasBJ==0){
@@ -572 +663 @@
           double F1e = veff.at(1);
           
-          double sWei  = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
+          double sWei  = m_sig2*F1e*WeiTot/(m_N1*F1e + m_N0*F0e);
           double sWeiN = sWei*(k+1) - k*WeiTot;
+
+          //cout << "sWei  = " << sWei << endl;
+          //cout << "sWeiN = " << sWeiN << endl;
           
           //For 1 fb-1
@@ -586 +680 @@
           if (i==1){
             h[1]->Fill(mt,WeiTot);
+            //h[4]->Fill(mt,sWei);
             twoDim_mt_HT_Sig->Fill(mt,HT,sWei);
             twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);    
@@ -609 +704 @@
   for (int i = 1; i < 151; i++){
     double sP = h[3]->GetBinContent(i);
-    double B  = (sig2 - N1)*h[2]->GetBinContent(i);
     double sB = h[5]->GetBinContent(i);
-    double W  = sP - B;
     double sW  = sP - sB;
-    double Tot1fb = h[0]->GetBinContent(i);
-    double BkgPure = Tot1fb - sP;
+
+    //double B  = (m_sig2 - m_N1)*h[2]->GetBinContent(i);
+    //double W  = sP - B;
+    //double Tot1fb = h[0]->GetBinContent(i);
+    //double BkgPure = Tot1fb - sP;
     
     double vb1 = 10*(i-1);
@@ -620 +716 @@
     double SumErrBin = 0;
     double SumErrNBin = 0;
-    for (int ic = 0; ic < icount; ic++){
+    for (int ic = 0; ic < m_icount; ic++){
       double mt  = mtV[ic];
       double errV  = err[ic];
@@ -630 +726 @@
     }
     h[4]->SetBinContent(i,sW);
-    h[4]->SetBinError(i,sqrt(SumErrBin));
-    h[6]->SetBinError(i,sqrt(SumErrNBin));
+    if (SumErrBin < 0) cout << "Warning negative error squared in bin " << i << " for sP" << endl; 
+    h[4]->SetBinError(i,sqrt(fabs(SumErrBin)));
+    if (SumErrNBin < 0) cout << "Warning negative error squared in bin " << i << " for esP" << endl; 
+    h[6]->SetBinError(i,sqrt(fabs(SumErrNBin)));
   }
   
@@ -642 +740 @@
   for (int i = 0; i < 7; i++) {
     //h[i]->Rebin(2);
-    h[i]->SetLineWidth(3);
+    h[i]->SetLineWidth(2);
     TAxis *ax = h[i]->GetXaxis();
     int b1 = ax->FindBin(x1);
@@ -650 +748 @@
 
   //plot for sPlot with M0 known
-  TCanvas *csP = new TCanvas("csP", "", 800, 400);
+  TCanvas *csP = new TCanvas("csP", "csP", 800, 400);
   csP->Divide(2);
   csP->cd(1);
@@ -659 +757 @@
   h[1]->Draw("same");
   h[4]->SetLineColor(1);
-  h[4]->SetLineWidth(3);
+  h[4]->SetLineWidth(2);
   h[4]->DrawCopy("hist,same");
   csP->cd(2);
@@ -665 +763 @@
 
   //plot for sPlot with M0 UNknown
-  TCanvas *cesP = new TCanvas("cesP", "", 800,400);
+  TCanvas *cesP = new TCanvas("cesP", "cesP", 800,400);
   cesP->Divide(2);
   cesP->cd(1);
@@ -672 +770 @@
   h[6]->DrawCopy("E1");
   h[6]->SetLineColor(1);
-  h[6]->SetLineWidth(3);
+  h[6]->SetLineWidth(2);
   h[6]->Draw("hist,same");
   cesP->cd(2);