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

// get FWHHM
vector<double> getFWHM(RooRealVar *mass, RooAbsPdf *pdf, RooDataSet *data, double wmin=110., double wmax=130., double step=0.0004) {
 
  cout << "Computing FWHM...." << endl;
  double nbins = (wmax-wmin)/step;
  TH1F *h = new TH1F("h","h",int(floor(nbins+0.5)),wmin,wmax);
  if (data){
    pdf->fillHistogram(h,RooArgList(*mass),data->sumEntries());
  }
  else {
    pdf->fillHistogram(h,RooArgList(*mass));
  }
  
  double hm = h->GetMaximum()*0.5;
  double low = h->GetBinCenter(h->FindFirstBinAbove(hm));
  double high = h->GetBinCenter(h->FindLastBinAbove(hm));

  cout << "FWHM: [" << low << "-" << high << "] Max = " << hm << endl;
  vector<double> result;
  result.push_back(low);
  result.push_back(high);
  result.push_back(hm);
  result.push_back(h->GetBinWidth(1));
  
  delete h;
  return result;
}

TH1F * DrawOverflow(TH1F *h)
{
      // This function paint the histogram h with an extra bin for overflows
   UInt_t nx    = h->GetNbinsX()+1;
   Double_t *xbins= new Double_t[nx+1];
   for (UInt_t i=0;i<nx;i++)
     xbins[i]=h->GetBinLowEdge(i+1);
   xbins[nx]=xbins[nx-1]+h->GetBinWidth(nx);
   char *tempName= new char[strlen(h->GetName())+10];
   sprintf(tempName,"%swtOverFlow",h->GetName());
   // Book a temporary histogram having ab extra bin for overflows
   TH1F *htmp = new TH1F(tempName, h->GetTitle(), nx, xbins);
   // Reset the axis labels
   htmp->SetXTitle(h->GetXaxis()->GetTitle());
   htmp->SetYTitle(h->GetYaxis()->GetTitle());
   // Fill the new hitogram including the extra bin for overflows
   for (UInt_t i=1; i<=nx; i++)
     htmp->Fill(htmp->GetBinCenter(i), h->GetBinContent(i));
   // Fill the underflows
   htmp->Fill(h->GetBinLowEdge(1)-1, h->GetBinContent(0));
   // Restore the number of entries
   htmp->SetEntries(h->GetEntries());
   // FillStyle and color
   htmp->SetFillStyle(h->GetFillStyle());
   htmp->SetFillColor(h->GetFillColor());
   return htmp;
}

void Difference(TH1* iH0,TH1 *iH1) {
  std::string lName = std::string(iH0->GetName());
  //TH1F *lHDiff  = new TH1F((lName+"Diff").c_str(),(lName+"Diff").c_str(),50,0,300); lHDiff->Sumw2();
  TH1F *lHDiff  = new TH1F((lName+"Diff").c_str(),(lName+"Diff").c_str(),iH0->GetNbinsX(),iH0->GetXaxis()->GetXmin(),iH0->GetXaxis()->GetXmax()); lHDiff->Sumw2();
  lHDiff->SetFillColor(kViolet); lHDiff->SetFillStyle(1001); lHDiff->SetLineWidth(1);
  TH1F *lXHDiff1 = new TH1F((lName+"XDiff1").c_str(),(lName+"XDiff1").c_str(),iH0->GetNbinsX(),iH0->GetXaxis()->GetXmin(),iH0->GetXaxis()->GetXmax());
  int i1 = 0;
  lXHDiff1->SetLineStyle(2); lXHDiff1->SetLineWidth(2); lXHDiff1->SetLineColor(kRed);

  lHDiff->GetYaxis()->SetRangeUser(0,2);
  lHDiff->GetYaxis()->SetTitleOffset(0.6);
  lHDiff->GetYaxis()->SetTitleSize(0.08);
  lHDiff->GetYaxis()->SetLabelSize(0.08);
  lHDiff->GetYaxis()->CenterTitle();
  lHDiff->GetXaxis()->SetTitleOffset(1.2);
  lHDiff->GetXaxis()->SetTitleSize(0.10);
  lHDiff->GetXaxis()->SetLabelSize(0.08);
  lHDiff->GetXaxis()->SetTitle("#slash{E}_{T} [GeV]");
  //lHDiff->GetXaxis()->CenterTitle();
  //lHDiff->GetYaxis()->SetTitle(YLabel);

  for(int i0 = 0; i0 < lHDiff->GetNbinsX()+1; i0++) {
    double lXCenter = lHDiff->GetBinCenter(i0);
    double lXVal     = iH0   ->GetBinContent(i0);
    lXHDiff1->SetBinContent(i0, 1.0);
    while(iH1->GetBinCenter(i1) < lXCenter) {i1++;}
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinContent(i0,(lXVal-iH1->GetBinContent(i0))/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinError  (i0,iH0->GetBinError(i0)/(iH1->GetBinContent(i0)));
  }
  lHDiff->SetMarkerStyle(kFullCircle); lHDiff->SetLineColor(kBlack); lHDiff->SetMarkerColor(kBlack);
  lHDiff->Draw("E1");
  lXHDiff1->Draw("hist sames");
}

TH1* AddOverflow(TH1* h) {

  ++overflowCounter;

  TString  name = h->GetName();
  Int_t    nx   = h->GetNbinsX()+1;
  Double_t bw   = h->GetBinWidth(nx);
  Double_t x1   = h->GetBinLowEdge(1);
  Double_t x2   = h->GetBinLowEdge(nx) + bw;
  
  // Book a new histogram having an extra bin for overflows
  TH1F* htmp = new TH1F(Form(name + "_overflow_%d", overflowCounter), "", nx, x1, x2);

  // Fill the new histogram including the extra bin for overflows
  for (Int_t i=1; i<=nx; i++) {
    htmp->Fill(htmp->GetBinCenter(i), h->GetBinContent(i));
    htmp->SetBinError(i, h->GetBinError(i));
  }

  // Fill the underflow
  htmp->Fill(x1-1, h->GetBinContent(0));

  // Restore the number of entries
  htmp->SetEntries(h->GetEntries());

  // Cosmetics
  htmp->SetLineColor(h->GetLineColor());
  htmp->SetLineWidth(h->GetLineWidth());
  htmp->GetXaxis()->SetTitleOffset(1.5);

  return htmp;
}

void PaintOverflow(TH1 *h)
{
  // This function paint the histogram h with an extra bin for overflows
  
  char* name  = h->GetName();
  char* title = h->GetTitle();
  Int_t nx    = h->GetNbinsX()+1;
  Double_t x1 = h->GetBinLowEdge(1);
  Double_t bw = h->GetBinWidth(nx);
  Double_t x2 = h->GetBinLowEdge(nx)+bw;
  
  // Book a temporary histogram having ab extra bin for overflows
  TH1F *htmp = new TH1F(name, title, nx, x1, x2);
  
  // Fill the new hitogram including the extra bin for overflows
  for (Int_t i=1; i<=nx; i++) {
    htmp->Fill(htmp->GetBinCenter(i), h->GetBinContent(i));
  }
  
  // Fill the underflows
  htmp->Fill(x1-1, h->GetBinContent(0));
  
  // Restore the number of entries
  htmp->SetEntries(h->GetEntries());
  
  // Draw the temporary histogram
  htmp->Draw();
  /*
  TText *t = new TText(x2-bw/2,h->GetBinContent(nx),"Overflow");
  t->SetTextAngle(90);
  t->SetTextAlign(12);
  t->SetTextSize(0.03);;
  t->Draw();
  */
}

//Difference plotting
void drawDifference(TH1* iH0,TH1 *iH1,TH1 *iHH=0,TH1 *iHL=0,TH1 *iHH1=0,TH1 *iHL1=0) {
  std::string lName = std::string(iH0->GetName());
  TH1F *lHDiff   = (TH1F*) iH0->Clone("Diff");
  TH1F *lHDiffH  = (TH1F*) iH0->Clone("DiffH");
  TH1F *lHDiffL  = (TH1F*) iH0->Clone("DiffL"); 
  TH1F *lHDiffH1 = (TH1F*) iH0->Clone("DiffH1");
  TH1F *lHDiffL1 = (TH1F*) iH0->Clone("DiffL1"); 
  lHDiff  ->SetFillColor(kViolet); lHDiff->SetFillStyle(1001); lHDiff->SetLineWidth(1);
  lHDiffL ->SetLineWidth(1); lHDiffL ->SetLineColor(iHL ->GetLineColor());
  lHDiffH ->SetLineWidth(1); lHDiffH ->SetLineColor(iHH ->GetLineColor());
  lHDiffL1->SetLineWidth(1); lHDiffL1->SetLineColor(iHL1->GetLineColor());
  lHDiffH1->SetLineWidth(1); lHDiffH1->SetLineColor(iHH1->GetLineColor());
  TH1F *lXHDiff1 = new TH1F((lName+"XDiff1").c_str(),(lName+"XDiff1").c_str(),iH0->GetNbinsX(),iH0->GetXaxis()->GetXmin(),iH0->GetXaxis()->GetXmax());
  TH1F *lXHDiff2 = new TH1F((lName+"XDiff2").c_str(),(lName+"XDiff2").c_str(),iH0->GetNbinsX(),iH0->GetXaxis()->GetXmin(),iH0->GetXaxis()->GetXmax());
  int i1 = 0;
  lXHDiff1->SetLineWidth(2); lXHDiff1->SetLineColor(kRed);
  lXHDiff2->SetLineWidth(2); lXHDiff2->SetLineColor(kRed);

  lXHDiff1->GetYaxis()->SetTitle("Ratio");
  lXHDiff1->GetYaxis()->SetRangeUser(0.2,1.8);
  lXHDiff1->GetYaxis()->SetTitleOffset(0.4);
  lXHDiff1->GetYaxis()->SetTitleSize(0.2);
  lXHDiff1->GetYaxis()->SetLabelSize(0.11);
  for(int i0 = 0; i0 < lHDiff->GetNbinsX()+1; i0++) {
    double lXCenter = lHDiff->GetBinCenter(i0);
    double lXVal     = iH0   ->GetBinContent(i0);
    double lXValH    = iHH   ->GetBinContent(i0);
    double lXValL    = iHL   ->GetBinContent(i0);
    double lXValH1   = iHH1  ->GetBinContent(i0);
    double lXValL1   = iHL1  ->GetBinContent(i0);
    lXHDiff1->SetBinContent(i0, 1.0);
    lXHDiff2->SetBinContent(i0, 1.0);
    while(iH1->GetBinCenter(i1) < lXCenter) {i1++;}
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinContent(i0,lXVal      /(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinError  (i0,sqrt(lXVal)/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiffL->SetBinContent(i0,lXValL/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiffH->SetBinContent(i0,lXValH/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiffL1->SetBinContent(i0,lXValL1/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiffH1->SetBinContent(i0,lXValH1/(iH1->GetBinContent(i0)));
   //if(iH1->GetBinContent(i0) > 0)  cout << "unc" << lXVal << " -- " << sqrt(lXVal)/(iH1->GetBinContent(i0)) << endl;
  }
  lHDiff->SetMarkerStyle(kFullCircle);
  //lHDiff->Draw("EP");
  
  lXHDiff1->SetStats(0);
  lXHDiff2->SetStats(0);
  lHDiff->SetStats(0);
  lHDiffH->SetStats(0);
  lHDiffL->SetStats(0);
  lHDiffH1->SetStats(0);
  lHDiffL1->SetStats(0);

  lXHDiff1->Draw("hist");
  lXHDiff2->Draw("hist sames");
  lHDiff->Draw("EP sames");
  lHDiffH ->Draw("hist sames");
  lHDiffL ->Draw("hist sames");
  lHDiffH1->Draw("hist sames");
  lHDiffL1->Draw("hist sames");
}

void FitterCBForSignificance(){



  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_gg_TuneD6T_Emine2013.root");

  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_qg_TuneD6T_Emine2013.root");



  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_Qstar_2012_D6T_ak5_fat30_save.root");
  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_RSGraviton_2012_D6T_ak5_GGtoGG_fat30_save.root");
  //TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_RSGraviton_2012_D6T_ak5_QQtoQQ_fat30_save.root");

  //TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_Qstar_ak5_fat30.root");
  for (int i = 1; i < 2; i++){
	   TH1F* shape = (TH1F*) _shapes->Get(Form("h_qg_%d;1", 1000+i*100));
	   TH1F* shape_largeBinning = new TH1F("shape_large_binning", "", NBINS, BOUNDARIES);

	   for (int j = 1; j < shape->GetNbinsX()+1; j++){
	     double binCenter = shape->GetBinCenter(j);
	     double weight = shape->GetBinContent(j);
	     double binWidth = shape->GetBinWidth(j);
	     shape_largeBinning->Fill(binCenter, weight);
	     //	     cout << "j = " << j << " binWidth = " << binWidth << " weight = " << weight << endl;
	     shape->SetBinContent(j, weight/binWidth);
	   }

	   for (int j = 1; j < shape_largeBinning->GetNbinsX()+1; j++){
	     double weight = shape_largeBinning->GetBinContent(j);
	     double binWidth = shape_largeBinning->GetBinWidth(j);
	     shape_largeBinning->SetBinContent(j, weight/binWidth*1000000);
	   }

	   shape->Draw();

	   shape_largeBinning->SetLineColor(kRed);
	   shape_largeBinning->Draw("SAME");

	   //   
	   FitHistWithCBShape(shape_largeBinning, 1000.+i*100);		
	}

	  cout << "mass\tsMean\tsSigma\tsAlphaHigh\tsAlphaLow\tsNHigh\tsNLow\tsFrac"<<endl;

	for (int i = 0; i < 40; i++){
	
	  cout << Masses[i] << "\t" << Means[i] << "\t" << Sigmas[i] << "\t"
	       << alphaHights[i] << "\t\t" << alphaLows[i] << "\t"
	       << nHighs[i] << "\t\t" << nLows[i] << "\t" << fracs[i] << endl;

	}
}

bool isgood(TH1F h, float sigma){
  if(sigma == -1) return true;
  bool Flag = true;
  int imaxbin = h.GetMaximumBin();
  float bincenter = h.GetBinCenter(imaxbin);
  for(int ibin=0;ibin<h.GetNbinsX();ibin++){
      if((h.GetBinCenter(ibin)<(bincenter-sigma)) || (h.GetBinCenter(ibin)>(bincenter+sigma))){
          if(h.GetBinContent(ibin)!=0){
              Flag = false;
              break;
          }
      }
  }

  if(Flag){
  h.SetBinContent(imaxbin,0);
  imaxbin = h.GetMaximumBin();
  bincenter = h.GetBinCenter(imaxbin);
  for(int ibin=0;ibin<h.GetNbinsX();ibin++){
      if((h.GetBinCenter(ibin)<(bincenter-sigma)) || (h.GetBinCenter(ibin)>(bincenter+sigma))){
          if(h.GetBinContent(ibin)!=0){
              Flag = false;
              break;
          }
      }
  }
  }
  return Flag;
}

/**
 * Saving Histogram Data
 */
void ExportHistogram(TFile* f,const char* histKey,const char* outputfile){
  // Root Setup
  TH1F* h = (TH1F*) f->Get(histKey);

  // File Setup
  ofstream out;
  out.open(outputfile);
  out<<"Bin,Value,Error"<<endl;
  for (int i = 1; i < h->GetNbinsX()-1; i++){
    out<<h->GetBinCenter(i)<<","<<h->GetBinContent(i)<<","<<h->GetBinError(i)<<endl;
  }
  out.close();
}

void drawPileup(const std::string& fileName)
{
  TFile* f = TFile::Open(((fileName)+".root").c_str(),"READ");
  TH1F* h = (TH1F*)( f->Get("pileup") );
  
  TCanvas* c = new TCanvas("c","pileup");
  c -> cd();
  c -> SetGridx();
  c -> SetGridy();
  
  h -> GetXaxis() -> SetTitle("n_{PU}");
  h -> GetXaxis() -> SetLabelSize(0.04);
  h -> GetYaxis() -> SetLabelSize(0.04);
  h -> GetYaxis() -> SetTitle("event fraction");
  h -> Scale(1./h->Integral());
  h -> SetMarkerStyle(20);
  h -> Draw("P");
  
  c -> Update();
  
  for(int bin = 1; bin <= h->GetNbinsX(); ++bin)
  {
    if( int(h->GetBinCenter(bin)) < 10 )
      std::cout << "distrPU_DATA[" << int(h->GetBinCenter(bin)) <<"]  = ";
    else
      std::cout << "distrPU_DATA[" << int(h->GetBinCenter(bin)) <<"] = ";
    std::cout << std::fixed << std::setprecision(9) << std::setw(11) << h -> GetBinContent(bin);
    std::cout << ";" << std::endl;
  }
  
  
  
  TFile* o = new TFile(((fileName)+"_plot.root").c_str(),"RECREATE");
  o -> cd();
  
  h -> Write();
  
  o -> Close();
}

// Called within makePlot when making the pu plot
// Basically adds in the +/- 5% systematic error bars
TGraphAsymmErrors makePUPlot( THStack stack, TString anaType, const float lumi ) {

  // Got the stack, now sum the =/- 5% histograms and make two new stacks
  // Make dummy legend etc.
  THStack p5Stack("Background MC p5","");
  THStack m5Stack("Background MC m5","");
  TLegend *legend= new TLegend(0.1,0.1,0.1,0.1);

  addBackgroundHistos( anaType, "nRecoPV_p5", p5Stack, legend, lumi );
  addBackgroundHistos( anaType, "nRecoPV_m5", m5Stack, legend, lumi );

  // Get last histogram in stacks i.e. sum of all backgrounds
  TH1F * hist = new TH1F( *(TH1F*)(stack.GetStack()->Last()) );
  TH1F * p5Hist = new TH1F( *(TH1F*)p5Stack.GetStack()->Last() );
  TH1F * m5Hist = new TH1F( *(TH1F*)m5Stack.GetStack()->Last() );

  // Create structures to make TGraphAssymErrors
  Int_t n=p5Hist->GetNbinsX();
  Double_t x[n]; // x values
  Double_t y[n]; // y values
  Double_t exl[n]; // low x error
  Double_t eyl[n]; // low y error
  Double_t exh[n]; // high x error
  Double_t eyh[n]; // high y error

  // Loop over histos and get values
  for ( int bin=1; bin < n+1; bin++ ) {
    x[bin-1]=hist->GetBinCenter( bin );
    y[bin-1]=hist->GetBinContent( bin );
    exl[bin-1]=hist->GetBinWidth( bin )/2;
    eyl[bin-1]=p5Hist->GetBinContent( bin )-y[bin-1];
    exh[bin-1]=hist->GetBinWidth( bin )/2;
    eyh[bin-1]=y[bin-1]-m5Hist->GetBinContent( bin );
  }


  // Make TGraphAssymErrors
  TGraphAsymmErrors gr(n,x,y,exl,exh,eyl,eyh);

  return gr;
}

//----------------------------------------//
// Get profile from all events mixed
//----------------------------------------//
TProfile* getMixProfile(TFile* file, int nbins,
			float xmin, float xmax,
			int color, int marker)
{

  // Histograms all have same base name
  string base = "h_evt";
  
  // String stream easy for manipulation
  stringstream ss;
  
  // Create Profile
  TProfile* p = new TProfile("mixed","",nbins,xmin,xmax);
  p->SetMarkerSize(0.75);
  p->SetMarkerStyle(marker);
  p->SetMarkerColor(color);
  p->SetLineColor(color);

  // Now loop over events and get result
  TH1F* h = NULL;
  for(int i=0; i<m_nEvent; ++i){
    ss.str("");
    ss << base << i;
    h = (TH1F*) file->Get(ss.str().c_str());

    // Loop over bin content
    for(int bin=1; bin<=h->GetNbinsX(); ++bin){
      float content = h->GetBinContent(bin);
      float center  = h->GetBinCenter(bin);
      p->Fill(center, content);
    }
    
    h = NULL;

  }// end loop events

  return p;

}

void Plot::FitSignal(int mode, int fitMode) {

	const int nPar = 6;
	TRandom ran;

	if(mode==0) {
		gStyle->SetOptLogy(1);
	} else {
		gStyle->SetOptLogy(0);
	}
	gStyle->SetOptStat(0);
	gStyle->SetOptTitle(0);
	const float limitBinSize = 2.0; // **** bin size here
	TCanvas* c = NewCanvas();
	c->Divide(1,2);

	gROOT->cd();
	TH1F* cc = new TH1F("CCSignal","CC Signal",500,0.0,1000.0);
	TH1F* ccBg = new TH1F("CCBgFit","CC Bg Fit",500,0.0,1000.0);
	TH1F* ccBgErr = new TH1F("CCBgErr","CC Bg Err",500,0.0,1000.0);
	TH1F* ccBgP = new TH1F("CCBgPlus","CC Bg Plus",500,0.0,1000.0);
	TH1F* ccBgM = new TH1F("CCBgMinus","CC Bg Minus",500,0.0,1000.0);
	TH1F* cp = new TH1F("CPSignal","CP Signal",500,0.0,1000.0);
	TH1F* cpBg = new TH1F("CPBgFit","CP Bg Fit",500,0.0,1000.0);
	TH1F* cpBgErr = new TH1F("CPBgErr","CP Bg Err",500,0.0,1000.0);
	TH1F* cpBgP = new TH1F("CPBgPlus","CP Bg Plus",500,0.0,1000.0);
	TH1F* cpBgM = new TH1F("CPBgMinus","CP Bg Minus",500,0.0,1000.0);
	TMatrixD matrix(nPar,nPar);

	fd->cd();

	TH1F* hInt,*hBgInt;
	char fitname[100];
	for(int ind=0; ind<2; ind++) {
		if(debug) printf("starting ind %i\n",ind);
		c->cd(ind+1);
		gStyle->SetOptLogy(1);
		printf("Starting %i ######################################\n",ind);

		TH1F* h;
		//char cind[20];
		//char handle[100];
		//sprintf(handle,"side_1exp_%02i_%02i_%02i",ind,mode,fitMode);
		TF1* fits[4];
		//TF1* dpx[4];
		if(debug) printf("looking for h %i\n",int(fd));
		if(ind==0) {
			h = (TH1F*)fd->FindObjectAny("pair_mass_2GeV1");
		} else if(ind==1) {
			h = (TH1F*)fd->FindObjectAny("pair_mass_2GeV3");
		}
		if(debug) printf("new h %i\n",int(h));

		if(debug) printf("new fit\n");
		sprintf(fitname,"hfit_%1i",ind);
		fits[ind] = new TF1(fitname,"([0]*pow((x-30.0),[1])+[3]*pow((x-30.0),0.2))*([2]*exp(-[2]*(x-30.0))+[4]*[5]*exp(-[5]*(x-30.0)))",30.0,500.0);
		//fits[ind] = new TF1(fitname,"([0]*((1-[3])*pow((x-30.0),[1])+[3]*pow((x-30.0),0.2)))*(exp(-[2]*(x-30.0))+[4]*exp(-[5]*(x-30.0)))",30.0,500.0);
		fits[ind]->SetParameter(0,0.0004);
		fits[ind]->SetParameter(1,2);
		fits[ind]->SetParameter(2,0.02);
		fits[ind]->SetParameter(3,0.005);
		//fits[ind]->SetParameter(3,0.5);
		fits[ind]->SetParameter(4,1.005);
		fits[ind]->SetParameter(5,0.05);

		float llim = 30.0;
		h->Fit(fits[ind],"LN","",llim,1000.0);

		double par[20],parMin[20],fval,fvalMin;
		for(int i=0; i<nPar; i++) parMin[i] = fits[ind]->GetParameter(i);

		gMinuit->Eval(nPar,0,fvalMin,parMin,0);
		//printf("got back %10.5f\n",fvalMin);

		// save the fit results in a histogram, for limit program
		for(int ibin=16; ibin<250; ibin++) {
			float xx = h->GetBinCenter(ibin);
			float yy = fits[ind]->Eval(xx);
			if(ind==0) {
				cc->SetBinContent(ibin,h->GetBinContent(ibin));
				ccBg->SetBinContent(ibin,yy);
				ccBgErr->SetBinContent(ibin,0.0);
				ccBgP->SetBinContent(ibin,0.0);
				ccBgM->SetBinContent(ibin,99999.0);
			} else {
				cp->SetBinContent(ibin,h->GetBinContent(ibin));
				cpBg->SetBinContent(ibin,yy);
				cpBgErr->SetBinContent(ibin,0.0);
				cpBgP->SetBinContent(ibin,0.0);
				cpBgM->SetBinContent(ibin,99999.0);
			}
		}

		//vary the parameters to find an error envelope
		double par2[20],fval2=1e10;
		int pslim = (ind==0?25000:150000);
		for(int ips=0; ips<pslim; ips++) {
			if(ips%10000==0) printf("Processing %d\n",ips);
			for(int i=0; i<nPar; i++) {
				par[i] = parMin[i];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }
  for(i=0;i<fo_BIN_NUM;i++)
    {
      hptAlice->SetBinContent(i+1,a11_Alice[i]*pdfrac);
      hminallAlice->SetBinContent(i+1,min_all_Alice[i]*pdfrac);
      hmaxallAlice->SetBinContent(i+1,max_all_Alice[i]*pdfrac);
      hpt11Alice->SetBinContent(i+1,a11_Alice[i]*pdfrac);
      hpt55Alice->SetBinContent(i+1,a55_Alice[i]*pdfrac);
      hpt22Alice->SetBinContent(i+1,a22_Alice[i]*pdfrac);
      hpt21Alice->SetBinContent(i+1,a21_Alice[i]*pdfrac);
      hpt12Alice->SetBinContent(i+1,a12_Alice[i]*pdfrac);
      hpt15Alice->SetBinContent(i+1,a15_Alice[i]*pdfrac);
      hpt51Alice->SetBinContent(i+1,a51_Alice[i]*pdfrac);
      hptAnly->SetBinContent(i+1,a11_Anly[i]*pdfrac);
      hpt11Anly->SetBinContent(i+1,a11_Anly[i]*pdfrac);
      hpt55Anly->SetBinContent(i+1,a55_Anly[i]*pdfrac);
      hpt22Anly->SetBinContent(i+1,a22_Anly[i]*pdfrac);
      hpt21Anly->SetBinContent(i+1,a21_Anly[i]*pdfrac);
      hpt12Anly->SetBinContent(i+1,a12_Anly[i]*pdfrac);
      hpt15Anly->SetBinContent(i+1,a15_Anly[i]*pdfrac);
      hpt51Anly->SetBinContent(i+1,a51_Anly[i]*pdfrac);
    }

  //Calculate data weighted ptbin centers
  cout<<"---- Calculate data weighted ptbin centers"<<endl;
  Double_t da_bincenter[da_BIN_NUM];
  Double_t summul=0,sum=0;
  for(i=0;i<da_BIN_NUM;i++)
    {
      summul=0;
      sum=0;
      for(j=da_ptbin[i]*4;j<da_ptbin[i+1]*4;j++)
	{
	  summul+=hptAlice->GetBinContent(j+1)*hptAlice->GetBinCenter(j+1);
	  sum+=hptAlice->GetBinContent(j+1);
	}
      da_bincenter[i] = summul/sum;
    }
  
  //Back up
  TH1F* hptAliceRef=(TH1F*)hptAlice->Clone();
  hptAliceRef->SetName("hptAliceRef");
  TH1F* hptAnlyRef=(TH1F*)hptAnly->Clone();
  hptAnlyRef->SetName("hptAnlyRef");

  Double_t apt[da_BIN_NUM],aptl[da_BIN_NUM],aptlsyst[da_BIN_NUM],aptlscaling[da_BIN_NUM];
  TH1F* hptAlice_rebinalice = (TH1F*)hptAliceRef->Rebin(da_BIN_NUM,"hptAlice_rebinalice",da_ptbin);
  for(i=0;i<da_BIN_NUM;i++)
    {
      apt[i] = hptAlice_rebinalice->GetBinCenter(i+1);
      aptl[i] = hptAlice_rebinalice->GetBinWidth(i+1)/2.;
      aptlsyst[i] = 1./6.;
      aptlscaling[i] = 1./4.;
    }
  
  //Test ALICE results
  if(testAliceData)
    {
      cout<<"---- Test Alice data points"<<endl;
      Double_t asigma[da_BIN_NUM],aerrorl[da_BIN_NUM],aerrorh[da_BIN_NUM],aminall[da_BIN_NUM],amaxall[da_BIN_NUM];
      TH1F* hminallAlice_rebinalice = (TH1F*)hminallAlice->Rebin(da_BIN_NUM,"hminallAlice_rebinalice",da_ptbin);
      TH1F* hmaxallAlice_rebinalice = (TH1F*)hmaxallAlice->Rebin(da_BIN_NUM,"hmaxallAlice_rebinalice",da_ptbin);
      for(i=0;i<da_BIN_NUM;i++)
	{
	  hptAlice_rebinalice->SetBinContent(i+1,hptAlice_rebinalice->GetBinContent(i+1)/(hptAlice_rebinalice->GetBinWidth(i+1)*4));
	  asigma[i]=hptAlice_rebinalice->GetBinContent(i+1);

//.........这里部分代码省略.........
         CLb+= (1.)/nToyMC;
      if(obsTSatObsUL <= toyTSatObsUL) // not sure about <= part yet
         CLbinclusive+= (1.)/nToyMC;


      // loop over points in belt to find upper limit for this toy data
      double thisUL = 0;
      for(Int_t i=0; i<parameterScan->numEntries(); ++i){
         tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
         double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
         firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
         //   double thisTS = profile->getVal();
         double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);

         //   cout << "poi = " << firstPOI->getVal()
         // << " max is " << arMax << " this profile = " << thisTS << endl;
         //      cout << "thisTS = " << thisTS<<endl;
         if(thisTS<=arMax){
            thisUL = firstPOI->getVal();
         } else{
            break;
         }
      }



      /*
      // loop over points in belt to find upper limit for this toy data
      double thisUL = 0;
      for(Int_t i=0; i<histOfThresholds->GetNbinsX(); ++i){
         tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
         cout <<"----------------  "<<i<<endl;
         tmpPoint->Print("v");
         cout << "from hist " << histOfThresholds->GetBinCenter(i+1) <<endl;
         double arMax = histOfThresholds->GetBinContent(i+1);
         // cout << " threhold from Hist = aMax " << arMax<<endl;
         // double arMax2 = belt->GetAcceptanceRegionMax(*tmpPoint);
         // cout << "from scan arMax2 = "<< arMax2 << endl; // not the same due to TH1F not TH1D
         // cout << "scan - hist" << arMax2-arMax << endl;
         firstPOI->setVal( histOfThresholds->GetBinCenter(i+1));
         //   double thisTS = profile->getVal();
         double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);

         //   cout << "poi = " << firstPOI->getVal()
         // << " max is " << arMax << " this profile = " << thisTS << endl;
         //      cout << "thisTS = " << thisTS<<endl;

         // NOTE: need to add a small epsilon term for single precision vs. double precision
         if(thisTS<=arMax + 1e-7){
            thisUL = firstPOI->getVal();
         } else{
            break;
         }
      }
      */

      histOfUL->Fill(thisUL);

      // for few events, data is often the same, and UL is often the same
      //    cout << "thisUL = " << thisUL<<endl;

      delete toyData;
   }
   histOfUL->Draw();
   c1->SaveAs("one-sided_upper_limit_output.pdf");