C++ TF1::SetParameters方法代码示例

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);
}

	  TF1 * tripleGaussFit(TH1 * hist, bool autorange, double inner=0.01, double outer=0.1){
		  /**
		   * @brief performes a double gaussian Fit
		   */
		  Double_t parameters[9] = {0,0,0,0,0,0,0,0,0};


		  double center = hist->GetMean();


		  if (autorange){
			  double inner = fabs(center - hist->GetXaxis()->GetXmax()) /10;
			  double outer = fabs(center - hist->GetXaxis()->GetXmax()) /10*8;
		  }
		  double middle = (outer - inner)/2;

		  TF1 * fit1 = new TF1("fit1", "gaus", center-inner, center+inner );

		  hist->Fit(fit1, "0R");
		  fit1->GetParameters(&parameters[0]);

		  TF1 * fitmid = new TF1("fitmid", "gaus", center-middle, center+middle );
		  fitmid->SetParameter(1, parameters[1]);
		  hist->Fit(fitmid, "0R");
		  fitmid->GetParameters(&parameters[3]);

		  TF1 * fit2 = new TF1("fit2", "gaus", center-outer, center+outer );
		  fit2->SetParameter(1, parameters[1]);
		  hist->Fit(fit2, "0R");
		  fit2->GetParameters(&parameters[6]);

		  TF1 * fitproper = new TF1("fitproper", "gaus(0)+gaus(3)+gaus(6)", center-outer, center+outer);
		  fitproper->SetParameters(parameters);
		  fitproper->SetParName(0, "Const.(inner)");
		  fitproper->SetParName(1, "Mean (inner)");
		  fitproper->SetParName(2, "Sigma (inner)");
		  fitproper->SetParName(3, "Const.(middle)");
		  fitproper->SetParName(4, "Mean (middle)");
		  fitproper->SetParName(5, "Sigma (middle)");
		  fitproper->SetParName(6, "Const.(outer)");
		  fitproper->SetParName(7, "Mean (outer)");
		  fitproper->SetParName(8, "Sigma (outer)");

		  hist->Fit(fitproper, "0R");

//		  fitproper->SetLineColor(hist->GetLineColor());
//		  fitproper->SetLineWidth(hist->GetLineWidth());
////		  fitproper->SetLineStyle(2);

		  return fitproper;

	  }

TF1 * residualFit(double pars[], const std::string & index)
{
  ResidualFit * fobj = new ResidualFit;
  TF1 * functionToFit = new TF1(("fp"+index).c_str(), fobj, 60, 120, fobj->parNum(), "ResidualFit");

//   double pars[] = { -2.30972e-02, 7.94253e-01, 7.40701e-05, 2.79889e-06,
//                     -2.04844e-08, -5.78370e-43, 6.84145e-02, -6.21316e+00,
//                     -6.90792e-03, 6.43131e-01, -2.03049e-03, 3.70415e-05,
//                     -1.69014e-07 };
  functionToFit->SetParameters(pars);

  return functionToFit;
}

float getBkgErr(float mass)
{
  static bool  isInit=false;
  static TF1 ebErr2Func("ebErrFunc","pol4",0,3000);
  static TF1 eeErr2Func("eeErrFunc","pol4",0,3000);
  static TF1 pdfFunc("pdfFunc","pol3",0,3000);
  static TF1 ewkFunc("ewkFunc","pol2",0,3000);
  static TF1 muonErrFunc("muonEffFunc","pol2",0,3000);
  if(!isInit){
    ebErr2Func.SetParameters(19.2,6.79E-3,4.32E-5,9.81E-9,7.18E-12);
    eeErr2Func.SetParameters(21.2,6.79E-3,4.32E-5,9.81E-9,7.18E-12);
    pdfFunc.SetParameters(4.15,1.83E-3,2.68E-6);
    ewkFunc.SetParameters(-1,4.2E-3);
    muonErrFunc.SetParameters(1.045,1.043E-5,3.378E-8);
  }
  
  //if(isMuon){
  //float pdfErr = pdfFunc.Eval(mass);
  //    float ewkErr = ewkFunc.Eval(mass);
  return muonErrFunc.Eval(mass)-1;
  //}
}

Double_t calc_dNdY(const char *data_file = "pt_RFE", Bool_t show = kFALSE)
{
  TGraphErrors *g = new TGraphErrors(data_file);
  if (g->IsZombie()) return;

  TF1 *flt = new TF1("flt", LevyTsallis, 0., 5., 4);
  flt->SetParameters(2.93483e-02, 2.80382e-01, 8.10224e+00, 1.01944e+00);
  flt->FixParameter(3, 1.019445);

  Double_t fitmin = 0.25, fitmax = 5.25; // !!!
  g->Fit(flt, "Q", "", fitmin, fitmax);
  g->Fit(flt, "Q", "", fitmin, fitmax);
  g->Fit(flt, "Q", "", fitmin, fitmax);

  Double_t first, last;
  // !!! bining sensibility !!!
  Int_t graph_range[2] = {2, 23}; // [2] = 0.05 - 0.01, [23] = 4.70 + 0.01

  // interval 1
  first = 0.0;   // !!!
  last  = g->GetX()[graph_range[0]] - g->GetErrorX(graph_range[0]);
  Double_t fi1 = flt->Integral(first, last);

  // interval 2
  first = last;
  last  = g->GetX()[graph_range[1]] + g->GetErrorX(graph_range[1]);
  Double_t fi2 = flt->Integral(first, last);
  Double_t gi2 = 0.0;
  for (Int_t ip = graph_range[0]; ip <= graph_range[1]; ip++)
    gi2 += g->GetY()[ip]*g->GetErrorX(ip)*2.0;

  // interval 3
  first = last;
  last  = 30.0;   // !!!
  Double_t fi3 = flt->Integral(first, last);

  Double_t result = fi1 + gi2 + fi3;

  if (!show) {
    delete g;
    delete flt;
    return result;
  }

  Printf("function: %f \t %f \t %f", fi1, fi2, fi3);
  Printf("graph: \t \t \t %f", gi2);
  Printf("result => %f", result);
  g->Draw("AP");
  gPad->SetLogy();
  return result;
}

void fitCauchy()
{

  TCanvas * c1 = new TCanvas("mycan","mycan",800,600);
  TFile * data = new TFile("~/unpacker/sort.root");
  fstream file1;
  file1.open("gaus3.par");


  TH1F * dist = (TH1F*)data->Get("/corr/Li7/Li7_theta_reactCoM_12C")->Clone("1");

  Double_t fparams[9] = {6.,1.,12.,1.,18.,1.,20000,15000,10000};
  TF1 * f = new TF1("f",Cauchy3,0,180,9);
  TF1 * f2 = new TF1("f2",Gaus3,0,180,9);
  f2->SetParameters(fparams);
  f->SetParameters(fparams);
  dist->Fit("f2","M0");
  dist->Draw();
  f2->Draw("same");

  if(file1.is_open()) cout <<"file1 is open" << endl;
  else cout << "failed to load file1" << endl;

   for(int i=0;i<9;i++)
     {
       file1 << f2->GetParameter(i) << endl;
     }

  TCanvas * c2 = new TCanvas("mycan2","mycan2",800,600);
  TH1F * check = new TH1F("check","",125,0,40);

  for(int i =0;i<300000;i++)
    {
      float  r = f2->GetRandom();
      check->Fill(r);
    }

  check->Draw();
    
  Double_t  par[9];

  //  for(int i=0;i<9;i++)
  //{
  // file1 >> par[i] >> endl;
  //  cout << "par[" << i << "] = " << par[i] << endl;
  //}
  //cout << "what" << endl;
  //file1.close();

  return;
}

TF1* fitTools::fitResponseGraph( TGraphErrors* graph, std::string funcType, std::string funcName, const std::string& option, float rangeMax, float rangeMin) {

  
  TF1* fitFunction;

  if( funcType=="rpf" ) {
    fitFunction = new TF1(funcName.c_str(), rpf, rangeMin, rangeMax, 6); //will have to fix the range issue!
    fitFunction->SetParameters(100,0.85,4.2,80,250,1.);
    fitFunction->SetParLimits(1,0.5,1.0);
    fitFunction->SetParLimits(2,1.,10.);
    fitFunction->SetParLimits(4, 100., 500.);
    fitFunction->SetParameter(0, 0.6);
    fitFunction->SetParameter(5, 1.);
    //fitFunction->SetParameter(1, -0.6);
  } else if( funcType=="powerlaw" ) {
    fitFunction = new TF1(funcName.c_str(), "[0] - [1]/(pow(x, [2]))");
    fitFunction->SetRange(rangeMin, rangeMax);
    fitFunction->SetParameters(1., 1., 0.3);
  } else if( funcType=="powerlawL2L3" ) {
    fitFunction = new TF1(funcName.c_str(), "1. + [0]/(pow(x, [1]))");
    fitFunction->SetRange(rangeMin, rangeMax);
    fitFunction->SetParameters(0.7, 0.7);
  //fitFunction->SetParLimits(0, 0.5, 2.);
  //fitFunction->SetParLimits(1, 0.5, 2.);
  } else if( funcType=="powerlaw_corr" ) {
    fitFunction = new TF1(funcName.c_str(), "[0] - [1]/(pow(x, [2])) + [3]/x");
    fitFunction->SetRange(rangeMin, rangeMax);
    fitFunction->SetParameters(1., 1., 0.3, 1.);
  } else {
    std::cout << "Function '" << funcType << "' not implemented yet for fitResponseGraph. Exiting." << std::endl;
    exit(119);
  }

  graph->Fit(fitFunction, option.c_str());

  return fitFunction;

}

double findDelPhiAngle(TH1D *input, int parorErr){
  
  TF1 *fitter = new TF1("fitter","gaus",-2,2);
  //input->Scale(1./(input->GetEntries()));
  fitter->SetParameters(1./input->GetMaximum(),-0.1,1.0);
  input->Fit("fitter","qN0","",-1.5,1.5);
  if(parorErr==1){
    return fitter->GetParameter(1);
  }
  else if(parorErr==2){
    return fitter->GetParError(1);
  }
  else return 0;
}

void qa2() {
   //Fill a 1-D histogram from a parametric function
   TCanvas *c1 = new TCanvas("c1","The FillRandom example",0,0,700,500);
   c1->SetFillColor(18);

   gBenchmark->Start("fillrandom");
   //
   // A function (any dimension) or a formula may reference
   // an already defined formula
   //
   TFormula *form1 = new TFormula("form1","abs(sin(x)/x)");
   TF1 *sqroot = new TF1("sqroot","x*gaus(0) + [3]*form1",0,10);
   sqroot->SetParameters(10,4,1,20);

   //
   // Create a one dimensional histogram (one float per bin)
   // and fill it following the distribution in function sqroot.
   //
   TH1F *h1f = new TH1F("h1f","Test random numbers",200,0,10);
   h1f->SetFillColor(45);
   h1f->FillRandom("sqroot",100000);
   h1f->Draw();
   TPaveLabel *lfunction = new TPaveLabel(5,39,9.8,46,"The sqroot function");
   lfunction->SetFillColor(41);

   c1->SetGridx();
   c1->SetGridy();
   c1->GetFrame()->SetFillColor(42);
   c1->GetFrame()->SetBorderMode(-1);
   c1->GetFrame()->SetBorderSize(5);

   h1f->SetDirectory(0);

   c1->Update();

   sqroot->SetParameters(200,4,1,20);
}

TF1* fitPixelBeam(TH1D* hist, double pixelWidth, double beamSigma, bool display)
{
  TF1* conv = new TF1("conv",
                      "[2] + [4] * 0.5 * ( TMath::Erf(([0]/2 + [3] - x)/(sqrt(2)*[1])) + TMath::Erf(([0]/2 + [3] + x)/(sqrt(2)*[1])) )");
  conv->SetParNames("Width", "Sigma", "Background", "Offset", "Scale");

  const unsigned int numBins = hist->GetNbinsX();
  const double limitLow = hist->GetXaxis()->GetBinLowEdge(1);
  const double limitHigh = hist->GetXaxis()->GetBinUpEdge(numBins);

  // Estimate the background using the +/- 10% edges of the histogram
  double background = 0;
  const unsigned int numBackgroundBins = numBins * 0.1 + 1;
  for (unsigned int n = 0; n < numBackgroundBins; n++)
  {
    background += hist->GetBinContent(n + 1);
    background += hist->GetBinContent(numBins - n);
  }
  background /= (double)(2 * numBackgroundBins);

  // Estimate the scale
  double scale = 0;
  for (unsigned int bin = 1; bin <= numBins; bin++)
    if (hist->GetBinContent(bin) > scale) scale = hist->GetBinContent(bin);

  double offset = 0;

  conv->SetRange(limitLow, limitHigh);
  conv->SetParameters(pixelWidth, beamSigma, background, offset, scale);
  conv->SetParLimits(0, 0.1 * pixelWidth, 100 * pixelWidth);
  conv->SetParLimits(1, 0, 100 * beamSigma);
  conv->SetParLimits(2, 0, scale);
  conv->SetParLimits(3, -pixelWidth, pixelWidth);
  conv->SetParLimits(4, 0.1 * scale, 10 * scale);
  hist->Fit("conv", "QR0B");

  if (display)
  {
    TCanvas* can = new TCanvas();
    conv->SetLineColor(46);
    conv->SetLineWidth(1);
    hist->Draw();
    conv->Draw("SAME");
    can->Update();
    can->WaitPrimitive();
  }

  return conv;
}

void birks()
{
 TF1 *bf4  = new TF1("birks4",birksf,0,100,2);
 bf4->SetParameters(4.6345e-2,0.0);
 bf4->SetLineWidth(2); 
 bf4->SetLineColor(7);
 bf4->GetXaxis()->SetTitle("-dE/dx  [MeV/g cm^{2}]");
 bf4->GetYaxis()->SetTitle("Birks Factor");
 bf4->Draw();
 TF1 *bf  = new TF1("birks",birksf,0,100,2);
 bf->SetParameters(1.29e-2,9.59e-6);
 bf->SetLineWidth(2); 
 bf->SetLineColor(3);
 bf->Draw("SAME");
 TF1 *bf1  = new TF1("birks1",birksf,0,100,2);
 bf1->SetParameters(8.35e-3,0.0);
 bf1->SetLineWidth(2); 
 bf1->SetLineColor(4);
 bf1->Draw("SAME");
 TF1 *bf2  = new TF1("birks2",birksf,0,100,2);
 bf2->SetParameters(1.31e-2,0.0);
 bf2->SetLineWidth(2); 
 bf2->SetLineColor(2);
 bf2->Draw("SAME");
 TF1 *bf3  = new TF1("birks3",birksf,0,100,2);
 bf3->SetParameters(1.59e-2,0.0);
 bf3->SetLineWidth(2); 
 bf3->SetLineColor(6);
 bf3->Draw("SAME");
 TF1 *bf5  = new TF1("birks5",birksf,0,100,2);
 bf5->SetParameters(6.8552e-4,0.0);
 bf5->SetLineWidth(2); 
 bf5->SetLineColor(21);
 bf5->Draw("SAME");
 TF1 *bf6  = new TF1("birks6",birksf,0,100,2);
 bf6->SetParameters(3.5175e-3,0.0);
 bf6->SetLineWidth(2); 
 bf6->SetLineColor(8);
 bf6->Draw("SAME");
 TLegend *legend = new TLegend(.44, .64, .99, .99);
 legend->AddEntry(bf,  "CMS/ATLAS");
 legend->AddEntry(bf2, "Standard Birks law.");
 legend->AddEntry(bf1, "Zeus SCSN38, lower limit (-35%) IEEE TNS Vol. 39 NO.4 (1992), 511");
 legend->AddEntry(bf3, "Zeus SCSN38, upper limit (+21%) IEEE TNS Vol. 39 NO.4 (1992), 511");
 legend->AddEntry(bf4, "BGO, NIM A439 (2000) 158-166");
 legend->AddEntry(bf5, "CsI(Tl), NIM A439 (2000) 158-166");
 legend->AddEntry(bf6, "GSO(Ge), NIM A439 (2000) 158-166");
 legend->Draw();
}

int main(int argc, char const *argv[]) {
	TFile *myFile = new TFile("output.root","RECREATE");

	Float_t E1,Px1,Py1,Pz1,E2,Px2,Py2,Pz2;
	Double_t mass = 0;

	TH1D *m12 = new TH1D("m12","m12",60,0.99,1.08);
	TF1 *myBW = new TF1("myBW","TMath::BreitWigner(x,[0],[1])", 1.01, 1.03);

	myBW->SetParName(0,"Mass");
	myBW->SetParName(1,"#Gamma");

	TFile *inputFile = new TFile("data.root");
	TTree *data_tree = (TTree*)inputFile->Get("data");

	data_tree->SetBranchAddress("E1", &E1);
	data_tree->SetBranchAddress("Px1", &Px1);
	data_tree->SetBranchAddress("Py1", &Py1);
	data_tree->SetBranchAddress("Pz1", &Pz1);
	data_tree->SetBranchAddress("E2", &E2);
	data_tree->SetBranchAddress("Px2", &Px2);
	data_tree->SetBranchAddress("Py2", &Py2);
	data_tree->SetBranchAddress("Pz2", &Pz2);

	TLorentzVector vec1(0.0,0.0,0.0,0.0);
	TLorentzVector vec2(0.0,0.0,0.0,0.0);

	int count = data_tree->GetEntries();

	for (int i = 0; i < count; i++) {	
		data_tree->GetEntry(i);
		vec1.SetPxPyPzE(Px1,Py1,Pz1,E1);
		vec2.SetPxPyPzE(Px2,Py2,Pz2,E2);
		mass = (vec1 + vec2).M();

		m12->Fill(mass);
	}
	myBW->SetParameters(0,1.02,0.05);
	myBW->SetParLimits(0,1.01,1.03);
	m12->Fit("myBW","","",0,2);
	m12->GetXaxis()->SetTitle("Mass (GeV)");
	gStyle->SetOptFit(1111);

	myFile->cd();
	myFile->Write();
	myFile->Close();

	return 0;
}

  void XZFit(){   
    if (xerr[0] == 0 && xerr[1] == 0 && xerr[2] == 0 && zerr[0] == 0 && zerr[1] == 0 && zerr[2] == 0) {cout << "ERRORE: Incertezze non specificate" << endl;}
    else {
      //  TGraphErrors *graph1 = new TGraphErrors(3,x,y,xerr,yerr);
    //  TF1* fitfunc1 = new TF1("fittingfunction","[0]+[1]*x",-100,100);
    //   xyfitfunc = fitfunc1;
    // xygraph = graph1;
      delete xzfitfunc;  //Dealloca prima di riallocarne uno nuovo
      delete xzgraph;
xzfitfunc = new TF1("xzfittingfunction", "[0]+[1]*x",-100,100);
        xzgraph = new TGraphErrors(3,xv,yv,xerr,zerr);
    xzfitfunc->SetParameters(30,-2);
    xzgraph->Fit(xzfitfunc,"0S");
    }
  };

void fitExclude() {
   //Create a source function
   TF1 *f1 = new TF1("f1","[0] +[1]*x +gaus(2)",0,5);
   f1->SetParameters(6,-1,5,3,0.2);
   // create and fill histogram according to the source function
   TH1F *h = new TH1F("h","background + signal",100,0,5);
   h->FillRandom("f1",2000);
   TF1 *fl = new TF1("fl",fline,0,5,2);
   fl->SetParameters(2,-1);
   //fit only the linear background excluding the signal area
   reject = kTRUE;
   h->Fit(fl,"0");
   reject = kFALSE;
   //store 2 separate functions for visualization
   TF1 *fleft = new TF1("fleft",fline,0,2.5,2);
   fleft->SetParameters(fl->GetParameters());
   h->GetListOfFunctions()->Add(fleft);
   gROOT->GetListOfFunctions()->Remove(fleft);
   TF1 *fright = new TF1("fright",fline,3.5,5,2);
   fright->SetParameters(fl->GetParameters());
   h->GetListOfFunctions()->Add(fright);
   gROOT->GetListOfFunctions()->Remove(fright);
   h->Draw();
}

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;
}