Changeset 22 in huonglan for SPLOT/NEWREALEASE/sPlots_NewRel.C


Ignore:
Timestamp:
Sep 2, 2011, 3:07:15 PM (13 years ago)
Author:
huonglan
Message:

sPlots_NewRel improved

File:
1 edited

Legend:

Unmodified
Added
Removed

  • SPLOT/NEWREALEASE/sPlots_NewRel.C

    r19 r22  
    1212#include <vector>
    1313
     14#include <iostream>
     15
     16using namespace std;
     17
     18void sPlots_NewRel_N0known();
     19
    1420//Set the cut values
    1521double Cut1 = 0.7;
    … …  
    1824double Cut4 = 250;
    1925double Cut5 = 620;
    20 double sigmaMin;
    21 int icount = 0;
    22 double iTot[2];
    23 double iCatS[6], iCatB[6];
    24 double F1[6];
    25 double F0[6];
    26 
    27 double N1 = 0., N0 = 0., sig2 = 0.;
     26
     27//global variables
     28int    m_icount = 0;
     29double m_iTot[2];              //for (int i = 0; i < 2; i++){ m_iTot[i] = 0;}
     30double m_iCatS[6], m_iCatB[6]; //for (int i = 0; i < 6; i++){ m_iCatS[i] = 0; m_iCatB[i] = 0;}
     31double m_F1[6];                //for (int i = 0; i < 6; i++){ m_F1[i] = 0; }
     32double m_F0[6];                //for (int i = 0; i < 6; i++){ m_F0[i] = 0; }
     33double m_N1 = 0., m_N0 = 0., m_sig2 = 0.;
     34
     35TH1F *h[7];
    2836
    2937//Define class Event
    … …  
    4654
    4755std::vector<EventC> m_EventList;
    48 
    49 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
     56bool m_ListAlreadyDone = false;
     57
     58void fcnN1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
    5059  double L = 0;
    51   for (int iev = 0; iev < m_EventList.size(); iev++) {
     60  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
     61    EventC evt = m_EventList[iev];
     62    double smt = m_N0*evt.fb() + par[0]*evt.fs();
     63    if (smt > 0) L += evt.w()*log(smt);
     64  }
     65  f = par[0] - L;
     66};
     67
     68void fcnN0N1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
     69  double L = 0;
     70  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
    5271    EventC evt = m_EventList[iev];
    5372    double smt = par[0]*evt.fb() + par[1]*evt.fs();
    … …  
    5675  f = par[0] + par[1] - L;
    5776};
    58 
    5977
    6078vector<double> InvMatrix(double a, double b, double c, double d){
    … …  
    7896  //######## COMPUTE the sWeights
    7997  TString filename = "";
    80   //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
    81   TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
     98  TString type[2] = {"AllBkgTopMCAtNLO","Signal115420"};
    8299  TFile *file[2];
    83100  TTree *t[2];
    84101 
    85102  for (int i = 0; i < 2; i++){
    86     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
     103    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
    87104    file[i] = TFile::Open(filename);
    88105    t[i]    = (TTree*)file[i]->Get("T");
    89106  }
    90107
    91   for (int i = 0; i < 2; i++){
    92     iTot[i] = 0;
    93   }
    94   for (int i = 0; i < 6; i++){
    95     iCatS[i] = 0;
    96     iCatB[i] = 0;
    97   }
    98  
    99108  for (int i = 0; i < 2; i++){
    100109    int nevt = t[i]->GetEntries();
    … …  
    111120      //cout << "WeiTot = " << WeiTot << endl;
    112121
    113       int hasBJ = t[i]->GetLeaf("hasBadJet")->GetValue();
    114       double hf = t[i]->GetLeaf("hfor")->GetValue();
     122      int hasBJ = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     123      int hf    = (int)t[i]->GetLeaf("hfor")->GetValue();
    115124
    116125      if (hf!=4 && hasBJ==0){
    … …  
    120129     
    121130        if (mt > Cut4 && HT > Cut5) {
    122           icount++;
    123           iTot[i] += WeiTot;
     131          m_iTot[i] += WeiTot;
    124132          int it = -1;
    125133          if (cos1<Cut1 && cos2<Cut2){
    … …  
    138146          }
    139147          if (i == 0)
    140             iCatB[it] += WeiTot;
     148            m_iCatB[it] += WeiTot;
    141149          else
    142             iCatS[it] += WeiTot;
     150            m_iCatS[it] += WeiTot;
    143151        }
    144152      }
    … …  
    146154  }//End of reading 2 files
    147155
    148   cout << "iTot0 = " << iTot[0] << endl;
    149   cout << "iTot1 = " << iTot[1] << endl;
    150 
    151   N1 = iTot[1];
    152   N0 = iTot[0];
     156  m_N0 = m_iTot[0];
     157  m_N1 = m_iTot[1];
     158
     159  cout << "To calculate efficiencies using " << type[1] << endl;
     160  cout << "iTot0 = " << m_iTot[0] << endl;
     161  cout << "iTot1 = " << m_iTot[1] << endl;
    153162
    154163  double sigmaInv = 0;
    155164  double sumeffS = 0, sumeffB = 0;
    156165  for (int i = 0; i < 6; i++){
    157     F1[i] = iCatS[i]/iTot[1];
    158     F0[i] = iCatB[i]/iTot[0];
    159     sumeffS += F1[i];
    160     sumeffB += F0[i];
    161     cout << "F1[" << i << "] = " << F1[i] << endl;
    162     cout << "F0[" << i << "] = " << F0[i] << endl;
    163 
    164     double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
    165     if (F1[i]==0 && F0[i]==0)  term = 0;
     166    m_F1[i] = m_iCatS[i]/m_iTot[1];
     167    m_F0[i] = m_iCatB[i]/m_iTot[0];
     168    sumeffS += m_F1[i];
     169    sumeffB += m_F0[i];
     170    cout << "m_F1[" << i << "] = " << m_F1[i] << endl;
     171    cout << "m_F0[" << i << "] = " << m_F0[i] << endl;
     172
     173    double term = m_F1[i]*m_F1[i]/(m_N1*m_F1[i] + m_N0*m_F0[i]);
     174    if (m_F1[i]==0 && m_F0[i]==0)  term = 0;
    166175    sigmaInv += term;
    167176  }
    168 
    169   double sigma2 = 1/N1/N1/sigmaInv;
    170   sigmaMin = sqrt(sigma2);
    171  
    172   cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
    173   cout << "sigmaMin = " << sigmaMin << endl;
    174177}
    175178
    … …  
    195198  }
    196199 
    197   f0e = F0[it];
    198   f1e = F1[it];
     200  f0e = m_F0[it];
     201  f1e = m_F1[it];
    199202 
    200203  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
    … …  
    207210
    208211void MakeEventListf0f1(){
     212  m_ListAlreadyDone = true;
    209213  //######## COMPUTE the sWeights
    210214  TString filename = "";
    211   //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
    212   TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
     215  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
    213216  TFile *file[2];
    214217  TTree *t[2];
    215218 
    216219  for (int i = 0; i < 2; i++){
    217     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
     220    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
    218221    file[i] = TFile::Open(filename);
    219222    t[i]    = (TTree*)file[i]->Get("T");
    … …  
    231234      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    232235      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     236
     237      // Yes some events have 0 weight due to PU reweighting
     238      // (just to save time, remove them)
     239      if (fabs(WeiTot) < 1e-9) continue;
     240      //cout << "Reading event " << iev << " of type " << i << " w = " << WeiTot << endl;
    233241     
    234242      //cout << "WeiTot = " << WeiTot << endl;
    235243     
    236       int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    237       double hf = t[i]->GetLeaf("hfor")->GetValue();
     244      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     245      int hf     = (int)t[i]->GetLeaf("hfor")->GetValue();
    238246
    239247      if (hf!=4 && hasBJ==0){
    … …  
    255263}
    256264
    257 Fit(){
     265void Fit(bool FitN0){
    258266  // Test
    259   MakeEventListf0f1();
    260 
    261   int ifl = 0;
    262   double step = 0.01;
    263   double fitPar0 = iTot[0];
    264   double fitPar1 = iTot[1];
    265   double nmax0   = fitPar0 + 4*sqrt(iTot[0]);
    266   double nmin0   = fitPar0 - 4*sqrt(iTot[0]);
    267   double nmax1   = fitPar1 + 4*sqrt(iTot[1]);
    268   double nmin1   = - 4*sqrt(iTot[1]);
    269   double n0best, err0best;
    270   double n1best, err1best;
    271   TMinuit* myMinuit = new TMinuit(2);
    272   //myMinuit->Command("SET PRI -1");
    273   myMinuit->Command("SET NOW");
    274   myMinuit->Command("SET ERR 0.5");
    275   myMinuit->Command("SET STR 2");   // try to improve mimimum search
    276   myMinuit->SetFCN(fcn);
    277   myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl);
    278   myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl);
    279   //myMinuit->FixParameter(0);
    280   myMinuit->Command("MIGRAD 500 1");
    281   myMinuit->GetParameter(0,n0best,err0best);
    282   myMinuit->GetParameter(1,n1best,err1best);
    283   std::cout << "Expected N0 = " << iTot[0] << " fitted = " << n0best << " +- " << err0best << std::endl;
    284   std::cout << "Expected N1 = " << iTot[1] << " fitted = " << n1best << " +- " << err1best << std::endl;
    285 
    286   N0 = n0best;
    287   N1 = n1best;
    288 
    289   cout << "N0 = " << N0 << " N1 = " << N1 << endl;
    290   sig2 = sigmaMin*sigmaMin*N1*N1;
     267  if (!m_ListAlreadyDone) MakeEventListf0f1();
     268  cout << "Number of events in the LH " << m_EventList.size() << endl;
     269
     270  if (FitN0){
     271
     272    double N0_exp = m_iTot[0];
     273    double N1_exp = m_iTot[1];
     274
     275    cout << "FitN0N1 " << " N0_exp " << N0_exp << " N1_exp " << N1_exp << endl;
     276
     277    int ifl = 0;
     278    double step = 0.01;
     279
     280    double fitPar0 = N0_exp;
     281    double nmax0   = fitPar0 + 4*sqrt(N0_exp);
     282    double nmin0   = fitPar0 - 4*sqrt(N0_exp);
     283
     284    double fitPar1 = N1_exp;
     285    double nmax1   = fitPar1 + 4*sqrt(N1_exp);
     286    double nmin1   =         - 4*sqrt(N1_exp);
     287
     288    double n0best, err0best;
     289    double n1best, err1best;
     290
     291    TMinuit* myMinuit = new TMinuit(2);
     292    //myMinuit->Command("SET PRI -1");
     293    myMinuit->Command("SET NOW");
     294    myMinuit->Command("SET ERR 0.5");
     295    myMinuit->Command("SET STR 2");   // try to improve mimimum search
     296    myMinuit->SetFCN(fcnN0N1);
     297    myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl);
     298    myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl);
     299    //myMinuit->FixParameter(0);
     300    myMinuit->Command("MIGRAD 500 1");
     301    myMinuit->GetParameter(0,n0best,err0best);
     302    myMinuit->GetParameter(1,n1best,err1best);
     303
     304    std::cout << "Expected N0 = " << N0_exp << " fitted = " << n0best << " +- " << err0best << std::endl;
     305    std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
     306   
     307    m_N0 = n0best;
     308    m_N1 = n1best;
     309  }
     310  else {
     311    double N0_exp = m_iTot[0];
     312    double N1_exp = m_iTot[1];
     313
     314    cout << "FitN1 " << " N0_exp " << N0_exp << " N1_exp " << N1_exp << endl;
     315
     316    int ifl = 0;
     317    double step = 0.01;
     318    double fitPar1 = N1_exp;
     319    double nmax1   = fitPar1 + 4*sqrt(N1_exp);
     320    double nmin1   = - 4*sqrt(N1_exp);
     321    double n1best, err1best;
     322    TMinuit* myMinuitN1 = new TMinuit(1);
     323    //myMinuitN1->Command("SET PRI -1");
     324    myMinuitN1->Command("SET NOW");
     325    myMinuitN1->Command("SET ERR 0.5");
     326    myMinuitN1->Command("SET STR 2");   // try to improve mimimum search
     327    myMinuitN1->SetFCN(fcnN1);
     328    myMinuitN1->mnparm(0,"N1",fitPar1,step,nmin1,nmax1,ifl);
     329    //myMinuitN1->FixParameter(0);
     330    myMinuitN1->Command("MIGRAD 500 1");
     331    myMinuitN1->GetParameter(0,n1best,err1best);
     332    std::cout << "Expected N0 = " << N0_exp << endl;
     333    std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
     334   
     335    m_N1 = n1best;
     336    m_sig2 = err1best*err1best;
     337  }
     338
     339  cout << "bool FitN0 " << FitN0 << " m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
     340  //cout << "m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
    291341}
    292342
    … …  
    296346
    297347void sPlots_NewRel(){
    298   //gROOT->SetStyle("Plain");
    299   //Make the table of sP values first
    300 
    301348  cout << "Cut1 = " << Cut1 << endl;
    302349  cout << "Cut2 = " << Cut2 << endl;
    … …  
    305352  cout << "Cut5 = " << Cut5 << endl;
    306353
     354  //Compute the efficiencies first
    307355  CalF0F1();
    308   cout << "icount = " << icount << endl;
     356
     357  //Reset m_iTot[0] and m_iTot[1]
     358  m_iTot[0] = 0.; m_iTot[1] = 0.;
     359  //double F0N[6]
     360
     361  //Now compute m_N0, m_N1
     362  TString filename = "";
     363  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
     364  TFile *file[2];
     365  TTree *t[2];
     366 
     367  for (int i = 0; i < 2; i++){
     368    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
     369    file[i] = TFile::Open(filename);
     370    t[i]    = (TTree*)file[i]->Get("T");
     371  }
     372
     373  for (int i = 0; i < 2; i++){
     374    int nevt = t[i]->GetEntries();
     375    for (int iev = 0; iev < nevt; iev++) {
     376      t[i]->GetEntry(iev);
     377           
     378      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
     379      double HT     = t[i]->GetLeaf("HT")->GetValue();
     380      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     381
     382      int hasBJ = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     383      int    hf = (int)t[i]->GetLeaf("hfor")->GetValue();
     384
     385      if (hf!=4 && hasBJ==0){
     386        if (mt > Cut4 && HT > Cut5) {
     387          m_icount++;
     388          m_iTot[i] += WeiTot;
     389        }
     390      }
     391    }//Read each file
     392  }//End of reading 2 files
     393
     394  cout << "In main macro" << endl;
     395  cout << "iTot0 = " << m_iTot[0] << endl;
     396  cout << "iTot1 = " << m_iTot[1] << endl;
     397
     398  m_N1 = m_iTot[1];
     399  m_N0 = m_iTot[0];
     400
     401  cout << "m_icount = " << m_icount << endl;
    309402
    310403  //sPlots N0 known
    … …  
    314407  sPlots_NewRel_N0known();
    315408
    316   //Get right N0 and N1 for standard sPlot
    317409  cout << "***********************************************" << endl;
    318410  cout << "************* sPlot N0 unknown ****************" << endl;
    319411  cout << "***********************************************" << endl;
    320412
    321   Fit();
     413  //Get right N0 and N1 for standard sPlot
     414  Fit(true);
    322415 
    323   TString filename = "";
    324   //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
    325   TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
    326   TFile *file[2];
    327   TTree *t[2];
    328  
    329   for (int i = 0; i < 2; i++){
    330     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    331     file[i] = TFile::Open(filename);
    332     t[i]    = (TTree*)file[i]->Get("T");
    333   }
    334 
    335416  double VI00 = 0, VI01 = 0, VI10 = 0, VI11 = 0;
    336417  //start reading 2 files
    337   //int icount = 0;
     418  //int m_icount = 0;
    338419  for (int i = 0; i < 2; i++){
    339420    int nevt = t[i]->GetEntries();
    … …  
    351432      double cos2 = fabs(y);
    352433     
    353       int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    354       double hf = t[i]->GetLeaf("hfor")->GetValue();
     434      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     435      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
    355436
    356437      if (hf!=4 && hasBJ==0){
    … …  
    368449          double vi11;
    369450         
    370           vi00 = WeiTot*F0e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
    371           vi01 = WeiTot*F0e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
    372           vi10 = WeiTot*F1e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
    373           vi11 = WeiTot*F1e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
     451          vi00 = WeiTot*F0e*F0e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
     452          vi01 = WeiTot*F0e*F1e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
     453          vi10 = WeiTot*F1e*F0e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
     454          vi11 = WeiTot*F1e*F1e/(m_N1*F1e + m_N0*F0e)/(m_N1*F1e + m_N0*F0e);
    374455
    375456          VI00 += vi00;
    … …  
    382463  }//End of reading 2 files
    383464
    384   cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
     465  cout << "Main macro : m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
    385466
    386467  cout << "VI00 = " << VI00 << endl;
    … …  
    402483  cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl;
    403484
    404   TH1F *hSigN  = new TH1F("hSigN",     "", 35, 250, 600);
    405   TH1F *hBkgN  = new TH1F("hBkgN",     "", 35, 250, 600);
    406   TH2F *h2DSig = new TH2F("h2DSig",    "", 35, 250, 600, 50, 500, 1500);
    407   TH2F *h2DBkg = new TH2F("h2DBkg",    "", 35, 250, 600, 50, 500, 1500);
     485  TH1F *hSigN  = new TH1F("hSigN",     "", 150, 0, 1500);
     486  TH1F *hBkgN  = new TH1F("hBkgN",     "", 150, 0, 1500);
     487  TH2F *h2DSig = new TH2F("h2DSig",    "", 150, 0, 1500, 50, 500, 1500);
     488  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 150, 0, 1500, 50, 500, 1500);
    408489
    409490  //Try standard sPlots
    410491  double SumsP1 = 0.;
    411492  double SumsP0 = 0.;
    412   double *sPN = new double[icount];
    413   double *mtN = new double[icount];
     493  double *sPN = new double[m_icount];
     494  double *mtN = new double[m_icount];
    414495  int icc = 0;
    415496  for (int i = 0; i < 2; i++){
    … …  
    424505      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    425506      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    426      
    427       int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    428       double hf = t[i]->GetLeaf("hfor")->GetValue();
     507      if (fabs(WeiTot) < 1e-9) continue;
     508     
     509      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     510      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
    429511
    430512      if (hf!=4 && hasBJ==0){
    … …  
    439521          double F1e = veff.at(1);
    440522         
    441           double sP1 = (V11*F1e + V10*F0e)/(N1*F1e + N0*F0e);
    442           double sP0 = (V00*F0e + V10*F1e)/(N1*F1e + N0*F0e);
     523          double sP1 = (V11*F1e + V10*F0e)/(m_N1*F1e + m_N0*F0e);
     524          double sP0 = (V00*F0e + V10*F1e)/(m_N1*F1e + m_N0*F0e);
    443525         
    444526          hSigN->Fill(mt,sP1*WeiTot);
    … …  
    447529          h2DBkg->Fill(mt,HT,sP0*WeiTot);
    448530          mtN[icc] = mt;
    449           if (WeiTot==0) sPN[icc] = 0;
    450           else sPN[icc] = sP1*sP1*WeiTot;
     531          sPN[icc] = sP1*sP1*WeiTot;
    451532          icc++;
    452533         
    … …  
    465546  cout << "hBkgN integral  " << hBkgN->Integral(0,1500)  << endl;
    466547 
    467   for (int i = 1; i < 36; i++){
    468     double vb1 = 250 + (i-1)*10;
    469     double vb2 = 250 + i*10;
     548  for (int i = 1; i < 151; i++){
     549    double vb1 = (i-1)*10;
     550    double vb2 = i*10;
    470551    double ErrBin = 0.;
    471     for (int j = 0; j < icount; j++){
     552    for (int j = 0; j < m_icount; j++){
    472553      double mt = mtN[j];
    473554      double sp = sPN[j];
    … …  
    476557      }
    477558    }
    478     hSigN->SetBinError(i,sqrt(ErrBin));
    479   }
    480 
    481   TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
    482 
    483   TCanvas *c = new TCanvas("c", "", 900, 400);
     559    if (ErrBin<0) cout << "ErrBin is negative in sPlot standard " << i << " " << ErrBin << endl;
     560    hSigN->SetBinError(i,sqrt(fabs(ErrBin)));
     561  }
     562
     563  //TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
     564
     565  TCanvas *c = new TCanvas("c", "cStandard", 900, 400);
    484566  c->Divide(2);
    485567  c->cd(1);
    486568  //frame->Draw();
     569  double x1 = 250;
     570  double x2 = 600 - 0.01;
     571  TAxis *ax = hSigN->GetXaxis();
     572  int b1 = ax->FindBin(x1);
     573  int b2 = ax->FindBin(x2);
     574  ax->SetRange(b1,b2);
    487575  hSigN->SetLineColor(2);
    488576  hSigN->SetLineWidth(1);
    489577  hSigN->DrawCopy("E1");
    490578  hSigN->SetLineColor(1);
    491   hSigN->SetLineWidth(2.5);
     579  hSigN->SetLineWidth(2);
    492580  hSigN->DrawCopy("hist,same");
    493581  c->cd(2);
     582  TAxis *axB = hBkgN->GetXaxis();
     583  int b1B = axB->FindBin(x1);
     584  int b2B = axB->FindBin(x2);
     585  axB->SetRange(b1B,b2B);
    494586  hBkgN->SetLineColor(4);
    495   hBkgN->SetLineWidth(2.5);
     587  hBkgN->SetLineWidth(2);
    496588  hBkgN->Draw();
    497589
    … …  
    510602
    511603void sPlots_NewRel_N0known(){
    512   cout << "in function N0 known " << icount << endl;
    513   N0 = iTot[0];
    514   N1 = iTot[1];
    515   sig2 = sigmaMin*sigmaMin*N1*N1;
    516   cout << "sigmaMin = " << sigmaMin << " sig2 = " << sig2 << endl;
    517   cout << "N0 = " << N0 << " N1 = " << N1 << endl;
    518 
    519   double k = (sig2 - N1)/(N0 - (sig2 - N1));
     604
     605  Fit(false);
     606
     607  cout << "in function N0 known " << m_icount << endl;
     608  cout << "m_sig2 = " << m_sig2 << endl;
     609  cout << "m_N0 = " << m_N0 << " m_N1 = " << m_N1 << endl;
     610
     611  double k = (m_sig2 - m_N1)/(m_N0 - (m_sig2 - m_N1));
    520612 
    521613  TH2F *twoDim_mt_HT_Sig     = new TH2F("mt_HT_Sig",     "", 75, 250, 1500, 50, 500, 2000);
    522614  TH2F *twoDim_mt_HT_Sig_esP = new TH2F("mt_HT_Sig_esP", "", 75, 250, 1500, 50, 500, 2000);
    523615 
    524   TH1F *h[7];
    525616  h[0]  = new TH1F("All1fb",       "",  150,  0,  1500);
    526617  h[1]  = new TH1F("Sig1fb",       "",  150,  0,  1500);
    … …  
    531622  h[6]  = new TH1F("esPlotSig",    "",  150,  0,  1500);
    532623 
    533   double *mtV  = new double[icount];
    534   double *err  = new double[icount];
    535   double *errN = new double[icount];
     624  double *mtV  = new double[m_icount];
     625  double *err  = new double[m_icount];
     626  double *errN = new double[m_icount];
    536627  int icc = 0;
    537628 
    538629  TString filename = "";
    539   TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
     630  TString type[2] = {"AllBkgTopMCAtNLO","Signal115426"};
    540631  TFile *file[2];
    541632  TTree *t[2];
    542633  for (int i = 0; i < 2; i++){
    543     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
     634    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3.root";
    544635    file[i] = TFile::Open(filename);
    545636    t[i]    = (TTree*)file[i]->Get("T");
    … …  
    558649      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    559650
    560       int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    561       double hf = t[i]->GetLeaf("hfor")->GetValue();
     651      int hasBJ  = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     652      int     hf = (int)t[i]->GetLeaf("hfor")->GetValue();
    562653     
    563654      if (hf!=4 && hasBJ==0){
    … …  
    572663          double F1e = veff.at(1);
    573664         
    574           double sWei  = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
     665          double sWei  = m_sig2*F1e*WeiTot/(m_N1*F1e + m_N0*F0e);
    575666          double sWeiN = sWei*(k+1) - k*WeiTot;
     667
     668          //cout << "sWei  = " << sWei << endl;
     669          //cout << "sWeiN = " << sWeiN << endl;
    576670         
    577671          //For 1 fb-1
    … …  
    586680          if (i==1){
    587681            h[1]->Fill(mt,WeiTot);
     682            //h[4]->Fill(mt,sWei);
    588683            twoDim_mt_HT_Sig->Fill(mt,HT,sWei);
    589684            twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);   
    … …  
    609704  for (int i = 1; i < 151; i++){
    610705    double sP = h[3]->GetBinContent(i);
    611     double B  = (sig2 - N1)*h[2]->GetBinContent(i);
    612706    double sB = h[5]->GetBinContent(i);
    613     double W  = sP - B;
    614707    double sW  = sP - sB;
    615     double Tot1fb = h[0]->GetBinContent(i);
    616     double BkgPure = Tot1fb - sP;
     708
     709    //double B  = (m_sig2 - m_N1)*h[2]->GetBinContent(i);
     710    //double W  = sP - B;
     711    //double Tot1fb = h[0]->GetBinContent(i);
     712    //double BkgPure = Tot1fb - sP;
    617713   
    618714    double vb1 = 10*(i-1);
    … …  
    620716    double SumErrBin = 0;
    621717    double SumErrNBin = 0;
    622     for (int ic = 0; ic < icount; ic++){
     718    for (int ic = 0; ic < m_icount; ic++){
    623719      double mt  = mtV[ic];
    624720      double errV  = err[ic];
    … …  
    630726    }
    631727    h[4]->SetBinContent(i,sW);
    632     h[4]->SetBinError(i,sqrt(SumErrBin));
    633     h[6]->SetBinError(i,sqrt(SumErrNBin));
     728    if (SumErrBin < 0) cout << "Warning negative error squared in bin " << i << " for sP" << endl;
     729    h[4]->SetBinError(i,sqrt(fabs(SumErrBin)));
     730    if (SumErrNBin < 0) cout << "Warning negative error squared in bin " << i << " for esP" << endl;
     731    h[6]->SetBinError(i,sqrt(fabs(SumErrNBin)));
    634732  }
    635733 
    … …  
    642740  for (int i = 0; i < 7; i++) {
    643741    //h[i]->Rebin(2);
    644     h[i]->SetLineWidth(3);
     742    h[i]->SetLineWidth(2);
    645743    TAxis *ax = h[i]->GetXaxis();
    646744    int b1 = ax->FindBin(x1);
    … …  
    650748
    651749  //plot for sPlot with M0 known
    652   TCanvas *csP = new TCanvas("csP", "", 800, 400);
     750  TCanvas *csP = new TCanvas("csP", "csP", 800, 400);
    653751  csP->Divide(2);
    654752  csP->cd(1);
    … …  
    659757  h[1]->Draw("same");
    660758  h[4]->SetLineColor(1);
    661   h[4]->SetLineWidth(3);
     759  h[4]->SetLineWidth(2);
    662760  h[4]->DrawCopy("hist,same");
    663761  csP->cd(2);
    … …  
    665763
    666764  //plot for sPlot with M0 UNknown
    667   TCanvas *cesP = new TCanvas("cesP", "", 800,400);
     765  TCanvas *cesP = new TCanvas("cesP", "cesP", 800,400);
    668766  cesP->Divide(2);
    669767  cesP->cd(1);
    … …  
    672770  h[6]->DrawCopy("E1");
    673771  h[6]->SetLineColor(1);
    674   h[6]->SetLineWidth(3);
     772  h[6]->SetLineWidth(2);
    675773  h[6]->Draw("hist,same");
    676774  cesP->cd(2);
Note: See TracChangeset for help on using the changeset viewer.