C++ TH1F::AddBinContent方法代码示例

void dNdypi0_auau(double c, double p0, double a, double b, double n)
{
  float pt, mt;
  float weight = 1.0;
  int nbins = 10000;
  float ptmin = 0.0;
  float ptmax = 10.;
  float binwidth = (ptmax-ptmin)/(double)nbins;
  TH1F * pthist = new TH1F("pthist","pthist",nbins,ptmin,ptmax);

  double pi    = acos(-1.0);

  double integral = 0.0;

  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    pt = ptmin+(double)(ibin-1)*binwidth+binwidth/2.0;

    // minimum bias
    weight = 2.0*pi*pt*c/pow(exp(-a*pt-b*pt*pt)+pt/p0,n);

    pthist->AddBinContent(ibin,weight);
    integral = integral + weight*binwidth;
  }

  cout << "Integral (dN/dy): " << integral << endl;
 
}

TH1F * InitializeWhiteY(double YMin, double YMax)
{
  int    nbins, ibin;
  float  binwidth;
  float  weight;
  
  nbins = 1000;
  binwidth = (YMax-YMin)/(double)nbins;
  TH1F * yhistogram = new TH1F("y","y",nbins,YMin,YMax);

  for ( ibin=1; ibin<=nbins; ibin++ )
  {
    weight = 1.0;
    yhistogram->AddBinContent(ibin,weight);
  }
  
  yhistogram = InitializeRandomHist(yhistogram);

  return yhistogram;

}

//.........这里部分代码省略.........
    // double pipower    = 51.28/pow(pt,8.21);
    // double otherhagedorn = 2.99e+06 / 
    //   pow(2.00+sqrt(mtother*mtother-pimass*pimass),12.45);
    // double otherpower = 51.28/pow(sqrt(mtother*mtother-pimass*pimass),8.21);
    // 40 - 60 % central
    // double trans      = 1./(1.+exp(5.*pt-17.5));
    // double pihagedorn = 3.00e+04/pow(1.56+pt,10.76);
    // double pipower    = 15.75/pow(pt,7.97);
    // double otherhagedorn = 3.00e+04 / 
    //   pow(1.56+sqrt(mtother*mtother-pimass*pimass),10.76);
    // double otherpower = 15.75/pow(sqrt(mtother*mtother-pimass*pimass),7.97);
    // 60 - 92 % central
    // double trans      = 1./(1.+exp(5.*pt-17.5));
    // double pihagedorn = 8.73e+02/pow(1.29+pt,9.94);
    // double pipower    = 5.28/pow(pt,8.28);
    // double otherhagedorn = 8.73e+02 / 
    //   pow(1.29+sqrt(mtother*mtother-pimass*pimass),9.94);
    // double otherpower = 5.28/pow(sqrt(mtother*mtother-pimass*pimass),8.28);
    // if ( pt<1.5 ) weight = 2.0 * pi * pt * pihagedorn;
    // if ( pt>=1.5 && pt<=5.5 ) 
    //   weight = 2.0 * pi * pt * (trans*pihagedorn + (1.-trans)*pipower);
    // if ( pt>5.5 ) weight = 2.0*pi*pt*pipower;

    // pp (Xinhua's Run-2 parameterization)
    // weight = 2.0*pi*pt * (485.0/42.) / pow(exp(-0.5187*pt) + pt/0.6631,8.163);
    // pp (Run-2 hard parameterization)
    weight = 2.0*pi*pt * (450.8/42.) / pow(exp(-0.5553*pt) + pt/0.6467,8.046);
    // pp (Run-3 parameterization)
    // weight = 2.0*pi*pt * (413.5/42.2) / pow(exp(-0.2614*pt) + pt/0.8144,8.807);
    // pp (average Run-2/3 parameterization)
    // weight = 2.0*pi*pt * (546.2/42.2) / pow(exp(-0.2664*pt) + pt/0.7703,8.718);
    //    weight = weight*(1+0.6362/exp(7.736*pt));

    ptpion->AddBinContent(ibin,weight);
    intpion = intpion + weight*binwidth;

    // if ( pt<1.5 ) weight = 2.0 * pi * pt * otherhagedorn;
    // if ( pt>=1.5 && pt<=5.5 ) 
    //   weight = 2.0 * pi * pt * (trans*otherhagedorn + (1.-trans)*otherpower);
    // if ( pt>5.5 ) weight = 2.0*pi*pt*otherpower;

    weight = 2.0*pi*pt * (450.8/42.2) / 
      pow(exp(-0.5553*sqrt(mtother*mtother-pimass*pimass)) + 
	  sqrt(mtother*mtother-pimass*pimass)/0.6467,8.046);

    ptothermt->AddBinContent(ibin,weight);
    intothermt = intothermt + weight*binwidth;

    weight = norm*InitializePtHydro(omass,t_fo,beta,pt);
    ptother->AddBinContent(ibin,weight);
    intother = intother + weight*binwidth;
  }

  cout << "Integral pion    (dN/dy): " << intpion << endl;
  cout << "Integral other   (dN/dy): " << intother << endl;
  cout << "Integral othermt (dN/dy): " << intothermt << endl;

  norm = intpion/intothermt;
  ptothermt->Scale(norm);
 
  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    weight = ptothermt->GetBinContent(ibin)/ptpion->GetBinContent(ibin);
    ratiomt->AddBinContent(ibin,weight);
  }

TH1F * InitializeYCERES(int ParticleID, ParticlePropertyList *PPList)
{
  int    nbins, ibin;
  float  binwidth;
  float  y, ymin, ymax, weight;
  double mass;
  double sqrts, y0;
  double gamma, sigma_landau, sig, sigma;
  double emax_pi0, pmax_pi0, ymax_pi0;
  double emax_part, pmax_part, ymax_part, ymin_part;
  double ymean;
  
  double mnuc  = 0.938;
  double ebeam = 450.0;
  double Aproj = 1.0;
  double Atarg = 9.0;
  double bpar  = 0.0;
  double m_pi0 = 0.1349743;

  sqrts        = sqrt(2.0*mnuc*ebeam+2.0*mnuc*mnuc);
  y0           = log(2.0*ebeam/mnuc)/2.0;
  gamma        = sqrts/(2.0*mnuc);
  sigma_landau = sqrt(log(gamma));
  emax_pi0     = (sqrts*sqrts-4.0*mnuc*mnuc+m_pi0*m_pi0)/(2.0*sqrts);
  pmax_pi0     = sqrt(emax_pi0*emax_pi0-m_pi0*m_pi0);
  ymax_pi0     = log(sqrts/m_pi0);

  ParticleProperty * PParticle = PPList->GetByID(ParticleID);
  mass      = PParticle->GetMass();
  emax_part = (sqrts*sqrts-4.0*mnuc*mnuc+mass*mass)/(2.0*sqrts);
  pmax_part = sqrt(emax_pi0*emax_pi0-mass*mass);
  ymax_part = log(sqrts/mass);
  ymin_part = -1.0*ymax_part;
  sig       = sigma_landau*(ymax_part/ymax_pi0);

  ymax_part = ymax_part*y0;
  ymin_part = ymin_part*y0;

  ymean = y0-y_shift(ebeam,Aproj,Atarg,bpar);

  nbins = 1000;
  ymin  = -2.0;
  ymax  = 8.0;
  binwidth = (ymax-ymin)/(double)nbins;
  TH1F * yhistogram = new TH1F("y","y",nbins,ymin,ymax);

  for ( ibin=1; ibin<=nbins; ibin++ )
  {
    y  = ymin+(double)(ibin-1)*binwidth+binwidth/2.0;
    if ( y<=ymean )
    {
      sigma = sig*(1.0-((y0-ymean)*(y0-ymean)/(ymax_part-ymin_part)));
    }
    else
    {
      sigma = sig*(1.0+((y0-ymean)*(y0-ymean)/(ymax_part-ymin_part)));
    }
    weight = exp(0.75*(y-ymean)/sigma)+exp(-0.75*(y-ymean)/sigma);
    weight = 1.0/(weight*weight);
    yhistogram->AddBinContent(ibin,weight);
  }
  
  yhistogram = InitializeRandomHist(yhistogram);

  return yhistogram;

}

void AddIAVCorrection(){

  gStyle->SetOptStat(0);
  Double_t xval;
  Double_t yval;
  double junk1, junk2, junk3, junk4, junk5, junk6;

  //Read in reactor flux
  TString fullname;
  fullname = Form("./InputIBD.dat");
  cout<<fullname<<endl;
  ifstream datafile;

  ////Read in IAV Correction
  TFile *file = new TFile("./IAVDistortion.root");
  TH1F *Correciton = file->Get("Correction");
  TH1F *Uncorr = file->Get("UncUncorr");
  
  ////Make histogram of reactor data
  TH1F *Reactor = new TH1F("Reactor","Reactor",240,0,12);
  datafile.open(fullname, ios::in);
  for (Int_t j=0; j<240 ; j++) {
    datafile >> junk1 >> junk2 >> junk3 >> yval >> xval;
    //cout<<xval+0.01<<" "<<yval<<endl;
    Reactor->Fill(xval+0.01,yval);
  }
  datafile.close();

  //////////////////////////////////////////////////////////
  ///////////Add IAV Correction to Histogram: IMPORTANT PART
  //////////////////////////////////////////////////////////
  double CorrectedPercent=0;
  double CorrectedPercentUnc=0;
  double TotalVal=0;
  double StatUnc=0;
  TH1F *ReactorCorr = new TH1F("ReactorCorr","ReactorCorr",240,0,12);
  TH1F *Corr = new TH1F("Corr","Corr",240,0,12);
  for(int i=21;i<241;i++){
    //Total events in input spectrum bin i
    TotalVal = Reactor->GetBinContent(i);
    for(int j=0;j<i+1;j++){
      //Percent of events in bin i moved to bin j
      CorrectedPercent = Correction->GetBinContent(i,j)/Correction->Integral(i,i,0,-1);
      CorrectedPercentUnc = TMath::Sqrt(Correction->GetBinContent(i,j))/Correction->Integral(i,i,0,-1);
      //Add percent of total of input spectrum bin i to corrected spectrum bin j
      ReactorCorr->AddBinContent(j,TotalVal*CorrectedPercent);
      
      //Add statistical uncertainty to bin j
      //Get Existing Uncertainty in bin j
      StatUnc = TMath::Sqrt(TMath::Power(ReactorCorr->GetBinError(j),2)+TMath::Power(TotalVal*CorrectedPercentUnc,2));      
      ReactorCorr->SetBinError(j,StatUnc);
    }    
  }
  
  //Make correlated uncertainty histogram
  for(int i=0;i<241;i++)
    Corr->SetBinContent(i,ReactorCorr->GetBinError(i)/ReactorCorr->GetBinContent(i));
  

  //Draw Histograms
  TCanvas *c1 = new TCanvas("c1","c1",0,0,600,500);
  TCanvas *c2 = new TCanvas("c2","c2",0,0,600,500);
  TCanvas *c3 = new TCanvas("c3","c3",0,0,600,500);
  TLegend *leg = new TLegend(0.5,0.5,0.8,0.7);
  c1->cd();
  c1->SetLogy();
  //ReactorCorr->Divide(Reactor);
  //for(int i=0;i<241;i++)
  //  ReactorCorr->SetBinError(i,0);
  Reactor->Draw("");  
  ReactorCorr->SetLineColor(kRed);
  ReactorCorr->Draw("Same");
  leg->AddEntry(Reactor,"Input IBD Spectrum","EL");
  leg->AddEntry(ReactorCorr,"After IAV Correction","EL");
  leg->Draw();
  c2->cd();
  Corr->Draw("");
  c3->cd();
  Uncorr->Draw("");
  
  //cout<<Reactor->Integral(0,-1)<<endl;
  //cout<<ReactorCorr->Integral(0,-1)<<endl;
  //cout<<(Reactor->Integral(0,-1)-ReactorCorr->Integral(0,-1))/Reactor->Integral(0,-1)<<endl;
}