C++ TH1F::IntegralAndError方法代码示例

void ootpu_fit(TString files, TString comments="") {

  TChain* chain = new TChain("reduced_tree");
  chain->Add(files);

  TH1::SetDefaultSumw2();

  TH1F* hrelIso = new TH1F("hrelIso","", 30, 0, 0.5);
  TH1F* hOOTPU = new TH1F("hOOTPU","", 75, 0, 75);
  TH1F* hOOTPUvsIso = new TH1F("hOOTPUvsIso",";Out-of-time ints.;Isolation cut efficiency", 15, 0, 75);

  hrelIso->StatOverflows(true);
  
  int n1(0), n2(0);
  if (files.Contains("20bx25")) {
    n1=15;
    n2=75;
  }
  else if (files.Contains("S14")) {
    n1=0;
    n2=120;
  }
  else { // default: 8 TeV scenario
    n1=0;
    n2=70;
  }


  TString mu("num_gen_muons==1&&muon_reco_match>=0");
  int low = 15; 
  for (int bin(1); bin<hOOTPUvsIso->GetNbinsX()+1; bin++) {
    if (bin<4) continue;
    if (bin>4) low=low+5;
    if (low>hOOTPUvsIso->GetBinLowEdge(hOOTPUvsIso->GetNbinsX())) break;
    TString cuts = Form("(%s)&&(oot_pu>=%d&&oot_pu<%d)",mu.Data(),low,low+4);
    cout << "Cuts: " << cuts.Data() << endl;
    chain->Project(hrelIso->GetName(), "muon_relIso", cuts);
    chain->Project(hOOTPU->GetName(), "oot_pu", cuts);
    hrelIso->Scale(1/hrelIso->GetEntries());
    Double_t left(0.), lerror(0.), right(0.), rerror(0.);
    left = hrelIso->IntegralAndError(1,12,lerror);
    right = hrelIso->IntegralAndError(13,31,rerror);
    float rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
    printf("bin1: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
    hOOTPUvsIso->SetBinContent(bin,left/(left+right));
    hOOTPUvsIso->SetBinError(bin,rat_error);
  }
  
  hOOTPUvsIso->SetLineColor(1);
  hOOTPUvsIso->SetMarkerColor(1);
  hOOTPUvsIso->SetMarkerStyle(20);

  TCanvas* c1 = new TCanvas();

  hOOTPUvsIso->Draw("e1");

  // do a linear fit
  TF1 *flin = new TF1("flin", "1 ++ x", 15, 75);
  hOOTPUvsIso->Fit(flin);
  flin=hOOTPUvsIso->GetFunction("flin");
  flin->SetLineColor(kBlue);
  flin->SetLineWidth(2);

  TString plotTitle ="OOTPU_fit"+comments+".pdf";
  c1->Print(plotTitle);

  // cout << "Rejection rates" << endl;
  // Double_t left(0.), lerror(0.), right(0.), rerror(0.);
  // left = hA->IntegralAndError(1,12,lerror);
  // right = hA->IntegralAndError(13,31,rerror);
  // float rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  // printf("bin1: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  // left = hB->IntegralAndError(1,12,lerror);
  // right = hB->IntegralAndError(13,31,rerror);
  // rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  // printf("bin2: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  // left = hC->IntegralAndError(1,12,lerror);
  // right = hC->IntegralAndError(13,31,rerror);
  // rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  // printf("bin3: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  // left = hD->IntegralAndError(1,12,lerror);
  // right = hD->IntegralAndError(13,31,rerror);
  // rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  // printf("bin4: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
}

//.........这里部分代码省略.........
  chain->Project("hC_l", var, cuts);
  cuts = Form("(%s)&&(eoot_pu>=%d)",mu.Data(),n4);
  cout << "Cuts: " << cuts.Data() << endl;
  chain->Project("hD", var, cuts);
  cuts = Form("(%s)&&(loot_pu>=%d)",mu.Data(),n4);
  chain->Project("hD_l", var, cuts);
  // }
  // else {
  //  }
  
  float avg1(hA->GetMean());
  float avg2(hB->GetMean());
  float avg3(hC->GetMean());
  float avg4(hD->GetMean());

  hA->Scale(1/hA->Integral());
  hB->Scale(1/hB->Integral());
  hC->Scale(1/hC->Integral());
  hD->Scale(1/hD->Integral());

  hA->SetLineColor(1);
  hB->SetLineColor(2);
  hC->SetLineColor(3);
  hD->SetLineColor(4);

  hA->SetLineWidth(2);
  hB->SetLineWidth(2);
  hC->SetLineWidth(2);
  hD->SetLineWidth(2);

  float avg1_l(hA_l->GetMean());
  float avg2_l(hB_l->GetMean());
  float avg3_l(hC_l->GetMean());
  float avg4_l(hD_l->GetMean());

  hA_l->Scale(1/hA_l->Integral());
  hB_l->Scale(1/hB_l->Integral());
  hC_l->Scale(1/hC_l->Integral());
  hD_l->Scale(1/hD_l->Integral());

  hA_l->SetLineColor(12);
  hB_l->SetLineColor(46);
  hC_l->SetLineColor(8);
  hD_l->SetLineColor(7);

  hA_l->SetLineWidth(2);
  hB_l->SetLineWidth(2);
  hC_l->SetLineWidth(2);
  hD_l->SetLineWidth(2);

  TCanvas* c1 = new TCanvas();
  float max = TMath::Max(hA->GetMaximum(), hB->GetMaximum());
  if (hC->GetMaximum()>max) max = hC->GetMaximum();
  if (hD->GetMaximum()>max) max = hD->GetMaximum();

  hA->SetMaximum(max*1.1);
  hA->Draw("hist");
  hB->Draw("hist,same");
  hC->Draw("hist,same");
  hD->Draw("hist,same");
  hA_l->Draw("hist,same");
  hB_l->Draw("hist,same");
  hC_l->Draw("hist,same");
  hD_l->Draw("hist,same");

  TLegend* leg = new TLegend(0.42,0.6,0.9,0.9);
  leg->SetFillStyle(0);
  char label[1000];
  sprintf(label,"%d#leqEarly Ints.<%d (#mu=%3.3f)",n1,n2,avg1);
  leg->AddEntry(hA,label,"lp");
  sprintf(label,"%d#leqEarly Ints.<%d (#mu=%3.3f)",n2,n3,avg2);
  leg->AddEntry(hB,label,"lp");
  sprintf(label,"%d#leqEarly Ints.<%d (#mu=%3.3f)",n3,n4,avg3);
  leg->AddEntry(hC,label,"lp");
  sprintf(label,"Early Ints.>%d (#mu=%3.3f)",n4,avg4);
  leg->AddEntry(hD,label,"lp");
  // leg->Draw();

  TString plotTitle ="relIso_vs_early_and_late_OOTPU_"+var+comments+".pdf";
  c1->Print(plotTitle);

  cout << "Rejection rates" << endl;
  Double_t left(0.), lerror(0.), right(0.), rerror(0.);
  left = hA->IntegralAndError(1,12,lerror);
  right = hA->IntegralAndError(13,31,rerror);
  float rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  printf("bin1: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  left = hB->IntegralAndError(1,12,lerror);
  right = hB->IntegralAndError(13,31,rerror);
  rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  printf("bin2: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  left = hC->IntegralAndError(1,12,lerror);
  right = hC->IntegralAndError(13,31,rerror);
  rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  printf("bin3: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
  left = hD->IntegralAndError(1,12,lerror);
  right = hD->IntegralAndError(13,31,rerror);
  rat_error=sqrt((left*left*rerror*rerror+right*right*lerror*lerror)/((left+right)*(left+right)));
  printf("bin4: %3.2f +/- %3.3f\n", left/(left+right),rat_error);
}

void bbtt_upg_em(std::string var,int nbins, double xmin, double xmax,std::string xtitle, std::string ytitle, double sigscale=1)
{

  TFile *outDC = new TFile("hh_em_inputs.root","RECREATE");

  SetStyle(); gStyle->SetLineStyleString(11,"20 10");
  TH1::SetDefaultSumw2(1);
 
  //std::string dir = "/data/blue/Bacon/029a/Upgrade/merged_talk_jun30/";
  std::string dir = "/afs/cern.ch/work/j/jlawhorn/public/ntuples/";
  
  std::stringstream scale; scale << sigscale;
  
  //Cut definitions
  double luminosity = 3000;
  std::stringstream lumi; lumi << luminosity;
  std::string objcut = "(tauCat1==3 && tauCat2==2 && ptTau1>20 && ptTau2>20 && tauIso2<0.4 && tauIso1<0.4 && (bTag1==2||bTag1==3||bTag1==6||bTag1==7) && (bTag2==1||bTag2==3||bTag2==6||bTag2==7) && ptB1>20 && ptB2>20 && sqrt( (etaTau1-etaTau2)**2 + (phiTau1-phiTau2)**2 )>0.4)";
  std::string jetcut = objcut+"*(mTT>20 && mTT<90)*(mBB1>70 && mBB1<140)*(mt2pileup>100)*(bdtVal>-0.05)";
  //signal region
  std::string mccut = jetcut+"*eventWeight*"+lumi.str();
  std::string sigcut = jetcut+"*eventWeight*"+lumi.str();
  std::string zjetcut = jetcut+"*eventWeight*(eventType==4)*"+lumi.str();
  std::string wjetcut = jetcut+"*eventWeight*(eventType==3)*"+lumi.str();
  std::string ewkcut = jetcut+"*eventWeight*(eventType!=1)*"+lumi.str();
  //--------------------------------------------------------------------------
  
  //Get the trees
  TTree *hhtree = load(dir+"HHToTTBB_14TeV.root"); 
  TTree *tttree = load(dir+"tt.root"); 
  //TTree *vbfhtree = load(dir+"VBFToTT_14TeV_phase2.root");
  //TTree *gfhtree = load(dir+"H.root");
  //TTree *vjettree = load(dir+"Vjets.root");
  //TTree *ewktree = load(dir+"diboson.root");
  
  //-------------------------------------------------------------------------
  
  //Get histograms
  TCanvas *canv0 = MakeCanvas("canv", "histograms", 600, 600);
  canv0->cd();
  std::string vardraw;
  /*  TH1F *Ztt = new TH1F("DY","",nbins,xmin,xmax);
  vardraw = var+">>"+"DY";
  vjettree->Draw(vardraw.c_str(),zjetcut.c_str());
  InitHist(Ztt  , xtitle.c_str(), ytitle.c_str(), TColor::GetColor(248,206,104), 1001);*/
  TH1F *ttbar = new TH1F("TTbar","",nbins,xmin,xmax);
  vardraw = var+">>"+"TTbar";
  tttree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(ttbar, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(155,152,204), 1001);
  /*  TH1F *wjets = new TH1F("Wjets","",nbins,xmin,xmax);
  vardraw = var+">>"+"Wjets";
  vjettree->Draw(vardraw.c_str(),wjetcut.c_str());
  InitHist(wjets, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(222,90,106), 1001);
  TH1F *ewk = new TH1F("Ewk","",nbins,xmin,xmax);
  vardraw = var+">>"+"Ewk";
  ewktree->Draw(vardraw.c_str(),ewkcut.c_str());
  InitHist(ewk, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  TH1F *vbfh = new TH1F("VBFH","",nbins,xmin,xmax);
  vardraw = var+">>"+"VBFH";
  vbfhtree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(vbfh, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(250,202,255), 1001);
  TH1F *ggh = new TH1F("GGH","",nbins,xmin,xmax);
  vardraw = var+">>"+"GGH";
  gfhtree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(ggh, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(250,202,255), 1001);*/
  TH1F *smhh = new TH1F("SMhh","",nbins,xmin,xmax);
  vardraw = var+">>"+"SMhh";
  hhtree->Draw(vardraw.c_str(),sigcut.c_str());
  InitSignal(smhh);
  smhh->SetLineColor(kBlack);
  delete canv0;
  //---------------------------------------------------------------------------
  //Print out the yields
  Double_t error=999;
  //ofstream outfile;
  //outfile.open("yields.txt");
  //outfile << "Yields for the signal region." << std::endl;
  cout << "SM hh   "  << smhh->IntegralAndError(0,smhh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  /*  outfile << "SM h   "  << smhh->IntegralAndError(0,smhh->GetNbinsX(),error) << "+/-" << error << endl;
      outfile << "Ztt    "  << Ztt->IntegralAndError(0,Ztt->GetNbinsX(),error) << "+/-" << error << endl;*/
  cout << "ttbar    "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  /*  outfile << "ewk    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
      outfile << "wjets    "  << wjets->IntegralAndError(0,wjets->GetNbinsX(),error) << "+/-" << error << endl;*/
  //--------------------------------------------------------------------------
  //outfile.close();
  outDC->cd();
  TDirectory* lTD = outDC->mkdir("emu");
  outDC->cd(lTD->GetPath());
  ttbar->SetName("data_obs");
  ttbar->SetTitle("data_obs");
  ttbar->Write();
  /*Ztt->SetName("ZTT");
  Ztt->SetTitle("ZTT");
  Ztt->Write();*/
  ttbar->SetName("TT");
  ttbar->SetTitle("TT");
  ttbar->Write();
  /*wjets->SetName("W");
  wjets->SetTitle("W");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    }
                    if (OSSFweight[1]>0){
                      totalOS *= OSSFweight[1];
                      totalOSerror += pow(OSSFweighterror[1]/OSSFweight[1], 2);
                    }
                    totalOSerror = sqrt(totalOSerror) * totalOS;
                    if (isocut==2){
                      sumOS = OSSFweight[0] + OSSFweight[1];
                      sumOSerror = sqrt(pow(OSSFweighterror[0], 2) + pow(OSSFweighterror[1], 2));
                    }
                    if (iCR==2){
                      totalOS*=MC_weight;
                      sumOS*=MC_weight;
                      totalOSerror*=MC_weight;
                      sumOSerror*=MC_weight;
                    }
                    wgt = totalOS;
                    wgterror = totalOSerror;
                  }
                  else{
                    if (iCR==2){
                      totalOS*=MC_weight;
                      sumOS*=MC_weight;
                    }
                  }

                  hmZZTxyOS->Fill(strdraw_mZZ, strdraw_Txy, totalOS);
                  hmZZmZ1mZ2OS->Fill(strdraw_mZZ, strdraw_mZ1, strdraw_mZ2, totalOS);

                  if (iwgt>0){
/*
                    int binx = hmZZTxyOS->GetXaxis()->FindBin(strdraw_mZZ);
                    int biny = hmZZTxyOS->GetYaxis()->FindBin(strdraw_Txy);
                    hmZZTxyOS->SetBinError(binx, biny, sqrt(pow(hmZZTxyOS->GetBinError(binx, biny), 2)+pow(totalOSerror, 2)));

                    binx = hmZZmZ1mZ2OS->GetXaxis()->FindBin(strdraw_mZZ);
                    biny = hmZZmZ1mZ2OS->GetYaxis()->FindBin(strdraw_mZ1);
                    int binz = hmZZmZ1mZ2OS->GetZaxis()->FindBin(strdraw_mZ2);
                    hmZZmZ1mZ2OS->SetBinError(binx, biny, binz, sqrt(pow(hmZZmZ1mZ2OS->GetBinError(binx, biny, binz), 2)+pow(totalOSerror, 2)));
*/
                    if (isocut==2){
                      hmZZTxyOS_extra->Fill(strdraw_mZZ, strdraw_Txy, sumOS);
                      hmZZmZ1mZ2OS_extra->Fill(strdraw_mZZ, strdraw_mZ1, strdraw_mZ2, sumOS);
/*
                      int binx = hmZZTxyOS_extra->GetXaxis()->FindBin(strdraw_mZZ);
                      int biny = hmZZTxyOS_extra->GetYaxis()->FindBin(strdraw_Txy);
                      hmZZTxyOS_extra->SetBinError(binx, biny, sqrt(pow(hmZZTxyOS_extra->GetBinError(binx, biny), 2)+pow(sumOSerror, 2)));

                      binx = hmZZmZ1mZ2OS_extra->GetXaxis()->FindBin(strdraw_mZZ);
                      biny = hmZZmZ1mZ2OS_extra->GetYaxis()->FindBin(strdraw_mZ1);
                      int binz = hmZZmZ1mZ2OS_extra->GetZaxis()->FindBin(strdraw_mZ2);
                      hmZZmZ1mZ2OS_extra->SetBinError(binx, biny, binz, sqrt(pow(hmZZmZ1mZ2OS_extra->GetBinError(binx, biny, binz), 2)+pow(sumOSerror, 2)));
*/
                    }
                  }

                }

                hsip->Fill(strdraw, wgt);
                if (iwgt>0){
                  int binx = hsip->GetXaxis()->FindBin(strdraw);
                  hsip->SetBinError(binx, sqrt(pow(hsip->GetBinError(binx), 2)+pow(wgterror, 2)));
                }

              }

              if (iCR!=1 && isocut>1 && catZ2<3){
                frecord->WriteTObject(hmZZTxyOS);
                delete hmZZTxyOS;
                frecord->WriteTObject(hmZZmZ1mZ2OS);
                delete hmZZmZ1mZ2OS;
                if (isocut==2 && iwgt>0){
                  frecord->WriteTObject(hmZZTxyOS_extra);
                  delete hmZZTxyOS_extra;
                  frecord->WriteTObject(hmZZmZ1mZ2OS_extra);
                  delete hmZZmZ1mZ2OS_extra;
                }
              }

              double totalerror = 0;
              double totalintegral = hsip->IntegralAndError(1,hsip->GetNbinsX(),totalerror);

              cout << '\t' << totalintegral << '\t' << totalerror;

              delete hsip;
            }
            cout << endl;
          }
        }
        cout << endl;
      }
      cout << endl;
    }
    cout << endl;
  }
  delete tc;

  for (int ww=0; ww<2; ww++) fweight[ww]->Close();
  frecord->Close();
}

//.........这里部分代码省略.........
  
    // Z -> ee observables 
  
    TH1F* ht = new TH1F("ht","ht",10,0,10000);
  
    ht -> Sumw2() ;

    TString cutszee = "cat==2&&minDelPhiNee>4&&nVLB>=1&&nJets>=";
    stringstream njcut ; njcut << nJetsCut ;
    cutszee += njcut.str();
    cutszee += "&&";

    for (int i = 0 ; i < nBinsMET ; i++) {
      for (int j = 0 ; j < nBinsHT ; j++) {
  
        TString obs_Zee = "N_Zee" ;
        obs_Zee = obs_Zee+sMbins[i]+sHbins[j] ;
        
        TString cut = "HT>";  
        cut += Hbins[j];
        cut += "&&HT<";
        cut += Hbins[j+1];
        cut += "&&METee>";
        cut += Mbins[i];
        cut += "&&METee<";
        cut += Mbins[i+1];
  
        TString allcutsZ = cutszee+cut ;
        allcutsZ += "&&" ;
        allcutsZ += leadJetPtCutString ;
  
        dyTree->Project("ht","HT",allcutsZ);
        double allerr(0.) ;
        double allval = ht->IntegralAndError(1,10,allerr) ;
        printf(" N_Zee -- HT,MET bins (%d,%d): events=%7.1f +/- %6.1f, cuts=%s\n", j,i,allval,allerr,allcutsZ.Data() ) ; cout << flush ;
          ht->Reset() ;
  
        ////// inFile << obs_Zee << "  \t" << (int)allval << endl;
        inFile << obs_Zee << "  \t" << allval << endl;

        //Z->ee counts, with 1 VLb and sig selection, so so MET>250, HT>400, mindelphi>4, 2e, 0mu, nJets >= 3
  
        int histbin = 1 + (nBinsHT+1)*i + j + 1 ;
  
        hmctruth_fit_zee_1b -> SetBinContent( histbin, allval ) ;
        hmctruth_fit_zee_1b -> SetBinError(   histbin, allerr ) ;
  
      }
    }
  
    
    printf("\n\n-----------------------------------------------------------------\n\n") ; cout << flush ;
  
    // Z -> mm observables
  
    TString cutszmm = "cat==1&&minDelPhiNmm>4&&nVLB>=1&&nJets>=";
    cutszmm += njcut.str();
    cutszmm += "&&";

    for (int i = 0 ; i < nBinsMET ; i++) {
      for (int j = 0 ; j < nBinsHT ; j++) {
  
        TString obs_Zmm = "N_Zmm" ;
        obs_Zmm = obs_Zmm+sMbins[i]+sHbins[j] ;
        
        TString cut = "HT>";  

void draw_tt(std::string var="mTT",int nbins=10, double xmin=0, double xmax=200,std::string xtitle="m_{#tau#tau}", std::string ytitle="Events")
{
  SetStyle(); gStyle->SetLineStyleString(11,"20 10");
  TH1::SetDefaultSumw2(1);

  std::string dir = "/afs/cern.ch/work/j/jlawhorn/public/vbf/ntuples/";
  double sigscale = 10;
  double sigscale1 = 10; 
  std::stringstream scale; scale << sigscale;
  std::stringstream scale1; scale1 << sigscale1;

  //Cut definitions
  double luminosity = 3000;
  std::stringstream lumi; lumi << luminosity;
  std::string objcut = "(ptTau1>45 && ptTau2>45 && abs(etaTau1) <4.0 && abs(etaTau1)<4.0 )";
  std::string jetcut = objcut+"*(ptJet1>30 && ptJet2>30 && abs(etaJet1) <4.7 && abs(etaJet2) <4.7 )";
  std::string ththcut = jetcut+"*(tauCat1==1 && tauCat2==1 && mJJ>500 && abs(dEta)>3.5)";

  //signal region
  std::string vbfcut = ththcut+"*eventWeight*(eventType==0)*"+lumi.str();
  std::string zttcut = ththcut+"*eventWeight*(eventType==1)*"+lumi.str();
  std::string ttbarcut = ththcut+"*eventWeight*(eventType==3)*"+lumi.str();
  //std::string ewkcut = ththcut+"*eventWeight*(eventType==2)*"+lumi.str();
  std::string othercut = ththcut+"*eventWeight*(eventType==4 || eventType==2)*"+lumi.str();
  
  //--------------------------------------------------------------------------
  
  //Get the trees
  TTree *tree = load(dir+"vbf_comb.root"); 

  //-------------------------------------------------------------------------
  
  //Get histograms
  TCanvas *canv0 = MakeCanvas("canv", "histograms", 600, 600);
  canv0->cd();
  std::string vardraw;
  TH1F *vbf = new TH1F("VBFH","",nbins,xmin,xmax);
  vardraw = var+">>"+"VBFH";
  tree->Draw(vardraw.c_str(),vbfcut.c_str());
  InitSignal(vbf);
  vbf->SetLineColor(kBlack);
  TH1F *ttbar = new TH1F("TTbar","",nbins,xmin,xmax);
  vardraw = var+">>"+"TTbar";
  tree->Draw(vardraw.c_str(),ttbarcut.c_str());
  InitHist(ttbar, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(155,152,204), 1001);
  TH1F *ztt = new TH1F("Ztt","",nbins,xmin,xmax);
  vardraw = var+">>"+"Ztt";
  tree->Draw(vardraw.c_str(),zttcut.c_str());
  InitHist(ztt, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(248,206,104), 1001);
  //TH1F *ewk = new TH1F("Ewk","",nbins,xmin,xmax);
  //vardraw = var+">>"+"Ewk";
  //tree->Draw(vardraw.c_str(),ewkcut.c_str());
  //InitHist(ewk, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  TH1F *other = new TH1F("Other","",nbins,xmin,xmax);
  vardraw = var+">>"+"Other";
  tree->Draw(vardraw.c_str(),othercut.c_str());
  InitHist(other, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  
  delete canv0;
  //----------------------------------------------------------------------------
  //Print out the yields
  Double_t error=0.0;
  ofstream outfile;
  outfile.open("yields_tt.txt");
  outfile << "Yields for the signal region." << std::endl;
  outfile << "VBF   "  << vbf->IntegralAndError(0,vbf->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "TTbar   "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(0,ztt->GetNbinsX(),error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "other   "  << other->IntegralAndError(0,other->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "S/sqrt(B)    "  << vbf->Integral()/(other->Integral()+ztt->Integral()+ttbar->Integral()/*+ewk->Integral()*/) << endl;
  //--------------------------------------------------------------------------
  //continue outputing
  //outfile << "Ewk total    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << endl << endl << endl;
  outfile << "In the signal region (100,150GeV)  " <<endl;
  outfile << "VBF   "  << vbf->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "TTbar    "  << ttbar->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(5,11,error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "other    "  << other->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile.close();
  //-----------------------------------------------------------------------
  //Draw the histograms
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);
  canv->cd();
  //ewk->Add(other); ttbar->Add(ewk); 
  //ztt->Add(ttbar); vbf->Add(ztt);
  //Error band stat
  TH1F* errorBand = (TH1F*)vbf ->Clone("errorBand");
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(13);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  //  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
  //     if(errorBand->GetBinContent(idx)>0){
  //       std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
  //       break;
  //     }
  //}
//.........这里部分代码省略.........

void MutualInfoSingleVar_reducedBias(TString theFormula = "fjet1MassSDb2",TString histFileName = "fjet1MassSDb2",
                                     int numBins = 320,Double_t histMin = -20,Double_t histMax = 300){

  TFile *signalFile = new TFile("signal_word.root");
  TTree *signalTree = (TTree*) signalFile->FindObjectAny("DMSTree");
  TFile *backgdFile = new TFile("background_word.root");
  TTree *backgdTree = (TTree*) backgdFile->FindObjectAny("DMSTree");

  TH1F *signalHist = new TH1F("Signal","Signal",numBins,histMin,histMax); signalHist->Sumw2();
  TH1F *backgdHist = new TH1F("Backgd","Backgd",numBins,histMin,histMax); backgdHist->Sumw2();
  TH1F *sumHist = new TH1F("Sum","Sum",numBins,histMin,histMax); sumHist->Sumw2();

  signalTree->Draw(theFormula + TString(">>Signal"),"weight*(abs(fjet1PartonId)==24)");
  backgdTree->Draw(theFormula + TString(">>Backgd"),"weight");
  *sumHist = *signalHist + *backgdHist;

  Double_t signalErr = 0.;
  Double_t signalWeights = signalHist->IntegralAndError(1,numBins,signalErr);
  Double_t backgdErr = 0.;
  Double_t backgdWeights = backgdHist->IntegralAndError(1,numBins,backgdErr);
  Double_t sumErr = 0.;
  Double_t sumWeights = sumHist->IntegralAndError(1,numBins,sumErr);

  Double_t signalFrac = signalWeights/sumWeights;
  Double_t signalFracErr = sqrt(TMath::Power(signalErr/sumWeights,2) + TMath::Power(signalWeights*sumErr/(TMath::Power(sumWeights,2)),2));

  Double_t truthEntropy = -1*signalFrac*TMath::Log2(signalFrac) - (1-signalFrac)*TMath::Log2(1-signalFrac);
  Double_t truthErr = (TMath::Abs(1+TMath::Log2(signalFrac)) + TMath::Abs(1+TMath::Log2(1-signalFrac)))*signalFracErr;

  Double_t varErr = 0.;
  Double_t varEntropy = histEntropy(sumHist,numBins,sumWeights,sumErr,varErr);
  Double_t unionErr = 0.;
  Double_t unionEntropy = histEntropy(signalHist,numBins,sumWeights,sumErr,unionErr) +
                          histEntropy(backgdHist,numBins,sumWeights,sumErr,unionErr);

  Double_t MutualInfo = truthEntropy + varEntropy - unionEntropy;
  Double_t MutualErr = sqrt(TMath::Power(truthErr,2) + TMath::Power(varErr,2) + TMath::Power(unionErr,2));

  cout << endl;
  cout << theFormula << endl;
  cout << endl;
  cout << "Old Way" << endl;
  cout << "H(T):   " << truthEntropy << " +- " << truthErr  << endl;
  cout << "H(A):   " << varEntropy   << " +- " << varErr    << endl;
  cout << "H(T,A): " << unionEntropy << " +- " << unionErr  << endl;
  cout << "I(T;A): " << MutualInfo   << " +- " << MutualErr << endl;

  TFile *outFile = new TFile(TString("scratch/") + histFileName + TString(".root"),"RECREATE");

  /////////////////////////////////////////////////////////////////////////////////////////////////////////
  //                            Attempts to reduce bias with A.20                                        //
  /////////////////////////////////////////////////////////////////////////////////////////////////////////

  Double_t signalEntropyErr = 0.;
  Double_t signalEntropy = histEntropy(signalHist,numBins,signalWeights,signalErr,signalEntropyErr);
  Double_t backgdEntropyErr = 0.;
  Double_t backgdEntropy = histEntropy(backgdHist,numBins,backgdWeights,backgdErr,backgdEntropyErr);

  TH1F *signalHist_0 = new TH1F("Signal_0","Signal_0",numBins,histMin,histMax); signalHist_0->Sumw2();
  TH1F *backgdHist_0 = new TH1F("Backgd_0","Backgd_0",numBins,histMin,histMax); backgdHist_0->Sumw2();
  TH1F *sumHist_0 = new TH1F("Sum_0","Sum_0",numBins,histMin,histMax); sumHist_0->Sumw2();

  signalTree->Draw(theFormula + TString(">>Signal_0"),"weight*(abs(fjet1PartonId)==24 && event % 2 == 0)");
  backgdTree->Draw(theFormula + TString(">>Backgd_0"),"weight*(event % 2 == 0)"); // This event % 2 term splits the sample in half
  *sumHist_0 = *signalHist_0 + *backgdHist_0;

  Double_t sumErr_0 = 0.;
  Double_t sumWeights_0 = sumHist_0->IntegralAndError(1,numBins,sumErr);

  Double_t mixedErr_0 = 0.;
  Double_t mixedEntropy_0 = histEntropy(sumHist_0,numBins,sumWeights_0,sumErr_0,mixedErr_0);

  Double_t MutualInfo_0 = mixedEntropy_0 - signalFrac*signalEntropy - (1-signalFrac)*backgdEntropy;
  Double_t MutualErr_0  = sqrt(TMath::Power(mixedErr_0,2) +
                               TMath::Power(signalFracErr*(signalEntropy-backgdEntropy),2) +
                               TMath::Power(signalFrac*signalEntropyErr,2) +
                               TMath::Power((1-signalFrac)*backgdEntropyErr,2));

  cout << endl;
  cout << "Reduced Bias Results (0)" << endl;
  cout << "f:         " << signalFrac     << " +- " << signalFracErr    << endl;
  cout << "H_sig+bck: " << mixedEntropy_0 << " +- " << mixedErr_0       << endl;
  cout << "H_sig:     " << signalEntropy  << " +- " << signalEntropyErr << endl;
  cout << "H_bck:     " << backgdEntropy  << " +- " << backgdEntropyErr << endl;
  cout << "I(T;A):    " << MutualInfo_0   << " +- " << MutualErr_0      << endl;

  TH1F *signalHist_1 = new TH1F("Signal_1","Signal_1",numBins,histMin,histMax); signalHist_1->Sumw2();
  TH1F *backgdHist_1 = new TH1F("Backgd_1","Backgd_1",numBins,histMin,histMax); backgdHist_1->Sumw2();
  TH1F *sumHist_1 = new TH1F("Sum_1","Sum_1",numBins,histMin,histMax); sumHist_1->Sumw2();

  signalTree->Draw(theFormula + TString(">>Signal_1"),"weight*(abs(fjet1PartonId)==24 && event % 2 == 1)");
  backgdTree->Draw(theFormula + TString(">>Backgd_1"),"weight*(event % 2 == 1)"); // This event % 2 term splits the sample in half
  *sumHist_1 = *signalHist_1 + *backgdHist_1;

  Double_t sumErr_1 = 0.;
  Double_t sumWeights_1 = sumHist_1->IntegralAndError(1,numBins,sumErr);

  Double_t mixedErr_1 = 0.;
  Double_t mixedEntropy_1 = histEntropy(sumHist_1,numBins,sumWeights_1,sumErr_1,mixedErr_1);

//.........这里部分代码省略.........

void listyieldsWZ_stat(std::string rootfile, int wzstep, std::string proc,int var,int lhe){
	string WZstepname [7]= {"Three leptons",
	"WZ candidate",
	"Z sel",
	"W sel",
	"M(3l) > 100 GeV",
	"1 or 0 b-jets",
	"0 b-jets"};
	
	string variation [5]= {"Central",
	"Up b-tag",
	"Down b-tag",
	"Up b-mistag",
	"Down b-mistag"};
	
	TH1F::SetDefaultSumw2();
	
	TH1F *histo = nullptr;
	TFile *f = new TFile(rootfile.c_str());
	if (var == 0) histo = (TH1F*)f->Get(Form("MET_WZSMstep%i/%s",wzstep, proc.c_str() ));
	if (var == 1) histo = (TH1F*)f->Get(Form("MET_WZSMstep%iUncBTAGUp/%s",wzstep, proc.c_str() ));
	if (var == 2) histo = (TH1F*)f->Get(Form("MET_WZSMstep%iUncBTAGDo/%s",wzstep, proc.c_str() ));
	if (var == 3) histo = (TH1F*)f->Get(Form("MET_WZSMstep%iUncBMISTAGUp/%s",wzstep, proc.c_str() ));
	if (var == 4) histo = (TH1F*)f->Get(Form("MET_WZSMstep%iUncBMISTAGDo/%s",wzstep, proc.c_str() ));
	
	if (var == 5) histo = (TH1F*)f->Get(Form("mZ1_WZaccden/%s", proc.c_str() ));
	if (var == 6) histo = (TH1F*)f->Get(Form("mZ1_WZaccnum/%s", proc.c_str() ));
	
	Double_t error;
	Int_t yield = histo->GetEntries();
	Double_t integral = histo->Integral(0,100000000);
	
	Double_t normfactor;
	
	if (proc.find("WZTo3LNu") != std::string::npos) normfactor=4.42965/1980800;
	else if (proc.find("DYJetsToLL_M10to50") != std::string::npos) normfactor=18610/2.967885e+07;
	else if (proc.find("DYJetsToLL_M50") != std::string::npos) normfactor=6025.2/2.874797e+07;
	else if (proc.find("_TW_") != std::string::npos) normfactor=35.6/995600;
	else if (proc.find("TbarW") != std::string::npos) normfactor=35.6/988500;
	else if (proc.find("TToLeptons_sch") != std::string::npos) normfactor=3.68064/984400;
	else if (proc.find("TToLeptons_tch") != std::string::npos) normfactor=44.071956/3299800;
	else if (proc.find("TbarToLeptons_tch") != std::string::npos) normfactor=26.2278/1680200;
	else if (proc.find("TTJets") != std::string::npos) normfactor=831.8/8.970451e+07;
	else if (proc.find("TTWToLNu") != std::string::npos) normfactor=0.2043/252908;
	else if (proc.find("TTZToLLNuNu") != std::string::npos) normfactor=0.2529/398000;
	else if (proc.find("TTHnobb") != std::string::npos) normfactor=0.2151/1562609;
	else if (proc.find("TTGJets") != std::string::npos) normfactor=3.697/4832230;
	else if (proc.find("WJetsToLNu") != std::string::npos) normfactor=61526.7/2.239467e+07;
	else if (proc.find("WWTo2L2Nu") != std::string::npos) normfactor=12.178/1965200;
	else if (proc.find("WZTo3LNu") != std::string::npos) normfactor=4.42965/1980800;
	else if (proc.find("_ZZ_") != std::string::npos) normfactor=16.523/996944;
	else if (proc.find("GGHZZ4L") != std::string::npos) normfactor=0.01212/498000;
	else if (proc.find("WGToLNuG") != std::string::npos) normfactor=585.8/6099599;
	else if (proc.find("_WWZ_") != std::string::npos) normfactor=0.1651/250000;
	else if (proc.find("_WZZ_") != std::string::npos) normfactor=0.05565/250000;
	else if (proc.find("_ZZZ_") != std::string::npos) normfactor=0.01398/250000;
	else if (proc.find("_ZGTo2LG_") != std::string::npos) normfactor=117.864/4451319;

	histo->Scale(normfactor*2260);
    //Double_t yield_scaled = integral*normfactor*2260;
    Double_t error_scaled;
    Double_t yield_scaled = histo->IntegralAndError(0,1000000,error_scaled);
    
    /*
    if (rootfile.find("all") != std::string::npos){
      cout << lhe << "\t";
      if (var==5) cout << "Den: " << integral;
      if (var==6) cout << "Num: " << integral << endl;
    }
    */
    
    
	if (rootfile.find("all") != std::string::npos){
		cout << WZstepname[wzstep] <<"\t" << yield_scaled << " $\\pm$ " << error_scaled << endl;
	} else if (rootfile.find("eee") != std::string::npos) {
		cout << "Selection step:	" << WZstepname[wzstep] << "	" << "eee	" << yield << endl;
	} else if (rootfile.find("eem") != std::string::npos) {
		cout << "Selection step:	" << WZstepname[wzstep] << "	" << "eem	" << yield << endl;
	} else if (rootfile.find("mme") != std::string::npos) {
		cout << "Selection step:	" << WZstepname[wzstep] << "	" << "mme	" << yield << endl;
	} else if (rootfile.find("mmm") != std::string::npos) {
		cout << "Selection step:	" << WZstepname[wzstep] << "	" << "mmm	" << yield << endl;
	}
	
}

//.........这里部分代码省略.........
  TH1F *hh_p5;
  if(hist)
    hh_p5= new TH1F("hh_p5","",nbins,xmin,xmax);
  else
    hh_p5= new TH1F("hh_p5","",13,massLEdges);
  vardraw = var+">>"+"hh_p5";
  hhtree_p5->Draw(vardraw.c_str(),sigcut.c_str());
  InitSignal(hh_p5);
  hh_p5->SetLineColor(kBlack);
  TH1F *hh_vbf;
  if(hist)
    hh_vbf= new TH1F("hh_vbf","",nbins,xmin,xmax);
  else
    hh_vbf= new TH1F("hh_vbf","",13,massLEdges);
  vardraw = var+">>"+"hh_vbf";
  hhtree_vbf->Draw(vardraw.c_str(),sigcut.c_str());
  InitSignal(hh_vbf);
  hh_vbf->SetLineColor(kBlue);
  delete canv0;
  //for(int i0 = 0; i0 < ewk->GetNbinsX()+1; i0++)
  //  if(ewk->GetBinContent(i0) == 0)
  //    ewk->SetBinError(i0,2.8);
  //---------------------------------------------------------------------------
  //Print out the yields
  Double_t error=999;
  Double_t sigN=0;
  Double_t sigSig=0;
  Double_t bgdN=0;
  Double_t bgdSig=0;
  //ofstream outfile;
  //outfile.open("yields.txt");
  //outfile << "Yields for the signal region." << std::endl;
  cout << jetcut << endl;
  cout << "SM hh   "  << smhh->IntegralAndError(0,smhh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  cout << " VBF HH "  << hh_vbf->IntegralAndError(0,hh_vbf->GetNbinsX(),error) << "+/-";
  sigN+=smhh->IntegralAndError(0,smhh->GetNbinsX(),error);
  sigSig+=error;
  cout << error << endl; error=999;
  cout << "ttbar   "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error);bgdSig+=error*error;
  cout << "Ztt     "  << Ztt->IntegralAndError(0,Ztt->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=Ztt->IntegralAndError(0,Ztt->GetNbinsX(),error);bgdSig+=error*error;
  cout << "ewk     "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=ewk->IntegralAndError(0,ewk->GetNbinsX(),error);bgdSig+=error*error;
  cout << "wjets   "  << wjets->IntegralAndError(0,wjets->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=wjets->IntegralAndError(0,wjets->GetNbinsX(),error);bgdSig+=error*error;
  cout << "ggH     "  << ggh->IntegralAndError(0,ggh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=ggh->IntegralAndError(0,ggh->GetNbinsX(),error);bgdSig+=error*error;
  cout << "vbfH    "  << vbfh->IntegralAndError(0,vbfh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=vbfh->IntegralAndError(0,vbfh->GetNbinsX(),error);bgdSig+=error*error;
  cout << "vH/ttH  "  << assoh->IntegralAndError(0,assoh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  bgdN+=assoh->IntegralAndError(0,assoh->GetNbinsX(),error); bgdSig+=error*error;
  cout << "S = " << sigN << "+/-" << sigSig << endl;
  cout << "B = " << bgdN << "+/-" << TMath::Sqrt(bgdSig) << endl;
  cout << "S/sqrt(B) = " << sigN/TMath::Sqrt(bgdN) << endl;

  //--------------------------------------------------------------------------
  //outfile.close();