C++ TH1::GetBinContent方法代码示例

/*============================================================================*/
void gaus1peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[5],epar[5],x[4],y[4];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  setcanvas(1);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x2;
  x[2] = x3;
  x[3] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x2);
  y[1] = hist->GetBinContent(bin2);
  Int_t bin3 = hist->FindBin(x3);
  y[2] = hist->GetBinContent(bin3);
  Int_t bin4 = hist->FindBin(x4);
  y[3] = hist->GetBinContent(bin4);
  TGraph *g = new TGraph(4,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[3]=fpol1->GetParameter(0);
  par[4]=fpol1->GetParameter(1);

  //--**-- Gaussian Peak estimation without background --**--//
  TF1 *fgaus = new TF1("GAUS","gaus",x2,x3);
  hist->Fit(fgaus,"RQN");
  fgaus->GetParameters(&par[0]);

  //--**-- Final Peak Fit with Background --**--//
  TF1 *func = new TF1("FGAUS","gaus(0)+pol1(3)",x1,x4);
  func->SetParameters(par);
  hist->Fit(func,"R+QN");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  Double_t fwhm = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N0 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r0 = epar[0]/par[0];
  Double_t r2 = epar[2]/par[2];
  Double_t eN0= N0*TMath::Sqrt(r0*r0+r2*r2);
  printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
          par[1],epar[1],fwhm,efwhm,N0,eN0);
  //printf("%11.4f %11.4f %11.0f %11.0f\n",
  //        par[1],epar[1],N0,eN0);
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
}

  void getPlotData() {
    TH1 * h = (TH1*) m_file->Get(m_direc.c_str());

    for (int i=0; i<h->GetXaxis()->GetNbins(); i++) {
      m_xs.push_back(h->GetXaxis()->GetBinCenter(i));
      m_ys.push_back(h->GetBinContent(i));
    }

    m_plot->m_xAxisTitle = std::string(h->GetXaxis()->GetTitle());
    m_plot->m_yAxisTitle = std::string(h->GetYaxis()->GetTitle());
    m_plot->m_title = std::string(h->GetTitle());

    std::stringstream ssN, ssMu, ssSig, ssUF, ssOF;
    ssN << std::setprecision(4) << h->GetEntries();
    ssMu << std::setprecision(4) << h->GetMean();
    ssSig << std::setprecision(4) << h->GetRMS();
    ssUF << std::setprecision(4) << h->GetBinContent(0);
    ssOF << std::setprecision(4) << h->GetBinContent(h->GetNbinsX() + 1);

    m_statsTitles.push_back("N:");
    m_statsTitles.push_back("mu:");
    m_statsTitles.push_back("sig:");
    m_statsTitles.push_back("UF:");
    m_statsTitles.push_back("OF:");

    m_statsValues.push_back(ssN.str());
    m_statsValues.push_back(ssMu.str());
    m_statsValues.push_back(ssSig.str());
    m_statsValues.push_back(ssUF.str());
    m_statsValues.push_back(ssOF.str());
  }

// ===  FUNCTION  ============================================================
//         Name:  TPlot::GetEffHist
//  Description:  
// ===========================================================================
TH1* TPlot::GetEffHist(std::string Eff, std::string det, std::string algo)
{

  TH1* hNum = GetHist1D(listEff[Eff].first, det, algo);

  TH1* hDem = GetHist1D(listEff[Eff].second, det, algo);

  TH1* temp = (TH1*)hNum->Clone(Eff.c_str());

  assert(hNum->GetNbinsX() == hDem->GetNbinsX());
  for (int i = 0; i < hNum->GetNbinsX(); ++i)
  {
    double val = hNum->GetBinContent(i) / hDem->GetBinContent(i);
    double valerr = val * sqrt( pow(hNum->GetBinError(i)/hNum->GetBinContent(i), 2) +
        pow(hDem->GetBinError(i)/hDem->GetBinContent(i), 2) );
    if (isnan(val)) { val = 0; valerr = 0; }
    temp->SetBinContent(i, val);
    temp->SetBinError(i, valerr);
    std::cout << " bin " <<i <<" val " << val << std::endl;
  }

  //temp->Divide(hDem);
  temp->GetYaxis()->SetTitle("Efficiency");

  return temp;
}       // -----  end of function TPlot::GetEffHist  -----

void KVCanvas::ProjectionX(TH2* hh)
{
   TString pname = Form("%s_px", hh->GetName());
   Int_t ip = 1;
   while (gROOT->FindObject(pname.Data())) {
      pname = Form("%s_px%d", hh->GetName(), ip);
      ip++;
   }

   TH1* px = hh->ProjectionX(pname.Data());
   if (!px) return;
   Double_t minY = (hh->GetYaxis()->GetXmin());
   Double_t maxY = (hh->GetYaxis()->GetXmax());
   Double_t dY = (maxY - minY) * 0.8;

   Double_t maxH = px->GetBinContent(px->GetMaximumBin());

   TGraph* gg = 0;
   if ((gg = (TGraph*)gROOT->FindObject(Form("%s_gjx", hh->GetName())))) gg->Delete();

   gg = new TGraph;
   for (int i = 0; i < px->GetNbinsX(); i++) {
      gg->SetPoint(i, px->GetBinCenter(i), minY + px->GetBinContent(i)*dY / maxH);
   }

   gg->SetName(Form("%s_gjx", hh->GetName()));
   gg->SetTitle(Form("%s_gjx", hh->GetName()));
   gg->SetLineColor(kBlack);
   gg->SetMarkerColor(kBlack);
   gg->SetMarkerStyle(8);
   gg->Draw("PL");

   Modified();
   Update();
}

void setRange(RooWorkspace& myws, RooPlot* frame, string dsName, bool setLogScale, double dMuonYmin)
{ 
  // Find maximum and minimum points of Plot to rescale Y axis
  TH1* h = myws.data(dsName.c_str())->createHistogram("hist", *myws.var("ctau"), Binning(frame->GetNbinsX(),frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax()));
  Double_t YMax = h->GetBinContent(h->GetMaximumBin());
  cout << YMax << endl;
  // Double_t YMin = min( h->GetBinContent(h->FindFirstBinAbove(0.0)), h->GetBinContent(h->FindLastBinAbove(0.0)) );
  Double_t YMin = 1e99;
  for (int i=1; i<=h->GetNbinsX(); i++) if (h->GetBinContent(i)>0) YMin = min(YMin, h->GetBinContent(i));
  
  bool isMC = false;
  if (dsName.find("MC")!=std::string::npos) isMC = true;
    
  Double_t Yup(0.),Ydown(0.);
  if(setLogScale)
  {
    if (isMC) Ydown = YMin*0.3;
    else Ydown = YMin/(TMath::Power((YMax/YMin), (0.1/(1.0-0.1-0.4))));
    Yup = YMax*TMath::Power((YMax/YMin), (0.4/(1.0-0.1-0.4)));
  }
  else
  {
    Ydown = max(YMin-(YMax-YMin)*(0.1/(1.0-0.1-0.4)),0.0);
    Yup = YMax+(YMax-YMin)*(0.4/(1.0-0.1-0.4));
  }
  cout << Ydown << " " << Yup << endl;
  frame->GetYaxis()->SetRangeUser(Ydown,Yup);
  delete h;
  
  // Create line to indicate upper fitting range for MC
  if (isMC)
  {
    if (dsName.find("JPSIP")!=std::string::npos)
    {
      TLine* line(0x0);
      if (dMuonYmin >= 1.6) line = new TLine(3.32,Ydown,3.32,Yup);
      else line = new TLine(3.26,Ydown,3.26,Yup);
      line->SetLineStyle(2);
      line->SetLineColor(1);
      line->SetLineWidth(3);
      
      frame->addObject(line);
    }
    else if (dsName.find("PSI2S")!=std::string::npos)
    {
      TLine* line(0x0);
      if (dMuonYmin >= 1.6) line = new TLine(3.95,Ydown,3.95,Yup);
      else line = new TLine(3.85,Ydown,3.85,Yup);
      line->SetLineStyle(2);
      line->SetLineColor(1);
      line->SetLineWidth(3);
      
      frame->addObject(line);
    }

  }
 
}

Double_t fitgp0( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){
    cout << hs << " does not exist\n";
    return 0;
  }

  h->SetMarkerStyle(21);
  h->SetMarkerSize(0.8);
  h->SetStats(1);
  gStyle->SetOptFit(101);

  gROOT->ForceStyle();

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1 = h->GetBinCenter(1);
  double x9 = h->GetBinCenter(nb);

  // create a TF1 with the range from x1 to x9 and 4 parameters
  TF1 *gp0Fcn = new TF1( "gp0Fcn", gp0Fit, x1, x9, 4 );

  gp0Fcn->SetParName( 0, "mean" );
  gp0Fcn->SetParName( 1, "sigma" );
  gp0Fcn->SetParName( 2, "area" );
  gp0Fcn->SetParName( 3, "BG" );

  gp0Fcn->SetNpx(500);
  gp0Fcn->SetLineWidth(4);
  gp0Fcn->SetLineColor(kMagenta);
  gp0Fcn->SetLineColor(kGreen);

  // set start values for some parameters:
  gp0Fcn->SetParameter( 0, xmax ); // peak position
  gp0Fcn->SetParameter( 1, 4*dx ); // width
  gp0Fcn->SetParameter( 2, nn ); // N
  gp0Fcn->SetParameter( 3, bg );

  // N: not drawing
  // Q: quiet
  // R: use specified range
  h->Fit( "gp0Fcn", "NQR", "ep" );

  return gp0Fcn->GetParameter(1);

}

// echo object at mouse position
void exec1()
{
   //example of macro called when a pad is redrawn
   //one must create a TExec object in the following way
   // TExec ex("ex",".x exec1.C");
   // ex.Draw();
   // this macro prints the bin number and the bin content when one clicks
   //on the histogram contour of any histogram in a pad
   //Author: Rene Brun
   
   if (!gPad) {
      Error("exec1", "gPad is null, you are not supposed to run this macro");
      return;
   }
   
   int event = gPad->GetEvent();
   if (event != 11) return;
   int px = gPad->GetEventX();
   TObject *select = gPad->GetSelected();
   if (!select) return;
   if (select->InheritsFrom(TH1::Class())) {
      TH1 *h = (TH1*)select;
      Float_t xx = gPad->AbsPixeltoX(px);
      Float_t x  = gPad->PadtoX(xx);
      Int_t binx = h->GetXaxis()->FindBin(x);
      printf("event=%d, hist:%s, bin=%d, content=%f\n",event,h->GetName(),binx,h->GetBinContent(binx));
   }
}

void dominik()
{
  TH1* matHistogramRoman = static_cast<TH1*>(extractObjectFromFile("lyRoman.root", "lightYieldProjectionY")->At(0));
  TList* objects = extractObjectFromFile("c.root", "chargeBins");
  TH1* matHistogramDominik = new TH1D("matHistogramDominik", ";channel;light yield / pixels", 512, -0.5, 512-0.5);
  int sipmIt = 0;
  for (int i = 0; i < objects->GetSize(); ++i) {
    TH1* h = static_cast<TH1*>(objects->At(i));
    if (h->GetLineColor() == 8) {
      for (int bin = 1; bin <= 128; ++bin) {
        matHistogramDominik->SetBinContent(512 - (sipmIt * 128 + bin - 1), h->GetBinContent(bin));
        if (h->GetBinError(bin) > 0)
          matHistogramDominik->SetBinError(512 - (sipmIt * 128 + bin - 1), h->GetBinError(bin));
      }
      ++sipmIt;
    }
  }
  TCanvas* c = new TCanvas;
  c->Divide(1, 2);
  c->cd(1);
  matHistogramDominik->Draw();
  matHistogramRoman->Draw("SAME");
  c->cd(2);
  TH1* h = static_cast<TH1*>(matHistogramDominik->Clone());
  h->Add(matHistogramRoman, -1);
  h->Draw();
}

  void initialize( const char* countdbfn = "counts.root" )
  {
    XSWeightTool::update( "Xsection8TeV_bkg_v4.txt" );
    XSWeightTool::update( "Xsection8TeV_sig_v4.txt" );
    // XSWeightTool::print();

    using namespace std;

    TFile* f = TFile::Open( countdbfn, "read" );
    if( f == NULL || f->IsZombie() )
    {
      return;
    }

    TIter next( f->GetListOfKeys() ); 
    TKey *key;
    int nhist = 0;
    while( (key = (TKey*)next()) ) {
      nhist++;
      // cout << "Key " << nhist << endl;
      // cout << " Classname " << key->GetClassName() << endl;
      // cout << " Title " <<key->GetTitle() << endl;
      if( key->ReadObj()->InheritsFrom( TH1::Class() ) )
      {
	TH1* h = (TH1*)( key->ReadObj() );
	long run_num = TString( h->GetName() ).Atoi();
	num_events[run_num] = h->GetBinContent( 1 );
	// cout << " Run, # of events : " << run_num << ", " <<  h->GetBinContent( 21 ) << endl;
      }
    }
    // cout << "Done" << endl;
  }

void makefinal(TFile* fp, TFile *fm,  TFile* of=0) {
  setTDRStyle();
  
  TH1 *hp = (TH1*) fp->Get("MYHA");
  TH1 *hm = (TH1*) fm->Get("MYHA");
  
  TH1 *h = hp->Clone("MYNA");
  for (int k=1; k<=h->GetNbinsX(); k++) {
    std::cout << hp->GetBinContent(k) <<  " " << hm->GetBinContent(k) << std::endl;
    h->SetBinContent(k,hp->GetBinContent(k)*.5 +hm->GetBinContent(k)*.5);
    std::cout << " NEW " << h->GetBinContent(k) << std::endl;
    //h->SetBinEntries(k,1);
  }
  
  if (of!=0) { of->cd(); h->Write();}

  
}

RooHistN::RooHistN(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the asymmetry between the specified TH1 objects
  // which may have fixed or variable bin widths, but which must both have
  // the same binning. The asymmetry is calculated as (1-2)/(1+2). Error bars are
  // calculated using Binomial statistics. Prints a warning and rounds
  // any bins with non-integer contents. Use the optional parameter to
  // specify the confidence level in units of sigma to use for
  // calculating error bars. The nominal bin width specifies the
  // default used by addAsymmetryBin(), and is used to set the relative
  // normalization of bins with different widths. If not set, the
  // nominal bin width is calculated as range/nbins.

  initialize();
  // copy the first input histogram's name and title
  SetName(data1.GetName());
  SetTitle(data1.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data1).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
		     data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
  // initialize our contents from the input histogram contents
  Int_t nbin= data1.GetNbinsX();
  if(data2.GetNbinsX() != nbin) {
    coutE(InputArguments) << "RooHistN::RooHistN: histograms have different number of bins" << endl;
    return;
  }
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data1.GetBinCenter(bin);
    if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
      coutW(InputArguments) << "RooHistN::RooHistN: histograms have different centers for bin " << bin << endl;
    }
    Stat_t y1= data1.GetBinContent(bin);
    Stat_t y2= data2.GetBinContent(bin);
    addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac);
  }
  // we do not have a meaningful number of entries
  _entries= -1;
}

  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

double fittp0sigma( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){ cout << hs << " does not exist\n"; return 0; }

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1 = h->GetBinCenter(1);
  double x9 = h->GetBinCenter(nb);

  // create a TF1 with the range from x1 to x9 and 5 parameters
  TF1 *tp0Fcn = new TF1( "tp0Fcn", tp0Fit, x1, x9, 5 );

  tp0Fcn->SetParName( 0, "mean" );
  tp0Fcn->SetParName( 1, "sigma" );
  tp0Fcn->SetParName( 2, "nu" );
  tp0Fcn->SetParName( 3, "area" );
  tp0Fcn->SetParName( 4, "BG" );
   
  // set start values for some parameters:
  tp0Fcn->SetParameter( 0, xmax ); // peak position
  tp0Fcn->SetParameter( 1, 4*dx ); // width
  tp0Fcn->SetParameter( 2, 2.2 ); // nu
  tp0Fcn->SetParameter( 3, nn ); // N
  tp0Fcn->SetParameter( 4, bg );
    
  h->Fit( "tp0Fcn", "Q R", "ep" );
  // h->Fit("tp0Fcn","V+","ep");
  TF1 *fit = h->GetFunction("tp0Fcn");
  return fit->GetParameter(1);
}

void TH2CB::FillElements(const TH1 &h)
{
    if( h.GetNbinsX() != GetNumberOfElements() ) {
        cerr << "WARNING: Number of bis don't match" << endl;
    }

    const Int_t n = min((Int_t)GetNumberOfElements(), h.GetNbinsX());

    for(Int_t i=1;i<=n; ++i) {
        SetElement(i-1, GetElement(i-1) + h.GetBinContent(i));
    }
    SetBinContentChanged(kTRUE);
}

TEST_F(EvalHistMethods, CreateHistogram1D) {
  evaluator->SetNormalizationBuffer(norm);
  evaluator->SetParameterBuffer(params);

  TH1* hist = evaluator->CreateHistogram();

  EXPECT_EQ(2, hist->GetNbinsX());
  ASSERT_FLOAT_EQ(1.0, hist->Integral("width"));
  ASSERT_FLOAT_EQ(1.6, hist->GetBinContent(hist->FindBin(0.25)));
  ASSERT_FLOAT_EQ(0.4, hist->GetBinContent(hist->FindBin(0.75)));

  delete hist;
}