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

void anaFragShape(TString infile="dj_HCPR-GoodTrkAndPixel_CleanEvt1130.root")
{
  TChain * djcalo = new TChain("djcalo/djTree");
  djcalo->Add(infile);
  aliases_dijet(djcalo);
  TChain * djcalopix = new TChain("djcalo_pxl/djTree");
  djcalopix->Add(infile);
  aliases_dijet(djcalopix);
  TString evtSel("(cent<10 && nljet>120 && abs(nljeta)<2 && aljet>50 && abs(aljeta)<2 && jdphi>2.5 && Aj>0.24)");
  TH1::SetDefaultSumw2();

  TCanvas * c0 = new TCanvas("c0","c0",500,500);
  djcalo->Draw("Aj>>hAj(20,0,1)",evtSel);
  djcalopix->Draw("Aj",evtSel,"Esame");

  TCanvas * c2 = new TCanvas("c2","c2",500,500);
  TH1D * hJDPhi = new TH1D("hJDPhi","hJDPhi",50,0,3.1416);
  TH1D * hJDPhi2 = new TH1D("hJDPhi2","hJDPhi",50,0,3.1416);
  Float_t numDJ = djcalo->Draw("jdphi>>hJDPhi",evtSel);
  Float_t numDJ2 = djcalopix->Draw("jdphi>>hJDPhi2",evtSel);
  cout << "num djs: " << numDJ << "  djs(pix)" << numDJ2 << endl;

  TH1D * hPNDRTrk = new TH1D("hPNDRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPNDRPix = new TH1D("hPNDRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRTrk = new TH1D("hPADRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRPix = new TH1D("hPADRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  djcalo->Draw("pndr>>hPNDRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalopix->Draw("pndr>>hPNDRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalo->Draw("padr>>hPADRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  djcalopix->Draw("padr>>hPADRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  hPNDRTrk->Scale(1./(numDJ*hPNDRTrk->GetBinWidth(1)));
  hPNDRPix->Scale(1./(numDJ*hPNDRPix->GetBinWidth(1)));
  hPADRTrk->Scale(1./(numDJ*hPADRTrk->GetBinWidth(1)));
  hPADRPix->Scale(1./(numDJ*hPADRPix->GetBinWidth(1)));
  hPNDRPix->SetMinimum(50);
  hPNDRPix->SetMaximum(100);
  hPNDRTrk->SetMinimum(10);
  hPNDRTrk->SetMaximum(2000);
  hPNDRPix->SetMarkerStyle(kFullCircle);
  hPNDRTrk->SetMarkerStyle(kFullSquare);
  hPNDRPix->SetMarkerColor(kRed);
  hPNDRTrk->SetMarkerColor(kRed);
  hPADRPix->SetMarkerStyle(kOpenCircle);
  hPADRTrk->SetMarkerStyle(kOpenSquare);
  hPADRPix->SetMarkerColor(kBlue);
  hPADRTrk->SetMarkerColor(kBlue);
  TCanvas * c3 = new TCanvas("c3","c3",500,500);
  hPNDRPix->Draw("E");
  hPADRPix->Draw("Esame");
  TCanvas * c3_2 = new TCanvas("c3_2","c3_2",500,500);
  c3_2->SetLogy();
  hPNDRTrk->Draw("E");
  hPADRTrk->Draw("Esame");
}

TH1D *
GetITSsaSpectrum(TFile *file, Int_t part, Int_t charge, Int_t cent, Bool_t cutSpectrum = kTRUE, Bool_t addSystematicError = kTRUE)
{
  /* pt limits for combined spectra */
  Double_t ptMin[AliPID::kSPECIES] = {0., 0., 0.1, 0.2, 0.3};
  Double_t ptMax[AliPID::kSPECIES] = {0., 0., 0.6, 0.5, 0.6};

  TList *list = (TList *)file->Get("output");
  TH1D *hin = (TH1D *)list->FindObject(Form("h_%s_%s_cen_%d", ITSsaPartName[part], ITSsaChargeName[charge], cent));
  if (!hin) return NULL;  

  /* get systematics */
  TFile *fsys = TFile::Open("SPECTRASYS_ITSsa.root");
  TH1 *hsys = fsys->Get(Form("hSystTot%s%s", ITSsaChargeName[charge], ITSsaPartName[part]));
			
  TH1D *h = new TH1D(Form("hITSsa_cent%d_%s_%s", cent, AliPID::ParticleName(part), chargeName[charge]), "ITSsa", NptBins, ptBin);
  Double_t pt, width, value, error, sys;
  Int_t bin;
  for (Int_t ipt = 0; ipt < NptBins; ipt++) {
    /* get input bin */
    pt = h->GetBinCenter(ipt + 1);
    width = h->GetBinWidth(ipt + 1);
    bin = hin->FindBin(pt);
    /* sanity check */
    if (TMath::Abs(hin->GetBinCenter(bin) - pt) > 0.001 ||
	TMath::Abs(hin->GetBinWidth(bin) - width) > 0.001)
      continue;
    /* check pt limits */
    if (cutSpectrum && (pt < ptMin[part] || pt > ptMax[part])) continue;
    /* copy bin */
    value = hin->GetBinContent(bin);
    error = hin->GetBinError(bin);
    /*** TEMP ADD SYS ***/
    if (addSystematicError) {
      sys = hsys->GetBinContent(bin) * value;
      error = TMath::Sqrt(error * error + sys * sys);
    }
    h->SetBinContent(ipt + 1, value);
    h->SetBinError(ipt + 1, error);
  }

  h->SetTitle("ITSsa");
  h->SetLineWidth(1);
  h->SetLineColor(1);
  h->SetMarkerStyle(20);
  h->SetMarkerColor(1);
  h->SetFillStyle(0);
  h->SetFillColor(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

void XS1S_Rap(){
  
  TFile *f = new TFile("XSection.root");
  TH1D *H;
  H = (TH1D*)gFile->Get("S1YieldEta");
  double tot(0.); double totE(0.); 
  
  ofstream OUT("XS_1S_rap.tex");
  OUT << "% ----------------------------------------------------------------------" << endl;
  OUT << "% -- Yields" << endl;
  for ( int x = 1; x <= H->GetNbinsX(); ++x ){
    OUT << Form("\\vdef{XS_rap%iS_bin%iContent} {\\ensuremath{ {%.3f } } }",1, x, H->GetBinContent(x)*H->GetBinWidth(x) ) << endl;
    OUT << Form("\\vdef{XS_rap%iS_bin%iError} {\\ensuremath{ {%.3f } } }",1, x, H->GetBinError(x)*H->GetBinWidth(x) ) << endl;
    tot += H->GetBinContent(x)*H->GetBinWidth(x);
    totE += (H->GetBinError(x)*H->GetBinWidth(x)*H->GetBinWidth(x)*H->GetBinError(x));
  }
		
  OUT << Form("\\vdef{sum%iS} {\\ensuremath{ {%.4f } } }",1, tot) << endl;
  OUT << Form("\\vdef{sum%iSError} {\\ensuremath{ {%.4f } } }",1 , TMath::Sqrt(totE) ) << endl;
  OUT.close();
  
}

void RRootHistRead_v1(){
	
	TFile *inputFile = new TFile("pPb_MBSpectra_Combine_-1_1.root");
	cout << "Is the file Open " << inputFile->IsOpen() << endl;

	//Copy the histogram
	TString hName = "Spectra_NtrkOffline0_inf";
	TH1D *h = (TH1D*)inputFile->Get(hName);

	cout << "Number of bins in the histogram: " << h->GetSize() << "\n";
	h->Draw();

	double nhBins = h->GetSize();
	for (int i=0; i < nhBins; i++) cout << h->GetBinWidth(i) << ", ";
	cout << endl;

	// Fit function
	double nParam = 2, maxValue = 5, minValue = 0;

	TF1 *fitFn = new TF1("fitFn", fitFunction, minValue, maxValue, nParam);

	// Set parameters for Fit
	fitFn->SetParameter(0, 0.008);

	// Get parameters from Fit
	double paramFromFit[2];
	fitFn->GetParameters(paramFromFit);

	// Generate data from fit
	double nBins = 27, pT[27], spectraFromFit[27]; // From the Xi 
	for (int i=0; i < nBins; i++){

		spectraFromFit[i] = paramFromFit[0]*exp(-i) + paramFromFit[1];
		cout << spectraFromFit[i] << ", ";
		pT[i] = i;
			
	}

	// Using Integral for the Fit
	h->Fit(fitFn, "I");

	TGraph *g = new TGraph(nBins, pT, spectraFromFit);
	g->SetMarkerStyle(7);
	g->SetMarkerSize(1.5);
	g->SetMarkerColor(kGreen);
	g->Draw("ALP");
	h->Draw("same");
	
}

TH1F* getTruthCrossSection(TFile *file, double x, TString name){

  TTree* genTree = (TTree*)file->Get("ttbarGenAna/tree");
  TH1D* h = new TH1D(Form("h_%s",name.Data()),"h",nGen, genBins);
  genTree->Project(Form("h_%s",name.Data()), "mTT");
  TH1F* hSigma = new TH1F(Form("hSigma_%s",name.Data()),"hSigma", nGen, genBins);

  for(int i=1; i <= nGen; i++){
    double dsigma = x * h->GetBinContent(i) / ( h->GetBinWidth(i) * h->Integral() );
    hSigma->SetBinContent(i,dsigma);
  }

  return hSigma;

}

TH1D* getFF_pp(double jetPt_low, double jetPt_high, const char* histTitle, int mode)
{
  TH1D * FF = new TH1D(histTitle,";p_{T trk};#frac{1}{N_{jet}} #frac{dN_{trk}}{dp_{t trk}}",trkBins ,yAxis);

  getSpectra(mode);

  //looping over bins in Trk pt
    for(int t = 1; t < trkBins+1; t++)
    {
      double sum   = 0;
      double error = 0;
      double norm  = 0;        

      //looping over jet pt
      for(int j = jet->FindBin(jetPt_low); j < jet->FindBin(jetPt_high); j++)
      {
        //reweighting factor to 5 TeV pPb
        //for (mode 2) 5TeV pPb jet_pPb = jet so this is by definition one
        double w = jet_pPb->GetBinContent(j)/jet->GetBinContent(j);

        //adding up contributions to the FF from each track and jet bin
        sum   += w*(trk->GetBinContent(j,t) - trkUE->GetBinContent(j,t));
        error += w*w*(TMath::Power(trk->GetBinError(j,t),2) + TMath::Power(trkUE->GetBinError(j,t),2));
        norm  += jet_pPb->GetBinContent(j);
      }

      if(norm!=0)
      {   
        sum   = sum/(norm*FF->GetBinWidth(t));
        error = TMath::Power(error,0.5)/(norm*FF->GetBinWidth(t));
        FF->SetBinContent(t,sum);
        FF->SetBinError(t,error);
      }
    }
  return FF;
}

void AdjustDensityForBinWidth(TH1D &h) {
    double entries = h.GetEntries();
    int nbins = h.GetNbinsX();
    double low = h.GetBinLowEdge(1);
    double high = h.GetBinLowEdge(nbins+1);
    double width = (high-low)/nbins;
    for(int bin = 1; bin <= nbins; ++bin) {
        double content = h.GetBinContent(bin);
        double error = h.GetBinError(bin);
        double this_width = h.GetBinWidth(bin);
        double scale = width/this_width;
        h.SetBinContent(bin, content*scale);
        h.SetBinError(bin, error*scale);
    }
    h.SetEntries(entries);
}

void writeDataBackgroundHistosForModel(const std::map<TString,TGraph*>& m_bkgds,
				       const std::vector<TH1D *>& vchans,
				       TFile  *allHistFile)
{
  for (std::map<TString,TGraph*>::const_iterator it = m_bkgds.begin();
       it != m_bkgds.end();
       it++) {
    const TString& name = it->first;

    // determine binning from the signal histogram for this channel
    // - (sigh) have to find it first...
    //
    TString channame = name.Copy().Remove(0,strlen("Bckgrndtot_"));
    TH1D *sigh=(TH1D*)0;
    for (int ichan=0; ichan<NUMCHAN; ichan++) {
      sigh = vchans.at(ichan);
      if (strstr(sigh->GetName(),channame.Data()))
	break;
    }
    assert (sigh);


    // for variable binning - all histos must have the same binning per channel
    TAxis *xax = sigh->GetXaxis();

    TVectorD xbins   = TVectorD(sigh->GetNbinsX(),xax->GetXbins()->GetArray());

    int lobin = xax->FindFixBin(sumwinmin);
    int hibin = xax->FindFixBin(sumwinmax)-1;
    int nbins = hibin-lobin+1;

    TVectorD xwindow = xbins.GetSub(lobin-1,hibin);
    
    printf("Booking TH1D(%s,%s,%d,xwindowarray)\n",
	   name.Data(),name.Data(),nbins);
    TH1D *h = new TH1D(name.Data(),name.Data(),nbins,xwindow.GetMatrixArray());

    for (int ibin=1; ibin <= nbins; ibin++)
      h->SetBinContent(ibin,
		       it->second->Eval(h->GetBinCenter(ibin))
		       * h->GetBinWidth(ibin)
		       );

    allHistFile->WriteTObject(h);
  }

}                                            //  writeDataBackgroundHistosForModel

//.........这里部分代码省略.........
    hbackground_1D[i] = (TH1D*)hbackground[i]->ProfileY(Form("background_1D_trg%d_ass0",i+1),1,-1,"e");
    hcorr_1D[i] = new TH1D(Form("Vn_1D_trg%d_ass0",i+1),";n;V_{n#Delta}",15,0.5,15.5);

    for(int n=1;n<=hcorr_1D[i]->GetNbinsX();n++)
    {
      if(n>=6 && hpt[i]->GetMean()>5.0) continue;
      if(n==5 && hpt[i]->GetMean()>8.0) continue;
      if(n==4 && hpt[i]->GetMean()>10.0) continue;

      double Vn=hsignal_1D[i]->GetBinContent(n)-hbackground_1D[i]->GetBinContent(n);
      double VnError=sqrt(hsignal_1D[i]->GetBinError(n)*hsignal_1D[i]->GetBinError(n)+hbackground_1D[i]->GetBinError(n)*hbackground_1D[i]->GetBinError(n));

      VnArray[n][i] = Vn;
      hcorr_1D[i]->SetBinContent(n,Vn);
      hcorr_1D[i]->SetBinError(n,VnError);
      
      double vn=Vn/vn0[n];
      if(i<=2) vn=Vn/vn01[n];
      double vnerr=vn*sqrt((VnError/Vn)*(VnError/Vn)+vn0err[n]*vn0err[n]);
      if(i<=2) vnerr=vn*sqrt((VnError/Vn)*(VnError/Vn)+vn01err[n]*vn01err[n]);
      gr[n-1]->SetPoint(i,hpt[i]->GetMean(),vn);
      gr[n-1]->SetPointError(i,0,vnerr);

/*
      double integral = hpt[i]->Integral();
      if(i==0) integral = hpt[1]->Integral();
      if(hpt[i]->GetMean()>3.0) continue;
      merit[n-1] =  merit[n-1] + integral/1000.*vn;
      merit_err[n-1] =  merit_err[n-1] + integral/1000.*vnerr;
      total[n-1] =  total[n-1] + integral/1000.;
      gr_merit[n-1]->SetPoint(i,hpt[i]->GetMean(),merit[n-1]/sqrt(total[n-1]));
if(n==2) cout<<"n="<<n<<" i="<<i<<" "<<merit[n-1]/sqrt(total[n-1])<<endl;
*/
      double integral = hptcorr->GetBinContent(hptcorr->FindBin(hpt[i]->GetMean()))*hptcorr->GetBinWidth(hptcorr->FindBin(hpt[i]->GetMean()));
      if(hpt[i]->GetMean()>3.0) continue;
      merit[n-1] =  merit[n-1] + integral/1000.*vn;
      merit_err[n-1] =  merit_err[n-1] + integral/1000.*vnerr;
      total[n-1] =  total[n-1] + integral/1000.;
      gr_merit[n-1]->SetPoint(i,hpt[i]->GetMean(),merit[n-1]/sqrt(total[n-1]));
if(n==2) cout<<"n="<<n<<" i="<<i<<" "<<merit[n-1]/sqrt(total[n-1])<<endl;
    }
  }

  TH1D* h1corr = (TH1D*)Get1DCFdPhiv3(filename.Data(),0,0,1.9,4.0);
  TH1D* h2corr = (TH1D*)Get1DCFdPhiv3(filename.Data(),11,0,1.9,4.0);
  TCanvas* ccorr = new TCanvas("ccorr","ccorr",900,400);
  ccorr->Divide(2,1);
  ccorr->cd(1);
  ccorr->GetPad(1)->SetLeftMargin(0.23);
  h1corr->SetAxisRange(50.9,52.,"Y");
  h1corr->GetYaxis()->SetTitleOffset(1.2);
  h1corr->GetXaxis()->SetLabelSize(h1corr->GetXaxis()->GetLabelSize()*0.98);
  h1corr->GetYaxis()->SetLabelSize(h1corr->GetXaxis()->GetLabelSize()*0.98);
  h1corr->GetXaxis()->SetTitleSize(h1corr->GetXaxis()->GetTitleSize()*0.98);
  h1corr->GetYaxis()->SetTitleSize(h1corr->GetXaxis()->GetTitleSize()*0.98);
  h1corr->GetYaxis()->SetTitleOffset(h1corr->GetYaxis()->GetTitleOffset()*1.4);
  h1corr->GetXaxis()->CenterTitle();
  h1corr->GetYaxis()->CenterTitle();
  h1corr->Draw("PE");
  ccorr->cd(2);
  ccorr->GetPad(2)->SetLeftMargin(0.23);
  h2corr->SetAxisRange(50.8,52.4,"Y");
  h2corr->GetYaxis()->SetTitleOffset(1.2);
  h2corr->GetXaxis()->SetLabelSize(h2corr->GetXaxis()->GetLabelSize()*0.98);
  h2corr->GetYaxis()->SetLabelSize(h2corr->GetXaxis()->GetLabelSize()*0.98);
  h2corr->GetXaxis()->SetTitleSize(h2corr->GetXaxis()->GetTitleSize()*0.98);

//.........这里部分代码省略.........
    W[8] = new TFile ("./root/TDAna_W41w_tk3.root");
    W[9] = new TFile ("./root/TDAna_W50w_tk3.root");
    W[10] = new TFile ("./root/TDAna_W51w_tk3.root");
    double Wxs[nwsamples] = { 45000., 9200., 250., 2500., 225., 590., 100., 125., 40., 85., 40. };
    double NW[nwsamples] = { 88000., 40000., 100530., 99523., 105255., 79000.,
                             88258., 83038., 30796., 59022., 41865.
                           };

    TFile * TTH = new TFile("./root/TDAna_ttH_120_tk3.root");
    double TTHxs = 0.667 ;
    double NTTH = 1000000.; // 62000.; // 1652000.; // 62000.;

    const int nttsamples=5;
    TFile * TT[nttsamples];
    TT[0] = new TFile("./root/TDAna_TT0_tk3.root");
    TT[1] = new TFile("./root/TDAna_TT1_tk3.root");
    TT[2] = new TFile("./root/TDAna_TT2_tk3.root");
    TT[3] = new TFile("./root/TDAna_TT3_tk3.root");
    TT[4] = new TFile("./root/TDAna_TT4_tk3.root");
    // double TTxs[5] = { 619., 176., 34.,  6., 1.5 };  // from web
    double TTxs[nttsamples] = { 434., 162., 43., 10., 1.9 };     // from note
    double NTT[nttsamples] = { 57900., 66000., 98159., 14768., 5304. };

    double Lumfactor = 100000; // 100/fb of luminosity assumed in histograms

    TH1D * Histo_TOT[nvars];
    TH1D * Histo_TTH[nvars];
    TH1D * Histo_TOTS[nvars];
    TH1D * Histo_TTHS[nvars];
    for ( int i=0; i<nvars; i++ ) {
        cout << i << endl;
        TH1D * H = dynamic_cast<TH1D*>(TTH->Get(pippo[i]));
        double minx=H->GetBinLowEdge(1);
        double maxx=nbins*H->GetBinWidth(1);
        Histo_TOT[i] = new TH1D ( pippotot[i],pippotot[i], nbins, minx, maxx );
        Histo_TTH[i] = new TH1D ( pippotth[i],pippotth[i], nbins, minx, maxx );
        Histo_TOTS[i] = new TH1D ( pippototS[i],pippototS[i], nbins, minx, maxx );
        Histo_TTHS[i] = new TH1D ( pippotthS[i],pippotthS[i], nbins, minx, maxx );
    }

    cout << "Starting loop on variables needing smoothing" << endl;

    // Loop on variables
    // -----------------
    for ( int ivar=0; ivar<nvars; ivar++ ) {

        // Extract sum histograms with the right normalization and errors
        // --------------------------------------------------------------
        double totWW[nwsamples][nbins]= {0.};
        double totW[nbins]= {0.};
        double s2_totW[nbins]= {0.};
        double totNW[nwsamples][nbins]= {0.};
        for ( int i=0; i<nwsamples; i++ ) {
            //cout << "Processing W file #" << i << " ..." << endl;
            TH1D * Histo = dynamic_cast<TH1D*>(W[i]->Get(pippo[ivar]));
            TH1D * HistoW = dynamic_cast<TH1D*>(W[i]->Get(pippo[ivar]+"W"));
            // For W, we need also total entries in histograms to add a
            // Poisson fluke contribution to total errors from matrix:
            // ----------------------------------------------------------
            TH1D * HistoN = dynamic_cast<TH1D*>(W[i]->Get(pippo[ivar]+"N"));
            for ( int ibin=1; ibin<=nbins; ibin++ ) {
                double t=Histo->GetBinContent(ibin);
                double s2t=HistoW->GetBinContent(ibin);
                double n=HistoN->GetBinContent(ibin);
                totWW[i][ibin-1]+=t*Wxs[i]/NW[i]*Lumfactor;
                s2_totW[ibin-1]+=s2t*pow(Wxs[i]/NW[i]*Lumfactor,2);

//.........这里部分代码省略.........
  if(iSpec == 2 &&  Weight==0  ) Z0Mass_1S = (TH2D*)fil2->Get("diMuonsGlobalInvMassVsY");
  if(iSpec == 3 &&  Weight==0  ) Z0Mass_1S = (TH2D*)fil2->Get("diMuonsGlobalInvMassVsCen");
  
  TH1D *service = (TH1D*)Z0Mass_1->ProjectionY("service");
  int pt_bin_bound[100];
  TH1D *dimuonsGlobalInvMassVsPt[100],  *dimuonsGlobalInvMassVsPtS[100];
  
  TCanvas *CanvPt_1 = new TCanvas("CanvPt_1"," Z0 Yield Vs. Pt ", 40,40,1000,700);
  if (Nptbin == 2)  CanvPt_1->Divide(2,1);
  if (Nptbin == 3 || Nptbin == 4)  CanvPt_1->Divide(2,2);
  if (Nptbin == 5 || Nptbin == 6)  CanvPt_1->Divide(3,2);
  if (Nptbin == 9 || Nptbin == 10)  CanvPt_1->Divide(5,2);
  cout << endl << Xname[iSpec] << "    Yield      Mass (GeV)    Width (GeV)    GauWidth    chi2/ndf " << endl << endl; 
  //ih loop
  for (Int_t ih = 0; ih < Nptbin; ih++) 
    {
      CanvPt_1->cd(ih+1);
      gPad->SetTickx();
      gPad->SetTicky();
      
      // Project 1 D 
      pt_bin_bound[ih] = Z0Mass_1->GetYaxis()->FindBin(pt_bound[ih]+0.0000001);
      pt_bin_bound[ih+1] = Z0Mass_1->GetYaxis()->FindBin(pt_bound[ih+1]-0.0000001);
      sprintf(namePt_1,"Z0_1_pt_%d",ih);
      sprintf(namePt_1S,"Z0_1S_pt_%d",ih);
      
      //printf(namePt_1,"Z0_1_pt_%d",ih);
      //cout<<endl<<endl;
      //continue;
      
      dimuonsGlobalInvMassVsPt[ih] = (TH1D*)Z0Mass_1->ProjectionX(namePt_1, pt_bin_bound[ih], pt_bin_bound[ih+1]-1+1, "e");
      dimuonsGlobalInvMassVsPtS[ih] = (TH1D*)Z0Mass_1S->ProjectionX(namePt_1S, pt_bin_bound[ih], pt_bin_bound[ih+1]-1+1,"e");
      if(iSpec == 1 || iSpec == 2) {sprintf(text," %s [%.1f, %.1f]", Xname[iSpec], service->GetBinLowEdge(pt_bin_bound[ih]), 
					    service->GetBinLowEdge(pt_bin_bound[ih+1])+service->GetBinWidth(pt_bin_bound[ih+1]));}
      if(iSpec == 3) {sprintf(text," %s [%.1f, %.1f] %s", Xname[iSpec], 2.5*service->GetBinLowEdge(pt_bin_bound[ih]), 
			      2.5*(service->GetBinLowEdge(pt_bin_bound[ih+1])+service->GetBinWidth(pt_bin_bound[ih+1])), "%");}
      
      dimuonsGlobalInvMassVsPt[ih]->Rebin(nrebin);
      dimuonsGlobalInvMassVsPtS[ih]->Rebin(nrebin);
      
      float m_low = 60.0;
      float m_high = 120.0;
      
      TAxis *axs   = dimuonsGlobalInvMassVsPt[ih]->GetXaxis();
      int binlow = axs->FindBin(m_low);
      int binhi  = axs->FindBin(m_high);
      //      Double_t bin_size = (1.0*dimuonsGlobalInvMassVsPt[ih]->GetNbinsX())/(axs->GetXmax() - axs->GetXmin());
      //    Float_t int_sig = 0.0;
      //for(Int_t bin = binlow; bin<=binhi;bin++) {
      //int_sig+ = dimuonsGlobalInvMassVsPt[ih]->GetBinContent(bin);
      //}   
      double recerror,recerrorS;
      //double yld;
      double yld_1 = dimuonsGlobalInvMassVsPt[ih]->IntegralAndError(binlow, binhi, recerror);
      double yldS_1 = dimuonsGlobalInvMassVsPtS[ih]->IntegralAndError(binlow, binhi, recerrorS); 
      eyld_cat_1[ih] =recerror;
      eyld_catS_1[ih] =recerrorS;
      cout<< "yld " << yld_1 <<"  eyld_cat_1[ih]  "<< eyld_cat_1[ih]<<" rec error "<<recerror<<endl;
      
      yld_cat_1[ih] = yld_1;
      yld_catS_1[ih] = yldS_1;
      
      if(isLog) gPad->SetLogy(1);
      TColor *pal = new TColor();
      Int_t kblue = pal->GetColor(9,0,200);
      //      Int_t korange  = pal->GetColor(101, 42,  0);

TH1D * getNiceHistogram(int nbins, std::vector<Double_t> binning, TFile * file, TString hist) {
  TH1D *sample = (TH1D*)(file->Get(hist)->Clone());
  TH1D *binned = dynamic_cast<TH1D*>(sample->Rebin(nbins,"madgraph",&binning.front()));
  for (Int_t bin=0; bin < nbins; bin++) {
    // Divide rebinned histogram's bin content by bin width factor (new/old):
    binned->SetBinError(bin+1,sqrt(binned->GetBinContent(bin+1))/((binning.at(bin+1)-binning.at(bin))/sample->GetBinWidth(1)));
    binned->SetBinContent(bin+1,binned->GetBinContent(bin+1)/((binning.at(bin+1)-binning.at(bin))/sample->GetBinWidth(1)));
  }
  binned->Scale(1/binned->Integral("width"));
  return binned;
}

void Smooth (TString sel) 
{

  const int nbins = 50;

  const int nvars=17;
  TString var[nvars] = { "C8", "M8", "C6", "M6", "MEt", "MEtSig", "CorrSumEt", "GoodHt", 
		      "M45bestall", "Chi2mass", "Chi2extall", "Mbbnoh", "DPbbnoh", 
		      "SumHED4", "SumHED6", "DP12", "MEtDP2" };
  TString pippo[nvars];
  TString pippotot[nvars];
  TString pippotth[nvars];
  TString pippototS[nvars];
  TString pippotthS[nvars];
  for ( int i=0; i<nvars; i++ ) { 
    pippo[i] = var[i]+sel; 
    pippotot[i]=var[i]+sel+"_bgr";
    pippotth[i]=var[i]+sel+"_sig";
    pippototS[i]=var[i]+sel+"_bgrS";
    pippotthS[i]=var[i]+sel+"_sigS";
  }
  
  const int nqcdsamples=8;
  TFile * QCD[nqcdsamples];     
  QCD[0] = new TFile("./root/et30_eta2.5/TDAna_QCD30-50_tk3.root");
  QCD[1] = new TFile("./root/et30_eta2.5/TDAna_QCD50-80_tk3.root");
  QCD[2] = new TFile("./root/et30_eta2.5/TDAna_QCD80-120_tk3.root");
  QCD[3] = new TFile("./root/et30_eta2.5/TDAna_QCD120-170_tk3.root");
  QCD[4] = new TFile("./root/et30_eta2.5/TDAna_QCD170-230_tk3.root");
  QCD[5] = new TFile("./root/et30_eta2.5/TDAna_QCD230-300_tk3.root");
  QCD[6] = new TFile("./root/et30_eta2.5/TDAna_QCD300-380_tk3.root");
  QCD[7] = new TFile("./root/et30_eta2.5/TDAna_QCD380incl_tk3.root");
  double QCDxs[nqcdsamples] = { 155929000., 20938850., 2949713., 499656., 100995.,  23855., 6391., 2821.};
  double NQCD[nqcdsamples] = { 86000., 78000., 104000., 96000., 100000., 102000., 112000., 102000.};


  const int nwsamples=11;
  TFile * W[nwsamples];
  W[0] = new TFile ("./root/et30_eta2.5/TDAna_W0w_tk3.root");
  W[1] = new TFile ("./root/et30_eta2.5/TDAna_W10w_tk3.root");
  W[2] = new TFile ("./root/et30_eta2.5/TDAna_W11w_tk3.root");
  W[3] = new TFile ("./root/et30_eta2.5/TDAna_W20w_tk3.root");
  W[4] = new TFile ("./root/et30_eta2.5/TDAna_W21w_tk3.root");
  W[5] = new TFile ("./root/et30_eta2.5/TDAna_W30w_tk3.root");
  W[6] = new TFile ("./root/et30_eta2.5/TDAna_W31w_tk3.root");
  W[7] = new TFile ("./root/et30_eta2.5/TDAna_W40w_tk3.root");
  W[8] = new TFile ("./root/et30_eta2.5/TDAna_W41w_tk3.root");
  W[9] = new TFile ("./root/et30_eta2.5/TDAna_W50w_tk3.root");
  W[10] = new TFile ("./root/et30_eta2.5/TDAna_W51w_tk3.root");
  double Wxs[nwsamples] = { 45000., 9200., 250., 2500., 225., 590., 100., 125., 40., 85., 40. };
  double NW[nwsamples] = { 88000., 40000., 100530., 99523., 105255., 79000., 
			   88258., 83038., 30796., 59022., 41865. };

  TFile * TTH = new TFile("./root/et30_eta2.5/TDAna_ttH_120_tk3.root");
  double TTHxs = 0.667 ;
  double NTTH = 62000.; // 1652000.; // 62000.;

  const int nttsamples=5;
  TFile * TT[nttsamples];
  TT[0] = new TFile("./root/et30_eta2.5/TDAna_TT0_tk3.root");
  TT[1] = new TFile("./root/et30_eta2.5/TDAna_TT1_tk3.root");
  TT[2] = new TFile("./root/et30_eta2.5/TDAna_TT2_tk3.root");
  TT[3] = new TFile("./root/et30_eta2.5/TDAna_TT3_tk3.root");
  TT[4] = new TFile("./root/et30_eta2.5/TDAna_TT4_tk3.root");
  // double TTxs[5] = { 619., 176., 34.,  6., 1.5 };  // from web
  double TTxs[nttsamples] = { 434., 162., 43., 10., 1.9 };     // from note
  double NTT[nttsamples] = { 57900., 66000., 98159., 14768., 5304. };

  double Lumfactor = 100000; // 100/fb of luminosity assumed in histograms

  TH1D * Histo_TOT[nvars];
  TH1D * Histo_TTH[nvars];
  TH1D * Histo_TOTS[nvars];
  TH1D * Histo_TTHS[nvars];
  for ( int i=0; i<nvars; i++ ) {
    cout << i << endl;
    TH1D * H = dynamic_cast<TH1D*>(TTH->Get(pippo[i]));
    double minx=H->GetBinLowEdge(1);
    double maxx=nbins*H->GetBinWidth(1);
    Histo_TOT[i] = new TH1D ( pippotot[i]," ", nbins, minx, maxx );
    Histo_TTH[i] = new TH1D ( pippotth[i]," ", nbins, minx, maxx );
    Histo_TOTS[i] = new TH1D ( pippototS[i]," ", nbins, minx, maxx );
    Histo_TTHS[i] = new TH1D ( pippotthS[i]," ", nbins, minx, maxx );
  }

  cout << "Starting loop on variables needing smoothing" << endl;

  // Loop on variables 
  // -----------------
  for ( int ivar=0; ivar<nvars; ivar++ ) {

    // Extract sum histograms with the right normalization and errors
    // --------------------------------------------------------------
    double totWW[nwsamples][nbins]={0.};
    double totW[nbins]={0.};
    double s2_totW[nbins]={0.};
    double totNW[nwsamples][nbins]={0.};
    for ( int i=0; i<nwsamples; i++ ) {
      cout << "Processing W file #" << i << " ..." << endl;
      TH1D * Histo = dynamic_cast<TH1D*>(W[i]->Get(pippo[ivar]));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if(chunkName[0]!='\0') fileIdentifier = chunkName;
    TFile *fin = TFile::Open(Form("%s/AnalysisResults_%s.root", inputDir, fileIdentifier.Data()), "READ");
    TTree *treeEvent = (TTree*)fin->Get("MultSelection/fTreeEvent");

    cout << "   - run number....................: " << runNo << endl;

    // define estimator histo for this run
    TH1D* hEstimator = new TH1D(Form("hEstimator_%d", runNo), "", lNDesiredBoundaries, lDesiredBoundaries);
    hEstimator->Sumw2();
    hEstimator->GetXaxis()->SetTitle("V0M Percentile");
    hEstimator->GetYaxis()->SetTitle("Counts");
    hEstimator->SetStats(0);
    hEstimator->SetLineColor(kRed);

    // get corresponding calibration histogram from OADB
    AliOADBMultSelection* lOADB = (AliOADBMultSelection*)lOADBcontainer->GetObject( runNo, "Default" );
    if( (Int_t)lOADBcontainer->GetIndexForRun( runNo )<0 ) {
      cout << "   ---> Warning: no calibration histo found for this run - skipping..." << endl;
      return;
    }

    // set the pointer to the calib histo for this run
    hCalib = (TH1D*)lOADB->GetCalibHisto( "hCalib_V0M" );;
    //
    Double_t nall = treeEvent->Draw(Form("get_percentile(fAmplitude_V0A+fAmplitude_V0C)>>hEstimator_%d", runNo), 
                                        Form("fRunNumber==%d && fEvSel_Triggered && fEvSel_IsNotPileupInMultBins && fEvSel_PassesTrackletVsCluster && fEvSel_INELgtZERO && fEvSel_HasNoInconsistentVertices && TMath::Abs(fEvSel_VtxZ)<=10.0 && isSelectedHM(fEvSel_TriggerMask)", runNo),
                                        "goff"); 

    hEstimator->Scale(1., "width");
    Double_t nevents = (Double_t)hEstimator->GetEntries();
    cout << "   - number of events (selected)...: " << nevents << endl;

    // draw histogram
    TCanvas *cEstimator = new TCanvas(Form("cEstimator_%d", runNo), "Estimator Distribution", 10, 10, 1000, 750);
    cEstimator->SetRightMargin(0.05);
    cEstimator->SetTopMargin(0.11);

    hEstimator->GetXaxis()->SetRangeUser(0., 0.2);
    hEstimator->Draw("hist e0");
    latex->SetNDC();
    latex->SetTextSize(0.06);
    latex->DrawLatex(0.1, 0.93, Form("Run: %d", runNo));

    // first, fit a pol0 in the flat region (usually up to 0.05)
    hEstimator->Fit(fpol0, "RQ0");
    Double_t flat_top = fpol0->GetParameter(0);

    // get standard deviantion of bin contents in the flat region
    Double_t flat_top_stdev = 0.;
    for(Int_t ibin=1; ibin<=hEstimator->FindBin(lMinimumAnchorPercentile); ++ibin) {
      Double_t content = hEstimator->GetBinContent(ibin);
      Double_t   width = hEstimator->GetBinWidth(ibin);
      flat_top_stdev += TMath::Power((content-flat_top), 2.)*width;
    }
    flat_top_stdev = TMath::Sqrt(flat_top_stdev/lMinimumAnchorPercentile) / 2.;
    fpol0_hi->SetParameter(0, flat_top+flat_top_stdev);
    fpol0_lo->SetParameter(0, flat_top-flat_top_stdev);
        
    // now, fix the constant parameter in the turnon function
    fturnon->SetParameters(1., 0.1, -1.);
    fturnon->FixParameter(0, flat_top);
        
    // get the maximum range to perform the fit
    Double_t range_max = (hEstimator->GetBinLowEdge(hEstimator->FindLastBinAbove())) / 1.8;
    fturnon->SetRange(0.005, (range_max>0.1) ? range_max : 0.1);

    // get anchor percentile
    Double_t anchor_percentile = -1.;
    TString fitstatus = "";
    if(nevents>0) {
      TFitResultPtr fitr = hEstimator->Fit(fturnon, "RQM");
      fturnon->Draw("lsame");
      fitstatus = gMinuit->fCstatu;
    }
    cEstimator->Flush();
    cEstimator->Update();
    cout << "   - fit status....................: " << fitstatus << endl;
    if( !fitstatus.Contains("OK") ) {
      if(gROOT->IsBatch()) {
         cout << "   ---> Warning: fit failed! -- skipping this run..." << endl;
         if(!automaticMode) {
            fap = fopen(Form("temp/anchors/Anchor_%s_%d_VHM.txt", lPeriodName.Data(), runNo), "w");
            fprintf(fap, "%d %d %.2lf %lf\n", runNo, runNo, -1., -1.);
         }
         return;
      }
      
      if(!automaticMode) {
         cout << "   - Please, provide an anchor percentile to continue: " << endl;
         cout << "     (entering a negative value will skip this run)" << endl;
         cout << "     >>>> anchor percentile: "; 
         cin >> anchor_percentile;
         if(anchor_percentile<0.) {
            cout << "   ---> Warning: percentile provided is negative -- skipping this run..." << endl;
            fap = fopen(Form("temp/anchors/Anchor_%s_%d_VHM.txt", lPeriodName.Data(), runNo), "w");
            fprintf(fap, "%d %d %.2lf %lf\n", runNo, runNo, -1., -1.);
            return;
         }
      }
      else return;      // in automatic mode we do not create an anchor file