Changeset 17 in huonglan


Ignore:
Timestamp:
Aug 29, 2011, 12:58:02 AM (13 years ago)
Author:
huonglan
Message:

upgrade sPlots

Location:
SPLOT/NEWREALEASE
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • SPLOT/NEWREALEASE/OptimizeNastyCut_2cosDelMinLep_NewRel.C

    r14 r17  
    1010
    1111  for (int i = 0; i < 2; i++){
    12     filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3Opt.root";
     12    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
    1313    file[i] = TFile::Open(filename);
    1414    t[i]    = (TTree*)file[i]->Get("T");

  • SPLOT/NEWREALEASE/sPlots_Data.C

    r14 r17  
    1818double Cut4 = 250;
    1919double Cut5 = 620;
    20 double sigmaMin = 0.185191;
     20double sigmaMin;
    2121double iTot[2];
    2222double iCatS[6], iCatB[6];
    … …  
    146146          else
    147147            iCatS[it] += WeiTot;
    148         }
    149       }
     148        }//Nasty cut
     149      }//hf and hadBJ cut
    150150    }//Read each file
    151151  }//End of reading 2 files
    … …  
    154154  cout << "iTot1 = " << iTot[1] << endl;
    155155
    156   //N1 = iTot[1];
    157   //N0 = iTot[0];
    158 
     156  N1 = iTot[1];
     157  N0 = iTot[0];
     158
     159  double sigmaInv = 0;
    159160  double sumeffS = 0, sumeffB = 0;
    160161  for (int i = 0; i < 6; i++){
    161    
    162162    F1[i] = iCatS[i]/iTot[1];
    163163    F0[i] = iCatB[i]/iTot[0];
    … …  
    166166    cout << "F1[" << i << "] = " << F1[i] << endl;
    167167    cout << "F0[" << i << "] = " << F0[i] << endl;
    168   }
     168
     169    double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
     170    if (F1[i]==0 && F0[i]==0)  term = 0;
     171    sigmaInv += term;
     172  }
     173
     174  double sigma2 = 1/N1/N1/sigmaInv;
     175  sigmaMin = sqrt(sigma2);
    169176 
    170177  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
     178  cout << "sigmaMin = " << sigmaMin << endl;
    171179}
    172180

  • SPLOT/NEWREALEASE/sPlots_Data_saved.C

    r14 r17  
    1818double Cut4 = 250;
    1919double Cut5 = 620;
    20 double sigmaMin = 0.171321;
     20double sigmaMin = 0.185191;
    2121double iTot[2];
    2222double iCatS[6], iCatB[6];
    … …  
    2626double N1 = 0., N0 = 0., sig2 = 0.;
    2727
     28//Define class Event
    2829class EventC {
    2930public:
    … …  
    5556};
    5657
     58
    5759vector<double> InvMatrix(double a, double b, double c, double d){
    5860  double detM = a*d - b*c;
    … …  
    107109      double y     = t[i]->GetLeaf("cosWqq")->GetValue();
    108110      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
    109       //double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
    110111      double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
    111112      double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    … …  
    116117      //cout << "WeiTot = " << WeiTot << endl;
    117118
    118       double cos1 = fabs(x);
    119       double cos2 = fabs(y);
    120      
    121       if (mt > Cut4 && HT > Cut5) {
    122         iTot[i] += WeiTot;
    123 
    124         if (cos1<Cut1 && cos2<Cut2){
    125           if ((dMLep<Cut3)){
    126             if (i==0) iCatB[0] += WeiTot;
    127             if (i==1) iCatS[0] += WeiTot;
     119      int hasBJ = t[i]->GetLeaf("hasBadJet")->GetValue();
     120      double hf = t[i]->GetLeaf("hfor")->GetValue();
     121
     122      if (hf!=4 && hasBJ==0){
     123       
     124        double cos1 = fabs(x);
     125        double cos2 = fabs(y);
     126     
     127        if (mt > Cut4 && HT > Cut5) {
     128          iTot[i] += WeiTot;
     129          int it = -1;
     130          if (cos1<Cut1 && cos2<Cut2){
     131            if (dMLep<Cut3) it = 0;
     132            if (dMLep>Cut3) it = 1;
    128133          }
    129           if ((dMLep>Cut3)){
    130             if (i==0) iCatB[1] += WeiTot;
    131             if (i==1) iCatS[1] += WeiTot;
     134          else if (cos1>Cut1 && cos2>Cut2){
     135            if (dMLep<Cut3) it = 2;
     136            if (dMLep>Cut3) it = 3;
    132137          }
    133         }
    134        
    135         if (cos1>Cut1 && cos2>Cut2){
    136           if ((dMLep<Cut3)){
    137             if (i==0) iCatB[2] += WeiTot;
    138             if (i==1) iCatS[2] += WeiTot;
     138          else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
     139            if (dMLep<Cut3) it = 4;
     140            if (dMLep>Cut3) it = 5;
     141          } else {
     142            cout << "Unknown category !!" << endl;
    139143          }
    140           if ((dMLep>Cut3)){
    141             if (i==0) iCatB[3] += WeiTot;
    142             if (i==1) iCatS[3] += WeiTot;
    143           }
    144         }
    145        
    146         if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
    147           if ((dMLep<Cut3)){
    148             if (i==0) iCatB[4] += WeiTot;
    149             if (i==1) iCatS[4] += WeiTot;
    150           }
    151           if ((dMLep>Cut3)){
    152             if (i==0) iCatB[5] += WeiTot;
    153             if (i==1) iCatS[5] += WeiTot;
    154           }
     144          if (i == 0)
     145            iCatB[it] += WeiTot;
     146          else
     147            iCatS[it] += WeiTot;
    155148        }
    156149      }
    … …  
    161154  cout << "iTot1 = " << iTot[1] << endl;
    162155
     156  //N1 = iTot[1];
     157  //N0 = iTot[0];
     158
    163159  double sumeffS = 0, sumeffB = 0;
    164160  for (int i = 0; i < 6; i++){
    165 
     161   
    166162    F1[i] = iCatS[i]/iTot[1];
    167163    F0[i] = iCatB[i]/iTot[0];
    … …  
    179175  double f0e = 0;
    180176  double f1e = 0;
    181 
     177 
    182178  //cout << "cos1 = " << cos1 << " cos2 = " << cos2 << " dMLep = " << dMLep << endl;
    183179 
     180  int it = -1;
    184181  if (cos1<Cut1 && cos2<Cut2){
    185     if ((dMLep<Cut3)){
    186       f0e = F0[0];
    187       f1e = F1[0];
    188     }
    189     if ((dMLep>Cut3)){
    190       f0e = F0[1];
    191       f1e = F1[1];
    192     }
    193   }
    194  
    195   if (cos1>Cut1 && cos2>Cut2){
    196     if ((dMLep<Cut3)){
    197       f0e = F0[2];
    198       f1e = F1[2];
    199     }
    200     if ((dMLep>Cut3)){
    201       f0e = F0[3];
    202       f1e = F1[3];
    203     }
    204   }
    205  
    206   if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
    207     if ((dMLep<Cut3)){
    208       f0e = F0[4];
    209       f1e = F1[4];
    210     }
    211     if ((dMLep>Cut3)){
    212       f0e = F0[5];
    213       f1e = F1[5];
    214     }
    215   }
    216 
     182    if (dMLep<Cut3) it = 0;
     183    if (dMLep>Cut3) it = 1;
     184  }
     185  else if (cos1>Cut1 && cos2>Cut2){
     186    if (dMLep<Cut3) it = 2;
     187    if (dMLep>Cut3) it = 3;
     188  }
     189  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
     190    if (dMLep<Cut3) it = 4;
     191    if (dMLep>Cut3) it = 5;
     192  }
     193 
     194  f0e = F0[it];
     195  f1e = F1[it];
     196 
     197  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
     198 
    217199  vector<double> vfe;
    218200  vfe.push_back(f0e);
    … …  
    221203}
    222204
     205int FillHisto(double cos1,double cos2,double dMLep){
     206  int it = -1;
     207  if (cos1<Cut1 && cos2<Cut2){
     208    if (dMLep<Cut3) it = 0;
     209    if (dMLep>Cut3) it = 1;
     210  }
     211  else if (cos1>Cut1 && cos2>Cut2){
     212    if (dMLep<Cut3) it = 2;
     213    if (dMLep>Cut3) it = 3;
     214  }
     215  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
     216    if (dMLep<Cut3) it = 4;
     217    if (dMLep>Cut3) it = 5;
     218  }
     219}
     220
    223221void MakeEventListf0f1(){
    224   //######## COMPUTE the sWeights
    225 
    226   TFile *fD = TFile::Open("resData.root");
     222  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
    227223  TTree *tD = (TTree*)fD->Get("T");
    228224  int nD = tD->GetEntries();
    229225  for (int iev = 0; iev < nD; iev++) {
    230226    tD->GetEntry(iev);
    231 
    232     double mt     = tD->GetLeaf("mtophad")->GetValue();
    233     double HT     = tD->GetLeaf("HT")->GetValue();
    234     double x      = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
    235     double y      = tD->GetLeaf("cosWqq")->GetValue();
    236     double dMLep  = tD->GetLeaf("MinLep1")->GetValue();
    237 
    238     //cout << "WeiEvt = " << WeiEvt << endl;
    239     //cout << "Wei1fb = " << Wei1fb << endl;
    240     //cout << "WeiTot = " << WeiTot << endl;
    241    
     227   
     228    double mt    = tD->GetLeaf("mtophad")->GetValue();
     229    double HT    = tD->GetLeaf("HT")->GetValue();
     230    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     231    double y     = tD->GetLeaf("cosWqq")->GetValue();
     232    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     233
     234    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
     235    int RunN  = tD->GetLeaf("RunNumber_ana")->GetValue();
     236   
     237    if (RunN >= 180614 && hasBJ==1) continue;
     238
    242239    double cos1 = fabs(x);
    243240    double cos2 = fabs(y);
    244      
    245     if (mt > Cut4 && HT > Cut5) {
     241   
     242    //sPlots total
     243    if (mt > Cut4 && HT > Cut5){
    246244      vector<double> fsfb =  vFe(cos1,cos2,dMLep);
    247245      double f0 = fsfb.at(0);
    … …  
    250248      EventC eventC(1,f1,f0);
    251249      m_EventList.push_back(eventC);
    252     }
     250    } 
    253251  }
    254252}
    … …  
    259257  MakeEventListf0f1();
    260258
    261   double Nexpt0 = iTot[0]/3;
    262   double Nexpt1 = iTot[1]/3;
    263 
    264259  int ifl = 0;
    265260  double step = 0.01;
    266   double fitPar0 = Nexpt0;
    267   double fitPar1 = Nexpt1;
    268   double nmax0   = fitPar0 + 4*sqrt(Nexpt0);
    269   double nmin0   = fitPar0 - 4*sqrt(Nexpt0);
    270   double nmax1   = fitPar1 + 4*sqrt(Nexpt1);
    271   double nmin1   = - 4*sqrt(Nexpt1);
     261  double fitPar0 = iTot[0];
     262  double fitPar1 = iTot[1];
     263  double nmax0   = fitPar0 + 4*sqrt(iTot[0]);
     264  double nmin0   = fitPar0 - 4*sqrt(iTot[0]);
     265  double nmax1   = fitPar1 + 4*sqrt(iTot[1]);
     266  double nmin1   = - 4*sqrt(iTot[1]);
    272267  double n0best, err0best;
    273268  double n1best, err1best;
    … …  
    284279  myMinuit->GetParameter(0,n0best,err0best);
    285280  myMinuit->GetParameter(1,n1best,err1best);
    286   std::cout << "Expected N0 = " << Nexpt0 << " fitted = " << n0best << " +- " << err0best << std::endl;
    287   std::cout << "Expected N1 = " << Nexpt1 << " fitted = " << n1best << " +- " << err1best << std::endl;
     281  std::cout << "Expected N0 = " << iTot[0] << " fitted = " << n0best << " +- " << err0best << std::endl;
     282  std::cout << "Expected N1 = " << iTot[1] << " fitted = " << n1best << " +- " << err1best << std::endl;
    288283
    289284  N0 = n0best;
    290285  N1 = n1best;
    291286
     287  cout << "N0 = " << N0 << " N1 = " << N1 << endl;
    292288  sig2 = sigmaMin*sigmaMin*N1*N1;
    293289}
    … …  
    309305  //Get right N0 and N1
    310306  Fit();
    311 
    312   cout << "After fit : " << endl;
    313   cout << "N0 = " << N0 << " N1 = " << N1 << endl;
    314307 
    315308  double VI00 = 0, VI01 = 0, VI10 = 0, VI11 = 0;
    316   //start reading 2 files
    317   int icount = 0;
    318 
    319   //Try standard sPlots
    320   TFile *fD = TFile::Open("resData.root");
     309
     310  int Nsel = 0;
     311  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
    321312  TTree *tD = (TTree*)fD->Get("T");
    322313  int nD = tD->GetEntries();
    323314  for (int iev = 0; iev < nD; iev++) {
    324315    tD->GetEntry(iev);
    325      
    326     double mt     = tD->GetLeaf("mtophad")->GetValue();
    327     double HT     = tD->GetLeaf("HT")->GetValue();
    328     double x      = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
    329     double y      = tD->GetLeaf("cosWqq")->GetValue();
    330     double dMLep  = tD->GetLeaf("MinLep1")->GetValue();
     316   
     317    double mt    = tD->GetLeaf("mtophad")->GetValue();
     318    double HT    = tD->GetLeaf("HT")->GetValue();
     319    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     320    double y     = tD->GetLeaf("cosWqq")->GetValue();
     321    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     322     
     323    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
     324    int RunN  = tD->GetLeaf("RunNumber_ana")->GetValue();
     325   
     326    if (RunN >= 180614 && hasBJ==1) continue;
     327
    331328    double cos1 = fabs(x);
    332329    double cos2 = fabs(y);
    333      
     330   
     331    //sPlots total
    334332    if (mt > Cut4 && HT > Cut5){
    335       icount++;
     333      Nsel++;
    336334      //Compute the matrix
    337335      vector<double> veff = vFe(cos1,cos2,dMLep);
    338336      double F0e = veff.at(0);
    339337      double F1e = veff.at(1);
     338     
    340339      double vi00;
    341340      double vi01;
    … …  
    347346      vi10 = F1e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
    348347      vi11 = F1e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
    349 
     348     
    350349      VI00 += vi00;
    351350      VI01 += vi01;
    352351      VI10 += vi10;
    353352      VI11 += vi11;
    354     }
    355   }
    356 
     353    }//Read each file
     354  }//End of reading 2 files
     355
     356  cout << "Nsel = " << Nsel << endl;
     357  cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
     358
     359  cout << "VI00 = " << VI00 << endl;
     360  cout << "VI01 = " << VI01 << endl;
     361  cout << "VI10 = " << VI10 << endl;
     362  cout << "VI11 = " << VI11 << endl;
    357363  vector<double> Vinv = InvMatrix(VI00,VI01,VI10,VI11);
    358364  double V00 = Vinv.at(0);
    … …  
    360366  double V10 = Vinv.at(2);
    361367  double V11 = Vinv.at(3);
     368  cout << "V00 = " << V00 << endl;
     369  cout << "V01 = " << V01 << endl;
     370  cout << "V10 = " << V10 << endl;
     371  cout << "V11 = " << V11 << endl;
    362372  double sigmaN0 = sqrt(V00);
    363373  double sigmaN1 = sqrt(V11);
    364374  double rho     = V10/sqrt(V00*V11);
    365 
    366 //   cout << "VI00 = " << VI00 << endl;
    367 //   cout << "VI01 = " << VI01 << endl;
    368 //   cout << "VI10 = " << VI10 << endl;
    369 //   cout << "VI11 = " << VI11 << endl;
    370 //   cout << "V00 = " << V00 << endl;
    371 //   cout << "V01 = " << V01 << endl;
    372 //   cout << "V10 = " << V10 << endl;
    373 //   cout << "V11 = " << V11 << endl;
    374 //   cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl;
    375 
    376   TH1F *hSigN  = new TH1F("hSigN",     "", 25, 0, 1500);
    377   TH1F *hBkgN  = new TH1F("hBkgN",     "", 25, 0, 1500);
    378   TH2F *h2DSig = new TH2F("h2DSig",    "", 50, 0, 1500, 50, 0, 1500);
    379   TH2F *h2DBkg = new TH2F("h2DBkg",    "", 50, 0, 1500, 50, 0, 1500);
    380 
    381   cout << "icount = " << icount << endl;
     375  cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl;
     376
     377  TH1F *hMTSig   = new TH1F("hMTSig",     "", 25, 0, 1500);
     378  TH1F *hMTBkg   = new TH1F("hMTBkg",     "", 25, 0, 1500);
     379  TH1F *hHTSig   = new TH1F("hHTSig",     "", 25, 0, 2000);
     380  TH1F *hHTBkg   = new TH1F("hHTBkg",     "", 25, 0, 2000);
     381  TH2F *hHTMTSig = new TH2F("hHTMTSig",   "", 25, 0, 2000, 25, 0, 1500);
     382  TH2F *hHTMTBkg = new TH2F("hHTMTBkg",   "", 25, 0, 2000, 25, 0, 1500);
     383
     384  TH1F *hmt[6];
     385  TH1F *hht[6];
     386  for (int i = 0; i < 6; i++){
     387    TString hname = "hmt";
     388    hname += "Cat";
     389    hname += (i+1);
     390    hmt[i] = new TH1F(hname, "", 25, 0, 1500);
     391
     392    hname = "hht";
     393    hname += "Cat";
     394    hname += (i+1);
     395    hht[i] = new TH1F(hname, "", 25, 0, 1500);
     396  }
    382397
    383398  double SumsP1 = 0.;
    384399  double SumsP0 = 0.;
    385   double sPN[28481];
    386   double mtN[28481];
     400  double *sPN = new double[Nsel];
     401  double *htN = new double[Nsel];
     402  double *mtN = new double[Nsel];
     403 
    387404  int icc = 0;
     405  int icountD = 0;
     406  int iCat[6];
     407  for (int i = 0; i < 6; i++) { iCat[i] = 0;}
    388408  //Try standard sPlots
    389   TFile *fD = TFile::Open("resData.root");
     409  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
    390410  TTree *tD = (TTree*)fD->Get("T");
    391411  int nD = tD->GetEntries();
    392412  for (int iev = 0; iev < nD; iev++) {
    393413    tD->GetEntry(iev);
    394      
    395     double mt     = tD->GetLeaf("mtophad")->GetValue();
    396     double HT     = tD->GetLeaf("HT")->GetValue();
    397     double x      = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
    398     double y      = tD->GetLeaf("cosWqq")->GetValue();
    399     double dMLep  = tD->GetLeaf("MinLep1")->GetValue();
    400      
     414   
     415    double mt    = tD->GetLeaf("mtophad")->GetValue();
     416    double HT    = tD->GetLeaf("HT")->GetValue();
     417    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
     418    double y     = tD->GetLeaf("cosWqq")->GetValue();
     419    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
     420     
     421    int hasBJ = tD->GetLeaf("hasBadJet")->GetValue();
     422    int RunN  = tD->GetLeaf("RunNumber_ana")->GetValue();
     423   
     424    if (RunN >= 180614 && hasBJ==1) continue;
     425
    401426    double cos1 = fabs(x);
    402427    double cos2 = fabs(y);
    403 
     428   
    404429    //sPlots total
    405430    if (mt > Cut4 && HT > Cut5){
     431      icountD++;
    406432      vector<double> veff = vFe(cos1,cos2,dMLep);
    407433      double F0e = veff.at(0);
    408434      double F1e = veff.at(1);
    409 
     435     
    410436      double sP1 = (V11*F1e + V10*F0e)/(N1*F1e + N0*F0e);
    411437      double sP0 = (V00*F0e + V10*F1e)/(N1*F1e + N0*F0e);
    412 
    413       hSigN->Fill(mt,sP1);
    414       hBkgN->Fill(mt,sP0);
     438     
     439      hMTSig->Fill(mt,sP1);
     440      hMTBkg->Fill(mt,sP0);
     441      hHTSig->Fill(HT,sP1);
     442      hHTBkg->Fill(HT,sP0);
     443      hHTMTSig->Fill(HT,mt,sP1);
     444      hHTMTBkg->Fill(HT,mt,sP0);
    415445      sPN[icc] = sP1*sP1;
    416446      mtN[icc] = mt;
     447      htN[icc] = HT;
     448      icc++;
     449     
     450      int it = FillHisto(cos1,cos2,dMLep);
     451      hmt[it]->Fill(mt,sP1);
     452      hht[it]->Fill(HT,sP1);
     453      iCat[it]++;
    417454
    418455      SumsP1 += sP1;
    419456      SumsP0 += sP0;
    420     }
    421   }
     457    }//Read each file
     458  }//End of reading 2 files
     459 
     460  cout << "icountD = " << icountD << endl;
     461  cout << "SumsP0 = " << SumsP0 << " SumsP1 = " << SumsP1 << endl;
    422462
    423463  cout << endl << "*******************" << endl;
    424464  cout << "STANDARD SPLOTS" << endl;
    425   cout << "hSigN integral  " << hSigN->Integral(0,1500)  << endl;
    426   cout << "hBkgN integral  " << hBkgN->Integral(0,1500)  << endl;
    427 
    428   /*
    429   for (int i = 0; i < 28481; i++){
    430     if (sPN[i]==0) continue;
    431     cout << "mtN[" << i << "] = " << mtN[i] << endl;
    432     cout << "sPN[" << i << "] = " << sPN[i] << endl;
    433   }
    434   */
    435 
    436   for (int i =1; i < 26; i++){
     465  cout << "hMTSig integral  " << hMTSig->Integral(0,1500)  << endl;
     466  cout << "hMTBkg integral  " << hMTBkg->Integral(0,1500)  << endl;
     467
     468//   for (int i = 0; i < Nsel; i++){
     469//     cout << "mt[" << i << "] = " << mtN[i] << " ";
     470//     cout << "sP[" << i << "] = " << sPN[i] << endl;
     471//   }
     472 
     473  for (int i = 1; i < 26; i++){
    437474    double vb1 = (i-1)*60;
    438475    double vb2 = i*60;
    439476    double ErrBin = 0.;
    440     for (int j = 0; j < 28481; j++){
     477    for (int j = 0; j < Nsel; j++){
    441478      double mt = mtN[j];
    442479      double sp = sPN[j];
    … …  
    445482      }
    446483    }
    447     hSigN->SetBinError(i,sqrt(ErrBin));
    448   }
    449 
    450   TCanvas *c = new TCanvas("c", "", 900, 400);
    451   c->Divide(2);
    452   c->cd(1);
    453   hSigN->Draw("E1");
    454   //hSigN->Draw("hist,same");
    455   c->cd(2);
    456   hBkgN->Draw();
    457 
    458   //TCanvas *c2 = new TCanvas("c2", "", 500, 400);
    459   //h2DSig->Draw("colz");
    460 
    461   //TCanvas *c3 = new TCanvas("c3", "", 500, 400);
    462   //h2DBkg->Draw("colz");
    463 }
     484    hMTSig->SetBinError(i,sqrt(ErrBin));
     485  }
     486
     487  for (int i = 1; i < 26; i++){
     488    double vb1 = (i-1)*80;
     489    double vb2 = i*80;
     490    double ErrBin = 0.;
     491    for (int j = 0; j < Nsel; j++){
     492      double ht = htN[j];
     493      double sp = sPN[j];
     494      if (ht >= vb1 && ht < vb2){
     495        ErrBin += sp;
     496      }
     497    }
     498    hHTSig->SetBinError(i,sqrt(ErrBin));
     499  }
     500
     501  TCanvas *c1 = new TCanvas("c1", "", 900, 400);
     502  c1->Divide(2);
     503  c1->cd(1);
     504  hMTSig->SetLineColor(2);
     505  hMTSig->DrawCopy("E1");
     506  hMTSig->SetLineColor(1);
     507  hMTSig->Draw("hist,same");
     508  c1->cd(2);
     509  hMTBkg->Draw();
     510
     511  TCanvas *c2 = new TCanvas("c2", "", 900, 400);
     512  c2->Divide(2);
     513  c2->cd(1);
     514  hHTSig->SetLineColor(2);
     515  hHTSig->DrawCopy("E1");
     516  hHTSig->SetLineColor(1);
     517  hHTSig->Draw("hist,same");
     518  c2->cd(2);
     519  hHTBkg->Draw();
     520
     521  TCanvas *c3 = new TCanvas("c3", "", 900, 400);
     522  c3->Divide(2);
     523  c3->cd(1);
     524  hHTMTSig->Draw("colz");
     525  c3->cd(2);
     526  hHTMTBkg->Draw("colz");
     527
     528  TCanvas *c4 = new TCanvas("c4", "", 1000, 600);
     529  c4->Divide(6,2);
     530  for (int i = 1; i < 7; i++){
     531    c4->cd(i);
     532    hmt[i-1]->Draw();
     533    c4->cd(i+6);
     534    hht[i-1]->Draw();
     535   
     536    cout << "iCat[" << i << "] = " << iCat[i-1] << endl;
     537    cout << "integral of " << hmt[i-1]->GetName() << " is " << hmt[i-1]->Integral(0,100) << endl;
     538    cout << "integral of " << hht[i-1]->GetName() << " is " << hht[i-1]->Integral(0,100) << endl;
     539  }
     540}

  • SPLOT/NEWREALEASE/sPlots_NewRel.C

    r14 r17  
    1818double Cut4 = 250;
    1919double Cut5 = 620;
    20 double sigmaMin = 0.163173;
     20double sigmaMin;
    2121double iTot[2];
    2222double iCatS[6], iCatB[6];
    … …  
    157157  N0 = iTot[0];
    158158
     159  double sigmaInv = 0;
    159160  double sumeffS = 0, sumeffB = 0;
    160161  for (int i = 0; i < 6; i++){
    161    
    162162    F1[i] = iCatS[i]/iTot[1];
    163163    F0[i] = iCatB[i]/iTot[0];
    … …  
    166166    cout << "F1[" << i << "] = " << F1[i] << endl;
    167167    cout << "F0[" << i << "] = " << F0[i] << endl;
    168   }
     168
     169    double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
     170    if (F1[i]==0 && F0[i]==0)  term = 0;
     171    sigmaInv += term;
     172  }
     173
     174  double sigma2 = 1/N1/N1/sigmaInv;
     175  sigmaMin = sqrt(sigma2);
    169176 
    170177  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
     178  cout << "sigmaMin = " << sigmaMin << endl;
    171179}
    172180
    … …  
    400408  cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl;
    401409
    402   TH1F *hSigN  = new TH1F("hSigN",     "", 25, 0, 1500);
    403   TH1F *hBkgN  = new TH1F("hBkgN",     "", 25, 0, 1500);
    404   TH2F *h2DSig = new TH2F("h2DSig",    "", 50, 0, 1500, 50, 0, 1500);
    405   TH2F *h2DBkg = new TH2F("h2DBkg",    "", 50, 0, 1500, 50, 0, 1500);
     410  TH1F *hSigN  = new TH1F("hSigN",     "", 35, 250, 600);
     411  TH1F *hBkgN  = new TH1F("hBkgN",     "", 35, 250, 600);
     412  TH2F *h2DSig = new TH2F("h2DSig",    "", 35, 250, 600, 50, 500, 1500);
     413  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 35, 250, 600, 50, 500, 1500);
    406414
    407415  cout << "icount = " << icount << endl;
    … …  
    446454          hSigN->Fill(mt,sP1*WeiTot);
    447455          hBkgN->Fill(mt,sP0*WeiTot);
    448           sPN[icc] = sP1*sP1*WeiTot;
     456          h2DSig->Fill(mt,HT,sP1*WeiTot);
     457          h2DBkg->Fill(mt,HT,sP0*WeiTot);
    449458          mtN[icc] = mt;
     459          if (WeiTot==0) sPN[icc] = 0;
     460          else sPN[icc] = sP1*sP1*WeiTot;
     461          icc++;
    450462         
    451463          SumsP1 += sP1*WeiTot;
    … …  
    463475  cout << "hBkgN integral  " << hBkgN->Integral(0,1500)  << endl;
    464476 
    465   for (int i =1; i < 26; i++){
    466     double vb1 = (i-1)*60;
    467     double vb2 = i*60;
     477  for (int i = 1; i < 36; i++){
     478    double vb1 = 250 + (i-1)*16;
     479    double vb2 = 250 + i*16;
    468480    double ErrBin = 0.;
    469     for (int j = 0; j < 23054; j++){
     481    for (int j = 0; j < icount; j++){
    470482      double mt = mtN[j];
    471483      double sp = sPN[j];
    … …  
    479491  TCanvas *c = new TCanvas("c", "", 900, 400);
    480492  c->Divide(2);
    481   c->cd(1);
    482   hSigN->Draw("E1");
     493  c->cd(1);
     494  hSigN->SetLineColor(2);
     495  hSigN->SetLineWidth(1);
     496  //hSigN->DrawCopy("E1");
     497  hSigN->SetLineColor(1);
     498  hSigN->SetLineWidth(2.5);
     499  hSigN->DrawCopy("hist");
    483500  c->cd(2);
     501  hBkgN->SetLineColor(4);
     502  hBkgN->SetLineWidth(2.5);
    484503  hBkgN->Draw();
     504
     505  TCanvas *c2D = new TCanvas("c2D", "", 900, 400);
     506  c2D->Divide(2);
     507  c2D->cd(1);
     508  h2DSig->Draw("colz");
     509  c2D->cd(2);
     510  h2DBkg->Draw("colz");
    485511}

  • SPLOT/NEWREALEASE/sPlots_NewRel_N0known.C

    r14 r17  
    1818double Cut4 = 250;
    1919double Cut5 = 620;
    20 double sigmaMin = 0.163173;
     20double sigmaMin;
    2121double iTot[2];
    2222double iCatS[6], iCatB[6];
    … …  
    157157  N0 = iTot[0];
    158158
     159  double sigmaInv = 0;
    159160  double sumeffS = 0, sumeffB = 0;
    160161  for (int i = 0; i < 6; i++){
    161    
    162162    F1[i] = iCatS[i]/iTot[1];
    163163    F0[i] = iCatB[i]/iTot[0];
    … …  
    166166    cout << "F1[" << i << "] = " << F1[i] << endl;
    167167    cout << "F0[" << i << "] = " << F0[i] << endl;
    168   }
     168
     169    double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
     170    if (F1[i]==0 && F0[i]==0)  term = 0;
     171    sigmaInv += term;
     172  }
     173
     174  double sigma2 = 1/N1/N1/sigmaInv;
     175  sigmaMin = sqrt(sigma2);
    169176 
    170177  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
     178  cout << "sigmaMin = " << sigmaMin << endl;
    171179}
    172180
    … …  
    359367  cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
    360368 
    361   TH2F *twoDim_mt_HT     = new TH2F("twoDim_mt_HT",     "", 150, 0, 1500, 100, 0, 2000);
    362   TH2F *twoDim_mt_HT_Bkg = new TH2F("twoDim_mt_HT_Bkg", "", 150, 0, 1500, 100, 0, 2000);
    363   TH2F *twoDim_mt_HT_Sig = new TH2F("twoDim_mt_HT_Sig", "", 150, 0, 1500, 100, 0, 2000);
    364   TH2F *twoDim_mt_HT_Sig_esP = new TH2F("twoDim_mt_HT_Sig_esP", "", 150, 0, 1500, 100, 0, 2000);
    365  
    366   TH1F *hSig        = new TH1F("hSig",    "", 150, 0, 1500);
    367   TH1F *hBkg        = new TH1F("hBkg",    "", 150, 0, 1500);
     369  TH2F *twoDim_mt_HT_Sig     = new TH2F("twoDim_mt_HT_Sig",     "", 75, 250, 1500, 50, 500, 2000);
     370  TH2F *twoDim_mt_HT_Bkg     = new TH2F("twoDim_mt_HT_Bkg",     "", 75, 250, 1500, 50, 500, 2000);
     371  TH2F *twoDim_mt_HT_Sig_esP = new TH2F("twoDim_mt_HT_Sig_esP", "", 75, 250, 1500, 50, 500, 2000);
     372  TH2F *twoDim_mt_HT_Bkg_esP = new TH2F("twoDim_mt_HT_Bkg_esP", "", 75, 250, 1500, 50, 500, 2000);
     373 
    368374  TH1F *h[11];
    369375  h[0]  = new TH1F("hAll1ifb",       "",  150,  0,  1500);
    … …  
    404410      double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
    405411      double WeiTot = WeiEvt*Wei1fb;
    406      
     412
    407413      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
    408414      double hf = t[i]->GetLeaf("hfor")->GetValue();
    … …  
    410416      if (hf!=4 && hasBJ==0){
    411417       
     418        //cout << "WeiEvt " << WeiEvt << endl;
     419        //cout << "Wei1fb " << Wei1fb << endl;
     420     
    412421        double cos1 = fabs(x);
    413422        double cos2 = fabs(y);
    … …  
    419428          double F1e = veff.at(1);
    420429         
    421           double sWei = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
     430          double sWei  = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
    422431          double sWeiN = sWei*(k+1) - k*WeiTot;
    423432         
    424           twoDim_mt_HT->Fill(mt,HT,sWei);       
    425           if (i==0)
     433          if (i==0){
    426434            twoDim_mt_HT_Bkg->Fill(mt,HT,sWei);
    427          
    428           twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);     
     435            twoDim_mt_HT_Bkg_esP->Fill(mt,HT,sWeiN);   
     436          }
     437          if (i==1){
     438            twoDim_mt_HT_Sig->Fill(mt,HT,sWei);
     439            twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);   
     440          }
    429441         
    430442          //1D plot
    … …  
    437449            h[8]->Fill(mt,sWei);
    438450          }
     451         
    439452          //sPlots
     453          //M0 distribution known
    440454          h[3]->Fill(mt,sWei);
     455          //M0 distribution unknown
    441456          h[7]->Fill(mt,sWeiN);
    442457          mtV[icc] = mt;
    443           err[icc] = sWei*sWei/WeiTot;
    444           errN[icc] = sWeiN*sWeiN*WeiTot;
     458          if (WeiTot==0){
     459            err[icc] = 0;
     460            errN[icc] = 0;
     461          } else {
     462            err[icc] = sWei*sWei/WeiTot;
     463            errN[icc] = sWeiN*sWeiN*WeiTot;
     464          }
    445465          icc++;
    446466          //if (i==0)
    447467          //esPlotBkg->Fill(mt,sWei - sWeiN);
    448           //if (i==1)
    449           //esPlotSig->Fill(mt,sWei);
     468          if (i==1)
     469            esPlotSig->Fill(mt,sWei);
    450470        }
    451471      }
    … …  
    459479  cout << "h7 integral is " << h[7]->Integral(0,1500) << endl;
    460480
    461   for (int i = 1; i < 151; i++){
    462     double WSig = 0., WBkg = 0.;
    463     for (int j = 1; j < 101; j++){
    464       double vtot = twoDim_mt_HT->GetBinContent(i,j);
    465       double vBkg = twoDim_mt_HT_Bkg->GetBinContent(i,j);
    466      
    467       double vSig = vtot - vBkg;
    468       twoDim_mt_HT_Sig->SetBinContent(i,j,vSig);
    469       WSig += vSig;
    470       WBkg += vBkg;
    471     }
    472     hSig->SetBinContent(i,WSig);
    473     hBkg->SetBinContent(i,WBkg);
    474   }
    475  
    476481  for (int i = 1; i < 151; i++){
    477482    double sP = h[3]->GetBinContent(i);
    … …  
    487492    double SumErrBin = 0;
    488493    double SumErrNBin = 0;
    489     for (int ic = 0; ic < 10616; ic++){
     494    for (int ic = 0; ic < icount; ic++){
    490495      double mt  = mtV[ic];
    491496      double errV  = err[ic];
    … …  
    496501    }
    497502   
     503    //cout << "SumErrBin " << SumErrBin << endl;
     504
    498505    h[4]->SetBinContent(i,W);
    499506    h[9]->SetBinContent(i,sW);
    … …  
    526533  double x2 = 600 - 0.01;
    527534  for (int i = 0; i < 11; i++) {
     535    //h[i]->Rebin(2);
    528536    h[i]->SetLineWidth(3);
    529537    TAxis *ax = h[i]->GetXaxis();
    … …  
    536544  cout << "h2 integral is " << h[2]->Integral(0,1500) << endl;
    537545  cout << "h7 integral is " << h[7]->Integral(0,1500) << endl;
    538 
    539   for (int i = 1; i < 151; i++){
    540     double WSig = 0., WBkg = 0.;
    541     for (int j = 1; j < 101; j++){
    542       double vtot = twoDim_mt_HT->GetBinContent(i,j);
    543       double vBkg = twoDim_mt_HT_Bkg->GetBinContent(i,j);
    544 
    545       double vSig = vtot - vBkg;
    546       twoDim_mt_HT_Sig->SetBinContent(i,j,vSig);
    547       WSig += vSig;
    548       WBkg += vBkg;
    549     }
    550     hSig->SetBinContent(i,WSig);
    551     hBkg->SetBinContent(i,WBkg);
    552   }
    553 
    554   TCanvas *cS = new TCanvas("cS", "", 800, 600);
     546  cout << "h[9] integral is " << h[9]->Integral(0,1000) << endl;
     547
     548  /*  TCanvas *cS = new TCanvas("cS", "", 800, 600);
    555549  cS->Divide(3);
    556550  cS->cd(1);
    … …  
    562556  cS->cd(3);
    563557  h[9]->SetLineColor(2);
     558  h[9]->Draw("E1");
    564559  h[1]->SetLineColor(kGreen+3);
    565   h[1]->Draw();
     560  h[1]->Draw("same");
     561  h[9]->SetLineWidth(1);
    566562  h[9]->DrawCopy("hist,same");
    567 
    568   //TCanvas *csPlot = new TCanvas("csPlot", "", 600,500);
    569   //csPlot->cd();
    570   //h[9]->SetLineWidth(1);
    571   //h[9]->DrawCopy("E1");
    572   //h[9]->SetLineWidth(3);
    573   //h[9]->Draw("hist,same");
    574   //h[1]->Draw("same");
    575563
    576564  TCanvas *cN = new TCanvas("cN", "", 1000, 400);
    … …  
    586574  h[7]->SetLineWidth(1);
    587575  h[7]->SetLineColor(2);
    588   h[7]->Draw("hist");
    589   h[7]->Draw("E1,same");
     576  h[7]->DrawCopy("E1");
     577  h[7]->Draw("hist,same");
    590578  //h[10]->SetLineColor(2);
    591579  //h[10]->Draw();
    592580  h[1]->Draw("same");
    593 
    594   TCanvas *cF = new TCanvas("cF", "", 500, 400);
    595   cF->cd();
    596   h[10]->Draw();
    597 
     581  */
     582
     583  //plot for sPlot with M0 known
     584  TCanvas *csP = new TCanvas("csP", "", 800, 400);
     585  csP->Divide(2);
     586  csP->cd(1);
     587  h[9]->SetLineColor(2);
     588  h[9]->SetLineWidth(1);
     589  h[9]->DrawCopy("E1");
     590  h[1]->SetLineColor(kGreen+3);
     591  h[1]->Draw("same");
     592  h[9]->SetLineColor(1);
     593  h[9]->SetLineWidth(3);
     594  h[9]->DrawCopy("hist,same");
     595  csP->cd(2);
     596  twoDim_mt_HT_Sig->Draw("colz,0");
     597
     598  //plot for sPlot with M0 UNknown
    598599  TCanvas *cesP = new TCanvas("cesP", "", 800,400);
    599600  cesP->Divide(2);
    600601  cesP->cd(1);
    601   h[7]->Draw("E1");
     602  h[7]->SetLineColor(2);
     603  h[7]->SetLineWidth(1);
     604  h[7]->DrawCopy("E1");
     605  h[7]->SetLineColor(1);
     606  h[7]->SetLineWidth(3);
    602607  h[7]->Draw("hist,same");
    603608  cesP->cd(2);
    604   twoDim_mt_HT_Sig_esP->Draw("colz,0");
     609  twoDim_mt_HT_Sig_esP->Draw("colz,0"); 
    605610}
Note: See TracChangeset for help on using the changeset viewer.