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

double Prediction::GetLimitedRatio(TH1D* h1, TH1D* h2, float min, float max, bool binomial, float &relerror){
	if(h1->GetNbinsX() != h2->GetNbinsX()){cout << "GetIntegral ERROR! " << endl; exit(-1);}
	TH1D* h1cp = h1->Clone("h1cp");
	TH1D* h2cp = h2->Clone("h2cp");
	for(int i=1; i<=h1cp->GetNbinsX(); ++i){
		if(h1cp->GetBinLowEdge(i) <min || h1cp->GetBinLowEdge(i)+h1cp->GetBinWidth(i) >max){
			h1cp->SetBinContent(i, 0);
			h1cp->SetBinError(i, 0);
			h2cp->SetBinContent(i, 0);
			h2cp->SetBinError(i, 0);
		}
	}
	TH1D* h1cp_2 = GetScaledHisto(h1cp, 1, 0, h1cp->GetNbinsX());
	TH1D* h2cp_2 = GetScaledHisto(h2cp, 1, 0, h1cp->GetNbinsX());

	TH1D* r      = h1cp_2->Clone("r");
	r->Divide(h1cp_2, h2cp_2, 1, 1, binomial? "B":"");

	*fLogStream << "Getting Limited Ratio for " << h1->GetName()  << " " << h2->GetName() << endl;
	*fLogStream << ">>> ratio = " << r->GetBinContent(1)  << " pm " << r->GetBinError(1) << endl;
	if(fDoRatioStatErr){
		relerror = r->GetBinError(1)/r->GetBinContent(1);
	}
	double  ratio = r->GetBinContent(1);
	delete h1cp;
	delete h2cp;
	delete h1cp_2;
	delete h2cp_2;
	delete r;
	return ratio;
}

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;
}

TH1D *
GetITSTPCSpectrum(TFile *file, Int_t part, Int_t charge, Int_t cent)
{
  TList *list = (TList *)file->Get("output");
  TH1D *hin = (TH1D *)list->FindObject(Form("h_%s_%s_cen_%d", ITSTPCPartName[part], ITSTPCChargeName[charge], cent + 1));
  if (!hin) return NULL;

  TH1D *h = new TH1D(Form("hITSTPC_cent%d_%s_%s", cent, AliPID::ParticleName(part), chargeName[charge]), "ITSTPC", NptBins, ptBin);
  Double_t pt, width, value, error;
  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;
    /* copy bin */
    value = hin->GetBinContent(bin);
    error = hin->GetBinError(bin);
    h->SetBinContent(ipt + 1, value);
    h->SetBinError(ipt + 1, error);
  }
  
#if 0
  /* add systematic error */
  Double_t sys;
  if (part == 2) sys = 0.5;
  else sys = 0.1;
  Double_t cont, conte;
  for (Int_t ipt = 0; ipt < h->GetNbinsX(); ipt++) {
    cont = h->GetBinContent(ipt + 1);
    conte = h->GetBinError(ipt + 1);
    conte = TMath::Sqrt(conte * conte + sys * sys * cont * cont);
    h->SetBinError(ipt + 1, conte);
  }
#endif
  
  h->SetTitle("ITSTPC");
  h->SetLineWidth(1);
  h->SetLineColor(1);
  h->SetMarkerStyle(21);
  h->SetMarkerColor(2);
  h->SetFillStyle(0);
  h->SetFillColor(0);
  return h;
}

void MergeOverflow(TH1D &h, bool merge_underflow, bool merge_overflow) {
    if(merge_underflow) {
        h.SetBinContent(1, h.GetBinContent(0)+h.GetBinContent(1));
        h.SetBinContent(0, 0.);
        h.SetBinError(1, hypot(h.GetBinError(0), h.GetBinError(1)));
        h.SetBinError(0, 0.);
    }
    int nbins = h.GetNbinsX();
    if(merge_overflow) {
        h.SetBinContent(nbins, h.GetBinContent(nbins)+h.GetBinContent(nbins+1));
        h.SetBinContent(nbins+1, 0.);
        h.SetBinError(nbins, hypot(h.GetBinError(nbins), h.GetBinError(nbins+1)));
        h.SetBinError(nbins+1, 0.);
    }
}

TH1D* ratio_hist_to_func(TH1D* num, double par0, double par1, double par2) {
    cout<<"[Ratio to fit used]"<<endl;

    TH1D *hRatio = (TH1D*) num->Clone("hRatio");

    TF1 *f3 = new TF1("f3","[0]*(1+(sqrt(0.1396**2+x**2)-0.1396)/([1]*[2]))**(-[2])",0,6.5);
    f3->SetParameters(par0,par1,par2);
    //f3->SetLineColor(2);

    int nbin = num->GetEntries();
    for(int i=0; i<nbin; i++) {

        double cms_value = (double) f3->Eval(hRatio->GetBinCenter(i+1));

        if(hRatio->GetBinCenter(i+1)>0.4 && hRatio->GetBinCenter(i+1)<6.0) {
            //double ratio = cms_value/hRatio->GetBinContent(i+1);
            //double ratio_err = cms_value/hRatio->GetBinError(i+1);
            double ratio = hRatio->GetBinContent(i+1)/cms_value;
            double ratio_err = hRatio->GetBinError(i+1)/cms_value;
        } else {
            double ratio = -999;
            double ratio_err = 0.0;
        }
        //double ratio = hRatio->GetBinContent(i+1)/cms_value;
        //double ratio_err = hRatio->GetBinError(i+1)/cms_value;
        hRatio->SetBinContent(i+1,ratio);
        hRatio->SetBinError(i+1,ratio_err);
    }

    return hRatio;

}

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 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();
}

CombineSpectra(TH1D *h, TObjArray *oa)
{

  /* combine with weighted mean */
  Double_t ptcen, w, sumw, mean, meane, cont, conte, syse, tote, minvalue, maxvalue, maxerror, maxmindiff;
  Int_t ptbin;
  TH1D *hin;
  for (Int_t ipt = 0; ipt < h->GetNbinsX(); ipt++) {
    mean = 0.;
    sumw = 0.;
    minvalue = kMaxInt;
    maxvalue = 0.;
    maxerror = 0.;
    ptcen = h->GetBinCenter(ipt + 1);
    TProfile prof(Form("prof_%d", ipt), "", 1, 0., 1.); /* to get RMS */
    for (Int_t ih = 0; ih < oa->GetEntries(); ih++) {
      hin = (TH1D *)oa->At(ih);
      ptbin = hin->FindBin(ptcen);
      if (ptbin <= 0.) continue;
      cont = hin->GetBinContent(ptbin);
      conte = hin->GetBinError(ptbin);
      if (cont == 0. || conte == 0.)
	continue;
      w = 1. / (conte * conte);
      mean += cont * w;
      sumw += w;
      if (cont < minvalue)
	minvalue = cont;
      if (cont > maxvalue)
	maxvalue = cont;
      if (conte > maxerror)
	maxerror = conte;
      prof.Fill(0., cont);
    }
    if (sumw == 0.) {
      mean = 0.;
      meane = 0.;	
      syse = 0.;
      tote = 0.;
    }
    else {
      mean /= sumw;
      meane = TMath::Sqrt(1. / sumw);
      syse = 0.5 * (maxvalue -  minvalue);
      tote = TMath::Sqrt(meane * meane + syse * syse);
      //      if (tote < maxerror)
      //	tote = maxerror;
    }

    //    printf("pt = %f, mean = %f, meane = %f, syse = %f, tote = %f\n", ptcen, mean, meane, syse, tote);
    
    //    syse = prof.GetBinError(1);
    h->SetBinContent(ipt + 1, mean);
    h->SetBinError(ipt + 1, tote);
    //    h->SetBinError(ipt + 1, TMath::Sqrt(meane * meane + syse * syse + mean * mean * 0.1 * 0.1)); /* add extra 10% systematics */
  }

}

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();
  
}

UpsilonYields getValFromTH1D(TString kineLabel) {
  
  UpsilonYields yields;
  TString fname = Form("fitResults/fixParam1MuPt4_Feb02/fitresults_upsilon_singleCB_%s.root",kineLabel.Data());
  cout << "opening " << fname << " ....." <<endl;
  TFile* f1 = new TFile(fname.Data() );
  if ( f1->IsZombie() )  {
    cout << " no such a file exists..." << endl << endl;
    yields.yield1s = -1 ;     yields.err1s = -1 ; 
    yields.yield2s = -1 ;     yields.err2s = -1 ; 
    yields.yield3s = -1 ;     yields.err3s = -1 ; 
  }
  else { 
    TH1D* h = (TH1D*)f1->Get("fitResults");
    
    yields.yield1s =    h->GetBinContent(1);
    yields.err1s = h->GetBinError (1);
    yields.yield2s =    h->GetBinContent(2);
    yields.err2s = h->GetBinError (2);
    yields.yield3s =    h->GetBinContent(3);
    yields.err3s = h->GetBinError (3);
    yields.r23o123 = h->GetBinContent(4);
    yields.err23o123 = h->GetBinError(4);
    yields.r23o1 = h->GetBinContent(5);
    yields.err23o1 = h->GetBinError(5);

    cout << " 1s = " << yields.yield1s << " +/- " << yields.err1s << endl;
    cout << " 2s = " << yields.yield2s << " +/- " << yields.err2s << endl;
    cout << " 3s = " << yields.yield3s << " +/- " << yields.err3s << endl;
    cout << " (2S+3S)/(1S+2S+3S) = " << yields.r23o123 << " +/- " <<  yields.err23o123 << endl;
    cout << " (2S+3S)/1S = " << yields.r23o1 << " +/- " <<  yields.err23o1 << endl;
    f1->Close();
  }
  cout << "*===*===*===*===*===*===*===*===*===*===*===*===*===*===*"<<endl;
  return yields;
    
}

void Yield1S_1D(){
  
  TFile *f = new TFile("SigmaDATA_RapInt_3SBin.root");
  TH1D *H;
  H = (TH1D*)gFile->Get("hYield1S");
  double tot(0.); double totE(0.); 
  
  ofstream OUT("Yield_1S_1D.tex");
  OUT << "% ----------------------------------------------------------------------" << endl;
  OUT << "% -- Yields" << endl;
  for ( int x = 1; x <= H->GetNbinsX(); ++x ){
    OUT << Form("\\vdef{%iS_bin%iLowEdge} {\\ensuremath{ {%.2f } } }",1, x, H->GetXaxis()->GetBinLowEdge(x) ) << endl;
    if ( x == H->GetNbinsX() ) OUT << Form("\\vdef{%iS_bin%iHighEdge} {\\ensuremath{ {%.2f } } }",1 , x, H->GetXaxis()->GetBinUpEdge(x) ) << endl;
    OUT << Form("\\vdef{%iS_bin%iContent} {\\ensuremath{ {%.0f } } }",1, x, H->GetBinContent(x) ) << endl;
    OUT << Form("\\vdef{%iS_bin%iError} {\\ensuremath{ {%.0f } } }",1, x, H->GetBinError(x) ) << endl;
    tot += H->GetBinContent(x);
    totE += (H->GetBinError(x)*H->GetBinError(x));
  }
		
  OUT << Form("\\vdef{sum%iS} {\\ensuremath{ {%.0f } } }",1, tot) << endl;
  OUT << Form("\\vdef{sum%iSError} {\\ensuremath{ {%.0f } } }",1 , TMath::Sqrt(totE) ) << endl;
  OUT.close();
  
}

double countDataError(TChain* chain, TCut cut)
{
  // Create temporary histogram
  TH1D* tmp = new TH1D("tmp", "tmp", 30, 0, 1000);

  double DataError = 0;
  tmp->Sumw2();
  chain->Draw("Njet>>tmp", cut, "goff");

  for(int i=0; i<=(tmp->GetNbinsX()+1); i++)
    DataError += tmp->GetBinError(i);

  delete tmp;
  return DataError;
}

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);
}

/** 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 ;
}

// 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