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

void RecoverCSVHisto (TH1D *his, TH1D *CSV1, TH1D *CSV2){

  for(int ireg=0;ireg<20;ireg++){
    TH1D *htempreg = new TH1D ("","",20,0.,1.);       
    for(int ibin=1;ibin<=20;ibin++){
      htempreg->SetBinContent(ibin,his->GetBinContent(ibin+20*ireg));	  
      htempreg->SetBinError(ibin,his->GetBinError(ibin+20*ireg));	  
    } // for(ibin)
    CSV2 -> Add(htempreg);
    double bYields, bError;
    bYields = htempreg->IntegralAndError(1,20,bError);
    CSV1 -> SetBinContent((ireg+1),bYields);
    CSV1 -> SetBinError((ireg+1),bError);
  } // for(ireg)  

}

void draw_from_trees(TString var, TCut other_cuts,
		     TString weights, TString title, int nbinsx, 
		     double xlow, double xup,
		     TString options="plotSig:plotLog:plotData",
		     double cut_low=-1, double cut_high=-1,
		     TString plot_title="default")
{

  bool plotSig = options.Contains("plotSig") && (!options.Contains("!plotSig"));
  bool plotLog = options.Contains("plotLog") && (!options.Contains("!plotLog"));
  bool plotData = options.Contains("plotData") && (!options.Contains("!plotData"));
  bool sigStack = options.Contains("sigStack") && (!options.Contains("!sigStack"));

  // Book histograms
  TH1D * httbar = new TH1D("ttbar" , title, nbinsx, xlow, xup);
  TH1D * hqcd = new TH1D("qcd" , title, nbinsx, xlow, xup);
  TH1D * hznn = new TH1D("znn" , title, nbinsx, xlow, xup);
  TH1D * hwjets = new TH1D("wjets" , title, nbinsx, xlow, xup);
  TH1D * hother = new TH1D("other" , title, nbinsx, xlow, xup);
  TH1D * hmc_exp = new TH1D("mc_exp" , title, nbinsx, xlow, xup);
  TH1D * hsingle_top = new TH1D("single_top" , title, nbinsx, xlow, xup);  
  TH1D * ht1bbbb_1500_100 = new TH1D("t1bbbb_1500_100" , title, nbinsx, xlow, xup);
  TH1D * ht1bbbb_1000_900 = new TH1D("t1bbbb_1000_900" , title, nbinsx, xlow, xup);
  TH1D * ht1tttt_1500_100 = new TH1D("t1tttt_1500_100" , title, nbinsx, xlow, xup);
  TH1D * ht1tttt_1200_800 = new TH1D("t1tttt_1200_800" , title, nbinsx, xlow, xup);
  TH1D * ht1qqqq_1400_100 = new TH1D("t1qqqq_1400_100" , title, nbinsx, xlow, xup);
  TH1D * ht1qqqq_1000_800 = new TH1D("t1qqqq_1000_800" , title, nbinsx, xlow, xup);

  // Format cuts
  TCut cut(other_cuts);

  // TCut ttbar_weight("(weightppb*4000)/top_pt_weight_official");

  TCut ttbar_weight("(3.17760399999999981e-05*4000)");

  cout << "Filling histograms for " << var.Data() << endl;
  ttbar_ch->Project("ttbar",var,(cut)*ttbar_weight);
  qcd_ch->Project("qcd",var,cut*weights);
  znn_ch->Project("znn",var,cut*weights);
  wjets_ch->Project("wjets",var,(cut)*weights);
  other_ch->Project("other",var,cut*weights);
  single_top_ch->Project("single_top",var,cut*weights);
  t1bbbb_1500_100_ch->Project("t1bbbb_1500_100",var,(cut)*weights);
  t1bbbb_1000_900_ch->Project("t1bbbb_1000_900",var,(cut)*weights);
  t1tttt_1500_100_ch->Project("t1tttt_1500_100",var,(cut)*weights);
  t1tttt_1200_800_ch->Project("t1tttt_1200_800",var,(cut)*weights);
  t1qqqq_1400_100_ch->Project("t1qqqq_1400_100",var,(cut)*weights);
  t1qqqq_1000_800_ch->Project("t1qqqq_1000_800",var,(cut)*weights);


  bool addOverflow(true);
  Double_t e_overflow(0.), i_overflow(0.);
  if (addOverflow) {
    i_overflow=httbar->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    httbar->SetBinContent(nbinsx, i_overflow);
    httbar->SetBinError(nbinsx, e_overflow);
    i_overflow=hqcd->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    hqcd->SetBinContent(nbinsx, i_overflow);
    hqcd->SetBinError(nbinsx, e_overflow);
    i_overflow=hznn->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    hznn->SetBinContent(nbinsx, i_overflow);
    hznn->SetBinError(nbinsx, e_overflow);
    i_overflow=hwjets->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    hwjets->SetBinContent(nbinsx, i_overflow);
    hwjets->SetBinError(nbinsx, e_overflow);
    i_overflow=hsingle_top->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    hsingle_top->SetBinContent(nbinsx, i_overflow);
    hsingle_top->SetBinError(nbinsx, e_overflow);
    i_overflow=hother->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    hother->SetBinContent(nbinsx, i_overflow);
    hother->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1bbbb_1500_100->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1bbbb_1500_100->SetBinContent(nbinsx, i_overflow);
    ht1bbbb_1500_100->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1bbbb_1000_900->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1bbbb_1000_900->SetBinContent(nbinsx, i_overflow);
    ht1bbbb_1000_900->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1tttt_1500_100->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1tttt_1500_100->SetBinContent(nbinsx, i_overflow);
    ht1tttt_1500_100->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1tttt_1200_800->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1tttt_1200_800->SetBinContent(nbinsx, i_overflow);
    ht1tttt_1200_800->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1qqqq_1400_100->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1qqqq_1400_100->SetBinContent(nbinsx, i_overflow);
    ht1qqqq_1400_100->SetBinError(nbinsx, e_overflow);
    i_overflow=ht1qqqq_1000_800->IntegralAndError(nbinsx,nbinsx+1,e_overflow);
    ht1qqqq_1000_800->SetBinContent(nbinsx, i_overflow);
    ht1qqqq_1000_800->SetBinError(nbinsx, e_overflow);
  }

  
  // Add up MC histograms
  hmc_exp->Add(httbar);
  hmc_exp->Add(hqcd);
  hmc_exp->Add(hznn);
  hmc_exp->Add(hwjets);
  hmc_exp->Add(hsingle_top);
  hmc_exp->Add(hother);

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

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

  


    // If there is no template for this region, go back to the standard hybrid
    if (doHybridSimple || (doHybridInclusiveTemplate && h_MT2Template == 0)) {
      // hybrid method: use nominal MC CR yield histogram to determine how many MT2 bins to use
      //  by default: use all MT2 bins integrated (no bin-by-bin).
      //  choose the last bin to try to have at least hybrid_nevent_threshold integrated events

      // Calculate last bin on local histogram
      for ( int ibin=1; ibin <= hDY->GetNbinsX(); ++ibin ) {
	float top = 0, integratedYield = 0;
	//if (hDataEM) top = hDataEM->Integral(ibin,-1)*rSFOF;
	integratedYield = hDY->Integral(ibin,-1) - top;
	if (integratedYield < hybrid_nevent_threshold) {
	  if (ibin == 1) lastbin_hybrid = 1;
	  else {
	    lastbin_hybrid = ibin-1;
	    lastmt2val_hybrid = hDY->GetBinLowEdge(ibin);
	  }
	  break;
	}
      }
      cout<<"lastbin_hybrid for doHybridSimple: "<<lastbin_hybrid<<endl;
    }
    
    

    TH1D* ratio = (TH1D*) hZinv->Clone("ratio");
    ratio->Divide(hDY);
    double errNum, errDen;
    float ratioValue = hZinv->IntegralAndError(1,-1,errNum) / hDY->IntegralAndError(1,-1,errDen);
    float ratioErr = ratioValue*sqrt(pow(errNum/hZinv->Integral(), 2) + pow(errDen/hDY->Integral(),2));

    TH1D* CRyield = (TH1D*) hData->Clone("h_mt2binsCRyield");

    TH1D* purityMC = (TH1D*) hDY->Clone("h_mt2binsPurityMC");
    if (hTop) purityMC->Add(hTop, -1);
    purityMC->Divide(hDY);

    TH1D* purityData = (TH1D*) hData->Clone("h_mt2binsPurityData");
    if (hDataEM) purityData->Add(hDataEM, -1*rSFOF); 
    purityData->Divide(purityData, hData, 1, 1, "B");
    
    TH1D* Stat = (TH1D*) CRyield->Clone("h_mt2binsStat");
    Stat->Multiply(purityData);
    Stat->Multiply(ratio);

    TH1D* Syst = (TH1D*) Stat->Clone("h_mt2binsSyst");
    TH1D* pred = (TH1D*) Stat->Clone("h_mt2bins");
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { 
      Syst->SetBinError(ibin, (1-purityData->GetBinContent(ibin))*0.2*Stat->GetBinContent(ibin));
      double quadrature = Stat->GetBinError(ibin)*Stat->GetBinError(ibin) + Syst->GetBinError(ibin)*Syst->GetBinError(ibin);
      pred->SetBinError(ibin, sqrt(quadrature));
    }
    //pred->Print("all");

    // Inputs to cardMaker
    TH1D* ratioCard  = (TH1D*) ratio->Clone("ratioCard");
    TH1D* purityCard  = (TH1D*) purityData->Clone("purityCard");   
    TH1D* CRyieldCard  = (TH1D*) CRyield->Clone("CRyieldCard");

    TH1D *CRyieldEM = 0, *CRyieldEMCard = 0;
    if (hDataEM){

//.........这里部分代码省略.........
//      cout<<"could not find histogram "<<fullhistnamePT<<" in file "<<fZinv->GetName()<<endl;
//      continue;
//    }
    cout<<"Looking at histo "<<fullhistname<<endl;      
    if (hGJet->GetNbinsX() != hZinv->GetNbinsX() ) {
      cout<<"different binning for histograms "<<fullhistname<<endl;
      continue;
    }
    //hGJet->Print("all");
    hGJet->Scale(kFactorGJetForRatio); // The goal is LO(Z) / LO(gamma)

    // Make directory and plot(s) in the output file
    TDirectory* dir = 0;
    dir = (TDirectory*)outfile->Get(directory.Data());
    if (dir == 0) {
      dir = outfile->mkdir(directory.Data());
    } 
    dir->cd();
    
    TH1D* Stat = (TH1D*) hZinv->Clone("h_mt2binsStat");
    
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { // "<=" to deal with overflow bin
      if (hGJet->GetBinContent(ibin) > 0)
	Stat->SetBinError(ibin, hZinv->GetBinContent(ibin)/sqrt( hGJet->GetBinContent(ibin) ));
      else Stat->SetBinError(ibin, hZinv->GetBinContent(ibin));
    }
    

    // Zgamma ratio in each MT2bin -> to get MC stat error on ratio
    TH1D* ratio = (TH1D*) hZinv->Clone("h_mt2binsRatio");
    //ratio->Print("all");
    //hGJet->Print("all");
    ratio->Divide(hGJet);
    //ratio->Print("all");

    TH1D* ratioInt = (TH1D*) hZinv->Clone("h_mt2binsRatioInt");
    double nGammaErr = 0;
    double nGamma = hGJet->IntegralAndError(0, -1, nGammaErr);
    for ( int ibin = 0; ibin <= ratioInt->GetNbinsX(); ++ibin) {
      if (nGamma>0) {
        ratioInt->SetBinContent(ibin, ratioInt->GetBinContent(ibin)/nGamma);
        float errOld = ratioInt->GetBinError(ibin)/nGamma;
        float errNew = (nGammaErr/nGamma) * ratioInt->GetBinContent(ibin);
        //cout<<nGamma<<" "<<nGammaErr<<" "<<errOld<<" "<<errNew<<endl;
        ratioInt->SetBinError(ibin, sqrt( errOld*errOld + errNew*errNew ) );
      }
      else {
        ratioInt->SetBinContent(ibin, 0);
        ratioInt->SetBinError(ibin, 0);
      }
    }
    //ratioInt->Print("all");

    // MCStat: use relative bin error from ratio hist, normalized to Zinv MC prediction
    TH1D* MCStat = (TH1D*) hZinv->Clone("h_mt2binsMCStat");
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { 
      MCStat->SetBinError(ibin, MCStat->GetBinContent(ibin) * ratio->GetBinError(ibin) / ratio->GetBinContent(ibin) );
    }

    TH1D* Syst = (TH1D*) Stat->Clone("h_mt2binsSyst");
    TH1D* pred = (TH1D*) Stat->Clone("h_mt2bins");
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { 
      Syst->SetBinError(ibin, 0.);
      double quadrature = Stat->GetBinError(ibin)*Stat->GetBinError(ibin) + Syst->GetBinError(ibin)*Syst->GetBinError(ibin);
      pred->SetBinError(ibin, sqrt(quadrature));
    }
    //pred->Print("all");
    
    TH1D* CRyield = (TH1D*) hGJet->Clone("h_mt2binsCRyield");

    // Extrapolation to next bin: just a ratio of GJet_i/GJet_i-1, so that we can later obtain bin i prediction from bin i-1 yield
    // Instead of : GJet_i * R(Zinv_i/GJet_i), we will do GJet_i-1 * R(GJet_i/GJet_i-1) * R(Zinv_i/GJet_i)
    TH1D* PreviousBinRatio = (TH1D*) hGJet->Clone("h_mt2binsPreviousBinRatio");
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { 
      if (ibin<=1) PreviousBinRatio->SetBinContent(ibin, 1.);
      else {
	PreviousBinRatio->SetBinContent(ibin, hGJet->GetBinContent(ibin)/hGJet->GetBinContent(ibin-1));
	PreviousBinRatio->SetBinContent(ibin, hGJet->GetBinContent(ibin)/hGJet->GetBinContent(ibin-1));
      }
      PreviousBinRatio->SetBinError(ibin, 0.); // Ignore uncertainty (just MC anyway)
    }
    

    pred->Write();
    Stat->Write();
    Syst->Write();
    CRyield->Write();
    MCStat->Write();
    PreviousBinRatio->Write();
    ratio->Write();
    ratioInt->Write();
    if(hZinvHT)  hZinvHT->Write();
    if(hZinvHT2) hZinvHT2->Write();
    if(hZinvNJ)  hZinvNJ->Write();
    if(hZinvNB)  hZinvNB->Write();
    //    hZinvPT->Write();
  } // loop over signal regions

  return;
}

int purityRandNorm(TH1D* h_template, TString name , TFile * fData, TFile * fZinv, TFile * fDY, int & lastmt2val_hybrid) {

  //cout<<"purityRandNorm for template "<<name<<endl;
  //h_template->Print("all");
  int lastbin_hybrid = 0;
  lastmt2val_hybrid = 200;
  TString name_emu = name + "emu";
  TString name_zinv = name;
  name_zinv.ReplaceAll("crdy", "sr");
  TH1D* hEMU   = (TH1D*) fData->Get(name_emu);    
  TH1D* hDY    = (TH1D*) fDY->Get(name);    
  TH1D* hZinv  = (TH1D*) fZinv->Get(name_zinv); 
  if (h_template == 0) {
    cout<<"ZinvMaker::purityAndRatio  :  could not find input template"<<endl;
    return lastbin_hybrid;
  }
  if (hDY == 0 || hZinv == 0) {
    cout<<"ZinvMaker::purityAndRatio  :  could not find DY or Zinv MC histogram"<<endl;
    return lastbin_hybrid;
  }
  if (hEMU) h_template->Add(hEMU, -1*rSFOF);

  // find the last bin
  for ( int ibin=1; ibin <= hDY->GetNbinsX(); ++ibin ) {
    float integratedYield = 0;
    integratedYield = hDY->Integral(ibin,-1);
    if (integratedYield < hybrid_nevent_threshold) {
      if (ibin == 1) lastbin_hybrid = 1;
      else {
	lastbin_hybrid = ibin-1;
	lastmt2val_hybrid = hDY->GetBinLowEdge(ibin);
      }
      break;
    }
  }

  // multiply R(Znn/Zll)
  TH1D* ratio =   (TH1D*) hZinv->Clone("ratio");
  ratio->Divide(hDY);

  h_template->Multiply(ratio);

  // Get the integrals to normalize the Zinv tails
  // and the uncertainties on the CR yield (dominated by data stats in the last N bins)
  double integratedYieldErrZinv = 0;
  float integratedYieldZinv = hZinv->IntegralAndError(lastbin_hybrid, -1., integratedYieldErrZinv);
  float relativeErrorZinv = integratedYieldErrZinv/integratedYieldZinv;
  double integratedYieldErr = 0;
  float integratedYield = h_template->IntegralAndError(lastbin_hybrid,-1,integratedYieldErr); 
  float relativeError = integratedYieldErr/integratedYield;

  // Hybridize the template: last N bins have a common stat uncertainty, and they follow the Zinv MC shape
  for ( int ibin=1; ibin <= hZinv->GetNbinsX(); ++ibin ) {

    if (ibin < lastbin_hybrid) continue;

    float kMT2 = hZinv->GetBinContent(ibin) / integratedYieldZinv;
    float err = sqrt( relativeError*relativeError + relativeErrorZinv*relativeErrorZinv);
    h_template->SetBinContent(ibin, integratedYield * kMT2);
    h_template->SetBinError(ibin, integratedYield * kMT2 * err );

  }

  // Normalize it: we just need a shape after all
  h_template->Scale(1./h_template->Integral());

  //h_template->Print("all");


  return lastbin_hybrid;
}

//.........这里部分代码省略.........
            double ptCut = 20 + l*5;

            //for(int j=0; j<9; j++){
            for(int j=0; j<9; j++) {


                double iasCut = 0.00 + j*0.05;

                cout<<"---------------------------------------------"<<""<<endl;
                cout<<"pt cut = "<<ptCut<<", Ias Cut = "<<iasCut<<endl;
                //TFile *fBkg        = new TFile(Form("totalBkg" + folder + "/TotalBkg_metCutEq%.0f_ptCutEq%0.f_ECaloCutEq%.0f_IasCutEq0p%02.0f_",metCut,ptCut,ecaloCut,iasCut*100) + region +".root" ,"READ");
                //TString filename   = "signal" + folder + "/Signal_" + samples[i]  + Form("_metCutEq%.0f_ptCutEq%.0f_ECaloCutEq%.0f_IasCutEq0p%02.0f.root",metCut,ptCut,ecaloCut,iasCut*100);
                TFile *fBkg        = new TFile(Form("totalBkg" + folder + "/TotalBkg_metCutEq%.0f_ptGt%.0f_Le50000_ECaloCutEq%.0f_IasGt0p%02.0f_Le0p99_",metCut,ptCut,ecaloCut,iasCut*100) + region +".root" ,"READ");
                TString filename   = "signal" + folder + "/Signal_" + samples[i]  + Form("_metCutEq%.0f_ptGt%.0f_Le50000_ECaloCutEq%.0f_IasGt0p%02.0f_Le0p99.root",metCut,ptCut,ecaloCut,iasCut*100);
                TFile *fSignal     = new TFile(filename,"READ");
                TH1D  *hBkg        = 0;
                TH1D  *hLepton     = 0;
                TH1D  *hFake       = 0;
                TH1D  *hSignal     = 0;
                fBkg       -> GetObject("totalBkg",hBkg);
                fBkg       -> GetObject("leptonBkg",hLepton);
                fBkg       -> GetObject("fakeBkg",hFake);
                fSignal    -> GetObject(samples[i],hSignal);

                cout.precision(3);
                double nAll         = 0;
                double nAllErrorUp  = 0;
                double nAllErrorLow = 0;
                double n, nErrorUp, nErrorLow;

                TString eq = "=";

                int bin         = hBkg->GetXaxis()->FindBin(region);
                n = hBkg->IntegralAndError(bin,bin,nErrorUp);
                nErrorLow = nErrorUp;
                if(nErrorUp==0) eq = "<";
                else            eq = "=";
                nAll         = n ;
                nAllErrorUp  = nErrorUp;
                nAllErrorLow = nErrorLow;

                cout<<"nAll        = "<<nAll<<endl;
                cout<<"nAllErrorUp = "<<nAllErrorUp<<endl;
                cout<<"One Sided Upper 68% limit of bkg = "<<getOneSidedUpperLimit(nAll,0.6827)-nAll<<endl;

                cout<<"Total bkg                          = "<<nAll<<" + "<<nAllErrorUp<<" - "<<nAllErrorLow<<""<<endl;
                cout<<samples[i]<<" = "<<hSignal->GetBinContent(1)<<" +/- "<<hSignal->GetBinError(1)<<endl;
                cout<<"Statistics of signal = "<<pow(hSignal->GetBinContent(1),2)/pow(hSignal->GetBinError(1),2)<<endl;

                hBkgYield    -> SetBinContent(l+1,j+1,n);
                hBkgUnc      -> SetBinContent(l+1,j+1,nErrorUp/n);
                hSignalYield -> SetBinContent(l+1,j+1,hSignal->GetBinContent(1));

                double sOverDeltabUp    = hSignal->GetBinContent(1)/(getOneSidedUpperLimit(nAll,0.6827)-nAll);
                double sOverDeltabError = 0;
                sig1->SetBinContent(l+1,j+1,sOverDeltabUp);
                sig1->SetBinError(l+1,j+1,sOverDeltabError);


                double leptonStatError = hLepton -> GetBinError(bin);
                double fakeStatError   = hFake   -> GetBinError(bin);
                cout<<"leptonStatError  = "<<leptonStatError<<endl;
                cout<<"fakeStatError    = "<<fakeStatError<<endl;

                double leptonSysError = hLepton -> GetBinContent(bin)*1.0;  // 100% error
                double fakeSysError   = hFake   -> GetBinContent(bin)*0.2;  // 20% error

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

    if (iSpec == 1 &&  Weight==1 ) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsPtW");
    if (iSpec == 2 &&  Weight==1) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsYW");
    if (iSpec == 3 &&  Weight==1) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsCenW");
    
    if (iSpec == 1 &&  Weight==1) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsPtW");
    if (iSpec == 2 &&  Weight==1) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsYW");
    if (iSpec == 3 &&  Weight==1) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsCenW");
    
    if (iSpec == 1 &&  Weight==0) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsPt");
    if (iSpec == 2 &&  Weight==0) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsY");
    if (iSpec == 3 &&  Weight==0) genMass_1 = (TH2D*)fil1->Get("diMuonsGenInvMassVsCen");


    if (iSpec == 1 &&  Weight==0) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsPt");
    if (iSpec == 2 &&  Weight==0) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsY");
    if (iSpec == 3 &&  Weight==0) genMassS_1 = (TH2D*)fil2->Get("diMuonsGenInvMassVsCen");

    TH1D *ptaxis = (TH1D*)genMass_1->ProjectionY("ptaxis");
    
    for (Int_t ih = 0; ih < Nptbin; ih++) {
      cout<<pt_bound[ih]<<"  "<<pt_bound[ih+1]<<endl;
      int pt_bin1 = ptaxis->FindBin(pt_bound[ih]+0.0000001);
      int pt_bin2 = ptaxis->FindBin(pt_bound[ih+1]+0.0000001);
      
      cout<< pt_bin1<<"  "<< pt_bin2<<endl; 
      TH1D * genMassVsPt = (TH1D*)genMass_1->ProjectionX("genMassVsPt", pt_bin1, pt_bin2-1);
      TH1D * genMassVsPtS = (TH1D*)genMassS_1->ProjectionX("genMassVsPtS", pt_bin1, pt_bin2-1);
      
      //does not work with weight
      //gen_pt[ih] = genMassVsPt->GetEntries();
      
      double genError,genErrorS;
      gen_pt[ih] = genMassVsPt->IntegralAndError(1,8000,genError);
      gen_ptError[ih]= genError;
      
      gen_ptS[ih] = genMassVsPtS->IntegralAndError(1,8000,genErrorS);
      gen_ptErrorS[ih]= genErrorS;
      
      cout<<" gen entries : "<< gen_pt[ih]<<endl;
      cout<<"genErro "<<genError<<"  "<< gen_ptError[ih]<<endl<<endl;
    }
  }
  
  // Efficiency
  float Eff_cat_1[100]; 
  float Eff_catS_1[100]; 
  
  float errEff_cat_1[100];
  float errEff_catS_1[100];
  
  float errEff_cat_S1[100], errEff_cat_S2[100];
  float errEff_catS_S1[100], errEff_catS_S2[100];
  
  char *Xname[100] = {" ", "p_{T}^{Dimuon} (GeV/c)", "rapidity", "centrality"};
  
  double yld_cat_1[100];
  double yld_catS_1[100];
  double cyld_cat_1[100];
  double cyld_catS_1[100];
  double eyld_cat_1[100];
  double eyld_catS_1[100];
  double ceyld_cat_1[100], ceyld_catS_1[100];

  char namePt_1[500];
  char namePt_1S[500];

void PrintLine(TCut cut, string description) {

  TH1::SetDefaultSumw2(1);


  // cout << "Cut: " << description << endl;
 
  TH1D *h = new TH1D("h","",2,0,2);
  TCut weight = Form("weightppb*%f",int_lumi);
  double n_SM(0.), n_t1tttt_1500_100(0.), n_t1tttt_1200_800(0.), n_t1bbbb_1500_100(0.), n_t1bbbb_1000_900(0.);
  double sig_t1tttt_1500_100(0.), sig_t1tttt_1200_800(0.), sig_t1bbbb_1500_100(0.), sig_t1bbbb_1000_900(0.);
  Double_t e_SM(0.), e_t1tttt_1500_100(0.), e_t1tttt_1200_800(0.), e_t1bbbb_1500_100(0.), e_t1bbbb_1000_900(0.);

  // Fill hist to get each weighted yield
  t1bbbb_1500_100_ch->Project("h","1",cut*weight);
  n_t1bbbb_1500_100=h->IntegralAndError(1,3,e_t1bbbb_1500_100);
  t1bbbb_1000_900_ch->Project("h","1",cut*weight);
  n_t1bbbb_1000_900=h->IntegralAndError(1,3,e_t1bbbb_1000_900);
  t1tttt_1500_100_ch->Project("h","1",cut*weight);
  n_t1tttt_1500_100=h->IntegralAndError(1,3,e_t1tttt_1500_100);
  t1tttt_1200_800_ch->Project("h","1",cut*weight);
  n_t1tttt_1200_800=h->IntegralAndError(1,3,e_t1tttt_1200_800);

  SM_ch->Project("h","1",cut*weight);
  n_SM=h->IntegralAndError(1,3,e_SM);

  double fstat_bg = e_SM/n_SM;
  double fsyst_bg=0.2;
  double ferr_bg = TMath::Sqrt(fstat_bg*fstat_bg+fsyst_bg*fsyst_bg);

  sig_t1tttt_1500_100 = RooStats::NumberCountingUtils::BinomialObsZ(n_t1tttt_1500_100+n_SM,n_SM,ferr_bg);
  if (sig_t1tttt_1500_100<0) sig_t1tttt_1500_100=0;
  sig_t1tttt_1200_800 = RooStats::NumberCountingUtils::BinomialObsZ(n_t1tttt_1200_800+n_SM,n_SM,ferr_bg);
  if (sig_t1tttt_1200_800<0) sig_t1tttt_1200_800=0;
  sig_t1bbbb_1500_100 = RooStats::NumberCountingUtils::BinomialObsZ(n_t1bbbb_1500_100+n_SM,n_SM,ferr_bg);
  if (sig_t1bbbb_1500_100<0) sig_t1bbbb_1500_100=0;
  sig_t1bbbb_1000_900 = RooStats::NumberCountingUtils::BinomialObsZ(n_t1bbbb_1000_900+n_SM,n_SM,ferr_bg);
  if (sig_t1bbbb_1000_900<0) sig_t1bbbb_1000_900=0;

  // sig_t1tttt_1500_100 = 2*(TMath::Sqrt(n_t1tttt_1500_100+n_SM)-TMath::Sqrt(n_SM)); 
  // sig_t1tttt_1200_800 = 2*(TMath::Sqrt(n_t1tttt_1200_800+n_SM)-TMath::Sqrt(n_SM)); 
  // sig_t1bbbb_1500_100 = 2*(TMath::Sqrt(n_t1bbbb_1500_100+n_SM)-TMath::Sqrt(n_SM)); 
  // sig_t1bbbb_1000_900 = 2*(TMath::Sqrt(n_t1bbbb_1000_900+n_SM)-TMath::Sqrt(n_SM)); 

  // sig_t1tttt_1500_100 = n_t1tttt_1500_100/TMath::Sqrt(n_SM); 
  // sig_t1tttt_1200_800 = n_t1tttt_1200_800/TMath::Sqrt(n_SM); 
  // sig_t1bbbb_1500_100 = n_t1bbbb_1500_100/TMath::Sqrt(n_SM); 
  // sig_t1bbbb_1000_900 = n_t1bbbb_1000_900/TMath::Sqrt(n_SM); 

  // printf("%s, %3.1f+-%3.1f, %3.1f+-%3.1f, %3.1f+-%3.1f, %3.1f+-%3.1f , %3.1f+-%3.1f\n",
  // 	  description.c_str(), n_t1bbbb_1500_100, e_t1bbbb_1500_100, n_t1bbbb_1000_900, e_t1bbbb_1000_900,
  // 	  n_t1tttt_1500_100, e_t1tttt_1500_100, n_t1tttt_1200_800, e_t1tttt_1200_800, n_SM, e_SM);  
  // if (errors) fprintf(file, "%s, %3.1f+-%3.1f, %3.1f+-%3.1f, %3.1f+-%3.1f, %3.1f+-%3.1f , %3.1f+-%3.1f\n",
  // 		      description.c_str(), n_t1bbbb_1500_100, e_t1bbbb_1500_100, n_t1bbbb_1000_900, e_t1bbbb_1000_900,
  // 		      n_t1tttt_1500_100, e_t1tttt_1500_100, n_t1tttt_1200_800, e_t1tttt_1200_800, n_SM, e_SM);
  // else fprintf(file, "%s, %3.1f, %3.1f, %3.1f, %3.1f , %3.1f\n",
  // 	       description.c_str(), n_t1bbbb_1500_100, n_t1bbbb_1000_900,
  // 		      n_t1tttt_1500_100, n_t1tttt_1200_800, n_SM);

  printf("%s, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f , %3.1f\n",
	 description.c_str(), n_t1bbbb_1500_100, sig_t1bbbb_1500_100, n_t1bbbb_1000_900, sig_t1bbbb_1000_900,
	 n_t1tttt_1500_100, sig_t1tttt_1500_100, n_t1tttt_1200_800, sig_t1tttt_1200_800, n_SM); 
  fprintf(file, "%s, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f , %3.1f\n",
	  description.c_str(), n_t1bbbb_1500_100, sig_t1bbbb_1500_100, n_t1bbbb_1000_900, sig_t1bbbb_1000_900,
	  n_t1tttt_1500_100, sig_t1tttt_1500_100, n_t1tttt_1200_800, sig_t1tttt_1200_800, n_SM); 


  delete h;

}

void SignificanceTestDeltaR(){

setTDRStyle();
 
 //MC
TH1D* ttgamma = getSample("TTGamma", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* tt = getSample("TTJet", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* wjets = getSample("WJetsToLNu", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* DY1 = getSample("DYJetsToLL_M-10To50", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* DY2 = getSample("DYJetsToLL_M-50", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* T_tW = getSample("T_tW-channel", 1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* Tbar_tW = getSample("Tbar_tW-channel",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* ZZ = getSample("ZZtoAnything",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* WW = getSample("WWtoAnything",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* WZ = getSample("WZtoAnything",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);

//QCD
TH1D* QCD_Pt_20_30_BCtoE = getSample("QCD_Pt_20_30_BCtoE",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_Pt_20_30_EMEnriched = getSample("QCD_Pt_20_30_EMEnriched",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
//TH1D* QCD_Pt_20_MuEnrichedPt_15 = getSample("QCD_Pt_20_MuEnrichedPt_15",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_Pt_30_80_BCtoE = getSample("QCD_Pt_30_80_BCtoE",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_Pt_30_80_EMEnriched = getSample("QCD_Pt_30_80_EMEnriched",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_Pt_80_170_BCtoE = getSample("QCD_Pt_80_170_BCtoE",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_Pt_80_170_EMEnriched = getSample("QCD_Pt_80_170_EMEnriched",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
TH1D* QCD_all = getSample("QCD_Pt_20_MuEnrichedPt_15",1, Obj, RefSelection, Type, Next, Variable, RebinFact, Systematic, Cut);
  QCD_all->Add(QCD_Pt_20_30_BCtoE);
  QCD_all->Add(QCD_Pt_20_30_EMEnriched);
  QCD_all->Add(QCD_Pt_30_80_BCtoE);
  QCD_all->Add(QCD_Pt_30_80_EMEnriched);
  QCD_all->Add(QCD_Pt_80_170_BCtoE);
  QCD_all->Add(QCD_Pt_80_170_EMEnriched);

TH1D* allMC = (TH1D*)ttgamma->Clone("all MC");
  allMC->Add(tt);
  allMC->Add(wjets);
  allMC->Add(DY1);
  allMC->Add(DY2);
  allMC->Add(T_tW);
  allMC->Add(Tbar_tW);
  allMC->Add(ZZ);
  allMC->Add(WW);
  allMC->Add(WZ);
  allMC->Add(QCD_all); 

 
  	
TH1D* sigTest = new TH1D("sigTest", "#DeltaR(#gamma, jets) Significance Test", 500/RebinFact, 0, 5);
 	sigTest->GetYaxis()->SetTitle("Significance");
 	//sigTest->SetFillColor(kAzure);
 	//sigTest->SetLineColor(kAzure);
	
	if(Variable == "Photon_deltaR_muons_")
 		sigTest->GetXaxis()->SetTitle("#DeltaR(#gamma, #mu)");
	if(Variable == "Photon_deltaR_jets_")
 		sigTest->GetXaxis()->SetTitle("#DeltaR(#gamma, jets)"); 	
	if(Variable == "Photon_deltaR_electrons_")
 		sigTest->GetXaxis()->SetTitle("#DeltaR(#gamma, e)");
			
	for(int i = 1; i<ttgamma->GetNbinsX(); i++){	
 	
	double ttgammaErr, allErr;
	
	double sigVal = ttgamma->IntegralAndError(i,ttgamma->GetNbinsX(),ttgammaErr) / sqrt(ttgamma->IntegralAndError(i,allMC->GetNbinsX(), allErr));
	double sigErr = sigVal*sqrt(pow(ttgammaErr/ttgamma->Integral(i,ttgamma->GetNbinsX()),2)+pow(allErr/allMC->Integral(i,ttgamma->GetNbinsX()),2));
	
	if(sigVal > 0){
	 
 	sigTest->SetBinContent(i, sigVal);
	sigTest->SetBinError(i, sigErr);
	}
	
std::cout << sigVal << " , pm: " << sigErr << std::endl; 
	
 	}

TCanvas *c1 = new TCanvas("Plot");

	

	sigTest->Draw("E");

	TString plotName("plots/Control/" + Obj + Cut + Type + Next + Variable+"ge1b_sig");
	c1->SaveAs(plotName+".pdf");
	c1->SaveAs(plotName+".png");
	delete c1;

}

void fakerate_systematics(int istart = 0, int iend = 999)
{
  // histogram name to be used for plotting
  const char * hname = "m_smuon";

  // reference values from standard analysis

  // get default values
  setup("../config/plot.cfg");
  TH1D * hReferenceFakes = fake_estimate_1d("fullrun76", hname);
  double N_ref, R_ref; // N and RMS
  N_ref = hReferenceFakes->IntegralAndError(1, hReferenceFakes->GetNbinsX(), R_ref);

  // list of systematics
  const int nMax = 13; // dont forget to set this when you chance the array below!
  const syst_struct sel[nMax] = {
    { "fullrun75_LooseMuonRelIso_0.2", "" },
    { "fullrun75_LooseMuonRelIso_0.4", "" },
    { "fullrun75_LooseMuonRelIso_0.8", "" },
    { "fullrun75_LooseMuonRelIso_1.0", "" },

    { "fullrun75_TL_jetpt_min_40.", "" },
    { "fullrun75_TL_jetpt_min_60.", "" },
    { "fullrun75_TL_jetpt_min_70.", "" },

    { "fullrun75_TL_met_max_40.", "" },
    { "fullrun75_TL_met_max_60.", "" },
    { "fullrun75_TL_met_max_70.", "" },

    { "fullrun75_TL_mt_max_30.", "" },
    { "fullrun75_TL_mt_max_50.", "" },
    { "fullrun75_TL_mt_max_60.", "" } //,

    // { "fullrun75_smoothsys_off", "" }
  };

  // get all numbers
  TH1D * hFakes[nMax] = { 0 };
  for (int i = istart; i < TMath::Min(iend, nMax); i++) {
    setup(Form("../config/plot%s.cfg", sel[i].cfg));
    hFakes[i] = fake_estimate_1d(sel[i].sel, hname);
  }

  // compute differences
  for (int i = istart; i < TMath::Min(iend, nMax); i++) {
    double N = 0, R = 0; // N and RMS
    if (hFakes[i]) {
      N = hFakes[i]->IntegralAndError(1, hFakes[i]->GetNbinsX(), R);
      delete hFakes[i];
    }
    else {
      ERROR("could not get fakes for selection " << sel[i].sel);
      continue;
    }
    N -= N_ref;
    R = TMath::Sqrt(TMath::Abs(q(R_ref)-q(R)));
    // output table line
    cout << sel[i].sel << ": " 
	 << "N = " << 100.*N/N_ref << " +- " << 100*R/N_ref << "%" << endl;
  }
}

TGraphAsymmErrors* GetCombinedEfficiency(const TString cut = "min_delta_phi_met_N>4.", const TString preselection="")
{ 
  
  TH1::SetDefaultSumw2();

  const int nbins=16;
  const double varbins[nbins+1]={0,20,40,60,80,100,125,150,175,200,250,300,400,500,600,800,1000};


    
  TCut all(preselection);
  TCut ccut(cut);
  TCut pass(all+ccut);


  TFile *fout = new TFile("macros/qcd_control/output.root", "recreate");
  fout->cd();

  cout << "Denom. cuts: " << (string)all << endl;
  cout << "Num. cuts: " << (string)pass << endl;

  TList* pList=new TList();
  vector<TH1D> v_all, v_pass;

  for (unsigned int sample(0); sample<nSamples; sample++) {
      
    TChain * ch = new TChain();
    ch->Add(qcd_files[sample]+"/reduced_tree");


    TH1D* hMETall = new TH1D(Form("hMETall_%d",sample),"hMETall",nbins,varbins);
    TH1D* hMETpass = new TH1D(Form("hMETpass_%d",sample),"hMETpass",nbins,varbins);

    ch->Project(hMETall->GetName(),"met",all);
    ch->Project(hMETpass->GetName(),"met",pass);

    Double_t e_overflow(0.), i_overflow(0.);
    i_overflow=hMETall->IntegralAndError(nbins,nbins+1,e_overflow);
    hMETall->SetBinContent(nbins, i_overflow);
    hMETall->SetBinError(nbins, e_overflow);
    i_overflow=hMETpass->IntegralAndError(nbins,nbins+1,e_overflow);
    hMETpass->SetBinContent(nbins, i_overflow);
    hMETpass->SetBinError(nbins, e_overflow);

    printf("Sample %d: hMETpass/hMETall=%.1f/%.1f\n",sample, hMETpass->Integral(),hMETall->Integral());

      
    delete ch;

    v_all.push_back(*hMETall);
    v_pass.push_back(*hMETpass);

    delete hMETall;
    delete hMETpass;

  }

  for (unsigned int sample(0); sample<nSamples; sample++) {
    cout << "Found hist " << v_all[sample].GetName() << endl;
    v_all[sample].Write();
    v_pass[sample].Write();
    TEfficiency* pEff = new TEfficiency(v_pass[sample],v_all[sample]);
    pList->Add(pEff);
  }

  TGraphAsymmErrors* gr = TEfficiency::Combine(pList,"v,mode",nSamples,weights);
  return gr;
}

//.........这里部分代码省略.........
    //    for (int ichan=0;ichan<nrChannels;ichan++)
    //  {
    //      for (int ijet=0;ijet<nrJets;ijet++)
    //     {
    //        for (int ireg=0;ireg<nrRegions;ireg++)
    //        {
	  

          LOG(logINFO) << "ireg = " << ireg << " / " << nrRegions;
	  string regionName = r->name;//Aaron
          //string regionName = (*regions)[ireg].name; //Aaron comments out
          //Region* r = &(*regions)[ireg];// Aaron comments this out? why?

	  bool skipLoop = skipRegion(r, channelName, jetName);
	  //bool skipLoop = skipRegion(r, channelNames[ichan], jetNames[ijet]); //Aaron comments out
          if (skipLoop) continue;

          LOG(logDEBUG) << "Made it through the gauntlet";

          stringstream histName;
	  histName << channelName << "_" << (*regions)[ireg].name << "_" << jetName;//Aaron
          //histName << channelNames[ichan] << "_" << (*regions)[ireg].name << "_" << jetNames[ijet];
          TH1D* hist = (TH1D*)f->Get(histName.str().c_str());
          LOG(logDEBUG) << "hist: " << hist;
          
          if (!hist)
          {
            LOG(logERROR) << "Couldn't find hist: " << fileName.str().c_str() << "/" << histName.str() << " -- " << sampleNames[isam];
            return 1;
          }

          string nomName = sampleNames[isam] + "_" + histName.str();
          double integral, error;
          integral = hist->IntegralAndError(0, hist->GetNbinsX()+1, error);
          nominals[nomName] = make_pair(integral, error);

          if (isam == 0) nominals_tot[histName.str()] = 0;
          nominals_tot[histName.str()] += integral;

          LOG(logINFO) << "Exiting loop";
	  }//end channels loop
      }//end regions loop
  }//end samples loop
  

  //big loop that does the normalization
  LOG(logINFO) << "Normalizing histograms / computing norms";
  vector<SysInfo> removed;
  vector<SysInfo> warn;
  stringstream outFileName;
  system(("mkdir -vp rev/"+version+"/normHists").c_str());
  ofstream outFile(("rev/"+version+"/normHists/norms_"+smass+salt+(interpSigOnly?"_sig":"")+".txt").c_str());
  for (int isys=0;isys<nrSys;isys++)
  {
    Sys* sys = &(*fileNames)[isys];
    bool first = true;
    LOG(logINFO) << "Processing sys: " << sys->folder;

    bool sym = sys->fileUp == sys->fileDown; // are the variations symmetric?
    for (int isam=0;isam<nrSamples;isam++)
    {
      if (sys->sampleNames.find(sampleNames[isam]) == sys->sampleNames.end()) continue;

      string baseFolder = sys->folder+"/";
      system(("mkdir -vp rev/"+version+"/normHists/"+baseFolder+sys->fileUp).c_str());
      if (!sym) system(("mkdir -vp rev/"+version+"/normHists/"+baseFolder+sys->fileDown).c_str());