C++ TAxis::GetXmin方法代码示例

void BinLog(TH2F *h) 
{

   TAxis *axis = h->GetXaxis(); 
   int bins = axis->GetNbins();

   Axis_t from = axis->GetXmin();
   Axis_t to = axis->GetXmax();
   Axis_t width = (to - from) / bins;
   Axis_t *new_bins = new Axis_t[bins + 1];

   for (int i = 0; i <= bins; i++) {
     new_bins[i] = TMath::Power(10, from + i * width);
   } 
   axis->Set(bins, new_bins); 

   TAxis *axis2 = h->GetYaxis(); 
   int bins2 = axis2->GetNbins();

   Axis_t from2 = axis2->GetXmin();
   Axis_t to2 = axis2->GetXmax();
   Axis_t width2 = (to2 - from2) / bins2;
   Axis_t *new_bins2 = new Axis_t[bins2 + 1];

   for (int i = 0; i <= bins2; i++) {
     new_bins2[i] = TMath::Power(10, from2 + i * width2);
   } 
   axis2->Set(bins2, new_bins2); 

   delete new_bins;
   delete new_bins2;
}

 /** comapre axis for equalitiy
  *
  * check whether the cass histogram axis and the root histogram axis are the
  * same. Test for number of bins, low and high ends and the title of the axis.
  *
  * @return true when both axis are the same
  * @param ca the axis of the cass histogram
  * @param ra the axis of the root histogram
  *
  * @author Lutz Foucar
  */
 bool operator== (const cass::Result<float>::axe_t &ca, const TAxis &ra)
 {
   return (static_cast<int>(ca.nBins) == ra.GetNbins() &&
           fabs(ca.low - ra.GetXmin()) <  sqrt(numeric_limits<double>::epsilon()) &&
           fabs(ca.up  - ra.GetXmax())  <  sqrt(numeric_limits<double>::epsilon()) &&
           ca.title == ra.GetTitle());
 }

//______________________________________________________________________________
void drawsparse_draw(THnSparse* h)
{
   // Draw a THnSparse using TParallelCoord, creating a temporary TTree.

   TTree* tree = toTree(h);

   TString whatToDraw;
   TIter iLeaf(tree->GetListOfLeaves());
   const TLeaf* leaf = 0;
   while ((leaf = (const TLeaf*)iLeaf())) {
      if (whatToDraw.Length())
         whatToDraw += ":";
      whatToDraw += leaf->GetName();
   }
   tree->Draw(whatToDraw, "", "para");
   TParallelCoord* parallelCoord = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject("ParaCoord");

   TIter iVar(parallelCoord->GetVarList());
   TParallelCoordVar* var = 0;
   for (Int_t d = 0;(var = (TParallelCoordVar*) iVar()) && d < h->GetNdimensions(); ++d) {
      TAxis* axis = h->GetAxis(d);
      var->SetHistogramBinning(axis->GetNbins());
      var->SetCurrentLimits(axis->GetXmin(), axis->GetXmax());
      var->SetTitle(axis->GetTitle());
   }
   var->SetTitle("bin content");
}

int getBin(TH1* theHist, double x) {

    // Find the bin number for the given x value
    if (!theHist) {return 0;}

    TAxis* xAxis = theHist->GetXaxis();
    double xMin = xAxis->GetXmin();
    double xMax = xAxis->GetXmax();
    int nX = theHist->GetNbinsX();
    double dx = (xMax - xMin)/(nX*1.0);
    
    int theBin = int(((x - xMin)/dx) + 0.5);
    return theBin;
    
}

void findFilterEnergyLoss::BinLogY(TH2 *h) {
   TAxis *axis = h->GetYaxis();
   int bins = axis->GetNbins();

   Axis_t from = axis->GetXmin();
   Axis_t to = axis->GetXmax();
   Axis_t width = (to - from) / bins;

   Axis_t *new_bins = new Axis_t[bins+1];

   for (int i=0; i <= bins; i++) {
      new_bins[i] = TMath::Power(10, from + i * width);
   }

   axis->Set(bins, new_bins);
   delete new_bins;
}

void BinLogAxis(const TH1 *h)
{
  //
  // Method for the correct logarithmic binning of histograms
  //
  TAxis *axis = const_cast<TAxis*>(h->GetXaxis());
  const Int_t bins = axis->GetNbins();
  
  const Double_t from = axis->GetXmin();
  const Double_t to = axis->GetXmax();
  Double_t *newBins = new Double_t[bins + 1];
  
  newBins[0] = from;
  Double_t factor = pow(to / from, 1. / bins);
  
  for (Int_t i = 1; i <= bins; i++) {
    newBins[i] = factor * newBins[i - 1];
  }
  axis->Set(bins, newBins);
  delete [] newBins;
  
}

Double_t effSigma(TH1 *hist )
{

  TAxis *xaxis = hist->GetXaxis();
  Int_t nb = xaxis->GetNbins();
  if(nb < 10) {
    cout << "effsigma: Not a valid histo. nbins = " << nb << endl;
    return 0.;
  }
  
  Double_t bwid = xaxis->GetBinWidth(1);
  if(bwid == 0) {
    cout << "effsigma: Not a valid histo. bwid = " << bwid << endl;
    return 0.;
  }
  Double_t xmin = xaxis->GetXmin();
  Double_t ave = hist->GetMean();
  Double_t rms = hist->GetRMS();

  Double_t total=0.;
  for(Int_t i=0; i<nb+2; i++) {
    total+=hist->GetBinContent(i);
  }
  if(total < 100.) {
    cout << "effsigma: Too few entries " << total << endl;
    return 0.;
  }
  Int_t ierr=0;
  Int_t ismin=999;
  
  Double_t rlim=0.683*total;
  Int_t nrms=rms/(bwid);    // Set scan size to +/- rms
  if(nrms > nb/10) nrms=nb/10; // Could be tuned...

  Double_t widmin=9999999.;
  for(Int_t iscan=-nrms;iscan<nrms+1;iscan++) { // Scan window centre
    Int_t ibm=(ave-xmin)/bwid+1+iscan;
    Double_t x=(ibm-0.5)*bwid+xmin;
    Double_t xj=x;
    Double_t xk=x;
    Int_t jbm=ibm;
    Int_t kbm=ibm;
    Double_t bin=hist->GetBinContent(ibm);
    total=bin;
    for(Int_t j=1;j<nb;j++){
      if(jbm < nb) {
        jbm++;
        xj+=bwid;
        bin=hist->GetBinContent(jbm);
        total+=bin;
        if(total > rlim) break;
      }
      else ierr=1;
      if(kbm > 0) {
        kbm--;
        xk-=bwid;
        bin=hist->GetBinContent(kbm);
        total+=bin;
        if(total > rlim) break;
      }
      else ierr=1;
    }
    Double_t dxf=(total-rlim)*bwid/bin;
    Double_t wid=(xj-xk+bwid-dxf)*0.5;
    if(wid < widmin) {
      widmin=wid;
      ismin=iscan;
    }   
  }
  if(ismin == nrms || ismin == -nrms) ierr=3;
  if(ierr != 0) cout << "effsigma: Error of type " << ierr << endl;
  
  return widmin;
  
}

void makePlot_wrapper(const TString& title, 
		      const TH1* histogram_numerator_Data_passed, const TH1* histogram_denominator_Data_passed, 
		      const TH1* histogram_numerator_mcSum_passed, const TH1* histogram_denominator_mcSum_passed, 
		      const TH1* histogram_numerator_Data_failed, const TH1* histogram_denominator_Data_failed, 
		      const TH1* histogram_numerator_mcSum_failed, const TH1* histogram_denominator_mcSum_failed, 
		      const TString& legendEntry_Data, const TString& legendEntry_mcSum, 
		      const TString& region_passed, const TString& region_failed,
		      const TString& outputFileName)
{
  TAxis* xAxis = histogram_numerator_Data_passed->GetXaxis();
  Double_t xMin = xAxis->GetXmin();
  Double_t xMax = xAxis->GetXmax();

  std::cout << "processing Data, region " << region_passed.Data() << "..." << std::endl;
  TGraphAsymmErrors* graph_efficiency_Data_passed = 
    getEfficiency(histogram_numerator_Data_passed, histogram_denominator_Data_passed);   
  TF1* fit_efficiency_Data_passed = new TF1("fit_efficiency_Data_passed", &integralCrystalBall_f, xMin, xMax, 5);
  initializeCrystalBall(fit_efficiency_Data_passed);
  graph_efficiency_Data_passed->Fit(fit_efficiency_Data_passed, "0");
  TString legendEntry_Data_passed = Form("%s, %s", legendEntry_Data.Data(), region_passed.Data());

  std::cout << "processing mcSum, region " << region_passed.Data() << "..." << std::endl;
  TGraphAsymmErrors* graph_efficiency_mcSum_passed = 
    getEfficiency(histogram_numerator_mcSum_passed, histogram_denominator_mcSum_passed);   
  TF1* fit_efficiency_mcSum_passed = new TF1("fit_efficiency_mcSum_passed", &integralCrystalBall_f, xMin, xMax, 5);
  initializeCrystalBall(fit_efficiency_mcSum_passed);
  graph_efficiency_mcSum_passed->Fit(fit_efficiency_mcSum_passed, "0");
  TString legendEntry_mcSum_passed = Form("%s, %s", legendEntry_mcSum.Data(), region_passed.Data());

  std::cout << "processing Data, region " << region_failed.Data() << "..." << std::endl;
  TGraphAsymmErrors* graph_efficiency_Data_failed = 
    getEfficiency(histogram_numerator_Data_failed, histogram_denominator_Data_failed);   
  TF1* fit_efficiency_Data_failed = new TF1("fit_efficiency_Data_failed", &integralCrystalBall_f, xMin, xMax, 5);
  initializeCrystalBall(fit_efficiency_Data_failed);
  graph_efficiency_Data_failed->Fit(fit_efficiency_Data_failed, "0");
  TString legendEntry_Data_failed = Form("%s, %s", legendEntry_Data.Data(), region_failed.Data());

  std::cout << "processing mcSum, region " << region_failed.Data() << "..." << std::endl;
  TGraphAsymmErrors* graph_efficiency_mcSum_failed = 
    getEfficiency(histogram_numerator_mcSum_failed, histogram_denominator_mcSum_failed);   
  TF1* fit_efficiency_mcSum_failed = new TF1("fit_efficiency_mcSum_failed", &integralCrystalBall_f, xMin, xMax, 5);
  initializeCrystalBall(fit_efficiency_mcSum_failed);
  graph_efficiency_mcSum_failed->Fit(fit_efficiency_mcSum_failed, "0");
  TString legendEntry_mcSum_failed = Form("%s, %s", legendEntry_mcSum.Data(), region_failed.Data());

  makePlot(title,
	   graph_efficiency_Data_passed,  fit_efficiency_Data_passed,  legendEntry_Data_passed,
	   graph_efficiency_mcSum_passed, fit_efficiency_mcSum_passed, legendEntry_mcSum_passed,
	   graph_efficiency_Data_failed,  fit_efficiency_Data_failed,  legendEntry_Data_failed,
	   graph_efficiency_mcSum_failed, fit_efficiency_mcSum_failed, legendEntry_mcSum_failed,
	   outputFileName);

  delete graph_efficiency_Data_passed;
  delete fit_efficiency_Data_passed;
  delete graph_efficiency_mcSum_passed;
  delete fit_efficiency_mcSum_passed;
  delete graph_efficiency_Data_failed;
  delete fit_efficiency_Data_failed;
  delete graph_efficiency_mcSum_failed;
  delete fit_efficiency_mcSum_failed;
}

void BkgDemo_Diplot(TCanvas* c, TH1F** h, TLegend* l1, TLegend* l2){
		//Make a split canvas, comparison of methods and tag above, bkg subtracted tag below
//	printf("\nDo BkgDemo_Diplot for canvas %s\n",c->GetName());
//	printf("integrals stage1: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());

	
	TLatex * TEX_CMSPrelim;
	if(preliminary) TEX_CMSPrelim = new TLatex(0.177136,0.953368,"CMS Preliminary"); 
	else TEX_CMSPrelim = new TLatex(0.177136,0.953368,"CMS"); 
	PrettyLatex(TEX_CMSPrelim,0.03);
	TLatex * TEX_E_TeV = new TLatex(0.800251,0.953368,"#sqrt{s} = 8 TeV");
	PrettyLatex(TEX_E_TeV,0.03);
	TLatex * TEX_lumi_fb = new TLatex(0.621859,0.953368,Form("#intL dt = %.1f fb^{-1}",Integrated_Luminosity_Data));
	PrettyLatex(TEX_lumi_fb,0.03);
	string sp1 = string(c->GetName())+"_p1";
	string sp2 = string(c->GetName())+"_p2";
    TPad *p1 = new TPad((char*)sp1.c_str(),"",0.,0.3,1.,1.);
    p1->SetBottomMargin(0);
    p1->cd();
//	c->Divide(1,2);
//	c->cd(1);
	PrettyHist(h[3],kRed);
	PrettyHist(h[5],kBlue);
	PrettyHist(h[1]);
	PrettyHist(h[6],kGreen);
	//PrettyHist(h[3],kRed);
	PrettyHist(h[7],kTeal);
PrettyBlock2(h[5],kBlue,3354,2);
	//PrettyBlock(h[5],kBlue,string("//thatch"));//PrettyMarker(h[5],kBlue,4);
	PrettyMarker(h[1]);
PrettyBlock2(h[6],kGreen,3345,2);
	//PrettyBlock(h[6],kGreen,string("\\thatch"));//PrettyMarker(h[6],kGreen,4);
	//PrettyMarker(h[3],kRed);
	PrettyMarker(h[7],kTeal);
//	h[5]->Scale(mHwidth/sidebandwidth);//lsb scaled
//	h[6]->Scale(mHwidth/sidebandwidth);//usb scaled
//	printf("integrals stage2: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());
//	printf("ranges maximi before setrange tag %f lsb %f, usb %f\n",h[1]->GetMaximum(),h[5]->GetMaximum(),h[6]->GetMaximum());
//	printf("integrals stage3: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());
//	printf("ranges maximi before after tag %f lsb %f usb %f\n",h[1]->GetMaximum(),h[5]->GetMaximum(),h[6]->GetMaximum());
	h[3]->SetFillStyle(0);//open rectangle
	h[3]->SetLineColor(kRed);
	h[3]->SetLineWidth(4);

	h[6]->SetMinimum(0.0);
	//float linmax = h[6]->GetMaximum();
	//float linmin = h[6]->GetMinimum();
	playNiceWithLegend(h[6],0.30,0.0);
	playNiceWithLegend(h[3],0.30,0.0);
	playNiceWithLegend(h[5],0.30,0.0);
	playNiceWithLegend(h[1],0.30,0.0);
	SameRange(h[3],h[1]);
	SameRange(h[3],h[5],h[6]);
	SameRange(h[3],h[1]);
	h[3]->Draw("e2p");
	h[5]->Draw("e2psame");
	h[6]->Draw("e2psame");
	h[3]->Draw("e2psame");
	//if(showTag) h[1]->Draw("e1psame");//tag
	TPad *p2 = new TPad((char*)sp2.c_str(),"",0.,0.125,1.,0.3);
    p2->SetTopMargin(0);
    p2->cd();
//	c->cd(2);
	TAxis * x = h[7]->GetXaxis();
	TAxis * y = h[7]->GetYaxis();
	float fontsize = 0.25;
	float fontsizeY = 0.10;
	x->SetTitleSize(fontsize);
	y->SetTitleSize(fontsizeY);
	x->SetLabelSize(fontsize);
	y->SetLabelSize(fontsizeY);

	h[7]->GetYaxis()->SetRangeUser(0.,2.);
	h[7]->Draw("e1p");
	TLine *OneLine = new TLine(x->GetXmin(),1.0,x->GetXmax(),1.0);
	OneLine->SetLineColor(kBlack);
	OneLine->SetLineWidth(2);
	OneLine->SetLineStyle(7);//dashed.
	OneLine->Draw("same");
	h[7]->Draw("e1psame");
	
	p1->cd();
//  c->cd(1);
	if(showTag) l1->AddEntry(h[1],"Higgs Mass Region");
        TH1F * box = new TH1F("box","asdf",1,0,1);
        box->SetMarkerColor(kRed);
        box->SetMarkerStyle(25);
        box->SetMarkerSize(2);
        box->SetLineColor(0);
	l1->AddEntry(box,"Data Driven Background");
	l1->AddEntry(h[5],"Lower Mass Sideband");
	l1->AddEntry(h[6], "Upper Mass Sideband");
	l1->Draw("same");
        TEX_CMSPrelim->Draw("same");
        TEX_E_TeV->Draw("same");
        TEX_lumi_fb->Draw("same");
	p2->cd();
		//c->cd(2);
	l2->SetTextSize(fontsize*0.8);
	l2->AddEntry(h[7],"Lower/Upper Sideband Ratio");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			playNiceWithLegend(h[3],0.30,0.0);
                        playNiceWithLegend(h[6],0.30,0.0);
                        playNiceWithLegend(h[5],0.30,0.0);
                        playNiceWithLegend(h[1],0.30,0.0);
			h[6]->SetMinimum(0.0);
			SameRange(h[1],h[3]);
			SameRange(h[5],h[6],h[1]);
			SameRange(h[1],h[3]);
			h[5]->Draw("e2p");
			h[6]->Draw("e2psame");
			h[3]->Draw("e2psame");
			//if(showTag) h[1]->Draw("e1psame");//tag
			TLegend* l1 = makeL1_v2(0.443467,0.720207,0.643216,0.919689);
			if(showTag) l1->AddEntry(h[1],"Higgs Mass Region");
			l1->AddEntry(h[5],"Lower Mass Sideband");
			l1->AddEntry(h[6],"Upper Mass Sideband");
			l1->AddEntry(box,"Data Driven Background");//h[3
			l1->Draw("same");
		        TEX_CMSPrelim->Draw("same");
			TEX_E_TeV->Draw("same");
			TEX_lumi_fb->Draw("same");
			tempCanv->Write();
			if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggif");//save as all types
			//if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggifpdfeps");//save as all types

			tempCanvlog->cd();
			h[5]->SetMinimum(0.05);
                        tempCanvlog->SetLogy();
			h[5]->Draw("e2p");
			h[6]->Draw("e2psame");
			h[3]->Draw("e2psame");
			l1->Draw("same");
			TEX_CMSPrelim->Draw("same");
			TEX_E_TeV->Draw("same");
			TEX_lumi_fb->Draw("same");
			if(saveImages)SaveCanvas(tempCanvlog,plotsdir+tempCanvlog->GetName(),"ggif");//save as all types
			
		}//edn for each kinematic varriable
	}//end for each topology
	}//end if


//     ___  __          ___       __  _        ___  __     __    
//    / _ )/ /_____ _  / _ \___ _/ /_(_)__    / _ \/ /__  / /____
//   / _  /  '_/ _ `/ / , _/ _ `/ __/ / _ \  / ___/ / _ \/ __(_-<
//  /____/_/\_\\_, / /_/|_|\_,_/\__/_/\___/ /_/  /_/\___/\__/___/
//            /___/                                              

		///Demostrate background--Ratio. Plot the ratio of the two side bands, and nothing else. 
	if(makeBkgRat){
	for (int iTop = 0; iTop<nEventTopologies; iTop++) {
		for (int kKinVar = 0; kKinVar<nKinemVars; kKinVar++) {
			/*Debug*/ if(printlevel >6) cout << "Data MHT L/R bin 0:" <<KinVarHistMap["Data"]["NULL"]["MHT"][7]->GetBinContent(1)<<endl;
			if (s_EventTopology[iTop].compare("metCut")==0 && s_KinemVars[kKinVar].compare("MET")==0 ) continue;
			string canvName = string("BkgRat_")+s_KinemVars[kKinVar]+s_EventTopology_v2[iTop];
			TCanvas* tempCanv = newTCanvas((char*)canvName.c_str());
			TH1F** h = KinVarHistMap[Data][s_EventTopology[iTop]][s_KinemVars[kKinVar]];
			/*Debug*/ if(printlevel >6)	printf("Scan2 Ratio bin1 for %s %s  = %f\n",s_EventTopology[iTop].c_str(),s_KinemVars[kKinVar].c_str(),h[7]->GetBinContent(1));
			tempCanv->cd();
			PrettyHist(h[7],kTeal);
			h[7]->GetYaxis()->SetRangeUser(0.,2.);
			if (s_KinemVars[kKinVar].compare("HT") == 0) {
				h[7]->GetXaxis()->SetLabelSize(0.04);}

			PrettyMarker(h[7],kTeal);
			h[7]->Draw("e1p");
			TAxis* xaxis = h[7]->GetXaxis();
			TLine *OneLine = new TLine(xaxis->GetXmin(),1.0,xaxis->GetXmax(),1.0);
			OneLine->SetLineColor(kBlack);
			OneLine->SetLineWidth(2);
			OneLine->SetLineStyle(7);//dashed.
			OneLine->Draw("same");
			h[7]->Draw("e1psame");
			TLegend* l1 = makeL1_v2(0.330402,0.836788,0.530151,0.919689);
			l1->AddEntry(h[7],"Lower/Upper Sideband Ratio");
			l1->Draw("same");
			TEX_CMSPrelim->Draw("same");
			TEX_E_TeV->Draw("same");
			TEX_lumi_fb->Draw("same");

			tempCanv->Write();
			if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggif");//save as all types
			//if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggifpdfeps");//save as all types
		}//edn for each kinematic varriable
	}//end for each topology
	}//end if makeBkgRat

		///END MAKE PLOTS
	/*Debug*/ if(printlevel > 0) cout << "Close Files" << endl;
	fplots->Close(); 
		//itterate over the map. close everything. 
	for( TFileMap::iterator i=PostAnaAnaFiles.begin(); i!=PostAnaAnaFiles.end(); ++i){
		(*i).second->Close();
	}
	for( TFileMap::iterator i=MainAnaFiles.begin(); i!=MainAnaFiles.end(); ++i){
		(*i).second->Close();
	}


}//end format plots

void fit(TH1* histogram, double xMin, double xMax, double massValue, std::vector<fitParameterType>& fitParameter)
{
  // create fit variable
  TAxis* xAxis = histogram->GetXaxis();
  if ( xMin < 0. ) xMin = xAxis->GetXmin(); 
  if ( xMax < 0. ) xMax = xAxis->GetXmax(); 
  std::string fitVariableName = Form("%s_fitVariable", histogram->GetName());
  RooRealVar fitVariable(fitVariableName.data(), fitVariableName.data(), xMin, xMax);
   
  // convert histogram to RooFit format
  std::string histogramName_data = Form("%s_data", histogram->GetName());
  RooDataHist histogram_data(histogramName_data.data(), histogramName_data.data(), fitVariable, histogram, 1.0);

  // create product of linear * Heaviside * Gaussian Error function
  std::string svFitLineShapeName_slope = Form("%s_svFitLineShape_slope", histogram->GetName());
  RooRealVar svFitLineShape_slope(svFitLineShapeName_slope.data(), svFitLineShapeName_slope.data(), 2./(massValue*massValue), 0., 1.);
  std::string svFitLineShapeName_offset = Form("%s_svFitLineShape_offset", histogram->GetName());
  RooRealVar svFitLineShape_offset(svFitLineShapeName_offset.data(), svFitLineShapeName_offset.data(), 0., -1., +1.);

  std::string heavisideName_threshold = Form("%s_heaviside_threshold", histogram->GetName());
  RooRealVar heaviside_threshold(heavisideName_threshold.data(), heavisideName_threshold.data(), 0.2*massValue, 1.e-2*massValue, 0.5*massValue);

  std::string sculptingName_bias = Form("%s_sculpting_bias", histogram->GetName());
  RooRealVar sculpting_bias(sculptingName_bias.data(), sculptingName_bias.data(), 0.25*massValue, 0.*massValue, 0.9*massValue);
  std::string sculptingName_width = Form("%s_sculpting_width", histogram->GetName());
  RooRealVar sculpting_width(sculptingName_width.data(), sculptingName_width.data(), 0.1*massValue, 1.e-2*massValue, 1.0*massValue);

  std::string svFitLineShapeName = Form("%s_svFitLineShape", histogram->GetName());
  std::string svFitLineShapeFormula = "(@0*@1 + @2)*0.5*(1.0 + TMath::Sign(+1, @0 - @3))*0.5*(1.0 + TMath::Erf((@0 - @4)/@5))";
  RooArgList svFitLineShapeArgs(fitVariable, 
				svFitLineShape_slope,
				svFitLineShape_offset,
				heaviside_threshold,
				sculpting_bias,
				sculpting_width);
  RooGenericPdf svFitLineShape(svFitLineShapeName.data(), svFitLineShapeName.data(), svFitLineShapeFormula.data(), svFitLineShapeArgs);

  // create Gaussian 
  std::string gaussianName_mean = Form("%s_gaussian_mean", histogram->GetName());
  RooConstVar gaussian_mean(gaussianName_mean.data(), gaussianName_mean.data(), 0.);
  std::string gaussianName_sigma = Form("%s_gaussian_sigma", histogram->GetName());
  RooRealVar gaussian_sigma(gaussianName_sigma.data(), gaussianName_sigma.data(), 0.2*massValue, 1.e-2*massValue, 0.5*massValue);
  std::string gaussianName = Form("%s_gaussian", histogram->GetName());
  RooGaussian gaussian(gaussianName.data(), gaussianName.data(), fitVariable, gaussian_mean, gaussian_sigma);

  // numerically convolute both PDFs using fast Fourier transform
  std::string fitFunctionName = Form("%s_fitFunction", histogram->GetName());
  RooFFTConvPdf fitFunction(fitFunctionName.data(), fitFunctionName.data(), fitVariable, svFitLineShape, gaussian);

  // fit histogram
  fitFunction.fitTo(histogram_data);

  // save fit parameter
  fitParameter.resize(6);
  fitParameter[0] = fitParameterType(&svFitLineShape_slope);
  fitParameter[1] = fitParameterType(&svFitLineShape_offset);
  fitParameter[2] = fitParameterType(&gaussian_sigma);
  fitParameter[3] = fitParameterType(&heaviside_threshold);
  fitParameter[4] = fitParameterType(&sculpting_bias);
  fitParameter[5] = fitParameterType(&sculpting_width);

  // create control plot
  std::string frameTitle = Form("%s_frame", histogram->GetName());
  RooPlot* frame = fitVariable.frame(RooFit::Title(frameTitle.data()));
  histogram_data.plotOn(frame);
  fitFunction.plotOn(frame, RooFit::LineColor(kRed));
  
  std::string canvasName = Form("%s_canvas", histogram->GetName());
  TCanvas* canvas = new TCanvas(canvasName.data(), canvasName.data(), 800, 600);
  gPad->SetLeftMargin(0.15); 
  frame->GetYaxis()->SetTitleOffset(1.4);
  frame->Draw();
    
  std::string outputFileName_plot = Form("svFitPerformance_%s_fit", histogram->GetName());
  canvas->Print(std::string(outputFileName_plot).append(".eps").data());
  canvas->Print(std::string(outputFileName_plot).append(".png").data());
  canvas->Print(std::string(outputFileName_plot).append(".pdf").data());
  
  delete canvas;
}

//.........这里部分代码省略.........
    if ( nbinsx > 50 && nbinsx <= 100 )  nrebin = 10;//10
    if ( nbinsx > 20 && nbinsx <= 50 )   nrebin = 2;
    if ( nbinsx <= 20 ) nrebin = 1;  
    
    printf(" Saw nbins =%4i, rebinning by nrebin =%2i to final %3i bins \n",nbinsx,nrebin,int(nbinsx/nrebin));	

    if ( nrebin != 0 ) {
	h0->Rebin(nrebin); // data
	for (Int_t i = 0; i<=nbkg; i++){
	    TH1F * h = (TH1F*) addedhistos[i];
	    h->Rebin(nrebin);
	}
	for (Int_t i = 0; i<nsig; i++){
	    TH1F * h = (TH1F*) signalhistos[i];
	    h->Rebin(nrebin);
	}
    }

    // default text size: 0.045
    // make it bigger for the paper
    float textSize = 0.045;
    if(params->GetValue("Histo.Preliminary","yes")==TString("paper")) textSize=0.07;
    if(params->Defined("Histo.TextSize")) textSize=params->GetValue("Histo.TextSize",0.07);
    
    // Now, check largest dimensions so that we can plot all histograms at once.
    Float_t xmin=9999., xmax=-9999., ymin=9999., ymax=-9999.;
    for(Int_t i = 0; i<=nbkg; i++){
	TH1F * h = (TH1F*) addedhistos[i];
	ostringstream baseSrcName;
	baseSrcName << "Files." << i+1 << ".";
	TString bSrcName(baseSrcName.str().c_str());
	
	TAxis *axis = h->GetXaxis();
	if( axis->GetXmin() < xmin ) xmin = axis->GetXmin();
	if( axis->GetXmax() > xmax ) xmax = axis->GetXmax();
	if( h->GetMinimum() < ymin ) ymin = h->GetMinimum();
	if( h->GetMaximum() > ymax ) ymax = h->GetMaximum();
    }
    ymax = TMath::Nint(ymax*1.25+1); // Make enough room for the big legend
    TString title = h0->GetTitle();
    
    //
    // now check if we should simply use the ranges that was passed to us.
    if(params->Defined("Histo.Xmin")) xmin = params->GetValue("Histo.Xmin",0.);
    if(params->Defined("Histo.Xmax")) xmax = params->GetValue("Histo.Xmax",0.);
    if(params->Defined("Histo.Ymin")) ymin = params->GetValue("Histo.Ymin",0.);
    if(params->Defined("Histo.Ymax")) ymax = params->GetValue("Histo.Ymax",0.);
    
    // Now make the frame:
    TH2F * frame = new TH2F("frame","",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
    cout<<" frame has xmin "<<xmin<<", xmax "<<xmax<<", ymax "<<ymax<<endl;
    
    // get the x- and y-axis titles
    TString ytitle=params->GetValue("Histo.YTitle","");
    if ( params->Defined("Histo.XTitle")) {
	frame->SetXTitle(params->GetValue("Histo.XTitle",""));
    } else {
	frame->SetXTitle(h0->GetTitle());
    }
    frame->SetYTitle(ytitle.Data()); 
    // also set the text size for the X and Y axis titles and numbers
    // do this globally for the style we are using
    float axisLabelSize=textSize;
    frame->GetXaxis()->SetLabelSize(axisLabelSize);
    frame->GetYaxis()->SetLabelSize(axisLabelSize);
    frame->GetXaxis()->SetTitleSize(axisLabelSize);

//.........这里部分代码省略.........
            //	if(r->IsValid()) r->Print();
            //else cout<<"Fit not valid!!!\n"<<endl;
        }

        //------  get fit parameters
        double par[20];
        RBWPOL->GetParameters(par);

        float GGphDimuMass = RBWPOL->GetParameter(1);
        float GGZ0Width    = RBWPOL->GetParameter(2);
        float GauWidth     =0;
        if(nFitFunction  == 3 || nFitFunction  == 4) GauWidth = RBWPOL->GetParameter(3);

        double chisq      = RBWPOL->GetChisquare();
        int ndf           = RBWPOL->GetNDF();
        double chisqdf    =1000;
        if(ndf!=0) chisqdf=chisq/ndf;

        // +++ set backgroudn fit
        sprintf(szBuf,"pt_1B_%d",ih);
        TF1 *bkgFit_1 = new TF1(szBuf, Pol2, massFit_low, massFit_high, 3);
        // if(nFitFunction  == 4) bkgFit_1 = new TF1(namePt_1B, Exp, massFit_low, massFit_high, 2);

        bkgFit_1->SetParameters(&par[3]);
        if(nFitFunction  == 3 || nFitFunction  == 4) bkgFit_1->SetParameters(&par[4]);

        // ----------  Integrated Yield
        //    float massCount_low =GGphDimuMass-(4.0*GGZ0Width);
        //    float massCount_high =GGphDimuMass+(4.0*GGZ0Width);

        TAxis *axs        = dimuonsGlobalInvMassVsPt[ih]->GetXaxis();
        int binlow        = axs->FindBin(massCount_low);
        int binhi         = axs->FindBin(massCount_high);
        Double_t bin_size = (1.0*dimuonsGlobalInvMassVsPt[ih]->GetNbinsX())/(axs->GetXmax() - axs->GetXmin());

        Float_t int_sig   = 0.0;
        Float_t int_sig_sqr   = 0.0;
        for(Int_t bin = binlow; bin<=binhi; bin++) {
            //	  cout << "	  int_sig += dimuonsGlobalInvMassVsPt[ih]->GetBinContent(bin);"<<int_sig <<"+="<< "bin" << bin << " content"<<dimuonsGlobalInvMassVsPt[ih]->GetBinContent(bin)<<endl;
            int_sig += dimuonsGlobalInvMassVsPt[ih]->GetBinContent(bin);
            int_sig_sqr += pow(dimuonsGlobalInvMassVsPt[ih]->GetBinContent(bin),2);
        }

        if(getYield == 2) {
            int_sig =  RBWPOL->Integral(massCount_low, massCount_high)*bin_size;
            yld_cat_1[ih]   = int_sig - bin_size*bkgFit_1->Integral(massCount_low, massCount_high);
            eyld_cat_1[ih] = TMath::Sqrt(int_sig + bin_size*bkgFit_1->Integral(massCount_low, massCount_high) );
        }

        else {
            yld_cat_1[ih]   = int_sig ;
            eyld_cat_1[ih] = int_sig_sqr;
        }
        cout << "int_sig - bin_size*bkgFit_1->Integral(massCount_low, massCount_high);" << int_sig<< "  -"<< bin_size<<"*"<<bkgFit_1->Integral(massCount_low, massCount_high)<< " with low"<< massCount_low<<" high "<< massCount_high<<endl;
        //// Printing /////
        cout <<  PT[ih] << "    " << yld_cat_1[ih] << " +- " << eyld_cat_1[ih] <<"     " << GGphDimuMass << "    " << GGZ0Width  << "    " << GauWidth <<"   "<< chisq << "/" << ndf  << endl;


        // -------------- Draw
        //    dimuonsGlobalInvMassVsPt[ih]->SetMinimum(-.05*dimuonsGlobalInvMassVsPt[ih]->GetMaximum());

        if(isLog) gPad->SetLogy(1);

        TColor *pal    = new TColor();
        Int_t kblue    = pal->GetColor(9,0,200);
        //    Int_t korange  = pal->GetColor(101, 42,  0);

//from Andrius
//******************************************//
// author: Chris Seez
//
// This function evaluates the effective
// sigma of a histogram.
//
//******************************************//
Double_t StudySel::calc_effSigma(TH1 * hist)
{
  using std::cout;
  using std::endl;
  TAxis *xaxis = hist->GetXaxis();
  Int_t nb = xaxis->GetNbins();
  if(nb < 10) {
    cout << "effsigma: Not a valid histo. nbins = " << nb << endl;
    return 0.;
  }
  Double_t bwid = hist->GetBinWidth(1);
  if(bwid == 0) {
    cout << "effsigma: Not a valid histo. bwid = " << bwid << endl;
    return 0.;
  }
  //Double_t xmax = xaxis->GetXmax(); // not used
  Double_t xmin = xaxis->GetXmin();
  Double_t ave = hist->GetMean();
  Double_t rms = hist->GetRMS();
  Double_t total=0.;
  for(Int_t i=0; i<nb+2; i++) {
    total+=hist->GetBinContent(i);
  }
  if(total < 10.) {
    cout << "effsigma: Too few entries " << total << endl;
    return 0.;
  }
  Int_t ierr=0;
  Int_t ismin=999;
  Double_t rlim=0.683*total;
  Int_t nrms= Int_t(2.0*rms/(bwid)); // Set scan size to +/- rms
  //if(nrms > nb/10) nrms=nb/10; // Could be tuned...
  Double_t widmin=9999999.;
  for(Int_t iscan=-nrms;iscan<nrms+1;iscan++) { // Scan window centre
    Int_t ibm=Int_t((ave-xmin)/bwid+1+iscan);
    Double_t x=(ibm-0.5)*bwid+xmin;
    Double_t xj=x;
    Double_t xk=x;
    Int_t jbm=ibm;
    Int_t kbm=ibm;
    // ibm -- curent bin value
    // x, xj, xk -- value of the x-axis
    // jbm, kbm -- bin values
    Double_t bin=hist->GetBinContent(ibm);
    total=bin;
    for(Int_t j=1;j<nb;j++){
      if(jbm < nb) {
	jbm++;
	xj+=bwid;
	bin=hist->GetBinContent(jbm);
	total+=bin;
	if(total > rlim) break;
      }
      else ierr=1;
      if(kbm > 0) {
	kbm--;
	xk-=bwid;
	bin=hist->GetBinContent(kbm);
	total+=bin;
	if(total > rlim) break;
      }
      else ierr=1;
    }
    Double_t dxf=(total-rlim)*bwid/bin;
    Double_t wid=(xj-xk+bwid-dxf)*0.5;
    if(wid < widmin) {
      widmin=wid;
      ismin=iscan;
    }
  }
  if(ismin == nrms || ismin == -nrms) {
    cout << endl << endl << "CRAP!!!" << endl;
    cout << ismin << '\t' << nrms << endl;
    ierr=3;
  }
  if(ierr != 0 && ierr != 1) cout << "effsigma: Error of type " << ierr <<endl;
  return widmin;
}