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

double  Getmean(TString sys1, TString sys2, float cent1, float cent2){

  gROOT->SetStyle("Plain");
  TFile *infile = TFile::Open(Form("glau_%s%s_ntuple_1M.root",sys1.Data(),sys2.Data()));
  infile->cd();
  TH2F *pEcc2B = new TH2F("pEcc2B","#epsilon_{2} vs impact paramter in 200 GeV Collisions from Glauber MC;B;#epsilon_{n}",25200,-0.5,251.5,100,0,1);
  if(sys2.Contains("smeared")){
  TNtuple *nt = (TNtuple*)infile->Get("nt");
  nt->Project("pEcc2B","Ecc3G:B");
  }
  else{
  TNtuple *nt = (TNtuple*)infile->Get(Form("nt_%s_%s",sys1.Data(),sys2.Data()));
  nt->Project("pEcc2B","Ecc2:B");
  }
  TH1D* hB = (TH1D*)pEcc2B->ProjectionX("hB",0,-1);
  for(int ibin=0;ibin<hB->GetNbinsX();ibin++){
      if(hB->Integral(0,ibin) >= cent1 /100. * hB->Integral()){ int bin1 = ibin; break;}
  }
  for(int ibin=0;ibin<hB->GetNbinsX();ibin++){
      if(hB->Integral(0,ibin) >= cent2 /100. * hB->Integral()){ int bin2 = ibin; break;}
  }
  cout<< bin1 << "\t" << bin2 << endl;
  TH1D* hEcc2 = (TH1D*)pEcc2B->ProjectionY("hEcc2",bin1,bin2);
//  hEcc2->Draw();
  cout << hEcc2 -> GetEntries() << endl;
  double mean = hEcc2->GetMean();
  return mean;
}

void createSysStatErrorHist(TH1D *hist_pre){

  TH1D *histErr = (TH1D*) hist_pre->Clone("histErr");
  hsyserr = (TH1D*) hist_pre->Clone("hsyserr"), hsyserr->Reset();
  hsyserr->Sumw2();

  cout<<""<<endl;
  for(int i=0;i<histErr->GetNbinsX();i++){ // last bin is meaningless
    double pt = histErr->GetBinCenter(i+1);
    double statferr = histErr->GetBinError(i+1)/histErr->GetBinContent(i+1);
    if(i==(histErr->GetNbinsX()-1)) statferr = 0.9;
    //double sysferr = fsyserr->Eval(pt);
    double sysferr = returnCombinedSysError(pt);
    double sumferr = sqrt(statferr*statferr + sysferr*sysferr);
    //double sumerr = 1.*(sumferr/100);
    double sumerr = 1.*sumferr;
    hsyserr->SetBinContent(i+1,1);
    hsyserr->SetBinError(i+1,sumerr);
    //hsyserr->SetBinError(i+1,0.5); 
    cout<<"[setSysStatError] stat err (frac): "<<statferr<<" syst err (fac): "<<sysferr<<" combined (frac): "<<sumferr<<endl;
  }

  //tsyserr = (TGraphErrors*) TgraphIt(hsyserr);

}

void loglikdistrib(Int_t ntrials = 10000, Bool_t print = kFALSE)
{
  // compute distribution of log likelihood value
  TH1D * hmc   = gStack[gPadNr][gOrder[gPadNr][0]];
  TH1D * hdata = gStack[gPadNr][gMaxProcess-1];
  Int_t nbins = hmc->GetNbinsX();
  Double_t loglik = loglikelihood(hmc, hdata, 1, nbins);
  TH1D * htest = new TH1D(*hdata);
  TH1D * lldistrib = new TH1D("lldistrib", "log(Likelihood) distribution", 
			      1000, loglik-200, loglik+200);
  setopt(lldistrib);
  for (Int_t n = 0; n < ntrials; n++) {
    // generate poisson around theorie
    for (Int_t i = 1; i <= nbins; i++) {
      htest->SetBinContent(i, gRandom->Poisson(hmc->GetBinContent(i)));
    }
    lldistrib->Fill(loglikelihood(hmc, htest, 1, nbins));
  }
  TCanvas * llcanvas = new TCanvas("llcanvas", "Log(Likelihood) distribution", 
				   40, 40, 800, 600);
  setopt(llcanvas);
  lldistrib->SetFillColor(kYellow);
  lldistrib->Draw();
  lldistrib->GetYaxis()->SetTitle("Anzahl Ereignisse");
  lldistrib->GetXaxis()->SetTitle("-ln L");
  // autozoom
  Int_t lowbin = 1;
  while (lldistrib->GetBinContent(lowbin) == 0)
    lowbin++;
  Int_t highbin = lldistrib->GetNbinsX();
  while (lldistrib->GetBinContent(highbin) == 0)
    highbin--;
  lldistrib->SetAxisRange(lldistrib->GetBinLowEdge(lowbin), 
			  lldistrib->GetBinLowEdge(highbin));
  TH1D * hworse = (TH1D *) lldistrib->Clone();
  for (Int_t nbin = 1; nbin < 501; nbin++) {
    hworse->SetBinContent(nbin, 0);
  }
  hworse->SetFillColor(95);
  hworse->Draw("same");
  Double_t pvalue = lldistrib->Integral(501,1000) / lldistrib->Integral();
  TLatex * tex = new TLatex(0.18, 0.96, Form("-ln L_{obs} = %5.2f", loglik));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  tex = new TLatex(0.18, 0.86, Form("CL_{obs} = %.3f", pvalue));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  TLine * l = new TLine(loglik, 0, loglik, lldistrib->GetMaximum());
  l->SetLineWidth(3);
  l->SetLineColor(kBlue);
  l->Draw();
  llcanvas->Modified();
  llcanvas->Update();
  if (print)
    llcanvas->Print("lldistrib.pdf");
  cd(gPadNr+1);
}

void takeDirectlyFromMC(TFile* fin, TFile* fout, TString gentype) {

  TList* listOfDirs = fin->GetListOfKeys();
  for (auto k : *listOfDirs) {
    TString srname = k->GetName();
    if (!srname.Contains("sr")) continue;
    if (srname.Contains("base") || srname.Contains("incl") || srname.Contains("sb")) continue;
    if (gentype == "_1lepW" && !srname.EndsWith("2") && !srname.EndsWith("3")) continue;

    auto indir = (TDirectoryFile*) fin->Get(srname);
    auto outdir = (TDirectory*) fout->mkdir(srname);
    auto hlist = indir->GetListOfKeys();

    for (auto h : *hlist) {
      TString hname = h->GetName();
      if (!hname.BeginsWith("h_metbins" + gentype)) continue;
      TH1D* hin = (TH1D*) indir->Get(hname);
      outdir->cd();
      TH1D* hout = (TH1D*) hin->Clone(TString(hname).ReplaceAll(gentype, ""));
      for (int i = 1; i <= hout->GetNbinsX(); ++i) {
        // zero out negative yields
        if (hout->GetBinContent(i) < 0) {
          hout->SetBinContent(i, 0);
          hout->SetBinError(i, 0);
        }
      }
      hout->Write();

      if (yearSeparateSyst && (hname.EndsWith("Up") || hname.EndsWith("Dn"))) {
        for (int i = 1; i < 4; ++i) {
          TH1D* hcen_yi = (TH1D*) fbkgs[i]->Get(srname+"/h_metbins"+gentype);
          TH1D* hsys_yi = (TH1D*) fbkgs[i]->Get(srname+"/"+hname);
          if (hsys_yi && !hcen_yi) {
            cout << "Find " << srname+"/"+hname << " from " << fbkgs[i]->GetName() << " but not hcen " << srname+"/h_metbins"+gentype << " Should not happen?" << endl;
          }
          TH1D* hout_yi = (TH1D*) fin->Get(srname+"/h_metbins"+gentype)->Clone(TString(hname).Insert(hname.Length()-2, Form("%d", 15+i)).ReplaceAll(gentype, ""));
          if (hcen_yi) hout_yi->Add(hcen_yi, -1);
          if (hsys_yi) hout_yi->Add(hsys_yi);
          hout_yi->Write();
        }
      }
    }
    if (!outdir->Get("h_metbins")) {
      cout << "Didn't find yield hist for " << gentype << " in " << fin->GetName() << ":" << srname << "/. Faking a 0 one!" << endl;
      outdir->cd();
      // Get the MET binning from h_metbins, which shall always exist, and then set all bins to 0
      TH1D* hout = (TH1D*) fin->Get(srname + "/h_metbins")->Clone("h_metbins");
      for (int i = 1; i <= hout->GetNbinsX(); ++i) {
        hout->SetBinContent(i, 0);
        hout->SetBinError(i, 0);
      }
      hout->Write();
    }
  }
}

TH1D GetSOverSqrtB(TH1D sig, TH1D back){
  if(sig.GetNbinsX()!=back.GetNbinsX()) return sig;
  sig.GetYaxis()->SetTitle("S/#sqrt{B}");
  for(int i(1); i<=sig.GetNbinsX(); ++i){
    if(back.GetBinContent(i)>0.0 && sig.GetBinContent(i)>0.0){
      sig.SetBinContent(i,sig.GetBinContent(i)/sqrt(back.GetBinContent(i)));
    }else{
      sig.SetBinContent(i,0.0);
    }
  }
  return sig;
}

TH1D* getQCD(int rebinFact){
	TString dir = "rootFiles/";
	
	TFile* file = new TFile(dir +"qcdest.root");		
	TH1D* plot = (TH1D*) file->Get("muon_AbsEta_0btag");

// 	for(int i = 1; i <= plot->GetNbinsX(); i++){
// 	plot->SetBinError(i, 0.0);
// 	}

	plot->SetFillColor(kYellow);
	plot->SetLineColor(kYellow);
	plot->SetMarkerStyle(1);
		
	TH1D* copyplot = new TH1D("qcd plot", "qcd plot", 30, 0.0, 3.0);
	
	for(int i = 1; i <= plot->GetNbinsX(); i++){
	copyplot->SetBinContent(i, plot->GetBinContent(i));
	//copyplot->SetBinError(i, plot->GetBinError(i));
	}
	copyplot->SetFillColor(kYellow);
	copyplot->SetLineColor(kYellow);
	copyplot->SetMarkerStyle(1);
	copyplot->Scale(1./copyplot->Integral());	
	copyplot->Rebin(rebinFact);
	
	//file->Close("R");
	return copyplot;
	//file->Close();
}

void QinvKstar (const char* filename, const char* histname) {

	TFile*  fileIn = new TFile(filename,"update");
	TH1D* hist = (TH1D*)fileIn->Get(histname);
	TH1D* histnew = new TH1D(Form("kstar_%s",histname),"",hist->GetNbinsX(),0,0.25);//(TH1D*)fileIn->Get(histname);

//   for (int ibins = 1; ibins <= hist->GetNbinsX()*2; ibins++ ) {
//     hist->SetBinContent(ibins, hist->GetBinContent(ibins));
//     //    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
//   }

	int bin=2;//1

	for (double kstar = 0.00375; kstar <= 0.25; kstar += 0.0025 ) {

//	for (double kstar = 0.00375; kstar <= 0.25; kstar += 0.0025 ) {
		//  for (double kstar = 0.005; kstar <= 0.25; kstar += 0.01 ) {
		//        hist->Fill(kstar,hist->GetBinContent(hist->FindBin(kstar*2)));

		histnew->SetBinContent(bin,hist->GetBinContent(hist->FindFixBin(kstar*2)));
		histnew->SetBinError(bin++,hist->GetBinError(hist->FindFixBin(kstar*2)));

		//    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
		//    hist->SetBinContent(ibins, hist->GetBinContent(ibins));
	}

	TFile* fileOut = new TFile(Form("kstar_%s",filename), "update");
	hist->Write();
	histnew->Write();
}

// Merge same modes
int dmMerge(TObject *a1, TObject *b1,TList *keys)
{
	TDecayMode *a=(TDecayMode*)a1;
	TDecayMode *b=(TDecayMode*)b1;
	TIter  nexthist(keys);
	TKey  *key_hist=0;
	while(key_hist=(TKey*)nexthist())
	{
		int cycle=-10;
		if(strcmp(key_hist->GetClassName(),"TH1D")==0)
		{
			TH1D  *h=0;
			int cycleh=key_hist->GetCycle();
			if(cycleh<cycle) { cout<<"Skipping..."<<endl; continue; }
			if(cycle<cycleh) cycle=cycleh;
			h=(TH1D*)key_hist->ReadObj();
			if(!h)
			{
				cout<<"Cannot read: "<<key_hist->GetName()<<endl;
				exit(-2);
			}
			TH1D *eh = (TH1D*) a->histograms->FindObject(h->GetName());
			if(!eh) continue;
			if(eh->GetNbinsX()!=h->GetNbinsX() ||
			   eh->GetXaxis()->GetXmax()!=h->GetXaxis()->GetXmax())
			   return -1;
			eh->Add(h);
		}
	}
        a->SetNEntries(a->GetNEntries()+b->GetNEntries());
        a->SetSumw(a->GetSumw()+b->GetSumw());
        a->SetSumw2(a->GetSumw2()+b->GetSumw2());
	return 0;
}

std::pair<float,float> getFitRange(TH1D& hist) {
  int maxBin = hist.GetMaximumBin();
  float targetVal = hist.GetBinContent( maxBin ) /2.;


  int nBins = hist.GetNbinsX();
  
  std::pair<float,float> results = {hist.GetBinLowEdge(1),hist.GetBinLowEdge(nBins+1)};

  int diff=0;
  while( maxBin-(++diff) >0) {
    if( hist.GetBinContent(maxBin-diff) < targetVal ){
      results.first = hist.GetBinCenter(maxBin-diff);
      break;
    }
  }

  diff=0;
  while( maxBin+(++diff) <=nBins) {
    if( hist.GetBinContent(maxBin+diff) < targetVal ){
      results.second = hist.GetBinCenter(maxBin+diff);
      break;
    }
  }
  
  
  return results;
}

void getStatPowerLoss() {

  TFile* mc = TFile::Open("MyMCPileupHistogram.root");
  TH1D* hmc = (TH1D*) mc->Get("pileup");
  
  TFile* wgt = TFile::Open("MyRatioPileupHistogram.root");
  TH1D* hwgt = (TH1D*) wgt->Get("pileup");
  
  double old_sumx = 0;
  double old_sumx2 = 0;
  double new_sumx = 0;
  double new_sumx2 = 0;
  for (int i=0; i<=hmc->GetNbinsX()+1; i++) {
    double val = hmc->GetBinContent(i);
    double factor = hwgt->GetBinContent(i);
    old_sumx += val;
    old_sumx2 += val*val;
    
    new_sumx += val*factor;
    new_sumx2 += val*factor*val*factor;
  }
  
  double old_stat = old_sumx*old_sumx/old_sumx2;
  double new_stat = new_sumx*new_sumx/new_sumx2;

  double loss = 1-new_stat/old_stat;
  
  cout << "The fractional loss in statistical power is " << std::setprecision(4) << loss*100 << "%" << endl;

  mc->Close();
  wgt->Close();
}

Double_t fTH1toTF1(Double_t *x, Double_t *par)
{
	Double_t g4  = x[0];          // -> g^4 value
	Int_t    mll = (Int_t)par[0]; // -> dilepton invariant mass bin number
	
	TH2D* h2 = (TH2D*)h2Template->Clone();
	TH1D* hTmp = (TH1D*)TH2toTH1(h2,mll); // TH1 histogram of g^4 in the dilepton mass bin number mll
	
	Int_t nbinsx = hTmp->GetNbinsX();
	Double_t xval = -999.;
	if(doInterpolation)
	{
		Int_t bin = hTmp->FindBin(g4);
		
		if(bin==0 || bin==1)   xval = hTmp->GetBinContent(1);     // cannot interpolate on the first bin or the underflow
		else if(bin==nbinsx+1) xval = hTmp->GetBinContent(bin-1); // cannot interpolate on the last bin or the overflow
		else if(bin==nbinsx)   xval = hTmp->GetBinContent(bin);   // cannot interpolate on the last bin or the overflow
		else                   xval = hTmp->Interpolate(g4);      // interpolated dN/dg^4 at g^4, between the 2 adjacent bin centers
	}
	else
	{
		Int_t bin = hTmp->FindBin(g4);
		if(bin==nbinsx+1) xval = hTmp->GetBinContent(bin-1); // overflow 
		else if(bin==0)   xval = hTmp->GetBinContent(1);     // underflow
		else              xval = hTmp->GetBinContent(bin);   // normal
	}
	delete h2;
	delete hTmp;
	return xval;
}

/** Personalized
build a TGraphError starting from a TH1F
*/
TGraphErrors buildGEfromH_Personalized (const TH1D & histo) 
{
 TVectorF xV(histo.GetNbinsX());
 TVectorF yV(histo.GetNbinsX());
 TVectorF errxV(histo.GetNbinsX());
 TVectorF erryV(histo.GetNbinsX());
 for (int iBin = 0; iBin<histo.GetNbinsX(); iBin++){
  xV[iBin] = histo.GetBinCenter(iBin);
  yV[iBin] = histo.GetBinContent(iBin);
  errxV[iBin] = 0;
//   histo.GetBinWidth(iBin);
  erryV[iBin] = histo.GetBinError(iBin);
 }  
 TGraphErrors g (xV, yV, errxV, erryV) ;
 return g ;
}

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

void Normalize(TH1D &h, double normalization, bool norm_per_avg_width) {
    int nbins = h.GetNbinsX();
    double low = h.GetBinLowEdge(1);
    double high = h.GetBinLowEdge(nbins+1);
    double width = (high-low)/nbins;
    if(norm_per_avg_width) normalization *= width;
    double integral = h.Integral("width");
    h.Scale(normalization/integral);
}

TH1D* Prediction::GetScaledHisto(TH1D* histo, float scale_fact, float scale_fact_err){
	// takes histo, scale factor and uncertainty on scale factor
	// and returns scaled and rebinned histo with 1 bin with uncertainty on scale factor propagated.
	TH1D *h        = (TH1D*) histo->Clone(histo->GetName());
	h->Rebin(h->GetNbinsX());
	h->SetBinError(1, sqrt(h->GetBinError(1)*  h->GetBinError(1)   *scale_fact    *scale_fact + 
			  h->GetBinContent(1)*h->GetBinContent(1) *scale_fact_err*scale_fact_err));
	h->SetBinContent(1, h->GetBinContent(1)*scale_fact);
	return h;
}