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

void erfFit(pair<int, double> bump, vector<pair<int, double> > shower, 
	    TH1D* energyResoX, TH1D* energyResoY, TH1D* posResoX,
	    TH1D* posResoY)
 {
   TF1* erfX = new TF1("erfX", erfPulse, -5, 5, 3);
   erfX->SetParameters(IDtoX(bump.first), 1., bump.second);

   TF1* erfY = new TF1("erfY", erfPulse, -5, 5, 3); 
   erfY->SetParameters(IDtoY(bump.first), 1., bump.second);
 
   TH2D* histo = new TH2D("histo", "ShowerFit",
			  11, -5, 5,
			  11, -5, 5);


   vector<pair<int, double> >::iterator b; 
   for (b=shower.begin(); b!=shower.end(); b++)
	 {
	   histo->Fill(IDtoX(b->first), IDtoY(b->first), b->second);
	 }
   TH1D* histoX = histo->ProjectionX();
   TH1D* histoY = histo->ProjectionY();
   histoX->Fit(erfX, "Q");
   energyResoX->Fill(erfX->Integral(-5, 5)-energyGamma);
   posResoX->Fill(erfX->GetParameter(0));
   histoY->Fit(erfY, "Q");
   energyResoY->Fill(erfY->Integral(-5, 5)-energyGamma);
   posResoY->Fill(erfY->GetParameter(0));
 `
			  
 }

  /*--------------------------------------------*/
 void FitCosTheta()
 /*--------------------------------------------*/
 {

//TH1F *brgauss = new TH1F("breitg","", 131, 0, 130);
TFile * file = new TFile("cosTheta_allnadjieh.root","READ");
 TH1D * h = (TH1D*)file->Get("cosThetaGen");
//h->Scale(1./h->Integral());

 h->Sumw2();
/*TFile * file= new TFile("WPol_SelectedTTBars.root","read");
TH1D* h = (TH1D*)file->Get("costheta_mm/hCosThetaPosLepton_Gen");*/
h->Rebin(10);
//h->Scale(1./h->Integral());
 TF1 *f = new TF1("f",cosTheta, -1, 1 ,3);
 Double_t par[3];
 par[0] = 1;
 par[1] = 1;
 par[2] = h->GetEntries();
 


 f->SetParameters(par);
 
 h->Fit(f, "RBO");
 h->Fit(f, "RBO");
 h->Fit(f, "RBO");
 h->Fit(f, "RBO");
 h->Fit(f, "RBO");

 }

void FindMCToDataScaleFactor() {

  std::string data_filename = "~/data/out/v92/total.root";
  TFile* data_file = new TFile(data_filename.c_str(), "READ");
  std::string data_plotname = "TME_Al50_EvdE/all_particles/SiL_EvdE";
  TH2F* hEvdE_data = (TH2F*) data_file->Get(data_plotname.c_str());
  hEvdE_data->SetDirectory(0);
  data_file->Close();

  std::string MC_filename = "plots.root";
  TFile* MC_file = new TFile(MC_filename.c_str(), "READ");
  std::string MC_plotname = "hAll_EvdE_SiL";
  TH2F* hEvdE_MC = (TH2F*) MC_file->Get(MC_plotname.c_str());
  hEvdE_MC->SetDirectory(0);
  MC_file->Close();

  double energy_slice = 3000;
  int energy_slice_bin = hEvdE_data->GetXaxis()->FindBin(energy_slice);
  TH1D* hDataProjection = hEvdE_data->ProjectionY("_py", energy_slice_bin, energy_slice_bin);  
  energy_slice_bin = hEvdE_MC->GetXaxis()->FindBin(energy_slice);
  TH1D* hMCProjection = hEvdE_MC->ProjectionY("_py", energy_slice_bin, energy_slice_bin);

  TFitResultPtr data_fit_result = hDataProjection->Fit("gaus", "S");
  TFitResultPtr MC_fit_result = hMCProjection->Fit("gaus", "S");

  hDataProjection->SetLineColor(kBlack);
  hDataProjection->GetXaxis()->SetRangeUser(0, energy_slice);
  hDataProjection->Draw();
  hMCProjection->SetLineColor(kRed);
  hMCProjection->GetXaxis()->SetRangeUser(0, energy_slice);
  hMCProjection->Draw("SAME");

  //EvdE_data->Draw("COLZ");
  double data_mean = data_fit_result->Parameter(1);
  double MC_mean = MC_fit_result->Parameter(1);
  double scale_factor = data_mean / MC_mean;
  std::cout << "data / MC = " <<  data_mean << " / " << MC_mean << " = " << scale_factor << std::endl;

  // Go through the MC projection and scale each bin by the scale factor
  int n_bins = hMCProjection->GetNbinsX();
  double min_x = hMCProjection->GetXaxis()->GetXmin();
  double max_x = hMCProjection->GetXaxis()->GetXmax();
  TH1F* hMCProjection_scaled = new TH1F("hMCProjection_scaled", "", n_bins,min_x,max_x);

  for (int i_bin = 1; i_bin <= n_bins; ++i_bin) {
    double old_energy = hMCProjection->GetBinLowEdge(i_bin);
    double new_energy = old_energy * scale_factor;

    double old_bin_content = hMCProjection->GetBinContent(i_bin);
    hMCProjection_scaled->Fill(new_energy, old_bin_content);
  }

  hMCProjection_scaled->Draw();
  hDataProjection->Draw("SAME");
}

TH2D* GetJetCorrFunc2D_ZYAM(int itrg, int jass)
{
                    TH2D* hcorr = (TH2D*)GetRawCorrFunc2D_ratio(itrg,jass);
                    TH2D* hcorr_clone = (TH2D*)hcorr->Clone(Form("corr_clone_itrg%d_jass%d",itrg,jass));
                    hcorr_clone->Reset();
		    for(int ietabin=1;ietabin<=hcorr->GetNbinsX();ietabin++)
		    {
                      TH1D* hcorrphi = (TH1D*)hcorr->ProjectionY(Form("corrphi_%d",ietabin),ietabin,ietabin,"e");  
                     // hcorrphi->SetAxisRange(0.35,1.57,"X");
                      double histminX = hcorrphi->GetBinCenter(hcorrphi->GetMinimumBin());
                      double histminY = hcorrphi->GetBinContent(hcorrphi->GetMinimumBin());
	   	      fitfunc = new TF1(Form("fitfunc_%d",ietabin),"[0]+[1]*(x-[2])*(x-[2])",0.6,1.2);		    
                      fitfunc->SetParameters(histminY,1,histminX);
                      fitfunc->SetParLimits(1,0,0.10000);
                  //    fitfunc->SetParLimits(2,0.35,1.57);
		      for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%d",ietabin),"RNO");                     		  
		      float level = fitfunc->GetParameter(0);
		      for(int iphibin=1;iphibin<=hcorr->GetNbinsY();iphibin++) 
                        hcorr_clone->SetBinContent(ietabin,iphibin,hcorr->GetBinContent(ietabin,iphibin)-level);
		    }
                    
                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");		    
                    return hcorr_clone;
}

void fitSlices(TH2* hCorr, TF1* func){

   int nBins = hCorr->GetNbinsX();

   TH1D* hMean = new TH1D(Form("%s_1",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());
   TH1D* hSigma = new TH1D(Form("%s_2",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());

   for(int i = 1; i < nBins+1; ++i){
      int bin = nBins - i;
      TH1D* h = hCorr->ProjectionY(Form("%s_bin%d",hCorr->GetName(),bin),i,i);

      func->SetParameter(0,h->GetMaximum());
      func->SetParameter(1,h->GetMean());
      func->SetParameter(2,h->GetRMS());

      if(useFits) h->Fit(func);

      hMean->SetBinContent(i,func->GetParameter(1));
      hMean->SetBinError(i,func->GetParError(1));
      hSigma->SetBinContent(i,func->GetParameter(2));
      hSigma->SetBinError(i,func->GetParError(2));
      
      if(onlySaveTable){
	 h->Delete();
      }
   }
}

//------------------------------------------------------------------------
void fitSlices(TH2* hCorr, TF1* func){

   int nBins = hCorr->GetNbinsX();

   TH1D* hMean = new TH1D(Form("%s_1",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());
   TH1D* hSigma = new TH1D(Form("%s_2",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());

   for(int i = 1; i < nBins+1; i++){
      int bin = nBins - i;
      TH1D* h = hCorr->ProjectionY(Form("%s_bin%d",hCorr->GetName(),bin),i,i);

      func->SetParameter(0,h->GetMaximum());
      func->SetParameter(1,h->GetMean());
      func->SetParameter(2,h->GetRMS());

      h->Fit(func,"Q");

      /*hMean->SetBinContent(i,func->GetParameter(1));
      hMean->SetBinError(i,func->GetParError(1));
      hSigma->SetBinContent(i,func->GetParameter(2));
      hSigma->SetBinError(i,func->GetParError(2));*/

      hMean->SetBinContent(i,h->GetMean());
      hMean->SetBinError(i,func->GetParError(1)); //errors are not use later for the actual table
      hSigma->SetBinContent(i,h->GetRMS());
      hSigma->SetBinError(i,func->GetParError(2));

   }
}

void tofQAMaker( string file, int xRange = 10 ){

	TFile * f = new TFile( file.c_str(), "READ" );


	TH1D * tpc = (TH1D*)f->Get( "zvertex" );
	TH1D * vpd = (TH1D*)f->Get( "zvertexVPD" );
	TH1D * delta = (TH1D*)f->Get( "zvertexDelta" );

	Reporter rp( "vertexQA.pdf" );

	TF1* gaus = new TF1( "g", "gaus" );
	gaus->SetRange( -xRange, xRange );

	rp.newPage();
	tpc->Draw();
	rp.savePage();
	rp.newPage();
	vpd->Draw();
	rp.savePage();
	rp.newPage();
	delta->GetXaxis()->SetRangeUser( -25, 25 );
	gStyle->SetOptFit( 111 );
	delta->Draw();
	delta->Fit( gaus, "R" );
	rp.savePage();

}

void Draw_KL_Test(){

  TChain* ch = new TChain("Tree");
  TChain* ch1 = new TChain("Tree");
  
  TH1D* his  = new TH1D("Klong6g","Klong6g",20,450,550);
  TH1D* his1 = new TH1D("Klong4g","Klong4g",20,450,550);
  TH1D* his2 = new TH1D("Klong4gAll","Klong4gAll",60,250,550);

  for( int i = 0; i< 68; i++){
    ch->Add(Form("klongRootFile/kl%d.root" ,4162+i));
    ch1->Add(Form("klongRootFile/ks%d.root",4162+i));
  }
  ch->Project(his->GetName()  ,"KlongMass[0]","CutCondition==0");
  ch1->Project(his1->GetName(),"KlongMass[0]","CutCondition==0");
  ch1->Project(his2->GetName(),"KlongMass[0]","CutCondition==0");
  
  TF1* func = new TF1("func","gaus(0)+expo(3)",0,550);
  func->SetParameter(1,498);
  func->SetParameter(2,5);
  TF1* func2 = new TF1("func2","gaus(0)",0,550);
  func2->SetParameter(1,498);
  func2->SetParameter(2,5);



  TCanvas* can = new TCanvas("can","",1200,600);
  can->Divide(2,1);
  can->cd(1);
  his2->Fit(func->GetName(),"","",450,550);
  his2->Draw();
  TF1* func1 = new TF1("Test","gaus",450,550);
  func1->SetParameter(0,func->GetParameter(0));
  func1->SetParameter(1,func->GetParameter(1));
  func1->SetParameter(2,func->GetParameter(2));

  can->cd(2);
  his1->SetLineColor(2);  
  his->Draw();
  his->Fit(func2->GetName(),"","",450,550);
  func->Draw("same");
  his1->Draw("same");

  std::cout<< func2->GetParameter(0) << " " 
	   << func->GetParameter(0)  << " " 
	   << func->GetParameter(0)/func2->GetParameter(0)<< std::endl;
  std::cout<< func2->Integral(450,550) << " " 
	   << func1->Integral(450,550)  << " " 
	   << func1->Integral(450,550)/func2->Integral(450,550)
	   << std::endl;
  //ch->Draw("KlongPt[0]:KlongMass[0]>>(400,200,600,50,0,20)","(CutCondition&(1|2|4|8))==0","colz");
  gPad->SetLogz();
  TText* text = new TText(0.5,0.5,"");
  TText* text1 = new TText(0.5,0.5,"");
  text->DrawTextNDC(0.5,0.5,Form("Integral:%2.3lf",func1->Integral(450,550)));
  text1->DrawTextNDC(0.5,0.6,Form("Integral:%2.3lf",func2->Integral(450,550)));
  
}

void pixelHistos(const char *filenames = "/export/data1/dambach/hardware/ntpls/gradeC/*.root", const char *cuts =  "defects == 0")
{

  Init();      // do not forget to load the Utilities.C

  TH1D *gainH = new TH1D("gainH", "Gain", 1000, 0., 10.);
  TH1D *pedH = new TH1D("pedH", "Pedestal", 1000, -1000., 1000.);
  TH1D *noiseH = new TH1D("noiseH", "Noise", 1100, -100., 1000.);
  TH1D *tthrH = new TH1D("tthrH", "Trimmed Threshold", 100, 0., 100.);
  
  TChain *c = new TChain("mod");
  c->Add(filenames);
  
  canvas->Clear();
  canvas->Divide(2,2);

  canvas->cd(1);
  InitPad("log");
  gainH->GetXaxis()->SetTitle("ADC units/DAC units");
  gainH->GetYaxis()->SetTitle("# pixels");
  gainH->GetXaxis()->SetTitleSize(0.055);
  gainH->GetYaxis()->SetTitleSize(0.055);
  c->Draw("gain>>gainH", cuts);

  TBox *box = new TBox();

  canvas->cd(2);
  InitPad("log");
  pedH->GetXaxis()->SetTitle("DAC units");
  pedH->GetYaxis()->SetTitle("# pixels");
  pedH->GetXaxis()->SetTitleSize(0.055);
  pedH->GetYaxis()->SetTitleSize(0.055);
  c->Draw("ped>>pedH", cuts);

  canvas->cd(3);
  InitPad("log");
  noiseH->GetXaxis()->SetTitle("e^{-}");
  noiseH->GetYaxis()->SetTitle("# pixels");
  noiseH->GetYaxis()->SetRangeUser(0.1,400000);
  noiseH->GetXaxis()->SetTitleSize(0.055);
  noiseH->GetYaxis()->SetTitleSize(0.055);
  c->Draw("noise>>noiseH", cuts);
  noiseH->Fit("gaus");
//   box->SetFillStyle(3003);
//   box->SetFillColor(kBlue);
//   box->DrawBox(0.1,0.,500.,400000.);
//   box->SetFillColor(kRed);
//   box->DrawBox(500.,0.,1000.,400000.);

  canvas->cd(4);
  tthrH->GetXaxis()->SetTitle("DAC units");
  tthrH->GetYaxis()->SetTitle("# pixels");
  tthrH->GetXaxis()->SetTitleSize(0.055);
  tthrH->GetYaxis()->SetTitleSize(0.055);
  InitPad("log");
  c->Draw("tthr>>tthrH", cuts);

}

void DynamicExec()
{
   // Example of function called when a mouse event occurs in a pad.
   // When moving the mouse in the canvas, a second canvas shows the
   // projection along X of the bin corresponding to the Y position
   // of the mouse. The resulting histogram is fitted with a gaussian.
   // A "dynamic" line shows the current bin position in Y.
   // This more elaborated example can be used as a starting point
   // to develop more powerful interactive applications exploiting CINT
   // as a development engine.
   //
   // Author:  Rene Brun
   
   TObject *select = gPad->GetSelected();
   if(!select) return;
   if (!select->InheritsFrom("TH2")) {gPad->SetUniqueID(0); return;}
   TH2 *h = (TH2*)select;
   gPad->GetCanvas()->FeedbackMode(kTRUE);

   //erase old position and draw a line at current position
   int pyold = gPad->GetUniqueID();
   int px = gPad->GetEventX();
   int py = gPad->GetEventY();
   float uxmin = gPad->GetUxmin();
   float uxmax = gPad->GetUxmax();
   int pxmin = gPad->XtoAbsPixel(uxmin);
   int pxmax = gPad->XtoAbsPixel(uxmax);
   if(pyold) gVirtualX->DrawLine(pxmin,pyold,pxmax,pyold);
   gVirtualX->DrawLine(pxmin,py,pxmax,py);
   gPad->SetUniqueID(py);
   Float_t upy = gPad->AbsPixeltoY(py);
   Float_t y = gPad->PadtoY(upy);

   //create or set the new canvas c2
   TVirtualPad *padsav = gPad;
   TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
   if(c2) delete c2->GetPrimitive("Projection");
   else   c2 = new TCanvas("c2","Projection Canvas",710,10,700,500);
   c2->SetGrid();
   c2->cd();

   //draw slice corresponding to mouse position
   Int_t biny = h->GetYaxis()->FindBin(y);
   TH1D *hp = h->ProjectionX("",biny,biny);
   hp->SetFillColor(38);
   char title[80];
   sprintf(title,"Projection of biny=%d",biny);
   hp->SetName("Projection");
   hp->SetTitle(title);
   hp->Fit("gaus","ql");
   hp->GetFunction("gaus")->SetLineColor(kRed);
   hp->GetFunction("gaus")->SetLineWidth(6);
   c2->Update();
   padsav->cd();
}

void ReFit(int lo, int hi) {
	TF1* floo = (TF1*) gROOT->FindObject("floo");
	if (floo == 0) return;
	TH1D* hCe = (TH1D*) gROOT->FindObject("hCe");
	if (hCe == 0) return;
	hCe->GetXaxis()->SetRangeUser(lo,hi);
	double max = hCe->GetMaximum();
	int binmax = hCe->GetMaximumBin();
	double xmax = hCe->GetXaxis()->GetBinCenter(binmax);
	floo->SetParameters(max,xmax,10,max*0.1,50,20,5);
	hCe->Fit(floo,"Q","",lo,hi);
}

void EfficiencyFit(TString inputdata="ROOTfiles/MCstudiesPP_Fine.root", TString label="PP_Fine"){

  gStyle->SetTextSize(0.05);
  gStyle->SetTextFont(42);
  gStyle->SetPadRightMargin(0.043);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadBottomMargin(0.145);
  gStyle->SetTitleX(.0f);

  TFile* inf = new TFile(inputdata.Data());
  TH1D* hEff = (TH1D*)inf->Get("hEff");
  TCanvas* c= new TCanvas(Form("c"),"",600,600);
  hEff->Draw();


  TF1 *myfit = new TF1("myfit","[0]+[1]*x+x*x*[2]+x*x*x*[3]+x*x*x*x*[4]", 7, 50);
  //TF1 *myfit = new TF1("myfit","[0]+[1]*x+x*x*[2]+x*x*x*[3]", 7, 50);
  //TF1 *myfit = new TF1("myfit","[0]+[1]*x+x*x*[2]", 7, 50);
  hEff->Fit("myfit","L q","",7,50);
  hEff->Fit("myfit","L q","",7,50);
  hEff->Fit("myfit","L q","",7,50);
  hEff->Fit("myfit","L m","",7,50);
  double par0=myfit->GetParameter(0);
  double par1=myfit->GetParameter(1);
  double par2=myfit->GetParameter(2);
  double par3=myfit->GetParameter(3);
  double par4=myfit->GetParameter(4);
  TString myEffWeight=Form("(%e+%e*Bpt+Bpt*Bpt*%e+Bpt*Bpt*Bpt*%e+Bpt*Bpt*Bpt*Bpt*%e)",par0,par1,par2,par3,par4);
  std::cout<<myEffWeight<<std::endl;

  c->SaveAs(Form("plotEffFit/%s.png", label.Data()));
  TFile *fout=new TFile(Form("plotEffFit/%s.root", label.Data()),"recreate");
  fout->cd();
  hEff->Write();
  myfit->Write();
  fout->Close();
}

void RRootHistRead_v1(){
	
	TFile *inputFile = new TFile("pPb_MBSpectra_Combine_-1_1.root");
	cout << "Is the file Open " << inputFile->IsOpen() << endl;

	//Copy the histogram
	TString hName = "Spectra_NtrkOffline0_inf";
	TH1D *h = (TH1D*)inputFile->Get(hName);

	cout << "Number of bins in the histogram: " << h->GetSize() << "\n";
	h->Draw();

	double nhBins = h->GetSize();
	for (int i=0; i < nhBins; i++) cout << h->GetBinWidth(i) << ", ";
	cout << endl;

	// Fit function
	double nParam = 2, maxValue = 5, minValue = 0;

	TF1 *fitFn = new TF1("fitFn", fitFunction, minValue, maxValue, nParam);

	// Set parameters for Fit
	fitFn->SetParameter(0, 0.008);

	// Get parameters from Fit
	double paramFromFit[2];
	fitFn->GetParameters(paramFromFit);

	// Generate data from fit
	double nBins = 27, pT[27], spectraFromFit[27]; // From the Xi 
	for (int i=0; i < nBins; i++){

		spectraFromFit[i] = paramFromFit[0]*exp(-i) + paramFromFit[1];
		cout << spectraFromFit[i] << ", ";
		pT[i] = i;
			
	}

	// Using Integral for the Fit
	h->Fit(fitFn, "I");

	TGraph *g = new TGraph(nBins, pT, spectraFromFit);
	g->SetMarkerStyle(7);
	g->SetMarkerSize(1.5);
	g->SetMarkerColor(kGreen);
	g->Draw("ALP");
	h->Draw("same");
	
}

TF1* GetFitFunc_ZYAM_MC(TH1D* h)
{
                    TH1D* hcorrphi = (TH1D*)h->Clone(h->GetName());
                    double histminY = hcorrphi->GetBinContent(10);
                    double histminX = 1.0;

                    //hcorrphi->SetAxisRange(-0.01,1.2,"X");
                    TF1* fitfunc = new TF1(Form("fitfunc_%s",h->GetName()),"[0]+[1]*(x-[2])*(x-[2])",0.8,2.8); //std 0.15-1.8
                    fitfunc->SetParameters(histminY,0.0002,histminX);
                    fitfunc->SetParLimits(1,0,1000);
                    fitfunc->SetParLimits(2,0.05,1000);

                    for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%s",h->GetName()),"RNO");
                    return fitfunc;
}

TF1* GetFitFunc_ZYAM_AllPhysics(TH1D* h)
{
                    TH1D* hcorrphi = (TH1D*)h->Clone(h->GetName());
                    double histminY = hcorrphi->GetBinContent(10);
                    double histminX = 1.0;

                    hcorrphi->SetAxisRange(-0.01,1.5,"X");
                    TF1* fitfunc = new TF1(Form("fitfunc_%s",h->GetName()),"[0]+[1]*(x-[2])*(x-[2])+[3]*(x-[2])*(x-[2])*(x-[2])",0.5,1.65); //std 0.6 1.55 vs pT ; 0.6 1.8 vs eta
                    fitfunc->SetParameters(histminY,0.0002,histminX,0.0001);
                    fitfunc->SetParLimits(1,0,1000);
                    fitfunc->SetParLimits(2,0.5,1000);
                  //  fitfunc->SetParLimits(3,0,1000);

                    for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%s",h->GetName()),"RNO");
                    return fitfunc;
}