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

void makePlots_InvariantMass_ElEl(){


  gROOT->SetStyle("Plain");
  TCanvas *c1 = new TCanvas("c1_fit1","The Fit Canvas",200,10,700,500);
  c1->SetGridx();
  c1->SetGridy();
  c1->GetFrame()->SetFillColor(21);
  c1->GetFrame()->SetBorderMode(-1);
  c1->GetFrame()->SetBorderSize(5); 
 
  TFile* file = TFile::Open("Output_LowPtSUSY_Tree_PPZ.root");
  TH1F *h_InvMass = (TH1F*) file->Get("h_InvariantMass_ElEl");
  h_InvMass->SetLineColor(kBlue);
  h_InvMass->SetLineWidth(2);
  h_InvMass->Rebin(60);
  h_InvMass->GetYaxis()->SetTitle("Events/2.0 GeV");
  h_InvMass->GetYaxis()->SetTitleOffset(1.4);
  h_InvMass->GetXaxis()->SetRangeUser(70, 110); 
  TF1 *fSignal = new TF1("fSignal","gaus",70.0,110.0);
  fSignal->SetParLimits(0, 11500.0, 15000.0);
  fSignal->SetParLimits(1, 80.0, 95.0);
  fSignal->SetParLimits(2, 0.0, 1.0);
  h_InvMass->Draw("HIST");
  h_InvMass->Fit("fSignal", "", "", 70.0, 110.0);
  c1->SaveAs("h_InvariantMass_ElEl.pdf");
  c1->SaveAs("h_InvariantMass_ElEl.png");

} 

TH2D* GetJetCorrFunc2D_doublegaussian(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");  
                      float min = hcorrphi->GetMinimum();
                      hcorrphi->SetAxisRange(-1,1,"X");
                      float nearmax = hcorrphi->GetMaximum();
                      hcorrphi->SetAxisRange(PI-1,PI+1,"X");
                      float awaymax = hcorrphi->GetMaximum();

		      TF1* fitfunc = new TF1("fitfunc",doubleGaussian,-PI/2.,3.*PI/2.,5);		    
                      fitfunc->SetParameters(min,nearmax-min,0.3,awaymax-min,0.5);  
                      fitfunc->SetParLimits(0,0,100000);
                      fitfunc->SetParLimits(1,0,100000);
                      fitfunc->SetParLimits(2,0,100000);
                      fitfunc->SetParLimits(3,0,100000);
                      fitfunc->SetParLimits(4,0,100000);		    
		      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);                      
                      delete fitfunc;
		    }

                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");
                    return hcorr_clone;
}

Bool_t misalignmentDependence(TCanvas *c1old,
                              Int_t nFiles,TString *names,TString misalignment,Double_t *values,Double_t *phases,TString xvar,TString yvar,
                              Bool_t drawfits = true,
                              Bool_t resolution = false,
                              TString saveas = "")
{
    TF1 *f = 0;
    TString functionname = "";

    //if only one parameter is of interest
    TString parametername = "";
    Int_t parameter = 9999;

    //if multiple parameters are of interest
    Int_t nParameters = -1;
    TString *parameternames = 0;
    Int_t *parameters = 0;

    if (misalignment == "breitwigner")
    {
        if (xvar == "" && yvar == "ZZMass")
        {
            f = new TF1("breitwigner","[2]*TMath::BreitWigner(x,[0],[1])");
            f->SetParameters(91,2);
            f->SetParLimits(0,70,105);
            f->SetParLimits(1,0,10);
            nParameters = 2;
            Int_t tempParameters[2] = {0,1};
            TString tempParameterNames[2] = {"M;GeV","#gamma;GeV"};
            parameters = tempParameters;
            parameternames = tempParameterNames;
            functionname = "Breit Wigner";
        }
    }

    if (xvar == "" && functionname == "")
    {
        if (c1old == 0 || misalignment == "" || values == 0) return false;
        misalignmentDependence(c1old,nFiles,names,misalignment,values,phases,xvar,yvar,(TF1*)0,0,"","",resolution,saveas);
        return true;
    }

    if (functionname == "") return false;
    if (drawfits)
    {
        parameter = -parameter-1;
        for (int i = 0; i < nParameters; i++)
            parameters[i] = -parameters[i]-1;
    }
    if (nParameters > 0)
        misalignmentDependence(c1old,nFiles,names,misalignment,values,phases,xvar,yvar,
                               f,nParameters,parameters,parameternames,functionname,resolution,saveas);
    else
        misalignmentDependence(c1old,nFiles,names,misalignment,values,phases,xvar,yvar,
                               f,parameter,parametername,functionname,resolution,saveas);
    delete f;
    return true;
    
}

/************************************************************
*  Fit ADC spectrum
************************************************************/
int FitAdcHist(TH1F *h, float data[NDATA], const int print ) {
  TF1 *fadc;
  TPaveText *adcPT;

  // Protect for null histograms
  if( h == 0 || h->GetEntries()<1000. ) return -1;
  printf("Fitting %s\n",h->GetName());

  gStyle->SetOptStat("emr"); 
  adcPT = new TPaveText(0.58,0.315,0.98,0.635,"NDC");

//  Refit ADC
  Double_t n    = h->GetEntries();
  Double_t rms  = h->GetRMS();
  Double_t mean = h->GetMean();
  Double_t ymin = mean-2*rms;
  if( ymin < 60. ) ymin=60.;
  Double_t ymax = mean+1.5*rms;
  fadc = new TF1("adc_fun",skewnormal,ymin,ymax,4);
  fadc->SetParameters(n,mean,rms,1.5);
  fadc->SetParLimits(1,mean-rms,mean+0.5*rms);
  fadc->SetParLimits(2,rms*0.5,rms*1.5);
  fadc->SetParLimits(3,1.,4.);
  h->Fit("adc_fun","r"); 
//  If skew is too high refit with fixed skew
  if( fadc->GetParameter(3) > 2.5 ) {
    printf("REFIT %s\n",h->GetName());
    fadc->SetParameter(3,1.6);
    fadc->SetParLimits(3,1.6,1.6);
    h->Fit("adc_fun","r") ;
  }
// if have enough hits use fit results  
  if( n > MINHITS ) {
    data[13] = fadc->GetParameter(1);
    data[14] = fadc->GetParError(1);
    data[15] = fadc->GetParameter(2);
    data[16] = fadc->GetParError(2);
    data[17] = fadc->GetParameter(3);
    data[18] = fadc->GetParError(3);
    double dof = fadc->GetNDF();
    if( dof > 0. ) data[19] = fadc->GetChisquare()/dof;
    adcPT->AddText(Form("Peak %.1lf #pm %.1lf",data[13],data[14]));
    adcPT->AddText(Form("Width %.1lf #pm %.1lf",data[15],data[16]));
    adcPT->AddText(Form("Skew %.2lf #pm %.2lf",	data[17],data[18]));
    if( dof > 0. ) adcPT->AddText(Form("#Chi^{2}/DoF     %.2lf",data[19]));
    adcPT->SetTextColor(2);
    adcPT->Draw();
// if have too few hits use fit stats instead of fit results  
  } else {
    data[12] = mean;
    data[13] = h->GetMeanError();
    data[14] = rms;
    data[15] = h->GetRMSError();
    data[16] = h->GetSkewness();
  }
  if( print ) c1->Print(Form("%s_%s.png",fname,h->GetName()));
  return 0;
}    //end FitAdcHist

void example(double E0 = 50, int nevents = 100000)
{
   TStopwatch timer;

   // 12C* -> 3(4He)
   // compound nucleus = carbon
   KVNucleus CN(6, 12);
   CN.SetExcitEnergy(E0);

   // decay products
   KVEvent decay;
   KVNucleus* n = decay.AddParticle();
   n->SetZandA(2, 4);
   n = decay.AddParticle();
   n->SetZandA(2, 4);
   n = decay.AddParticle();
   n->SetZandA(2, 4);

   MicroStat::mdweight gps;
   Double_t etot = E0 + decay.GetChannelQValue();

   if (etot <= 0) {
      printf("Break-up channel is not allowed\n");
      return;
   }
   gps.SetWeight(&decay, etot);
   gps.initGenerateEvent(&decay);

   TH1F* h1 = new TH1F("h1", "Kinetic energy of alpha particle 3", 200, 0, etot * 2. / 3.);
   h1->Sumw2();

   KVEvent event;

   while (nevents--) {
      gps.GenerateEvent(&decay, &event);
      h1->Fill(event.GetParticle(3)->GetKE());
      gps.resetGenerateEvent();
   }

   h1->Draw();
   TF1* EDis = new TF1("EDis", edist, 0., etot, 3);
   EDis->SetNpx(500);
   EDis->SetParLimits(0, 0, 1.e+08);
   EDis->SetParLimits(1, 0, 2 * etot);
   EDis->FixParameter(2, 3);
   gStyle->SetOptFit(1);
   h1->Fit(EDis, "EM");

   timer.Print();
}

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

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

TF1* linLorentzianFit(TH1* histoY){
  TString name = histoY->GetName() + TString("Fit");

  // Fit slices projected along Y from bins in X 
  TF1 *fit = new TF1(name,lorentzianPlusLinear,70,110,5);
  fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean(),10,-0.1);
  fit->SetParNames("norm","width","mean","offset","slope");
  fit->SetParLimits(1,-10,10);
  //fit->SetParLimits(2,85,95);
  //fit->SetParLimits(2,90,93);
  fit->SetParLimits(2,2,4);
  fit->SetParLimits(3,0,100);
  fit->SetParLimits(4,-1,0);
  fit->SetLineWidth(2);
  //histoY -> Fit(name,"0","",85,97);
  histoY -> Fit(name,"0","",2,4);
  return fit;
}

TF1* GetFit(const char *address, int lyr, double xinit) {
  const char *bgr = Form("[6]*TMath::Exp([7]*(x-%f))",xinit);
  const char *la1 = "(1-[3]-[4]-[5])*TMath::Landau(x,[1],[2],1)";
  const char *la2 = "[3]*TMath::Landau(x,2*[1]+1.4*[2],2.0*[2],1)";
  const char *la3 = "[4]*TMath::Landau(x,3*[1]+3.3*[2],3.0*[2],1)";
  const char *la4 = "[5]*TMath::Landau(x,4*[1]+5.6*[2],4.0*[2],1)";
  TF1 *ret = new TF1("fit_H1",
                     Form("[0]*(%s+%s+%s+%s)+%s",la1,la2,la3,la4,bgr) );
  ret->SetParName(0,"A");
  ret->SetParName(1,"lambda");
  ret->SetParName(2,"sigma");
  ret->SetParName(3,"f2");
  ret->SetParameter(0,1e4);  ret->SetParLimits(0,1e2,1e7);
  ret->SetParameter(1,21);   ret->SetParLimits(1,13,25.0);
  ret->SetParameter(2,3.0);  ret->SetParLimits(2,1.5,5.0);
  if(lyr>1) {
    ret->SetParameter(1,15);   ret->SetParLimits(1,2,30);
    ret->SetParameter(2,5.0);  ret->SetParLimits(2,0.8,10.0);
  }
  ret->SetParameter(3,0.20); ret->SetParLimits(3,0,0.40);
  ret->SetParameter(4,0.0);
  ret->SetParameter(5,0.0);
  ret->SetParameter(6,500); ret->SetParLimits(6,1e2,1e7);
  ret->SetParameter(7,-1);  ret->SetParLimits(7,-10,-0.1);

  ret->SetLineColor(kRed-3);

  ifstream infit;
  infit.open( address );
  double tmp;
  bool found=false;
  for(int n=0; n!=ret->GetNpar(); ++n) {
    infit >> tmp;
    if(!infit.good()) break;
    ret->SetParameter(n,tmp);
    infit >> tmp;
    found = true;
  }
  infit.close();
  if(found) {
    cout << " Previous fit results found" << endl;
    ret->SetParLimits(4,0,0.40);
    ret->SetParLimits(5,0,0.40);
  } else {
    ret->SetParLimits(4,+1,-1);
    ret->SetParLimits(5,+1,-1);
  }
  return ret;
}

TF1* GetFitFunc_ZYAM_pp(TH1D* h)
{
                    TH1D* hcorrphi = (TH1D*)h->Clone(h->GetName());
//                    double histminY = hcorrphi->GetBinContent(hcorrphi->GetMinimumBin());
//                    double histminX = hcorrphi->GetBinCenter(hcorrphi->GetMinimumBin());
                    double histminY = hcorrphi->GetBinContent(hcorrphi->FindBin(1.2));
                    double histminX = hcorrphi->GetBinCenter(hcorrphi->FindBin(1.2));

                    //hcorrphi->SetAxisRange(-0.01,1.2,"X");
                    TF1* fitfunc = new TF1(Form("fitfunc_%s",h->GetName()),"[0]+[1]*(x-[2])*(x-[2])+[3]*abs((x-[2])*(x-[2])*(x-[2]))",0.5,2.5); //std 0.6 1.55 vs pT ; 0.6 1.8 vs eta
                    fitfunc->SetParameters(histminY,0.002,histminX,0.0);
                    fitfunc->SetParLimits(1,0,1000);
                    fitfunc->SetParLimits(2,0.2,2.0);
                   // fitfunc->SetParLimits(3,0,100);
                    fitfunc->FixParameter(3,0);

                    for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%s",h->GetName()),"NO","",histminX-0.6,histminX+0.6);
                    fitfunc->SetName(Form("fitfunc_%s_%.3f",h->GetName(),fitfunc->GetParameter(0)));
                    return fitfunc;
}

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

  
  TF1* fitFunction;

  if( funcType=="NSC" ) {
    fitFunction = new TF1(funcName.c_str(), NSC, rangeMin, rangeMax, 3); 
    fitFunction->SetParameters( 0., 1., 0.05);
    fitFunction->SetParLimits( 0, -2., 11.);
    fitFunction->SetParLimits( 1, -2., 2.);
    fitFunction->SetParLimits( 2, -0.3, 0.3 );
  } else if( funcType=="NSCPF" ) {
    fitFunction = new TF1(funcName.c_str(), NSCPF, rangeMin, rangeMax, 4);
    fitFunction->SetParameters( 0., 1., 0.05, 0.);
    fitFunction->SetParLimits( 0, -2., 5.);
    fitFunction->SetParLimits( 1, 0., 1.1);
    fitFunction->SetParLimits( 2, 0., 0.12 );
    fitFunction->SetParLimits( 3, -1., 1.);
  } else {
    std::cout << "Function '" << funcType << "' not implemented yet for fitResolutionGraph. Exiting." << std::endl;
    exit(119);
  }

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

  return fitFunction;

}

//______________________________________________________________________________
void Fit(TH1* h, Double_t* outPos, Double_t* outFWHM)
{
    // Perform fit.
    
    Char_t tmp[256];

    Double_t x1 = 50;
    Double_t x2 = 230;

    TF1* func = new TF1("func", "gaus(0)+pol5(3)", x1, x2);
    func->SetParameters(h->GetMaximum(), 135, 8, 1, 1, 0.1, 0.1);
    func->SetLineColor(kBlue);
    func->SetParLimits(1, 130, 140);
    func->SetParLimits(2, 1, 15);
    h->Fit(func, "RBQO");
    
    // get position and FWHM
    fPi0Pos = func->GetParameter(1);
    *outPos = fPi0Pos;
    *outFWHM = 2.35*func->GetParameter(2);

    // indicator line
    TLine* line = new TLine();
    line->SetLineWidth(2);
    line->SetX1(fPi0Pos);
    line->SetX2(fPi0Pos);
    line->SetY1(0);
    line->SetY2(h->GetMaximum());

    TF1* fBG = new TF1("funcBG", "pol5", x1, x2);
    for (Int_t i = 0; i < 6; i++) fBG->SetParameter(i, func->GetParameter(3+i));
    fBG->SetLineColor(kRed);

    // draw 
    h->GetXaxis()->SetRangeUser(0, 300);
    h->Draw();
    func->Draw("same");
    fBG->Draw("same");
    line->Draw("same");
}

TF1 * relativisticBWintFit(const std::string & index)
{
  RelativisticBWwithZGammaInterference * fobj = new RelativisticBWwithZGammaInterference;
  TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), fobj, 60, 120, fobj->parNum(), "RelativisticBWwithZGammaInterference");
  functionToFit->SetParameter(0, 2.);
  functionToFit->SetParameter(1, 90.);
  functionToFit->SetParameter(2, 1.);
  functionToFit->SetParameter(3, 1.);

  functionToFit->SetParLimits(3, 0., 1.);

  return functionToFit;
}

//--------------------------------------------------------------------------------------------------
// Fit histogram with function
//--------------------------------------------------------------------------------------------------
TF1* fitFakeRatePtHist(TH1F *FakeRatePt) {

  TF1 *fitFunction = new TF1("fitFunction", "[0]*(1 - exp(-1*[1]*(x - [2])))",-10,10);
  fitFunction->SetParameter(0,0.05);
  fitFunction->SetParameter(1,0.07);
  fitFunction->SetParameter(2,10.0);
  fitFunction->SetParLimits(2,9.5,10.5);
  
  FakeRatePt->Fit("fitFunction");
  
  return fitFunction;

}