Changeset 19 in huonglan


Ignore:
Timestamp:
Aug 30, 2011, 10:24:25 AM (13 years ago)
Author:
huonglan
Message:

sPlots for MC in the same file now, MakeFakeData with weight for trees contains scale factor for flavour fractions

Location:
SPLOT/NEWREALEASE
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • SPLOT/NEWREALEASE/sPlots_Data.C

    r17 r19  
    8181  TFile *file[2];
    8282  TTree *t[2];
    83   TTree *tw[2];
    8483 
    8584  for (int i = 0; i < 2; i++){
    86     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     85    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    8786    file[i] = TFile::Open(filename);
    8887    t[i]    = (TTree*)file[i]->Get("T");
    89     tw[i]   = (TTree*)file[i]->Get("Tw");
    90     t[i]->AddFriend(tw[i]);
    9188  }
    9289
    … …  
    104101      t[i]->GetEntry(iev);
    105102           
    106       double mt    = t[i]->GetLeaf("mtophad")->GetValue();
    107       double HT    = t[i]->GetLeaf("HT")->GetValue();
    108       double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    109       double y     = t[i]->GetLeaf("cosWqq")->GetValue();
    110       double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    111       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    112       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    113       double WeiTot = WeiEvt*Wei1fb;
    114 
    115       //cout << "WeiEvt = " << WeiEvt << endl;
    116       //cout << "Wei1fb = " << Wei1fb << endl;
     103      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
     104      double HT     = t[i]->GetLeaf("HT")->GetValue();
     105      double x      = t[i]->GetLeaf("cos1")->GetValue();
     106      double y      = t[i]->GetLeaf("cos2")->GetValue();
     107      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
     108      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     109
    117110      //cout << "WeiTot = " << WeiTot << endl;
    118111
    … …  
    327320    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
    328321    double y     = tD->GetLeaf("cosWqq")->GetValue();
    329     double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     322    double dMLep = tD->GetLeaf("MinLep1")->GetValue();           
    330323     
    331324    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
    … …  
    383376  cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl;
    384377
     378  //TH1F *hMTSig   = new TH1F("hMTSig",     "", 35, 250, 600);
     379  //TH1F *hMTBkg   = new TH1F("hMTBkg",     "", 35, 250, 600);
    385380  TH1F *hMTSig   = new TH1F("hMTSig",     "", 25, 0, 1500);
    386381  TH1F *hMTBkg   = new TH1F("hMTBkg",     "", 25, 0, 1500);
    387   TH1F *hHTSig   = new TH1F("hHTSig",     "", 25, 0, 2000);
    388   TH1F *hHTBkg   = new TH1F("hHTBkg",     "", 25, 0, 2000);
    389   TH2F *hHTMTSig = new TH2F("hHTMTSig",   "", 25, 0, 2000, 25, 0, 1500);
    390   TH2F *hHTMTBkg = new TH2F("hHTMTBkg",   "", 25, 0, 2000, 25, 0, 1500);
     382
     383  TH1F *hHTSig   = new TH1F("hHTSig",     "", 50, 500, 1500);
     384  TH1F *hHTBkg   = new TH1F("hHTBkg",     "", 50, 500, 1500);
     385
     386  TH2F *hMTHTSig = new TH2F("hMTHTSig",   "", 35, 250, 600, 50, 500, 1500);
     387  TH2F *hMTHTBkg = new TH2F("hMTHTBkg",   "", 35, 250, 600, 50, 500, 1500);
    391388
    392389  TH1F *hmt[6];
    … …  
    449446      hHTSig->Fill(HT,sP1);
    450447      hHTBkg->Fill(HT,sP0);
    451       hHTMTSig->Fill(HT,mt,sP1);
    452       hHTMTBkg->Fill(HT,mt,sP0);
     448      hMTHTSig->Fill(mt,HT,sP1);
     449      hMTHTBkg->Fill(mt,HT,sP0);
    453450      sPN[icc] = sP1*sP1;
    454451      mtN[icc] = mt;
    … …  
    474471  cout << "hMTBkg integral  " << hMTBkg->Integral(0,1500)  << endl;
    475472
    476 //   for (int i = 0; i < Nsel; i++){
    477 //     cout << "mt[" << i << "] = " << mtN[i] << " ";
    478 //     cout << "sP[" << i << "] = " << sPN[i] << endl;
     473//   for (int i = 1; i < 36; i++){
     474//     double vb1 = 250 +(i-1)*10;
     475//     double vb2 = 250 + i*10;
     476//     double ErrBin = 0.;
     477//     for (int j = 0; j < Nsel; j++){
     478//       double mt = mtN[j];
     479//       double sp = sPN[j];
     480//       if (mt >= vb1 && mt < vb2){
     481//      ErrBin += sp;
     482//       }
     483//     }
     484//     hMTSig->SetBinError(i,sqrt(ErrBin));
    479485//   }
    480  
     486
    481487  for (int i = 1; i < 26; i++){
    482488    double vb1 = (i-1)*60;
    … …  
    493499  }
    494500
    495   for (int i = 1; i < 26; i++){
    496     double vb1 = (i-1)*80;
    497     double vb2 = i*80;
     501  for (int i = 1; i < 51; i++){
     502    double vb1 = 500 + (i-1)*20;
     503    double vb2 = 500 + i*20;
    498504    double ErrBin = 0.;
    499505    for (int j = 0; j < Nsel; j++){
    … …  
    510516  c1->Divide(2);
    511517  c1->cd(1);
     518  //hMTSig->Rebin(5);
    512519  hMTSig->SetLineColor(2);
    513520  hMTSig->DrawCopy("E1");
    514521  hMTSig->SetLineColor(1);
     522  hMTSig->SetLineWidth(2);
    515523  hMTSig->Draw("hist,same");
    516524  c1->cd(2);
     525  //hMTBkg->Rebin(5);
     526  hMTBkg->SetLineWidth(2);
    517527  hMTBkg->Draw();
    518528
    … …  
    530540  c3->Divide(2);
    531541  c3->cd(1);
    532   hHTMTSig->Draw("colz");
     542  hMTHTSig->Draw("colz");
    533543  c3->cd(2);
    534   hHTMTBkg->Draw("colz");
     544  hMTHTBkg->Draw("colz");
    535545
    536546  TCanvas *c4 = new TCanvas("c4", "", 1000, 600);
    … …  
    547557  }
    548558}
     559
     560//##############################################################
     561//########### N0 known #########################################
     562//##############################################################
     563
     564void sPlots_Data_N0known(){
     565  //gROOT->SetStyle("Plain");
     566  //Make the table of sP values first
     567 
     568  //Get right N0 and N1
     569  //Fit();
     570  CalF0F1();
     571  MakeEventListf0f1();
     572  sig2 = sigmaMin*sigmaMin*N1*N1;
     573 
     574  TString filename = "";
     575  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
     576  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
     577  TFile *file[2];
     578  TTree *t[2];
     579 
     580  for (int i = 0; i < 2; i++){
     581    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
     582    file[i] = TFile::Open(filename);
     583    t[i]    = (TTree*)file[i]->Get("T");
     584  }
     585 
     586  for (int i = 0; i < 2; i++){
     587    iTot[i] = 0;
     588  }
     589
     590  TH1F *hBkg1fb = new TH1F("hBkg1fb", "", 35, 250, 600);
     591
     592  //start reading 2 files
     593  for (int i = 0; i < 2; i++){
     594    int nevt = t[i]->GetEntries();
     595    for (int iev = 0; iev < nevt; iev++) {
     596      t[i]->GetEntry(iev);
     597     
     598      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
     599      double HT     = t[i]->GetLeaf("HT")->GetValue();
     600      double x      = t[i]->GetLeaf("cos1")->GetValue();
     601      double y      = t[i]->GetLeaf("cos2")->GetValue();
     602      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
     603      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     604     
     605      double cos1 = fabs(x);
     606      double cos2 = fabs(y);
     607     
     608      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
     609      double hf = t[i]->GetLeaf("hfor")->GetValue();
     610     
     611      if (hf!=4 && hasBJ==0){
     612        if (mt > Cut4 && HT > Cut5){
     613          iTot[i] += WeiTot;
     614          if (i==0) hBkg1fb->Fill(mt,WeiTot);
     615        }
     616      }
     617    }//Read each file
     618  }//End of reading 2 files
     619
     620  hBkg1fb->Scale(1./hBkg1fb->Integral(0,1000)); 
     621
     622  N1 = iTot[1];
     623  N0 = iTot[0];
     624
     625  cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
     626 
     627  TH1F *hsPTot  = new TH1F("hsPTot" , "",  35, 250,  600);
     628  TH1F *hsP     = new TH1F("hsP" ,    "",  35, 250,  600);
     629  TH1F *hesP    = new TH1F("hesP",    "",  35, 250,  600);
     630 
     631  int icount = 0;
     632  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
     633  TTree *tD = (TTree*)fD->Get("T");
     634  int nD = tD->GetEntries();
     635  for (int iev = 0; iev < nD; iev++) {
     636    tD->GetEntry(iev);
     637   
     638    double mt    = tD->GetLeaf("mtophad")->GetValue();
     639    double HT    = tD->GetLeaf("HT")->GetValue();
     640    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     641    double y     = tD->GetLeaf("cosWqq")->GetValue();
     642    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     643     
     644    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
     645    int RunN  = tD->GetLeaf("RunNumber_ana")->GetValue();
     646   
     647    if (RunN >= 180614 && hasBJ==1) continue;
     648
     649    double cos1 = fabs(x);
     650    double cos2 = fabs(y);
     651   
     652    //sPlots total
     653    if (mt > Cut4 && HT > Cut5){
     654      icount++;
     655    }
     656  }
     657
     658  cout << "icount = " << icount << endl;
     659  double k = (sig2 - N1)/(N0 - (sig2 - N1));
     660 
     661  double *mtV  = new double[icount];
     662  double *err  = new double[icount];
     663  double *errN = new double[icount];
     664  int icc = 0;
     665 
     666  //Try standard sPlots
     667  for (int iev = 0; iev < nD; iev++) {
     668    tD->GetEntry(iev);
     669   
     670    double mt    = tD->GetLeaf("mtophad")->GetValue();
     671    double HT    = tD->GetLeaf("HT")->GetValue();
     672    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     673    double y     = tD->GetLeaf("cosWqq")->GetValue();
     674    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     675     
     676    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
     677    int RunN  = tD->GetLeaf("RunNumber_ana")->GetValue();
     678   
     679    if (RunN >= 180614 && hasBJ==1) continue;
     680
     681    double cos1 = fabs(x);
     682    double cos2 = fabs(y);
     683   
     684    //sPlots total
     685    if (mt > Cut4 && HT > Cut5){
     686      vector<double> veff = vFe(cos1,cos2,dMLep);
     687      double F0e = veff.at(0);
     688      double F1e = veff.at(1);
     689         
     690      double sWei  = sig2*F1e/(N1*F1e + N0*F0e);
     691      double sWeiN = sWei*(k+1) - k;
     692
     693      //M0 distribution unknown
     694      hsPTot->Fill(mt,sWei);
     695      hesP->Fill(mt,sWeiN);
     696      mtV[icc] = mt;
     697      err[icc] = sWei*sWei;
     698      errN[icc] = sWeiN*sWeiN;
     699      icc++;
     700    }//Read each file
     701  }//End of reading 2 files
     702 
     703  for (int i = 1; i < 36; i++){
     704    double sPTot = hsPTot->GetBinContent(i);
     705    double B     = (sig2-N1)*hBkg1fb->GetBinContent(i);
     706    double sPSig = sPTot - B;
     707   
     708    double vb1 = 250 + 10*(i-1);
     709    double vb2 = 250 + 10*i;
     710    double SumErrBin  = 0;
     711    double SumErrNBin = 0;
     712    for (int ic = 0; ic < icount; ic++){
     713      double mt    = mtV[ic];
     714      double errV  = err[ic];
     715      double errNV = errN[ic];
     716      if (mt>=vb1 && mt<vb2){
     717        SumErrBin += errV;
     718        SumErrNBin += errNV;
     719      }
     720    }
     721    hsP->SetBinContent(i,sPSig);
     722    hsP->SetBinError(i,sqrt(SumErrBin));
     723    hesP->SetBinError(i,sqrt(SumErrNBin));
     724  }
     725
     726  //plot for sPlot with M0 UNknown
     727  TCanvas *c = new TCanvas("c", "", 800,400);
     728  c->Divide(2);
     729  c->cd(1);
     730  hsP->SetLineColor(2);
     731  hsP->SetLineWidth(1);
     732  hsP->DrawCopy("E1");
     733  hsP->SetLineColor(1);
     734  hsP->SetLineWidth(2);
     735  hsP->Draw("hist,same");
     736  c->cd(2);
     737  hesP->SetLineColor(2);
     738  hesP->SetLineWidth(1);
     739  hesP->DrawCopy("E1");
     740  hesP->SetLineColor(1);
     741  hesP->SetLineWidth(2);
     742  hesP->Draw("hist,same");
     743}

  • SPLOT/NEWREALEASE/sPlots_NewRel.C

    r17 r19  
    1919double Cut5 = 620;
    2020double sigmaMin;
     21int icount = 0;
    2122double iTot[2];
    2223double iCatS[6], iCatB[6];
    … …  
    8182  TFile *file[2];
    8283  TTree *t[2];
    83   TTree *tw[2];
    84  
    85   for (int i = 0; i < 2; i++){
    86     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     84 
     85  for (int i = 0; i < 2; i++){
     86    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    8787    file[i] = TFile::Open(filename);
    8888    t[i]    = (TTree*)file[i]->Get("T");
    89     tw[i]   = (TTree*)file[i]->Get("Tw");
    90     t[i]->AddFriend(tw[i]);
    9189  }
    9290
    … …  
    104102      t[i]->GetEntry(iev);
    105103           
    106       double mt    = t[i]->GetLeaf("mtophad")->GetValue();
    107       double HT    = t[i]->GetLeaf("HT")->GetValue();
    108       double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    109       double y     = t[i]->GetLeaf("cosWqq")->GetValue();
    110       double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    111       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    112       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    113       double WeiTot = WeiEvt*Wei1fb;
    114 
    115       //cout << "WeiEvt = " << WeiEvt << endl;
    116       //cout << "Wei1fb = " << Wei1fb << endl;
     104      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
     105      double HT     = t[i]->GetLeaf("HT")->GetValue();
     106      double x      = t[i]->GetLeaf("cos1")->GetValue();
     107      double y      = t[i]->GetLeaf("cos2")->GetValue();
     108      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
     109      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     110
    117111      //cout << "WeiTot = " << WeiTot << endl;
    118112
    … …  
    126120     
    127121        if (mt > Cut4 && HT > Cut5) {
     122          icount++;
    128123          iTot[i] += WeiTot;
    129124          int it = -1;
    … …  
    218213  TFile *file[2];
    219214  TTree *t[2];
    220   TTree *tw[2];
    221  
    222   for (int i = 0; i < 2; i++){
    223     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     215 
     216  for (int i = 0; i < 2; i++){
     217    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    224218    file[i] = TFile::Open(filename);
    225219    t[i]    = (TTree*)file[i]->Get("T");
    226     tw[i]   = (TTree*)file[i]->Get("Tw");
    227     t[i]->AddFriend(tw[i]);
    228220  }
    229221 
    … …  
    235227      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
    236228      double HT     = t[i]->GetLeaf("HT")->GetValue();
    237       double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    238       double y      = t[i]->GetLeaf("cosWqq")->GetValue();
     229      double x      = t[i]->GetLeaf("cos1")->GetValue();
     230      double y      = t[i]->GetLeaf("cos2")->GetValue();
    239231      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    240       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    241       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    242       double WeiTot = WeiEvt*Wei1fb;
    243      
    244       //cout << "WeiEvt = " << WeiEvt << endl;
    245       //cout << "Wei1fb = " << Wei1fb << endl;
     232      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     233     
    246234      //cout << "WeiTot = " << WeiTot << endl;
    247235     
    … …  
    269257Fit(){
    270258  // Test
    271   CalF0F1();
    272259  MakeEventListf0f1();
    273260
    … …  
    318305  cout << "Cut5 = " << Cut5 << endl;
    319306
    320   //Get right N0 and N1
     307  CalF0F1();
     308  cout << "icount = " << icount << endl;
     309
     310  //sPlots N0 known
     311  cout << "***********************************************" << endl;
     312  cout << "************** sPlot N0 known *****************" << endl;
     313  cout << "***********************************************" << endl;
     314  sPlots_NewRel_N0known();
     315
     316  //Get right N0 and N1 for standard sPlot
     317  cout << "***********************************************" << endl;
     318  cout << "************* sPlot N0 unknown ****************" << endl;
     319  cout << "***********************************************" << endl;
     320
    321321  Fit();
    322322 
    … …  
    326326  TFile *file[2];
    327327  TTree *t[2];
    328   TTree *tw[2];
    329  
    330   for (int i = 0; i < 2; i++){
    331     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     328 
     329  for (int i = 0; i < 2; i++){
     330    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    332331    file[i] = TFile::Open(filename);
    333332    t[i]    = (TTree*)file[i]->Get("T");
    334     tw[i]   = (TTree*)file[i]->Get("Tw");
    335     t[i]->AddFriend(tw[i]);
    336333  }
    337334
    338335  double VI00 = 0, VI01 = 0, VI10 = 0, VI11 = 0;
    339336  //start reading 2 files
    340   int icount = 0;
     337  //int icount = 0;
    341338  for (int i = 0; i < 2; i++){
    342339    int nevt = t[i]->GetEntries();
    … …  
    346343      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
    347344      double HT     = t[i]->GetLeaf("HT")->GetValue();
    348       double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    349       double y      = t[i]->GetLeaf("cosWqq")->GetValue();
     345      double x      = t[i]->GetLeaf("cos1")->GetValue();
     346      double y      = t[i]->GetLeaf("cos2")->GetValue();
    350347      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    351       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    352       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    353       double WeiTot = WeiEvt*Wei1fb;
     348      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    354349
    355350      double cos1 = fabs(x);
    … …  
    361356      if (hf!=4 && hasBJ==0){
    362357        if (mt > Cut4 && HT > Cut5){
    363           icount++;
    364358          //Compute the matrix
    365359          vector<double> veff = vFe(cos1,cos2,dMLep);
    … …  
    413407  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 35, 250, 600, 50, 500, 1500);
    414408
    415   cout << "icount = " << icount << endl;
    416 
     409  //Try standard sPlots
    417410  double SumsP1 = 0.;
    418411  double SumsP0 = 0.;
    … …  
    420413  double *mtN = new double[icount];
    421414  int icc = 0;
    422   //Try standard sPlots
    423415  for (int i = 0; i < 2; i++){
    424416    int nevt = t[i]->GetEntries();
    … …  
    428420      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
    429421      double HT    = t[i]->GetLeaf("HT")->GetValue();
    430       double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    431       double y     = t[i]->GetLeaf("cosWqq")->GetValue();
     422      double x     = t[i]->GetLeaf("cos1")->GetValue();
     423      double y     = t[i]->GetLeaf("cos2")->GetValue();
    432424      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    433       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    434       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    435       double WeiTot = WeiEvt*Wei1fb;
     425      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    436426     
    437427      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    … …  
    476466 
    477467  for (int i = 1; i < 36; i++){
    478     double vb1 = 250 + (i-1)*16;
    479     double vb2 = 250 + i*16;
     468    double vb1 = 250 + (i-1)*10;
     469    double vb2 = 250 + i*10;
    480470    double ErrBin = 0.;
    481471    for (int j = 0; j < icount; j++){
    … …  
    489479  }
    490480
     481  TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
     482
    491483  TCanvas *c = new TCanvas("c", "", 900, 400);
    492484  c->Divide(2);
    493485  c->cd(1);
     486  //frame->Draw();
    494487  hSigN->SetLineColor(2);
    495488  hSigN->SetLineWidth(1);
    496   //hSigN->DrawCopy("E1");
     489  hSigN->DrawCopy("E1");
    497490  hSigN->SetLineColor(1);
    498491  hSigN->SetLineWidth(2.5);
    499   hSigN->DrawCopy("hist");
     492  hSigN->DrawCopy("hist,same");
    500493  c->cd(2);
    501494  hBkgN->SetLineColor(4);
    … …  
    510503  h2DBkg->Draw("colz");
    511504}
     505
     506
     507//##################################################################
     508//########### FUNCTION FOR BACKGROUND YIELD KNOWN ##################
     509//##################################################################
     510
     511void 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));
     520 
     521  TH2F *twoDim_mt_HT_Sig     = new TH2F("mt_HT_Sig",     "", 75, 250, 1500, 50, 500, 2000);
     522  TH2F *twoDim_mt_HT_Sig_esP = new TH2F("mt_HT_Sig_esP", "", 75, 250, 1500, 50, 500, 2000);
     523 
     524  TH1F *h[7];
     525  h[0]  = new TH1F("All1fb",       "",  150,  0,  1500);
     526  h[1]  = new TH1F("Sig1fb",       "",  150,  0,  1500);
     527  h[2]  = new TH1F("Bkg1fb",       "",  150,  0,  1500);
     528  h[3]  = new TH1F("sPlotAll",     "",  150,  0,  1500);
     529  h[4]  = new TH1F("sPlotSig",     "",  150,  0,  1500);
     530  h[5]  = new TH1F("sPlotBkg",     "",  150,  0,  1500);
     531  h[6]  = new TH1F("esPlotSig",    "",  150,  0,  1500);
     532 
     533  double *mtV  = new double[icount];
     534  double *err  = new double[icount];
     535  double *errN = new double[icount];
     536  int icc = 0;
     537 
     538  TString filename = "";
     539  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
     540  TFile *file[2];
     541  TTree *t[2];
     542  for (int i = 0; i < 2; i++){
     543    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
     544    file[i] = TFile::Open(filename);
     545    t[i]    = (TTree*)file[i]->Get("T");
     546  }
     547
     548  for (int i = 0; i < 2; i++){
     549    int nevt = t[i]->GetEntries();
     550    for (int iev = 0; iev < nevt; iev++) {
     551      t[i]->GetEntry(iev);
     552     
     553      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
     554      double HT     = t[i]->GetLeaf("HT")->GetValue();
     555      double x      = t[i]->GetLeaf("cos1")->GetValue();
     556      double y      = t[i]->GetLeaf("cos2")->GetValue();
     557      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
     558      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     559
     560      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
     561      double hf = t[i]->GetLeaf("hfor")->GetValue();
     562     
     563      if (hf!=4 && hasBJ==0){
     564       
     565        double cos1 = fabs(x);
     566        double cos2 = fabs(y);
     567       
     568        //sPlots total
     569        if (mt > Cut4 && HT > Cut5){
     570          vector<double> veff = vFe(cos1,cos2,dMLep);
     571          double F0e = veff.at(0);
     572          double F1e = veff.at(1);
     573         
     574          double sWei  = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
     575          double sWeiN = sWei*(k+1) - k*WeiTot;
     576         
     577          //For 1 fb-1
     578          h[0]->Fill(mt,WeiTot);
     579          h[3]->Fill(mt,sWei);//sPlot
     580          h[6]->Fill(mt,sWeiN);//esPlot
     581
     582          if (i==0) {
     583            h[2]->Fill(mt,WeiTot);
     584            h[5]->Fill(mt,sWei);
     585          }
     586          if (i==1){
     587            h[1]->Fill(mt,WeiTot);
     588            twoDim_mt_HT_Sig->Fill(mt,HT,sWei);
     589            twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);   
     590          }
     591         
     592          //uncertainty
     593          mtV[icc] = mt;
     594          if (WeiTot==0){
     595            err[icc] = 0;
     596            errN[icc] = 0;
     597          } else {
     598            err[icc] = sWei*sWei/WeiTot;
     599            errN[icc] = sWeiN*sWeiN/WeiTot;
     600          }
     601          icc++;
     602        }
     603      }
     604    }//Read each file
     605  }//End of reading 2 files
     606 
     607  h[2]->Scale(1./h[2]->Integral(0,1500));
     608
     609  for (int i = 1; i < 151; i++){
     610    double sP = h[3]->GetBinContent(i);
     611    double B  = (sig2 - N1)*h[2]->GetBinContent(i);
     612    double sB = h[5]->GetBinContent(i);
     613    double W  = sP - B;
     614    double sW  = sP - sB;
     615    double Tot1fb = h[0]->GetBinContent(i);
     616    double BkgPure = Tot1fb - sP;
     617   
     618    double vb1 = 10*(i-1);
     619    double vb2 = 10*i;
     620    double SumErrBin = 0;
     621    double SumErrNBin = 0;
     622    for (int ic = 0; ic < icount; ic++){
     623      double mt  = mtV[ic];
     624      double errV  = err[ic];
     625      double errNV = errN[ic];
     626      if (mt>=vb1 && mt<vb2){
     627        SumErrBin += errV;
     628        SumErrNBin += errNV;
     629      }
     630    }
     631    h[4]->SetBinContent(i,sW);
     632    h[4]->SetBinError(i,sqrt(SumErrBin));
     633    h[6]->SetBinError(i,sqrt(SumErrNBin));
     634  }
     635 
     636  cout << "h2 integral is " << h[2]->Integral(0,1500) << endl;
     637  cout << "h4 integral is " << h[4]->Integral(0,1000) << endl;
     638  cout << "h6 integral is " << h[6]->Integral(0,1500) << endl;
     639
     640  double x1 = 250;
     641  double x2 = 600 - 0.01;
     642  for (int i = 0; i < 7; i++) {
     643    //h[i]->Rebin(2);
     644    h[i]->SetLineWidth(3);
     645    TAxis *ax = h[i]->GetXaxis();
     646    int b1 = ax->FindBin(x1);
     647    int b2 = ax->FindBin(x2);
     648    ax->SetRange(b1,b2);
     649  }
     650
     651  //plot for sPlot with M0 known
     652  TCanvas *csP = new TCanvas("csP", "", 800, 400);
     653  csP->Divide(2);
     654  csP->cd(1);
     655  h[4]->SetLineColor(2);
     656  h[4]->SetLineWidth(1);
     657  h[4]->DrawCopy("E1");
     658  h[1]->SetLineColor(kGreen+1);
     659  h[1]->Draw("same");
     660  h[4]->SetLineColor(1);
     661  h[4]->SetLineWidth(3);
     662  h[4]->DrawCopy("hist,same");
     663  csP->cd(2);
     664  twoDim_mt_HT_Sig->Draw("colz,0");
     665
     666  //plot for sPlot with M0 UNknown
     667  TCanvas *cesP = new TCanvas("cesP", "", 800,400);
     668  cesP->Divide(2);
     669  cesP->cd(1);
     670  h[6]->SetLineColor(2);
     671  h[6]->SetLineWidth(1);
     672  h[6]->DrawCopy("E1");
     673  h[6]->SetLineColor(1);
     674  h[6]->SetLineWidth(3);
     675  h[6]->Draw("hist,same");
     676  cesP->cd(2);
     677  twoDim_mt_HT_Sig_esP->Draw("colz,0"); 
     678}

  • SPLOT/NEWREALEASE/sPlots_NewRel_N0known.C

    r17 r19  
    8181  TFile *file[2];
    8282  TTree *t[2];
    83   TTree *tw[2];
    84  
    85   for (int i = 0; i < 2; i++){
    86     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     83 
     84  for (int i = 0; i < 2; i++){
     85    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    8786    file[i] = TFile::Open(filename);
    8887    t[i]    = (TTree*)file[i]->Get("T");
    89     tw[i]   = (TTree*)file[i]->Get("Tw");
    90     t[i]->AddFriend(tw[i]);
    9188  }
    9289
    … …  
    106103      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
    107104      double HT    = t[i]->GetLeaf("HT")->GetValue();
    108       double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    109       double y     = t[i]->GetLeaf("cosWqq")->GetValue();
     105      double x     = t[i]->GetLeaf("cos1")->GetValue();
     106      double y     = t[i]->GetLeaf("cos2")->GetValue();
    110107      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    111       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    112       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    113       double WeiTot = WeiEvt*Wei1fb;
    114 
    115       //cout << "WeiEvt = " << WeiEvt << endl;
    116       //cout << "Wei1fb = " << Wei1fb << endl;
     108      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     109
    117110      //cout << "WeiTot = " << WeiTot << endl;
    118111
    … …  
    218211  TFile *file[2];
    219212  TTree *t[2];
    220   TTree *tw[2];
    221  
    222   for (int i = 0; i < 2; i++){
    223     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     213 
     214  for (int i = 0; i < 2; i++){
     215    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    224216    file[i] = TFile::Open(filename);
    225217    t[i]    = (TTree*)file[i]->Get("T");
    226     tw[i]   = (TTree*)file[i]->Get("Tw");
    227     t[i]->AddFriend(tw[i]);
    228218  }
    229219 
    … …  
    235225      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
    236226      double HT     = t[i]->GetLeaf("HT")->GetValue();
    237       double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    238       double y      = t[i]->GetLeaf("cosWqq")->GetValue();
     227      double x      = t[i]->GetLeaf("cos1")->GetValue();
     228      double y      = t[i]->GetLeaf("cos2")->GetValue();
    239229      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    240       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    241       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    242       double WeiTot = WeiEvt*Wei1fb;
    243      
    244       //cout << "WeiEvt = " << WeiEvt << endl;
    245       //cout << "Wei1fb = " << Wei1fb << endl;
     230      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
     231     
    246232      //cout << "WeiTot = " << WeiTot << endl;
    247233     
    … …  
    318304  TFile *file[2];
    319305  TTree *t[2];
    320   TTree *tw[2];
    321  
    322   for (int i = 0; i < 2; i++){
    323     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
     306 
     307  for (int i = 0; i < 2; i++){
     308    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
    324309    file[i] = TFile::Open(filename);
    325310    t[i]    = (TTree*)file[i]->Get("T");
    326     tw[i]   = (TTree*)file[i]->Get("Tw");
    327     t[i]->AddFriend(tw[i]);
    328311  }
    329312 
    … …  
    340323      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
    341324      double HT     = t[i]->GetLeaf("HT")->GetValue();
    342       double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    343       double y      = t[i]->GetLeaf("cosWqq")->GetValue();
     325      double x      = t[i]->GetLeaf("cos1")->GetValue();
     326      double y      = t[i]->GetLeaf("cos2")->GetValue();
    344327      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    345       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    346       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    347       double WeiTot = WeiEvt*Wei1fb;
     328      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    348329     
    349330      double cos1 = fabs(x);
    … …  
    404385      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
    405386      double HT     = t[i]->GetLeaf("HT")->GetValue();
    406       double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
    407       double y      = t[i]->GetLeaf("cosWqq")->GetValue();
     387      double x      = t[i]->GetLeaf("cos1")->GetValue();
     388      double y      = t[i]->GetLeaf("cos2")->GetValue();
    408389      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
    409       double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    410       double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    411       double WeiTot = WeiEvt*Wei1fb;
     390      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
    412391
    413392      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    … …  
    416395      if (hf!=4 && hasBJ==0){
    417396       
    418         //cout << "WeiEvt " << WeiEvt << endl;
    419         //cout << "Wei1fb " << Wei1fb << endl;
    420      
    421397        double cos1 = fabs(x);
    422398        double cos2 = fabs(y);
    … …  
    461437          } else {
    462438            err[icc] = sWei*sWei/WeiTot;
    463             errN[icc] = sWeiN*sWeiN*WeiTot;
     439            errN[icc] = sWeiN*sWeiN/WeiTot;
    464440          }
    465441          icc++;
    … …  
    496472      double errV  = err[ic];
    497473      double errNV = errN[ic];
    498       if (mt>=vb1 && mt<vb2)
     474      if (mt>=vb1 && mt<vb2){
    499475        SumErrBin += errV;
    500476        SumErrNBin += errNV;
     477      }
    501478    }
    502479   
    … …  
    588565  h[9]->SetLineWidth(1);
    589566  h[9]->DrawCopy("E1");
    590   h[1]->SetLineColor(kGreen+3);
     567  h[1]->SetLineColor(kGreen+1);
    591568  h[1]->Draw("same");
    592569  h[9]->SetLineColor(1);
Note: See TracChangeset for help on using the changeset viewer.