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

void zeeAnaDiff(TH1D* hsee, TCut addCut,TH1D* hseeMc) {
  
  TFile *f1=new TFile(fname1.Data());
  TTree *zp = (TTree*)f1->Get("tz");
  zp->AddFriend("yEvt=yongsunHiEvt"       ,fname1.Data());
  zp->AddFriend("ySkim=yongsunSkimTree"   ,fname1.Data());
  zp->AddFriend("yHlt=yongsunHltTree"     ,fname1.Data());
  zp->AddFriend("tgj",                     fname1.Data());
  
  TH1D* hmass = new TH1D("hmass",";inv. mass (GeV);event fraction",34,30,200);
  hmass->Sumw2();

  zp->Draw("invm>>hmass",addCut && "isEle1==1 && isEle2==1");
  TH1D* htemp = (TH1D*)hsee->Clone("htemp");   htemp->Reset();
  zp->Draw("see1>>htemp",addCut && "isEle1==1 && isEle2==1" && "invm>80 && invm<100" );
  hsee->Reset();
  hsee->Add(htemp);
  handsomeTH1(hmass,2);
  scaleInt(hmass);
  
  TH1D* htempMC;
  TFile *f1mc;
  TTree *zpmc ;
  f1mc=new TFile(fname1MC.Data());
  zpmc = (TTree*)f1mc->Get("tz");
  zpmc->AddFriend("yEvt=yongsunHiEvt"       ,fname1MC.Data());
  zpmc->AddFriend("ySkim=yongsunSkimTree"   ,fname1MC.Data());
  zpmc->AddFriend("yHlt=yongsunHltTree"     ,fname1MC.Data());
  zpmc->AddFriend("yPho=yongsunPhotonTree"     ,fname1MC.Data());
  zpmc->AddFriend("tgj",                     fname1MC.Data());
  hseeMc->Sumw2();
  htempMC = (TH1D*)hseeMc->Clone("htempMc");   htempMC->Reset();
  TH1D* hmassMC = new TH1D("hmassMC11",";inv. mass (GeV);event fraction",34,30,200);
  hmassMC->Sumw2();
  zpmc->Draw("invm>>hmassMC11",Form("(%s) *yPho.ncoll",addCut.GetTitle()));
  zpmc->Draw("see1>>htempMc",Form("(%s) *yPho.ncoll", (addCut && "invm>80 && invm<100").GetTitle()) );
  hseeMc->Reset();
  hseeMc->Add(htempMC);
  handsomeTH1(hmassMC,1);
  scaleInt(hmassMC);
  
  
  hmassMC->SetAxisRange(0,0.6,"Y");
  hmassMC->DrawCopy();
  hmass->DrawCopy("same");                                                                                             
  TLegend* leg1 =  new TLegend(0.5580645,0.5381356,0.8923387,0.7648305,NULL,"brNDC");
  easyLeg(leg1,"Z->ee mass");
  leg1->AddEntry(hmass,"DATA","pl");
  leg1->AddEntry(hmassMC,"MC","pl");
  leg1->Draw();
  
}

// ============================================================================
StatusCode Gaudi::Utils::Histos::fromXml
( TH1D& result , const std::string& input )
{
  //
  result.Reset() ;                                 // RESET old histogram
  //
  _Xml<TH1D> _xml ;
  std::auto_ptr<TH1D> histo =  _xml ( input ) ;
  if ( 0 == histo.get() ) { return StatusCode::FAILURE ; }        // RETURN
  //
  result.Reset() ;
  histo->Copy ( result ) ;
  //
  return StatusCode::SUCCESS ;
}

void  makeGaussianSignals(SigData_t& m_sigdata)
{
  //std::vector<TH1D *> vgsh(NUMCHAN);
  std::vector<TH1D *> vcdfh(NUMCHAN);

  if( m_sigdata.find("gs") == m_sigdata.end() ) {
    cerr << "Gaussian signal data not found, not making CDF signal!" << endl;
    return;
  }

  for (int ichan=0; ichan<NUMCHAN; ichan++) {
    TH1D *cdfh;

    TString channame(channames[ichan]);
    TString name;

    TH1D * gsh = m_sigdata["gs"].at(ichan);

    assert(gsh) ;

#if 0
    name = "Signalgs_"+channame;
    gsh = (TH1D *)tch->Clone(name.Data());

    assert(gsh);

    gsh->SetTitle("Gaussian signal");
    
    gsh->Reset();

    TF1 *g = (TF1 *)gROOT->GetFunction("gaus");
    g->SetParameters(1,gaussian_mean_mass_gev,gaussian_mass_sigma_gev);
    gsh->FillRandom("gaus",100000);

    // norm to 1pb signal with 1/fb integrated luminosity
    double norm = 1000 * gseffxacc[ichan]/gsh->Integral(0,gsh->GetNbinsX()+1,"width");

    //gsh->Scale(norm/eff_fudgefactor); // kludge: pre-undo the fudge in the next module
    gsh->Scale(norm);

    vgsh[ichan] = gsh;
#endif

    // New CDF bump, same as Gauss but set to CDF (obs/theor)*(LHC theor) = 3.43pb
    cdfh = (TH1D *)gsh->Clone("CDFbump");

    cdfh->Scale(3.43);

    vcdfh[ichan] = cdfh;

    cdfh->Draw();

    gsh->Draw("same");

  } // channel loop

  //m_sigdata["gs"]  = vgsh;
  m_sigdata["cdf"] = vcdfh;
}                                                           // makeGaussianSignals

///////////////////////////////////////////////////////////////////
//////// Go4 GUI example script fft.C
//          J.Adamczewski, gsi, 30 May 2012
// NOTE: to be run in Go4 GUI local command line only!
//       NEVER call this script in remote analysis process!!!
/////// Functionality:
// perfroms fft on histogram of name1 using the option as explained in root TVirtualFFT:FFT
/////// Usage:
// The draw flag switches if the results are displayed each time this macro is called
// if display is switched off, the result histogram is just updated in browser and existing displays
///////
Bool_t fft(const char* name1, Option_t*  opt = "R2C M", Bool_t draw=kTRUE)
{
   if(TGo4AbstractInterface::Instance()==0 || go4!=TGo4AbstractInterface::Instance()) {
      std::cout <<"FATAL: Go4 gui macro executed outside Go4 GUI!! returning." << std::endl;
      return kFALSE;
   }
   TString newname;
   TString fullname1 = go4->FindItem(name1);
   TObject* ob1 = go4->GetObject(fullname1,1000); // 1000=timeout to get object from analysis in ms

   if ((ob1==0) || !ob1->InheritsFrom("TH1")) {
     std::cout <<"fft could not get histogram "<<fullname1 << std::endl;
     return kFALSE;
   }

   if(ob1->InheritsFrom("TH2") || ob1->InheritsFrom("TH3")){  // 2d
      std::cout <<"fft does not support 2d/3d histogram "<<fullname1 << std::endl;
      return kFALSE;
   }

   TH1* his1=(TH1*)ob1;
   TString n1=his1->GetName();
   TString t1=his1->GetTitle();
   newname.Form("_fft_%s",opt);
   TString finalname = n1+newname;
   TString finaltitle = t1+newname;


   // do fft here:
   Int_t N = his1->GetNbinsX();
   TH1D* result = new TH1D(finalname, finaltitle,N,0,N);
   result->SetName(finalname);
   result->SetTitle(finaltitle);
   result->Reset("");
   Double_t *in = new Double_t[N];
   // since we do not know type of input histo, we copy contents to Double array:
   for(Int_t ix=0; ix<N;++ix)
   {
      in[ix]=his1->GetBinContent(ix+1);
   }
   TVirtualFFT *thefft = TVirtualFFT::FFT(1, &N, opt);
   thefft->SetPoints(in);
   thefft->Transform();
   Double_t re, im;
   for (Int_t i=0; i<N; i++) {
      thefft->GetPointComplex(i, re, im);
      result->SetBinContent(i+1,TMath::Sqrt(re*re + im*im));
   }
   result->SetDirectory(0);

   TString rname = go4->SaveToMemory("FFT", result, kTRUE);
   std::cout<< "Saved result histogram to " << rname.Data() <<std::endl;
   if(draw){
      ViewPanelHandle vpanel = go4->StartViewPanel();
      go4->DrawItem(rname, vpanel);
   }
   return kTRUE;
}

TH1D* merge(TH1D* bg1, TH1D* bg2, TH1D* signal)
{
  TList *list = new TList;
  list->Add(bg1);
  list->Add(bg2);
  list->Add(signal);
  ostringstream newname; newname << signal->GetName() << "_merged";
  TH1D* newsignal = (TH1D*)signal->Clone(newname.str().c_str());
  newsignal->Reset();
  newsignal->Merge(list);
  return newsignal;
}

void getseePhoMC(TH1D* hsee, TCut addCut){

  TFile *f1=new TFile(fPho.Data());
  TTree *photon1 = (TTree*)f1->Get("yongsunPhotonTree");
  photon1->AddFriend("yEvt=yongsunHiEvt"       ,fPho.Data());
  photon1->AddFriend("ySkim=yongsunSkimTree"   ,fPho.Data());
  photon1->AddFriend("yHlt=yongsunHltTree"     ,fPho.Data());
  photon1->AddFriend("tgj",                     fPho.Data());
 
  TH1D* htemp = (TH1D*)hsee->Clone("htemp");  htemp->Reset();
  photon1->Draw("sigmaIetaIeta>>htemp",Form("(%s) *ncoll",addCut.GetTitle()));
  hsee->Reset();
  hsee->Add(htemp);
}

TH1D* mergeSys( TH1D* h1, TH1D* h2) {

  TH1D* hres = (TH1D*)h1->Clone(Form("%s_merged",h1->GetName()) );
  hres->Reset();
  Int_t nBins = h1->GetNbinsX();
  for ( Int_t j=1; j<=nBins ;j++)
    {
      float y1 =  h1->GetBinContent(j);
      float y2 =  h2->GetBinContent(j);
      hres->SetBinContent(j,  sqrt(  y1*y1 + y2*y2) );
      hres->SetBinError(j,  0.00001);
    }
  return hres;
}

//_____________________________________________________________________________
TH1D *SetRMS(TH2D *h2, const TString s)
{
  TH1D *h = h2->ProjectionX(s.Data()); h->Reset();
  for (Int_t kx=1; kx<=h2->GetNbinsX(); kx++) {
    TH1D *hTmp = h2->ProjectionY(Form("hTmp_%d",kx), kx, kx);

    if (hTmp->Integral()<=0.)
      h->SetBinContent(kx, 0.);
    else
      h->SetBinContent(kx, hTmp->GetRMS());

      cout << h->GetBinContent(kx) << endl;
    delete hTmp; hTmp = 0;
  }

  return h;
}

// Rebin first histogram to match the second
TH1D *tools::Rebin(const TH1D *h, const TH1D* href) {

  //assert(href->GetNbinsX()<=h->GetNbinsX());
  if (!(href->GetNbinsX()<=h->GetNbinsX())) {
    cout << "Histo has less bins than ref: "
	 << h->GetNbinsX() << " vs " << href->GetNbinsX()
	 << " for " << h->GetName() << endl;
  }
  
  // First, we need to rebin inclusive jets to match b-tagged jets
  TH1D *hre = (TH1D*)href->Clone(Form("%s_rebin",h->GetName()));
  hre->Reset();

  for (int i = 1; i != h->GetNbinsX()+1; ++i) {

    double x = h->GetBinLowEdge(i);
    int j = hre->FindBin(x);
    // Check that h is fully contained within href bin
    if (h->GetBinContent(i)!=0) {

      if (!(h->GetBinLowEdge(i)>=hre->GetBinLowEdge(j) - 1e-5 &&
	    h->GetBinLowEdge(i+1)<=hre->GetBinLowEdge(j+1) + 1e-5)) {
	cerr << Form("Warning, bin edges overlapping: h=[%1.0f,%1.0f],"
		     " hre=[%1.0f,%1.0f] (%s)",
		     h->GetBinLowEdge(i), h->GetBinLowEdge(i+1),
		     hre->GetBinLowEdge(j), hre->GetBinLowEdge(j+1),
		     h->GetName()) << endl;
      }

      double y = ( hre->GetBinContent(j)*hre->GetBinWidth(j)
		   + h->GetBinContent(i)*h->GetBinWidth(i) )
	/ hre->GetBinWidth(j);
      //double ey = ( hre->GetBinError(j)*hre->GetBinWidth(j)
      //	    + h->GetBinError(i)*h->GetBinWidth(i) )
      // / hre->GetBinWidth(j);
      double ey = sqrt( pow(hre->GetBinError(j)*hre->GetBinWidth(j),2)
			+ pow(h->GetBinError(i)*h->GetBinWidth(i),2) )
	/ hre->GetBinWidth(j);
      hre->SetBinContent(j, y);
      hre->SetBinError(j, ey);
    }
  } // for i

  return hre;
} // Rebin

TH1D *LPFilter( TH1D *hin , Double_t Cend  ) {

    if (Cend==-1.0) Cend = CPF ; 
    Double_t RCns = 50. * Cend * 1.e-3 ;
    Double_t At = hin->GetBinCenter(2) - hin->GetBinCenter(1) ;
    Double_t alfa = At/(RCns + At) ;
    TH1D *hout = (TH1D *) hin->Clone(); hout->Reset();
    hout->SetTitle( Form("C=%d pF",TMath::Nint(Cend)) ) ;
    hout->SetName( Form("C=%d pF",TMath::Nint(Cend)) ) ;

    hout->SetBinContent(1, hin->GetBinContent(1) );
    for ( Int_t i=2 ; i<=hin->GetNbinsX() ; i++ ) {
      Double_t val = (1.0-alfa)*hout->GetBinContent(i-1) + alfa*hin->GetBinContent(i) ; 
      hout->SetBinContent(i,val);
    }
    
    return hout ;
    
}

void make(TDirectory & out, TObject * o) {
  TDirectory * dir;
  TH1F * th1f;
  TH1D * th1d;
  TH2F * th2f;
  TH2D * th2d;
  out.cd();
  if((dir = dynamic_cast<TDirectory*>(o)) != 0) {
    TDirectory * outDir = out.mkdir(dir->GetName(), dir->GetTitle());
    TIter next(dir->GetListOfKeys());
    TKey *key;
    while( (key = dynamic_cast<TKey*>(next())) ) {
      string className(key->GetClassName());
      string name(key->GetName());
      TObject * obj = dir->Get(name.c_str());
      if(obj == 0) {
	cerr <<"error: key " << name << " not found in directory " << dir->GetName() << endl;
	exit(-1);
      }
      make(*outDir, obj);
    }
  } else if((th1f = dynamic_cast<TH1F*>(o)) != 0) {
    TH1F *h = (TH1F*) th1f->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th1d = dynamic_cast<TH1D*>(o)) != 0) {
    TH1D *h = (TH1D*) th1d->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2f = dynamic_cast<TH2F*>(o)) != 0) {
    TH2F *h = (TH2F*) th2f->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2d = dynamic_cast<TH2D*>(o)) != 0) {
    TH2D *h = (TH2D*) th2d->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  }
}

void TrkClosure(
    int anaMode=0, // 0=pt, 1=eta, 2=dphi, 3=xi
    TString outdir="fig/trkcorrv14_closure"
)
{
    TH1::SetDefaultSumw2();
    gSystem->mkdir(outdir,kTRUE);

    bool doCorr=true;
    float maxEta=2.4;
    float ptmin=0.5,ptmax=200;
    TString mod="Forest2_MergedGeneral";
    TrackingCorrections trkCorr("Forest2STAv14",mod);
    trkCorr.AddSample("trkcorr/IterTrkCorrv14XSec/IterTrkCorrv14XSec_hy18dj80to100_akPu3Calo_100_-1_-1000_genJetMode0.root",80);
    trkCorr.Init();

//   TFile * inf = new TFile("output_hy18dj100_fv27_ntv1_corrv14_akPu3PF_100_-1_-1000_saveTrk1_jmin100_tmin1_genJetMode0.root");
    TFile * inf = new TFile("output_hy18dj80_fv27_ntv1_corrv14_akPu3PF_100_-1_-1000_saveTrk1_jmin100_tmin1_genJetMode0.root");
//    TFile * inf = new TFile("output_sigdj80_fv27_ntv1_akPu3PF_100_-1_-1000_saveTrk1_jmin100_tmin1_genJetMode0_corrhi.root");
//    TFile * inf = new TFile("output_sigdj80_fv27_ntv1_akPu3PF_100_-1_-1000_saveTrk1_jmin100_tmin1_genJetMode0_corrppv14jetfinexsec.root");
//     TFile * inf = new TFile("output_sigdj80_fv27_ntv1_corrv14_akPu3PF_100_-1_-1000_saveTrk1_jmin100_tmin1_genJetMode0.root");
    TTree * t = (TTree*)inf->Get("tgj");

    //////////////////////////////////////////
    // Analysis Setup
    //////////////////////////////////////////
    TString tag=Form("trkClos%d_0to12_jpt100_pthat80",anaMode);
    TCut sel = "cBin<12&&pt1>100";
//     TString tag=Form("trkClos%d_pp_jpt100_pthat80",anaMode);
//     TCut sel = "pt1>100";
    TCut genpSel = "";
    TCut trkSel = "";
//     TCut genpSel = "simAsso==1||simAsso==2";
//     TCut trkSel = "trkAsso==1||trkAsso==2";
    TString trkWt = "trkWt";

    float nEvt = t->GetEntries(sel);
    cout << sel << ": " << nEvt << " events" << endl;

    Compare cmp("cmp","");
    TH1D * hGenp, *hTrk, *hTrkCorr;
    float xmin=1, xmax=119.9;
    bool doLogx=true, doLogy=true;
    TString genVar, trkVar;
    TCut finalGenSel,finalTrkSel;
    if (anaMode==0) {
        TH1D * hPtBin = (TH1D*)trkCorr.ptBin_->Clone("hPtBin");
        hPtBin->Reset();
        hPtBin->SetAxisRange(0.5,179.9,"X");
        hGenp = (TH1D*)trkCorr.ptBin_->Clone("hGenp");
        hTrk = (TH1D*)trkCorr.ptBin_->Clone("hTrk");
        hTrkCorr = (TH1D*)trkCorr.ptBin_->Clone("hTrkCorr");
        cmp.Legend(0.24,0.23,0.7,0.4);
        cmp.leg->AddEntry(hTrk,Form("|#eta| < %.1f",maxEta),"");
        tag += Form("_eta%.0f",maxEta*10);
        genVar="simPt";
        trkVar="trkPt";
        finalGenSel = sel&&genpSel&&Form("abs(simEta)<%.2f",maxEta);
        finalTrkSel = sel&&trkSel&&Form("abs(trkEta)<%.2f",maxEta);
    } else if (anaMode==1) {
        TH1D * hEtaBin = (TH1D*)trkCorr.etaBin_->Clone("hEtaBin");
        hEtaBin->Reset();
        hGenp = (TH1D*)hEtaBin->Clone("hGenp");
        hTrk = (TH1D*)hEtaBin->Clone("hTrk");
        hTrkCorr = (TH1D*)hEtaBin->Clone("hTrkCorr");
        doLogx=false;
        doLogy=false;
        cmp.Legend(0.19,0.76,0.65,0.93);
        cmp.leg->AddEntry(hTrk,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax),"");
        tag += Form("_pt%.0f",ptmin);
        genVar="simEta";
        trkVar="trkEta";
        finalGenSel = sel&&genpSel&&Form("simPt>=%.2f&&simPt<%.2f",ptmin,ptmax);
        finalTrkSel = sel&&trkSel&&Form("trkPt>=%.2f&&trkPt<%.2f",ptmin,ptmax);
    } else if (anaMode==2) {
        TH1D * hBin = new TH1D("hBin",";#Delta#phi(trk,jet1)",25,0,3.14159);
        hGenp = (TH1D*)hBin->Clone("hGenp");
        hTrk = (TH1D*)hBin->Clone("hTrk");
        hTrkCorr = (TH1D*)hBin->Clone("hTrkCorr");
        doLogx=false;
        doLogy=true;
        cmp.Legend(0.19,0.76,0.65,0.93);
        cmp.leg->AddEntry(hTrk,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax),"");
        tag += Form("_pt%.0f",ptmin);
        genVar="acos(cos(simPhi-phi1))";
        trkVar="acos(cos(trkPhi-phi1))";
        finalGenSel = sel&&genpSel&&Form("simPt>=%.2f&&simPt<%.2f",ptmin,ptmax);
        finalTrkSel = sel&&trkSel&&Form("trkPt>=%.2f&&trkPt<%.2f",ptmin,ptmax);
    } else if (anaMode==3) {
        TH1D * hBin = new TH1D("hBin",";#xi=1/ln(z)",12,0,6);
        hGenp = (TH1D*)hBin->Clone("hGenp");
        hTrk = (TH1D*)hBin->Clone("hTrk");
        hTrkCorr = (TH1D*)hBin->Clone("hTrkCorr");
        doLogx=false;
        doLogy=false;
        cmp.Legend(0.19,0.76,0.65,0.93);
        cmp.leg->AddEntry(hTrk,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax),"");
        tag += Form("_pt%.0f",ptmin);
        genVar="log(pt1/simPt)";
        trkVar="log(pt1/trkPt)";
//.........这里部分代码省略.........

void drawSys_merged() {
  //  TH1::SetDefaultSumw2();
    gStyle -> SetOptStat(0);
    const int kPPcentral = 1;
    const int kPPperipheral =2 ;
    const int kHIcentral =  3;
    const int kHIperipheral = 4;
    const int kPADATA = 5;
    const int kPAMC = 6;

  TH1D* dNdJetPt[7][5]; // [collision] [ ptbin]  [Before/After variation]
  TH1D* dNdJetPtSys[7][5]; // [collision] [ ptbin]  [Before/After variation]
  TH1D* dNdJetPtSysIaa[7][5]; // [collision] [ ptbin]  [Before/After variation]
  TH1D* meanJetPt[7];   // [collision] [Before/After variation]
  TH1D* meanJetPtSys[7];   // [collision] [Before/After variation]
  TH1D* meanRjg[7];
  TH1D* meanRjgSys[7];
  
  TH1D* bmeanJetPt[7];   // [collision] [Before/After variation]
  TH1D* bmeanJetPtSys[7];   // [collision] [Before/After variation]
  TH1D* bmeanRjg[7];
  TH1D* bmeanRjgSys[7];
  
  TFile * f = new TFile("../resultHistograms_ppSmeared10030.root");
  bmeanJetPt[1]  =(TH1D*)f->Get(Form("meanJetPt_pp"));
  bmeanRjg[1]  =(TH1D*)f->Get(Form("meanRjg_pp"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[1][ipt] = (TH1D*)f->Get(Form("dNdJetPt_pp_ptBin%d", ipt ));
  }  

  f   = new TFile("../resultHistograms_ppSmeared13099.root");
  bmeanJetPt[2]  =(TH1D*)f->Get(Form("meanJetPt_pp"));
  bmeanRjg[2]  =(TH1D*)f->Get(Form("meanRjg_pp"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[2][ipt] = (TH1D*)f->Get(Form("dNdJetPt_pp_ptBin%d", ipt ));
  }  
  
  // pbpb
  f =   new TFile("../resultHistograms_jetResCorrection.root");
  bmeanJetPt[3]  =(TH1D*)f->Get(Form("meanJetPt_pbpb_centralityBin1"));
  bmeanRjg[3]  =(TH1D*)f->Get(Form("meanRjg_pbpb_centralityBin1"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[3][ipt] = (TH1D*)f->Get(Form("dNdJetPt_pbpb_centralityBin1_ptBin%d", ipt ));
  }
  bmeanJetPt[4]  =(TH1D*)f->Get(Form("meanJetPt_pbpb_centralityBin2"));
  bmeanRjg[4]  =(TH1D*)f->Get(Form("meanRjg_pbpb_centralityBin2"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[4][ipt] = (TH1D*)f->Get(Form("dNdJetPt_pbpb_centralityBin2_ptBin%d", ipt ));
  }
  // pPb
  f =   new TFile("../resultHistograms_jetResCorrection.root");
  bmeanJetPt[5]  =(TH1D*)f->Get(Form("meanJetPt_ppb"));
  bmeanRjg[5]  =(TH1D*)f->Get(Form("meanRjg_ppb"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[5][ipt] = (TH1D*)f->Get(Form("dNdJetPt_ppb_ptBin%d", ipt ));
  }
  
  f =   new TFile("../resultHistograms_MC_recoIsoPhoton.root");
  bmeanJetPt[6]  =(TH1D*)f->Get(Form("meanJetPt_ppb"));
  bmeanRjg[6]  =(TH1D*)f->Get(Form("meanRjg_ppb"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPt[6][ipt] = (TH1D*)f->Get(Form("dNdJetPt_ppb_ptBin%d", ipt ));
  }  
  
  
  TFile* fSys = new TFile("../merged/relativeSys_merged.root");
  
  bmeanJetPtSys[1]  =(TH1D*)fSys->Get("meanJetPt_pp_uncertainty_merged");
  bmeanRjgSys[1]  =(TH1D*)fSys->Get(Form("meanRjg_pp_uncertainty_merged"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPtSysIaa[1][ipt] = (TH1D*)fSys->Get(Form("dNdJetPt_IaaBin_pp_ptBin%d_uncertainty_merged",  ipt ));
  }
  bmeanJetPtSys[2] = (TH1D*)bmeanJetPtSys[1]->Clone(Form("%s_2",bmeanJetPtSys[1]->GetName()));
  bmeanRjgSys[2] = (TH1D*)bmeanRjgSys[1]->Clone(Form("%s_2",bmeanRjgSys[1]->GetName()));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPtSysIaa[2][ipt] = (TH1D*)dNdJetPtSysIaa[1][ipt]->Clone(Form("%s_2",dNdJetPtSysIaa[1][ipt]->GetName()));
  }
  bmeanJetPtSys[3]  =(TH1D*)fSys->Get("meanJetPt_pbpb_centralityBin1_uncertainty_merged");
  bmeanRjgSys[3]  =(TH1D*)fSys->Get(Form("meanRjg_pbpb_centralityBin1_uncertainty_merged"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPtSysIaa[3][ipt] = (TH1D*)fSys->Get(Form("dNdJetPt_IaaBin_pbpb_centralityBin1_ptBin%d_uncertainty_merged",  ipt ));
  }
  bmeanJetPtSys[4]  =(TH1D*)fSys->Get("meanJetPt_pbpb_centralityBin2_uncertainty_merged");
  bmeanRjgSys[4]  =(TH1D*)fSys->Get(Form("meanRjg_pbpb_centralityBin2_uncertainty_merged"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPtSysIaa[4][ipt] = (TH1D*)fSys->Get(Form("dNdJetPt_IaaBin_pbpb_centralityBin2_ptBin%d_uncertainty_merged",  ipt ));
  }
  bmeanJetPtSys[5]  =(TH1D*)fSys->Get("meanJetPt_ppb_uncertainty_merged");
  bmeanRjgSys[5]  =(TH1D*)fSys->Get(Form("meanRjg_ppb_uncertainty_merged"));
  for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
    dNdJetPtSysIaa[5][ipt] = (TH1D*)fSys->Get(Form("dNdJetPt_IaaBin_ppb_ptBin%d_uncertainty_merged",  ipt ));
  }


  
  // transplant the uncertainty  
  for ( int coll = 1 ; coll<=5 ; coll++) {
    for ( int ipt = 1 ; ipt <=nPtBin ; ipt++) {
      dNdJetPtSys[coll][ipt] = (TH1D*) dNdJetPt[coll][ipt]->Clone(Form("dNdJetPtSys_coll%d_ipt%d",coll,ipt));
      dNdJetPtSys[coll][ipt]->Reset();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    std::cout<<"\nScale Histograms..."<<std::endl;
    th1Top->Scale(ttbarA_scale);
    th1stops->Scale(STopS_scale);
    th1stopt->Scale(STopT_scale);
    th1stoptw->Scale(STopTW_scale);
    th1stopps->Scale(SToppS_scale);
    th1stoppt->Scale(SToppT_scale);
    th1stopptw->Scale(SToppTW_scale);
    th1wajets->Scale(WAJets_scale);
    th1wajetsPT100->Scale(WAJetsPT100_scale);
    th1wwa->Scale(WWA_scale);
    th1wwa2->Scale(WWA2_scale);
    th1wza->Scale(WZA_scale);
    th1zz->Scale(ZZ_scale);
    th1zajets->Scale(ZAJets_scale);


    ///////////////////////////
    // Combine certain samples:
    th1stop->Add(th1stopptw,1);
    th1stop->Add(th1stoppt,1);
    th1stop->Add(th1stopps,1);
    th1stop->Add(th1stoptw,1);
    th1stop->Add(th1stopt,1);
    th1stop->Add(th1stops,1);

    th1wwa->Add(th1wwa2,1);
    th1wajets->Add(th1wajetsPT100,1);


    ///////////////////////////////////////////////////////////////
    // Combine all background + SM WWA samples for QCD Estimations:
    TH1D *background = (TH1D*)th1wajets->Clone("background");
    background->Reset();
    background->Add(th1fkdata,1);
    background->Add(th1Top,1);
    background->Add(th1stop,1);
    background->Add(th1wajets,1);
    background->Add(th1zajets,1);
    background->Add(th1zz,1);
    background->Add(th1wwa,1);
    background->Add(th1wza,1);

    float qcd_scale=0;
    if(th1qcd->Integral()>0 ) {
        qcd_scale = 0.0637*background->Integral()/th1qcd->Integral();
        th1qcd->Scale(qcd_scale);
        background->Add(th1qcd,1);
    }

    /*
    std::cout << "\nSample Contribution:" << std::endl;
    std::cout << "WWA: " << th1wwa->Integral() << std::endl;
    std::cout << "Top: " << th1Top->Integral() << std::endl;
    std::cout << "Single Top: " << th1stop->Integral() << std::endl;
    std::cout << "WA+Jets: " << th1wajets->Integral() << std::endl;
    std::cout << "ZA+Jets: " << th1zajets->Integral() << std::endl;
    std::cout << "WZA: " << th1wza->Integral() << std::endl;
    std::cout << "ZZ: " << th1zz->Integral() << std::endl;
    std::cout << "Fake Photons: " << th1fkdata->Integral() << std::endl;
    std::cout << "Data: " << data_obs->Integral() << std::endl;
    std::cout << "QCD: " << th1qcd->Integral() << std::endl;

    std::cout << "\nSample Contribution:" << std::endl;
    std::cout << "Top: " << th1Top->GetBinError(1) << std::endl;
    std::cout << "WA+Jets: " << th1wajets->GetBinError(1) << std::endl;

void toyMC_bkgTemp(int runIeta=-1, int runIpt=-1) {

    cout << "ROOT version = " << gROOT->GetVersion() << endl;
    //   gSystem->mkdir("toysPlot");
    char tmp[1000];


    TH1D* htoyResult_pull[nEtaBin][nPtBin];
    TH1D* htoyResult_bias[nEtaBin][nPtBin];

    TFile *fsumFile = new TFile("/afs/cern.ch/user/s/syu/scratch0/LxplusArea/proj_comb_comb3Iso_template.root");
    TFile* finFile = new TFile("/afs/cern.ch/user/s/syu/scratch0/LxplusArea/template_comb3Iso_template.root");
    TFile* zeeFile = new TFile("/afs/cern.ch/user/s/syu/scratch0/LxplusArea/anadipho_Zee_Vg_3pb.root");
    TH1D* hTemplate = (TH1D*)finFile->FindObjectAny("h_EB_comb3Iso_EGdata_pt21");
    hTemplate->Reset();


    char* dec[2] = {"EB","EE"};
    for(int ieta=0; ieta<nEtaBin; ieta++) {
        for(int ipt=0; ipt < nPtBin; ipt++) {

            htoyResult_pull[ieta][ipt] = new TH1D(Form("hpull_%s_pt_%d",
                                                  dec[ieta],
                                                  (int)fBinsPt[ipt]),
                                                  "",50,-5.0,5.0);


            htoyResult_bias[ieta][ipt] = new TH1D(Form("hbias_%s_pt_%d",
                                                  dec[ieta],
                                                  (int)fBinsPt[ipt]),
                                                  "",100,-0.5,0.5);


        }
    }



    TH1D* hfit_sig;
    TH1D* hfit_bkg;

    TH1D* hTemplate_S[nEtaBin][nPtBin];
    TH1D* hTemplate_B[nEtaBin][nPtBin];
    TH1D* hZeeTemplate_S[nEtaBin];
    TH1D* hdata_data[nEtaBin][nPtBin];
    TH1D* htemp;
    TH1D* hTemplate_ShiftB[nEtaBin][nPtBin];

    for(int ieta=0; ieta< nEtaBin; ieta++) {

        // getting a different signal template
        if(ieta==0) {
            sprintf(tmp,"h_%s_combIso",dec[ieta]);
            cout << "looking for histogram " << tmp << " in file " <<
                 zeeFile->GetName() << endl;
            hZeeTemplate_S[ieta]= (TH1D*)zeeFile->FindObjectAny(tmp);
            hZeeTemplate_S[ieta]->Rebin(REBINNINGS_TEMP);
        }
        else {
            sprintf(tmp,"h_%s_comb3Iso_sig_sum_SIG",dec[ieta]); //no pt dep.
            cout << "looking for histogram " << tmp << " in file " <<
                 fsumFile->GetName() << endl;
            hZeeTemplate_S[ieta]= (TH1D*)fsumFile->FindObjectAny(tmp);
            hZeeTemplate_S[ieta]->Rebin(REBINNINGS_TEMP);
        }

        for(int ipt=0; ipt < nPtBin; ipt++) {

            if(runIeta>=0 && ieta!=runIeta)continue;
            if(runIpt>=0 && ipt!=runIpt)continue;

            // getting histograms from data root file
            sprintf(tmp,"h_%s_comb3Iso_EGdata_pt%d",dec[ieta],(int)fBinsPt[ipt]);
            cout << "looking for histogram " << tmp << " in file " <<
                 finFile->GetName() << endl;
            hdata_data[ieta][ipt] = (TH1D*)finFile->FindObjectAny(tmp);
            hdata_data[ieta][ipt]->Rebin(REBINNINGS_DATA);

            // filling unbinned data
            htemp = (TH1D*)hdata_data[ieta][ipt]->Clone("htemp");
            htemp->Reset();


            sprintf(tmp,"h_%s_comb3Iso_sig_pt%d",dec[ieta],(int)fBinsPt[ipt]);
            cout << "looking for histogram " << tmp << " in file " <<
                 finFile->GetName() << endl;
            hTemplate_S[ieta][ipt] = (TH1D*)finFile->FindObjectAny(tmp);
            hTemplate_S[ieta][ipt]->Rebin(REBINNINGS_TEMP);

            if(ieta==0 && fBinsPt[ipt]>=50)
                sprintf(tmp,"h_%s_comb3Iso_bkg_pt%d",dec[ieta],50);
            else if(ieta==0 )
                sprintf(tmp,"h_%s_comb3Iso_bkg_pt%d",dec[ieta],(int)fBinsPt[ipt]);
            else if(ieta==1 && fBinsPt[ipt]>=60)
                sprintf(tmp,"h_%s_comb3IsoSB_EGdata_pt%d",dec[ieta],60);
            else if(ieta==1)
                sprintf(tmp,"h_%s_comb3IsoSB_EGdata_pt%d",dec[ieta],(int)fBinsPt[ipt]);
            cout << "looking for histogram " << tmp << " in file " <<
                 finFile->GetName() << endl;
            hTemplate_B[ieta][ipt] = (TH1D*)finFile->FindObjectAny(tmp);
//.........这里部分代码省略.........