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/EFFICIENCY/EfficiencyDependOnCut.C

    r21 r22  
    1010#include <TLegend.h>
    1111#include <TROOT.h>
     12#include <TMinuit.h>
    1213
    1314#include <iostream>
    … …  
    1516using namespace std;
    1617
     18//Set the cut values
     19double Cut1 = 0.7;
     20double Cut2 = 0.7;
     21double Cut3 = 48;
     22double Cut4 = 250;
     23double Cut5 = 620;
     24
     25//Define class Event
     26class EventC {
     27public:
     28  EventC() { m_w = 1; m_fs = 1; m_fb = 0; m_cat = 1; }
     29  EventC(double w0, double fs0, double fb0) { m_w = w0; m_fs = fs0; m_fb = fb0; }
     30  EventC(double w0, int cat0) { m_w = w0; m_cat = cat0; m_fs = 1; m_fb = 0; }
     31  double w() { return m_w; }
     32  double fs() { return m_fs; }
     33  double fb() { return m_fb; }
     34  int cat() { return m_cat; }
     35  void fs(double fs0) { m_fs = fs0; }
     36  void fb(double fb0) { m_fb = fb0; }
     37
     38private:
     39  double m_w, m_fs, m_fb;
     40  int m_cat;
     41};
     42
     43//global variables
     44std::vector<EventC> m_EventList;
     45double m_F1[3][31][6];
     46double m_F0[3][31][6];
     47double m_N1exp[3][31];
     48double m_N0exp[3][31];
     49
     50void fcnN0N1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
     51  double L = 0;
     52  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
     53    EventC evt = m_EventList[iev];
     54    double smt = par[0]*evt.fb() + par[1]*evt.fs();
     55    if (smt > 0) L += evt.w()*log(smt);
     56  }
     57  f = par[0] + par[1] - L;
     58};
     59
     60vector<double> vFe(double cos1,double cos2,double dMLep, int i, int j)
     61{
     62  double f0e = 0;
     63  double f1e = 0;
     64 
     65  //cout << "cos1 = " << cos1 << " cos2 = " << cos2 << " dMLep = " << dMLep << endl;
     66 
     67  int it = -1;
     68  if (cos1<Cut1 && cos2<Cut2){
     69    if (dMLep<Cut3) it = 0;
     70    if (dMLep>Cut3) it = 1;
     71  }
     72  else if (cos1>Cut1 && cos2>Cut2){
     73    if (dMLep<Cut3) it = 2;
     74    if (dMLep>Cut3) it = 3;
     75  }
     76  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
     77    if (dMLep<Cut3) it = 4;
     78    if (dMLep>Cut3) it = 5;
     79  }
     80 
     81  f0e = m_F0[i][j][it];
     82  f1e = m_F1[i][j][it];
     83 
     84  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
     85 
     86  vector<double> vfe;
     87  vfe.push_back(f0e);
     88  vfe.push_back(f1e);
     89  return vfe;
     90}
     91
     92//Make list of events in Data file
     93void MakeEventListf0f1(int i, int j){
     94  m_EventList.clear(); 
     95
     96  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
     97  TTree *tD = (TTree*)fD->Get("T");
     98  int nD = tD->GetEntries();
     99  for (int iev = 0; iev < nD; iev++) {
     100    tD->GetEntry(iev);
     101   
     102    double mt    = tD->GetLeaf("mtophad")->GetValue();
     103    double HT    = tD->GetLeaf("HT")->GetValue();
     104    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     105    double y     = tD->GetLeaf("cosWqq")->GetValue();
     106    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     107
     108    int hasBJ = (int)tD->GetLeaf("hasBadJet")->GetValue();
     109    int RunN  = (int)tD->GetLeaf("RunNumber_ana")->GetValue();
     110   
     111    if (RunN >= 180614 && hasBJ==1) continue;
     112
     113    double cos1 = fabs(x);
     114    double cos2 = fabs(y);
     115   
     116    //sPlots total
     117    if (mt > Cut4 && HT > Cut5){
     118      vector<double> fsfb =  vFe(cos1,cos2,dMLep,i,j);
     119      double f0 = fsfb.at(0);
     120      double f1 = fsfb.at(1);
     121      //cout << "f0 = " << f0 << " f1 = " << f1 << endl;
     122      EventC eventC(1,f1,f0);
     123      m_EventList.push_back(eventC);
     124    }
     125  }
     126}
     127
     128//Fit the number N0 and N1 using the events from data files
     129std::pair<pair<double,double>,double> Fit(double N0_exp, double N1_exp, int i, int j){
     130  // Test
     131  MakeEventListf0f1(i,j);
     132  cout << "Number of events in the LH " << m_EventList.size() << endl;
     133
     134  int ifl = 0;
     135  double step = 0.01;
     136 
     137  double fitPar0 = N0_exp;
     138  double nmax0   = fitPar0 + 4*sqrt(N0_exp);
     139  double nmin0   = fitPar0 - 4*sqrt(N0_exp);
     140 
     141  double fitPar1 = N1_exp;
     142  double nmax1   = fitPar1 + 4*sqrt(N1_exp);
     143  double nmin1   =         - 4*sqrt(N1_exp);
     144 
     145  double n0best, err0best;
     146  double n1best, err1best;
     147 
     148  TMinuit* myMinuit = new TMinuit(2);
     149  //myMinuit->Command("SET PRI -1");
     150  myMinuit->Command("SET NOW");
     151  myMinuit->Command("SET ERR 0.5");
     152  myMinuit->Command("SET STR 2");   // try to improve mimimum search
     153  myMinuit->SetFCN(fcnN0N1);
     154  myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl);
     155  myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl);
     156  //myMinuit->FixParameter(0);
     157  myMinuit->Command("MIGRAD 500 1");
     158  myMinuit->GetParameter(0,n0best,err0best);
     159  myMinuit->GetParameter(1,n1best,err1best);
     160 
     161  std::cout << "Expected N0 = " << N0_exp << " fitted = " << n0best << " +- " << err0best << std::endl;
     162  std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
     163 
     164  double N0 = n0best;
     165  double N1 = n1best;
     166  double sig2 = err1best*err1best;
     167
     168  pair<double,double> N0N1 = make_pair<double,double> (N0,N1);
     169 
     170  return make_pair<pair<double,double>,double> (N0N1,sig2);
     171};
     172
     173
     174//##############################################################
     175//####################### MAIN PROGRAM #########################
     176//##############################################################
     177
    17178void EfficiencyDependOnCut(){
    18   TTree *t[3];
     179  TTree *t[4];
    19180
    20181  TFile *fB = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_AllBkg.root");
    … …  
    24185  t[1] = (TTree*)fS->Get("T");
    25186
     187  TFile *fBS = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_FakeData.root");
     188  t[2] = (TTree*)fBS->Get("T");
     189
    26190  TFile *fD = TFile::Open("/exp/atlas/huonglan/DATA2011/resData.root");
    27   t[2] = (TTree*)fD->Get("T");
    28 
    29   double Cut1 = 0.7;
    30   double Cut2 = 0.7;
    31   double Cut3 = 48;
    32 
    33   double   e1[3][3][31];
    34   double   e2[3][3][31];   
    35   double   e3[3][3][31];   
    36   double   e4[3][3][31];   
    37   double   e5[3][3][31];   
    38   double   e6[3][3][31];   
    39   double iTot[3][3][31];
     191  t[3] = (TTree*)fD->Get("T");
     192
     193  double   e1[3][4][31];
     194  double   e2[3][4][31];   
     195  double   e3[3][4][31];   
     196  double   e4[3][4][31];   
     197  double   e5[3][4][31];   
     198  double   e6[3][4][31];   
     199  double iTot[3][4][31];
    40200
    41201  for (int k = 0; k < 3; k++){
    42     for (int i = 0; i < 3; i++){
     202    for (int i = 0; i < 4; i++){
    43203      for (int j = 0; j < 31; j++){
    44204        iTot[k][i][j] =0 ;
    … …  
    56216  double mjjjErr[31]; for (int i = 0; i < 31; i++) { mjjjErr[i] = 0;}
    57217
    58   for (int i = 0; i < 3; i++){
     218  for (int i = 0; i < 4; i++){
    59219    int nevt = t[i]->GetEntries();
    60220    for (int iev = 0; iev < nevt; iev++) {
    … …  
    64224      double wei = 1;
    65225      double cos1 = 1, cos2 = 1;
     226      int hfor = -1;
    66227
    67228      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
    68229      double ht   = t[i]->GetLeaf("HT")->GetValue();
    69       if (i==0||i==1) {
     230      int hasBJ   = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
     231      if (i==0||i==1||i==2) {
    70232        wei  = t[i]->GetLeaf("wei")->GetValue();
    71233        cos1 = fabs(t[i]->GetLeaf("cos1")->GetValue());
    72234        cos2 = fabs(t[i]->GetLeaf("cos2")->GetValue());
    73       } else if (i==2){
     235        hfor = (int)t[i]->GetLeaf("hfor")->GetValue();
     236      } else if (i==3){
    74237        cos1 = fabs(t[i]->GetLeaf("cos_tlepCMtpair")->GetValue());
    75238        cos2 = fabs(t[i]->GetLeaf("cosWqq")->GetValue());
    76239      }
    77240      double dM   = t[i]->GetLeaf("MinLep1")->GetValue();
     241
     242      if (hasBJ!=0 || hfor==4) continue;
    78243
    79244      int it = -1;
    … …  
    110275  }
    111276
    112   double    eB[3][31][6],    eS[3][31][6],    eD[3][31][6];
    113   double eBErr[3][31][6], eSErr[3][31][6], eDErr[3][31][6];
     277  double    eB[3][31][6],    eS[3][31][6],    eBS[3][31][6],    eD[3][31][6];
     278  double eBErr[3][31][6], eSErr[3][31][6], eBSErr[3][31][6], eDErr[3][31][6];
    114279
    115280  for (int i = 0; i < 3; i++){
    … …  
    136301      if (iTot[i][1][j]==0) { for (int k = 0; k < 6; k++) { eS[i][j][k] = 0; } }
    137302     
    138       eD[i][j][0] = e1[i][2][j]/iTot[i][2][j];
    139       eD[i][j][1] = e2[i][2][j]/iTot[i][2][j];
    140       eD[i][j][2] = e3[i][2][j]/iTot[i][2][j];
    141       eD[i][j][3] = e4[i][2][j]/iTot[i][2][j];
    142       eD[i][j][4] = e5[i][2][j]/iTot[i][2][j];
    143       eD[i][j][5] = e6[i][2][j]/iTot[i][2][j];
    144 
    145       if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
     303      eBS[i][j][0] = e1[i][2][j]/iTot[i][2][j];
     304      eBS[i][j][1] = e2[i][2][j]/iTot[i][2][j];
     305      eBS[i][j][2] = e3[i][2][j]/iTot[i][2][j];
     306      eBS[i][j][3] = e4[i][2][j]/iTot[i][2][j];
     307      eBS[i][j][4] = e5[i][2][j]/iTot[i][2][j];
     308      eBS[i][j][5] = e6[i][2][j]/iTot[i][2][j];
     309
     310      if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eBS[i][j][k] = 0; } }
     311     
     312      eD[i][j][0] = e1[i][3][j]/iTot[i][3][j];
     313      eD[i][j][1] = e2[i][3][j]/iTot[i][3][j];
     314      eD[i][j][2] = e3[i][3][j]/iTot[i][3][j];
     315      eD[i][j][3] = e4[i][3][j]/iTot[i][3][j];
     316      eD[i][j][4] = e5[i][3][j]/iTot[i][3][j];
     317      eD[i][j][5] = e6[i][3][j]/iTot[i][3][j];
     318
     319      if (iTot[i][3][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
    146320     
    147321      //cout << "e1S[" << i << "][" << j << "] = " << e1S[i][j];
    … …  
    150324  }
    151325 
     326  //3 is the number of regions
     327  //31 is the number of the cut for mjjj
     328  //6 is the number of categories
    152329  for (int i = 0; i < 3; i++){
    153330    for (int j = 0; j < 31; j++){
    154331      for (int k = 0; k < 6; k++){
    155         eBErr[i][j][k] = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
    156         eSErr[i][j][k] = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
    157         eDErr[i][j][k] = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][2][j]);
    158       }
     332        eBErr[i][j][k]  = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
     333        eSErr[i][j][k]  = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
     334        eBSErr[i][j][k] = sqrt(eBS[i][j][k]*(1-eBS[i][j][k])/iTot[i][2][j]);
     335        eDErr[i][j][k]  = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][3][j]);
     336
     337        m_F1[i][j][k] = eS[i][j][k];
     338        m_F0[i][j][k] = eB[i][j][k];
     339        m_N1exp[i][j] = iTot[i][1][j];
     340        m_N0exp[i][j] = iTot[i][0][j];
     341
     342        cout << "region " << i << " cut mjjj " << mjjj[j] << " cate " << k+1 ;
     343        cout << " iTotB = " << iTot[i][0][j] ;
     344        cout << " iTotS = " << iTot[i][1][j] << endl;
     345      }
     346    }
     347  }
     348
     349  //3 is the number of regions
     350  //31 is the number of the cut for mjjj
     351  double N0fit[3][31], N1fit[3][31], N1err[3][31];
     352  for (int i = 0; i < 3; i++){
     353    for (int j = 0; j < 31; j++){
     354      pair<pair<double,double>,double> N0N1sig2 = Fit(m_N0exp[i][j],m_N1exp[i][j],i,j);
     355      pair<double,double> N0N1 = N0N1sig2.first;
     356      double N0   = N0N1.first;
     357      double N1   = N0N1.second;
     358      double sig2 = N0N1sig2.second;
     359      N0fit[i][j] = N0;
     360      N1fit[i][j] = N1;
     361      N1err[i][j] = sig2;
    159362    }
    160363  }
    … …  
    175378  TString cat[6] = {"cat1","cat2","cat3","cat4","cat5","cat6"};
    176379
    177   TGraphErrors *heffB[3][6];
    178   TGraphErrors *heffS[3][6];
    179   TGraphErrors *heffD[3][6];
    180 
     380//   TGraphErrors *heffB[3][6];
     381//   TGraphErrors *heffS[3][6];
     382//   TGraphErrors *heffBS[3][6];
     383//   TGraphErrors *heffD[3][6];
     384
     385  TGraph *heffB[3][6];
     386  TGraph *heffS[3][6];
     387  TGraph *heffBS[3][6];
     388  TGraph *heffD[3][6];
     389
     390  //*************Define histograms on N1************
     391  TGraphErrors *hN1[3];
     392  for (int i = 0; i < 3; i++){
     393    double N1fittemp[31], N1errtemp[31];
     394    for (int j = 0; j < 31; j++){
     395      N1fittemp[j] = N1fit[i][j];
     396      N1errtemp[j] = N1err[i][j];
     397      cout << "mjjj = " << mjjj[j] << endl;
     398      cout << " m_N1exp[" << i << "][" << j << "]= " << m_N1exp[i][j] ;
     399      cout << " m_N0exp[" << i << "][" << j << "]= " << m_N0exp[i][j] ;
     400      cout << " N1fit[" << i << "][" << j << "]= " << N1fit[i][j] ;
     401      cout << " N0fit[" << i << "][" << j << "]= " << N0fit[i][j] ;
     402      cout << " N1err[" << i << "][" << j << "]= " << N1err[i][j] << endl;
     403    }
     404    hN1[i] = new TGraphErrors(31,mjjj,N1fittemp,mjjjErr,N1errtemp);
     405  }
     406
     407  //*************Define histograms on efficiencies****************
    181408  for (int i = 0; i < 3; i++){
    182409    for (int j = 0; j < 6; j++){
    … …  
    189416        effErr[ic] = eBErr[i][ic][j];
    190417      }
    191       heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
    192       heffB[i][j]->SetLineColor(4);
     418      //heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     419      heffB[i][j] = new TGraph(31, mjjj, eff);
     420      heffB[i][j]->SetLineColor(8);
     421      heffB[i][j]->SetLineWidth(2);
    193422     
    194423      name = "heffS"; name += region[i]; name += cat[j];
    … …  
    197426        effErr[ic] = eSErr[i][ic][j];
    198427      }
    199       heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     428      //heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     429      heffS[i][j] = new TGraph(31, mjjj, eff);
    200430      heffS[i][j]->SetLineColor(2);
     431      heffS[i][j]->SetLineWidth(2);
     432     
     433      name = "heffBS"; name += region[i]; name += cat[j];
     434      for (int ic = 0; ic < 31; ic++) {
     435        eff[ic]    = eBS[i][ic][j];
     436        effErr[ic] = eBSErr[i][ic][j];
     437      }
     438      //heffBS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     439      heffBS[i][j] = new TGraph(31, mjjj, eff);
     440      heffBS[i][j]->SetLineColor(kMagenta);
     441      heffBS[i][j]->SetLineWidth(2);
    201442     
    202443      name = "heffD"; name += region[i]; name += cat[j];
    … …  
    205446        effErr[ic] = eDErr[i][ic][j];
    206447      }
    207       heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     448      //heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
     449      heffD[i][j] = new TGraph(31, mjjj, eff);
     450      heffD[i][j]->SetLineWidth(2);
    208451    }
    209452  }
    … …  
    211454  legend->AddEntry(heffB[0][1],"BkgMC");
    212455  legend->AddEntry(heffS[0][1],"SigMC");
     456  legend->AddEntry(heffBS[0][1],"Sig+BkgMC");
    213457  legend->AddEntry(heffD[0][1],"Data");
     458
     459  //****************** DRAWING *********************
     460
     461  TCanvas *cN1 = new TCanvas("cN1", "cN1", 1000, 400);
     462  cN1->Divide(3);
     463 
     464  cN1->cd(1);
     465  hN1[0]->Draw("ALP");
     466  cN1->cd(2);
     467  hN1[1]->Draw("ALP");
     468  cN1->cd(3);
     469  hN1[2]->Draw("ALP");
     470
     471  cN1->Print("N1fit.eps");
    214472
    215473  TCanvas *c[6];
    … …  
    223481    frame1->SetTitle(name);
    224482    frame1->Draw();
    225     heffB[0][i]->Draw("LP");
    226     heffS[0][i]->Draw();
     483    heffS[0][i]->Draw("LP");
    227484    heffD[0][i]->Draw();
     485    heffBS[0][i]->Draw();
     486    heffB[0][i]->Draw();
    228487    legend->Draw("same");
    229488
    … …  
    232491    frame2->SetTitle(name);
    233492    frame2->Draw();
    234     heffB[1][i]->Draw("LP");
    235     heffS[1][i]->Draw();
     493    heffS[1][i]->Draw("LP");
    236494    heffD[1][i]->Draw();
     495    heffBS[1][i]->Draw();
     496    heffB[1][i]->Draw();
    237497    legend->Draw("same");
    238498
    … …  
    241501    frame3->SetTitle(name);
    242502    frame3->Draw();
    243     heffB[2][i]->Draw("LP");
    244     heffS[2][i]->Draw();
     503    heffS[2][i]->Draw("LP");
    245504    heffD[2][i]->Draw();
     505    heffBS[2][i]->Draw();
     506    heffB[2][i]->Draw();
    246507    legend->Draw("same");
    247508
    … …  
    250511  }
    251512
    252   //********************************************
    253 
    254   //TCanvas *cN = new TCanvas("cN", "", 500, 400);
    255   //TGraph *hB = new TGraph(31, mjjj, N1);
    256   //hB->Draw("ALP");
    257513}
Note: See TracChangeset for help on using the changeset viewer.