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

TH1D* GetFinalCorrFunc1D(int itrg, int jass)
{
                    TH1D* hbackground_phi = (TH1D*)GetRawBackground1D(itrg,jass);
		    TH1D* hcorrphi = (TH1D*)GetJetCorrFunc1D_ZYAM(itrg,jass);
                    hcorrphi->Multiply(hbackground_phi);
                    float max = hcorrphi->GetBinContent(hcorrphi->GetMaximumBin());
		    hcorrphi->SetAxisRange(ymin,max*1.5,"Y");
                    hcorrphi->SetYTitle("#frac{1}{N^{trg}}#frac{d^{2}N}{d#Delta#phi}");
                    return hcorrphi;
}

//.........这里部分代码省略.........
  divideBinWidth(hPtMCrecoonly);
  //ntGen->Project("hPtGen","hiBin",(TCut(selmcgen.Data())));
  ntGen->Project("hPtGen","hiBin",TCut(weighpthat)*TCut(weightGpt)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);
  //ntGen->Project("hPtGenAcc","hiBin",(TCut(selmcgenacceptance.Data())));
  ntGen->Project("hPtGenAcc","hiBin",TCut(weighpthat)*TCut(weightGpt)*(TCut(selmcgenacceptance.Data())));
  divideBinWidth(hPtGenAcc);
  //ntGen->Project("hPtGenAccWeighted","hiBin",TCut(weightfunctiongen)*(TCut(selmcgenacceptance.Data())));
  ntGen->Project("hPtGenAccWeighted","hiBin",TCut(weighpthat)*TCut(weightGpt)*TCut(weightHiBin)*(TCut(selmcgenacceptance.Data())));
  divideBinWidth(hPtGenAccWeighted);

  ntMC->Project("hPthat","pthat","1");
  ntMC->Project("hPthatweight","pthat",TCut("1"));

  hPtMC->Sumw2();
  hPtGenAcc->Sumw2();
  hPtMCrecoonly->Sumw2();
  //Acceptance
  TH1D* hEffAcc = (TH1D*)hPtGenAcc->Clone("hEffAcc");
  hEffAcc->Sumw2();
  hEffAcc->Divide(hEffAcc,hPtGen,1,1,"b");
  //Selection
  TH1D* hEffSelection = (TH1D*)hPtMC->Clone("hEffSelection");
  hEffSelection->Sumw2();
  hEffSelection->Divide(hEffSelection,hPtGenAccWeighted,1,1,"b");
  //Acc * Eff (one shot)
  TH1D* hEffReco = (TH1D*)hPtMCrecoonly->Clone("hEffReco");
  hEffReco->Sumw2();
  hEffReco->Divide(hEffReco,hPtGen,1,1,"b");
  //Acc * Eff
  TH1D* hEff = (TH1D*)hEffSelection->Clone("hEff");
  hEff->Sumw2();
  //hEff->Divide(hPtMC,hPtGen,1,1,"");
  hEff->Multiply(hEff,hEffAcc,1,1);

  TH2F* hemptyEff=new TH2F("hemptyEff","",50,_ptBins[0]-5.,_ptBins[_nBins]+5.,10.,0,1.0);  
  hemptyEff->GetXaxis()->CenterTitle();
  hemptyEff->GetYaxis()->CenterTitle();
  //hemptyEff->GetYaxis()->SetTitle("acceptance x #epsilon_{reco} x #epsilon_{sel} ");
  hemptyEff->GetYaxis()->SetTitle("#alpha x #epsilon");
  hemptyEff->GetXaxis()->SetTitle("Centrality");
  hemptyEff->GetXaxis()->SetTitleOffset(0.9);
  hemptyEff->GetYaxis()->SetTitleOffset(0.95);
  hemptyEff->GetXaxis()->SetTitleSize(0.05);
  hemptyEff->GetYaxis()->SetTitleSize(0.05);
  hemptyEff->GetXaxis()->SetTitleFont(42);
  hemptyEff->GetYaxis()->SetTitleFont(42);
  hemptyEff->GetXaxis()->SetLabelFont(42);
  hemptyEff->GetYaxis()->SetLabelFont(42);
  hemptyEff->GetXaxis()->SetLabelSize(0.035);
  hemptyEff->GetYaxis()->SetLabelSize(0.035);  
  hemptyEff->SetMaximum(2);
  hemptyEff->SetMinimum(0.);
  hemptyEff->Draw();

  TH2F* hemptyEffAcc=(TH2F*)hemptyEff->Clone("hemptyEffAcc");
  TH2F* hemptyEffReco=(TH2F*)hemptyEff->Clone("hemptyEffReco");
  TH2F* hemptyEffSelection=(TH2F*)hemptyEff->Clone("hemptyEffSelection");
 

  TCanvas*canvasEff=new TCanvas("canvasEff","canvasEff",1000.,500);
  canvasEff->Divide(2,1);
  canvasEff->cd(1);
  
  hemptyEffAcc->SetYTitle("#alpha");
  hemptyEffAcc->Draw();

//.........这里部分代码省略.........
	//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){
        CRyieldEM = (TH1D*) hDataEM->Clone("h_mt2binsCRyieldEM");
        CRyieldEMCard = (TH1D*) CRyieldEM->Clone("CRyieldEMCard");
    }

    if (  doHybridSimple || (doHybridInclusiveTemplate && h_MT2Template==0) ) { 
      // purity needs to describe the integrated purity of the CR
      // ratio needs to be modified so that the last N bins include kMT2
      // CRyield needs to be modified so that the last N bins have the same yield (which is the integral over those N bins)
      if (verbose) cout<<" Implementing simple hybrid "<<endl;

      for ( int ibin=1; ibin <= hZinv->GetNbinsX(); ++ibin ) {

	if (ibin < lastbin_hybrid) continue;
	

int drawProgressivePtCut(std::string fileName) {

  // Create file on which histogram(s) can be saved.
  TFile *inFile = new TFile(fileName.c_str(), "READ");
  TH1D* ptProfile = ((TProfile*) inFile->Get( "pT bins" ))->ProjectionX(""); 
  TH1D* hcalProfile = ((TProfile*) inFile->Get( "hcal bins" ))->ProjectionX("");
  TH1D* gev3Profile = ((TProfile*) inFile->Get( "3gev bins" ))->ProjectionX("");
  TH1D* combProfile = (TH1D*) hcalProfile->Clone("hnew");
  combProfile->Multiply( gev3Profile );
  TH1D* hiEProfile = ((TProfile*) inFile->Get( "hie bins" ))->ProjectionX("");

  TCanvas *canv = new TCanvas("c1","c1",1200,1200);
  setTDRStyle();
  canv->UseCurrentStyle();
  ptProfile->GetYaxis()->SetTitleOffset(1.6);
  canv->SetSelected(canv);
  canv->cd();
  // Show histogram  
  TPad *pad1 = new TPad("pad1","",0,0,1,1);
  TPad *pad2 = new TPad("pad2","",0,0,1,1);
  TPad *pad3 = new TPad("pad3","",0,0,1,1);
  TPad *pad4 = new TPad("pad4","",0,0,1,1);
  TPad *pad5 = new TPad("pad5","",0,0,1,1);

  pad2->SetFillStyle(4000);
  pad3->SetFillStyle(4000);
  pad4->SetFillStyle(4000);
  pad5->SetFillStyle(4000);
  pad1->Draw();
  pad1->cd();
  ptProfile->GetXaxis()->SetTitle("p_{T} (GeV)");
  ptProfile->GetYaxis()->SetTitle("simulation for detected pt/full pt");
  ptProfile->GetYaxis()->SetTitleOffset(1.6);
  ptProfile->SetStats(0);
  ptProfile->GetXaxis()->SetMoreLogLabels();
  ptProfile->GetXaxis()->SetNoExponent();
  ptProfile->Draw();
  TLegend *leg = tdrLeg(0.4,0.73,0.6,0.93);
  leg->AddEntry( ptProfile, "All effects" );
  leg->AddEntry( combProfile, "HCal+3GeV" );
  leg->AddEntry( hcalProfile, "HCal" );
  leg->AddEntry( gev3Profile, "3GeV" );
  leg->AddEntry( hiEProfile, "High pt ineff." );
  leg->Draw();
  pythiaFinal();
  canv->Modified();
  pad1->UseCurrentStyle();
  ptProfile->SetMarkerStyle(kCircle);
  ptProfile->SetMaximum(1.04);
  ptProfile->SetMinimum(0.89);
  ptProfile->SetMarkerColor(7);
  pad1->Update();
  
  canv->cd();
  pad2->Draw();
  pad2->cd();
  combProfile->SetStats(0);
  combProfile->GetXaxis()->SetMoreLogLabels();
  combProfile->GetXaxis()->SetNoExponent();
  combProfile->Draw();
  pad2->UseCurrentStyle();
  combProfile->SetMaximum(1.04);
  combProfile->SetMinimum(0.89);
  combProfile->SetMarkerColor(3);
  pad2->Update();

  canv->cd();
  pad3->Draw();
  pad3->cd();
  hcalProfile->SetStats(0);
  hcalProfile->GetXaxis()->SetMoreLogLabels();
  hcalProfile->GetXaxis()->SetNoExponent();
  hcalProfile->Draw();
  pad3->UseCurrentStyle();
  hcalProfile->SetMaximum(1.04);
  hcalProfile->SetMinimum(0.89);
  hcalProfile->SetMarkerStyle(kStar);
  hcalProfile->SetMarkerColor(5);
  pad3->Update();
  
  canv->cd();
  pad4->Draw();
  pad4->cd();
  gev3Profile->SetStats(0);
  gev3Profile->GetXaxis()->SetMoreLogLabels();
  gev3Profile->GetXaxis()->SetNoExponent();
  gev3Profile->Draw();
  pad4->UseCurrentStyle();
  gev3Profile->SetMaximum(1.04);
  gev3Profile->SetMinimum(0.89);
  gev3Profile->SetMarkerStyle(kMultiply);
  pad4->Update();

  canv->cd();
  pad5->Draw();
  pad5->cd();
  hiEProfile->SetStats(0);
  hiEProfile->GetXaxis()->SetMoreLogLabels();
  hiEProfile->GetXaxis()->SetNoExponent();
  hiEProfile->Draw();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      auto hist_MC_SR = (TH1D*) fmc->Get(srname + "/" + hnameSR);

      if (!hist_MC_SR) {
        if (!hnameSR.Contains("cr2lTriggerSF"))
          cout << "Couldn't find yield hist for " << (srname + "/" + hnameSR) << " in " << fmc->GetName() << ". Use centralHist!" << endl;
        hist_MC_SR = (TH1D*) fmc->Get(hname_MC_SR)->Clone(hnameSR);
      }

      auto alphaHist = (TH1D*) hist_MC_SR->Clone(hname+"_alpha");
      if (useMetExtrapolation && extr_start_bin != lastbin) {
        // To take the MET distribution from the CR
        double cerr_SR = 0;
        double cyld_SR = alphaHist->IntegralAndError(extr_start_bin, -1, cerr_SR);
        double cyld_CR = hist_MC_CR->Integral(extr_start_bin, -1);
        for (int ibin = extr_start_bin; ibin <= lastbin; ++ibin) {
          double metfrac = hist_MC_CR->GetBinContent(ibin) / cyld_CR;
          alphaHist->SetBinContent(ibin, metfrac * cyld_SR);
          alphaHist->SetBinError(ibin, metfrac * cerr_SR);
        }
        combineYieldsInExtrBins(hist_MC_CR);
      }
      alphaHist->Divide(hist_MC_CR);

      for (int i = 1; i <= alphaHist->GetNbinsX(); ++i) {
        // zero out negative yields
        if (alphaHist->GetBinContent(i) < 0) {
          alphaHist->SetBinContent(i, 0);
          alphaHist->SetBinError(i, 0);
        }
      }

      outdir->cd();
      TH1D* hout = (TH1D*) alphaHist->Clone(hname);
      hout->Multiply(hist_data_CR);
      hout->Write();

      if (yearSeparateSyst && (hname.EndsWith("Up") || hname.EndsWith("Dn"))) {
        for (int i = 1; i < 4; ++i) {

          auto hcen_MC_CR = (TH1D*) fbkgs[i]->Get(crname+"/h_metbins");
          auto hcen_MC_SR = (TH1D*) fbkgs[i]->Get(srname+"/h_metbins"+gentype);
          auto hsys_MC_CR = (TH1D*) fbkgs[i]->Get(crname+"/"+hname);
          auto hsys_MC_SR = (TH1D*) fbkgs[i]->Get(srname+"/"+hnameSR);

          if (!hist_MC_SR) {
            if (!hnameSR.Contains("cr2lTriggerSF"))
              cout << "Couldn't find yield hist for " << (srname + "/" + hnameSR) << " in " << fmc->GetName() << ". Use centralHist!" << endl;
            hist_MC_SR = (TH1D*) fmc->Get(hname_MC_SR)->Clone(hnameSR);
          }

          auto alphaHist_yi = (TH1D*) fmc->Get(hname_MC_SR)->Clone(TString(hname).Insert(hname.Length()-2, Form("%d", 15+i)));
          auto h_MC_CR_yi = (TH1D*) fmc->Get(hname_MC_CR)->Clone(Form("%s_den_%d", hname.Data(), 15+i));
          if (hcen_MC_SR) alphaHist_yi->Add(hcen_MC_SR, -1);
          if (hsys_MC_SR) alphaHist_yi->Add(hsys_MC_SR);
          if (hcen_MC_CR) h_MC_CR_yi->Add(hcen_MC_CR, -1);
          if (hsys_MC_CR) h_MC_CR_yi->Add(hsys_MC_CR);

          if (useMetExtrapolation && extr_start_bin != lastbin) {
            // To take the MET distribution from the CR
            double cerr_SR = 0;
            double cyld_SR = alphaHist_yi->IntegralAndError(extr_start_bin, -1, cerr_SR);
            double cyld_CR = h_MC_CR_yi->Integral(extr_start_bin, -1);
            for (int ibin = extr_start_bin; ibin <= lastbin; ++ibin) {
              double metfrac = h_MC_CR_yi->GetBinContent(ibin) / cyld_CR;
              alphaHist_yi->SetBinContent(ibin, metfrac * cyld_SR);
              alphaHist_yi->SetBinError(ibin, metfrac * cerr_SR);

TH1D *H1DConvolution( TH1D *htf , TH1D *htct , Double_t Cend , int tid) { 
      
   //------------>Both input histograms should have the same bin width<-----------------
   
   //Here you can apply an extra LPFiltering
   //if (Cend!=0) htct = LPFilter( htct , Cend ); 

   //Convolute (commutative)
   //C(t) = Int[ tct(x) transferfunction(t-x) dx ]
   Double_t bw = htct->GetBinCenter(2) - htct->GetBinCenter(1);
   //Double_t bw = htct->GetXaxis()->GetBinCenter(2) - htct->GetXaxis()->GetBinCenter(1);
   //TF1 *f1 = new TF1("f1","abs(sin(x)/x)*sqrt(x)",0,10);
   //   float r = f1->GetRandom();
    TString tftit, tfname;
    tftit.Form("hConv_%d_%d", tid, count);
    tfname.Form("conv_%d_%d", tid, count);
   TH1D *hConv = new TH1D(tftit,tfname,2*htct->GetNbinsX(),-htct->GetNbinsX()*bw,htct->GetNbinsX()*bw);
   
   
   //The convoluted response to the TCT signal is going to be another histogram sized similar to htct
   Int_t Ntf = htf->GetNbinsX() , Ntct = htct->GetNbinsX();
   
   //Create the reverse histogram of the transfer function
   TH1D *hinv = (TH1D *) htf->Clone(); hinv->Reset();
   for (Int_t j=1; j<= Ntf ; j++) hinv->SetBinContent( j , htf->GetBinContent(Ntf-j+1) );  
  // std::basic_string<char> t= std::to_string(tid);
   TH1D *hg = (TH1D *) htct->Clone(); 
   
   hinv->Draw();
   #if CTRLPLOT == 1
   gPad->Print( "hgcontrol.pdf[" ) ; gPad->Print( "hgcontrol.pdf" );
   #endif
   
   for ( Int_t i=1 ; i<=2*Ntct ; i++ ) { 
          
     //Create a shifted histogram version of the inverse
     hg->Reset();
     //for (Int_t j=TMath::Nint(-0.5*Ntf); j<= TMath::Nint(0.5*Ntf) ; j++) hg->SetBinContent( i-j , hinv->GetBinContent( j+TMath::Nint(0.5*Ntf) ) );  
     if ( i<=Ntf ) {
       //Histogram is shifting in from the left
       for (Int_t j=1; j<=i ; j++) hg->SetBinContent( j , hinv->GetBinContent( Ntf-i+j ) );
     } else {
       //Histogram is shifting out. Leaving from the right
       Int_t cont=1 ;
       for (Int_t j=i-Ntf+1; j<=2*Ntf ; j++) {
         hg->SetBinContent( j , hinv->GetBinContent( cont ) );
	 cont++;
       }
  
     }
     
     
     //Multiply f(tau)*g(t-tau)
     hg->Multiply( htct );
     
     //Double_t fxg = hg->Integral("width");
     Double_t fxg = hg->Integral();
     
     hConv->SetBinContent(i,fxg);
     
     #if CTRLPLOT==1
       THStack *hst=new THStack("hst","conv");
       hConv->SetLineColor(2);hConv->SetLineWidth(2);
       hst->Add(hg) ; hst->Add(hConv);
       hst->Draw("nostack");
       if   (i==2*Ntct) { 
         gPad->Print( "hgcontrol.pdf" )  ; gPad->Print( "hgcontrol.pdf]" ) ; 
       } else if (i%10==0)  gPad->Print( "hgcontrol.pdf" )  ;
     #endif

   }

   gStyle->SetOptStat(0);
   gStyle->SetHistLineWidth(2);
   //hConv->SetLineColor(kRed) ;htct->SetLineColor(kBlack) ;htf->SetLineColor(kBlue) ;
   
   THStack *hs = new THStack();
   //hs->Add(htf);
   hs->Add(htct);
   hs->Add(hConv);
   TCanvas *c1=new TCanvas(); c1->cd();
   hs->Draw("nostack");
   //hs->GetXaxis()->SetRangeUser(-2.,10.);
   hs->GetXaxis()->SetTitle("Time [ns]") ;

   c1->SetGrid(1);

   TLegend* legend = c1->BuildLegend();
   legend->Draw();
   c1->Update();
   //   NOPDF for the moment
//   c1->Print( "convolution.pdf" );
 //too many files...
/*
    tftit.Form("conv_%d_%d.root", tid, count);
   TFile *f=new TFile(tftit,"UPDATE");
   hConv->Write();
   //f->Close();
   delete f;

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

//.........这里部分代码省略.........
          if(!htrigeffmc && dmc->cd(dir->GetName())) {

            TDirectory *dir1 = dmc->GetDirectory(dir->GetName()); assert(dir1);
            TH1D *hpty = (TH1D*)dir1->Get("hpt"); assert(hpty);
            assert(dmc->cd("mc"));
            dmc->cd("mc");
            TDirectory *dir2 = dmc->GetDirectory("mc"); assert(dir2);
            TH1D *hptx = (TH1D*)dir2->Get(Form("hpt_%s",dir->GetName()));

            outdir->cd();
            if (hpty && hptx) htrigeffmc = (TH1D*)hpty->Clone("htrigeffmc");
            if (hpty && hptx) htrigeffmc->Divide(hpty,hptx,1,1,"B");
          }

          // Add data/MC scale factor for trigger efficiency
          if (_nh_dt && !htrigeffsf) {

            assert(dmc->cd(dir->GetName()));
            dmc->cd(dir->GetName());
            TDirectory *dirmc = dmc->GetDirectory(dir->GetName()); assert(dirmc);
            TProfile *pm = (TProfile*)dirmc->Get("ptrigefftp");
            TProfile *pd = (TProfile*)dir->Get("ptrigefftp");

            outdir->cd();
            if (pm && pd) htrigeffsf = pm->ProjectionX("htrigeffsf");
            if (pm && pd) htrigeffsf->Divide(pd,pm,1);
          }

          // Combine MC trigger efficiency and scalefactor
          if (htrigeffmc) { // not available for 'mc' directory
            outdir->cd();
            htrigeff = (TH1D*)htrigeffmc->Clone("htrigeff");
            assert(!_nh_dt || htrigeffsf);
            if (_nh_dt) htrigeff->Multiply(htrigeffsf);

            TH1D *h = (TH1D*)dir->Get("hpt");
            assert(outdir->FindObject("hpt_notrigeff")==0);
            outdir->cd();
            hpt_notrigeff = (TH1D*)h->Clone("hpt_notrigeff");
          }

          fmc->Close();
        } // dotrigeff


        // Scale data to account for time dependence
        bool dotimedep = ((string(obj2->GetName())=="hpt") || isjk || isjet);
        TH1D *htimedep = (TH1D*)outdir->FindObject("htimedep");
        TH1D *htimefit = (TH1D*)outdir->FindObject("htimefit");
        TH1D *hpt_notimedep = 0, *hpt_withtimedep = 0;
        double ktime = 1.;

        if (!htimedep) {

          TH1D *h = (TH1D*)dir->Get("hpt");
          TH1D *hsel = (TH1D*)dir->Get("hselpt");
          TH1D *hpre = (TH1D*)dir->Get("hpt_pre");
          //TH1D *hlumi0 = (TH1D*)dir->Get("../jt450/hlumi");

          // Fix luminosity for unprescaled trigger
          //string strg = dir->GetName();
          //double lum0 = triglumi["jt450"];
          //for (int i = 1; i != hlumi0->GetNbinsX()+1; ++i) {
          //hlumi0->SetBinContent(i, lum0);
          //}

//.........这里部分代码省略.........
		 int title_id=iProf12%4;
		
		 if(title_id!=2){
			 profileHist1[iProf12][0]->Rebin(5);
			 profileHist1[iProf12][1]->Rebin(5);
		 }
	
		 TH1D* baseHist = (TH1D*)profileHist1[iProf12][0]->Clone();
		 TH1D* overHist = (TH1D*)profileHist1[iProf12][1]->Clone();
		 
		 if(MakeTotal && iProf12>7){ 
			 double binLow = baseHist->GetXaxis()->GetBinLowEdge(baseHist->GetXaxis()->GetFirst());
			 double binUp = baseHist->GetXaxis()->GetBinUpEdge(baseHist->GetXaxis()->GetLast());
			 int nBins = baseHist->GetXaxis()->GetNbins();
			 
			 TH1D* newBaseHist = new TH1D(TString(baseHist->GetName())+"TOT",baseHist->GetTitle(),nBins,binLow,binUp);
			 TH1D* newDMHist = new TH1D("newDMHist","",nBins,binLow,binUp);
			 TH1D* newOverHist = new TH1D(TString(overHist->GetName())+"TOTover",baseHist->GetTitle(),nBins,binLow,binUp);
			 TH1D* newDMHist2 = new TH1D("newDMHist2","",nBins,binLow,binUp);
			 
			 
			 for(int iBin=0; iBin <= baseHist->GetNbinsX(); iBin++)
			 {
				newBaseHist->SetBinContent(iBin,baseHist->GetBinContent(iBin));
				newBaseHist->SetBinError(iBin,baseHist->GetBinError(iBin));
				newDMHist->SetBinContent(iBin,profileHist1[title_id][0]->GetBinContent(iBin));
				newDMHist->SetBinError(iBin,profileHist1[title_id][0]->GetBinError(iBin));
				newOverHist->SetBinContent(iBin,overHist->GetBinContent(iBin));
				newOverHist->SetBinError(iBin,overHist->GetBinError(iBin));
				newDMHist2->SetBinContent(iBin,profileHist1[title_id][1]->GetBinContent(iBin));
				newDMHist2->SetBinError(iBin,profileHist1[title_id][1]->GetBinError(iBin));	
				
			 }
			 newBaseHist->Multiply(newDMHist);
			 newOverHist->Multiply(newDMHist2);
			 baseHist=newBaseHist;
			 overHist=newOverHist;
			 delete newDMHist;
			 delete newDMHist2;
		 }
	
		 if(title_id==0) baseHist->GetXaxis()->SetRangeUser(0,100);
		 baseHist->GetXaxis()->SetLabelOffset(0.1);		
		 baseHist->GetXaxis()->SetTitle(XTitle[title_id]);
		 baseHist->GetYaxis()->SetTitle("Fake rate");
		 baseHist->SetMarkerStyle(20);
		 baseHist->SetMarkerSize(2);
		 baseHist->SetLineWidth(2);
		 
		 overHist->SetMarkerColor(kRed);
		 overHist->SetMarkerStyle(25);
		 overHist->SetMarkerSize(2);
		 overHist->SetLineWidth(2);
		 overHist->SetLineStyle(2);
		 overHist->SetLineColor(kRed);
		
		 int binmax = baseHist->GetMaximumBin();
	     double max = 1.2*(baseHist->GetBinContent(binmax)+baseHist->GetBinError(binmax));
	     int binmax2 = overHist->GetMaximumBin();
	     double max2 = 1.2*(overHist->GetBinContent(binmax2)+overHist->GetBinError(binmax2));
	     max = max2 > max ? max2 : max;
	 
	     baseHist->SetMaximum(10.0);
	     baseHist->SetMinimum(1e-3);
	     
		 

//.........这里部分代码省略.........
	peakPosth[iSample]->SetBinError(1, error);
	
	integral = mth[iSample]->IntegralAndError( n3, n4, error);
	tailh[iSample]->SetBinContent(1, integral);
	tailh[iSample]->SetBinError(1, error);
      }

      ///////////////////////////////////////////
      // Mt peak weight
      ///////////////////////////////////////////

      // peak SF Pre

      peakSFPreh[1]->Add(peakPreh[0]);
      peakSFPreh[1]->Add(peakPreh[3], -1.); 
      peakSFPreh[1]->Add(peakPreh[4], -1.);

      denh->Clear(); denh->Reset();
      denh->Add(peakPreh[1]);
      denh->Add(peakPreh[2]);

      peakSFPreh[1]->Divide(denh);

      peakSFPre = peakSFPreh[1]->GetBinContent(1);
      peakSFPreErr = peakSFPreh[1]->GetBinError(1);

      peakSFPreh[2]->SetBinContent( 1, peakSFPre);
      peakSFPreh[2]->SetBinError( 1, peakSFPreErr);

      // peak SF Post

      peakSFPosth[2]->Add(peakPosth[0]);
      tmph->Clear(); tmph->Reset();
      peakSFPosth[2]->Add(peakPosth[3], -1.);
      peakSFPosth[2]->Add(peakPosth[4], -1.);

      
      tmph->Add(peakPosth[1]);
      tmph->Multiply(peakSFPreh[1]);
      tmph->Add(peakPosth[2]);
      peakSFPosth[2]->Divide(tmph);
      
      peakSFPost = peakSFPosth[2]->GetBinContent(1);
      peakSFPostErr = peakSFPosth[2]->GetBinError(1);

      ///////////////////////////////////////////
      // Tail and SF
      ///////////////////////////////////////////

      tailData = tailh[0]->GetBinContent(1);
      tailDataErr = tailh[0]->GetBinError(1);
      
      tmph->Clear(); tmph->Reset();
      denh->Clear(); denh->Reset();
      denh->Add(tailh[1]);
      denh->Multiply(peakSFPreh[1]);
      tmph->Add(tailh[2]);
      tmph->Multiply(peakSFPosth[2]);
      denh->Add(tmph);
      denh->Add(tailh[3]);
      denh->Add(tailh[4]);
      tailMC = denh->GetBinContent(1);
      tailMCErr = denh->GetBinError(1);
      
      tailDiLep = tailh[2]->GetBinContent(1);
      tailDiLepErr = tailh[2]->GetBinError(1);

      SFh->Add(tailh[0]);
      SFh->Divide(denh);
      SF = SFh->GetBinContent(1);
      SFErr = SFh->GetBinError(1);

      ///////////////////////////////////////////////
      // Fill outTree
      ///////////////////////////////////////////////

      tree[ilep]->Fill();
      
      //////////////////////////////////
      // Pritn some output
      //////////////////////////////////

      cout<<ID<<"\t";
      cout<<iSR<<" "<<ilep<<"\t";
      cout<<peakSFPre<<"+-"<<peakSFPreErr<<"\t";
      cout<<peakSFPost<<"+-"<<peakSFPostErr<<"\t";
      cout<<SF<<"+-"<<SFErr<<endl;      
    }      
  }

  for (int iSample = 0; iSample < NSamples; iSample++)
    inFile[iSample]->Close();

  for (int ilep = 0; ilep < NLep; ilep++){
    outFile->cd();
    tree[ilep]->Write();
  }

  return 0;
}