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

void multifit() {
   const Int_t np = 49;
   Float_t x[np] = {1.913521, 1.953769, 2.347435, 2.883654, 3.493567,
                    4.047560, 4.337210, 4.364347, 4.563004, 5.054247,
                    5.194183, 5.380521, 5.303213, 5.384578, 5.563983,
                    5.728500, 5.685752, 5.080029, 4.251809, 3.372246,
                    2.207432, 1.227541, 0.8597788,0.8220503,0.8046592,
                    0.7684097,0.7469761,0.8019787,0.8362375,0.8744895,
                    0.9143721,0.9462768,0.9285364,0.8954604,0.8410891,
                    0.7853871,0.7100883,0.6938808,0.7363682,0.7032954,
                    0.6029015,0.5600163,0.7477068,1.188785, 1.938228,
                    2.602717, 3.472962, 4.465014, 5.177035};

   TH1F *h = new TH1F("h","Example of several fits in subranges",np,85,134);
   h->SetMaximum(7);

   for (int i=0;i<np;i++) {
      h->SetBinContent(i+1,x[i]);
   }

   Double_t par[9];
   TF1 *g1    = new TF1("g1","gaus",85,95);
   TF1 *g2    = new TF1("g2","gaus",98,108);
   TF1 *g3    = new TF1("g3","gaus",110,121);
   TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",85,125);
   total->SetLineColor(2);
   h->Fit(g1,"R");
   h->Fit(g2,"R+");
   h->Fit(g3,"R+");
   g1->GetParameters(&par[0]);
   g2->GetParameters(&par[3]);
   g3->GetParameters(&par[6]);
   total->SetParameters(par);
   h->Fit(total,"R+");
}

void EtaMinv(Char_t* Spect)
{
  TH1F* Hist = gROOT->FindObject(Spect);
  Double_t par[5];

  TF1* f2 = new TF1("f2", "expo", 200, 700);
  TF1* f3 = new TF1("f3", "gaus", 500, 600);

  f2->SetLineWidth(1);
  f2->SetLineColor(8);
  Hist->Fit(f2,"R0");
  f3->SetLineWidth(1);
  f3->SetLineColor(6);
  Hist->Fit(f3,"R0+");

  TF1 *to= new TF1("to", "expo(0)+gaus(2)", 200, 700);

  f2->GetParameters(&par[0]);
  f3->GetParameters(&par[2]);

  to->SetParameters(par);
  to->SetLineWidth(2);
  to->SetLineColor(2);
  Hist->Fit(to,"R+");
}

double Alignment::Fit(TH1F &beta){
  double xmin=beta.GetXaxis()->GetXmin();
  double xmax=beta.GetXaxis()->GetXmax();

  double prev_norma=-1;
  double norma,mean,sigma;


  for(int i=0;i<10;i++){
    if(xmin>beta.GetXaxis()->GetXmax() || xmax<beta.GetXaxis()->GetXmin() || xmin>=xmax) return 0;  // Fail

    //    beta.Fit("gaus","0q","",xmin,xmax);
    beta.Fit("gaus","0q","",xmin,xmax);
    TF1 *func=beta.GetFunction("gaus");
    norma=func->GetParameter(0);
    mean=func->GetParameter(1);
    sigma=fabs(func->GetParameter(2));

    xmin=mean-1.0*sigma;
    xmax=mean+2.0*sigma;
    if(fabs(norma-prev_norma)<1e-3) break;

    prev_norma=norma;

  }

  return -norma;  // This have to be minimized
}

TF1 * iterateFitter(TH1F * histo, int i, TF1 * previousResiduals = 0, TCanvas * inputCanvas = 0)
{
  TCanvas * canvas = inputCanvas;
  std::stringstream ss;
  int iCanvas = i;
  ss << i;

  int nBins = histo->GetNbinsX();

  TF1 * functionToFit;
  if( previousResiduals == 0 ) {
    // functionToFit = relativisticBWFit(ss.str());
    // functionToFit = relativisticBWintFit(ss.str());
    // functionToFit = relativisticBWintPhotFit(ss.str());
    functionToFit = expRelativisticBWintPhotFit(ss.str());
    // functionToFit = crystalBallFit();
    // functionToFit = reversedCrystalBallFit();
    // functionToFit = lorentzianFit();
    // functionToFit = relativisticBWFit();
  }
  else {
    functionToFit = combinedFit(previousResiduals->GetParameters(), ss.str());
  }
  histo->Fit(functionToFit, "MN", "", 60, 120);


  double xMin = histo->GetXaxis()->GetXmin();
  double xMax = histo->GetXaxis()->GetXmax();
  double step = (xMax-xMin)/(double(nBins));
  TH1F * functionHisto = new TH1F(("functionHisto"+ss.str()).c_str(), "functionHisto", nBins, xMin, xMax);
  for( int i=0; i<nBins; ++i ) {
    functionHisto->SetBinContent( i+1, functionToFit->Eval(xMin + (i+0.5)*step) );
  }

  if( canvas == 0 ) {
    canvas = new TCanvas(("canvasResiduals"+ss.str()).c_str(), ("canvasResiduals"+ss.str()).c_str(), 1000, 800);
    canvas->Divide(2,1);
    canvas->Draw();
    iCanvas = 0;
  }
  canvas->cd(1+2*iCanvas);
  histo->Draw();
  functionToFit->SetLineColor(kRed);
  functionToFit->Draw("same");
  // functionHisto->Draw("same");
  // functionHisto->SetLineColor(kGreen);
  canvas->cd(2+2*iCanvas);
  TH1F * residuals = (TH1F*)histo->Clone();
  residuals->Add(functionHisto, -1);
  residuals->SetName("Residuals");

  // TF1 * residualFit = new TF1("residualFit", "-[0] + [1]*x+sqrt( ([1]-1)*([1]-1)*x*x + [0]*[0] )", 0., 1000. );
  // TF1 * residualFit = new TF1("residualFit", "[0]*(x - [1])/([2]*x*x + [3]*x + [4])", 0., 1000. );
  TF1 * residualFitFunction = residualFit(ss.str());
  residuals->Fit(residualFitFunction, "ME", "", 90.56, 120);

  residuals->Draw();

  return residualFitFunction;
}

void fitcont()
{
   //be sure default is Minuit since we will use gMinuit
   TVirtualFitter::SetDefaultFitter("Minuit");

   TCanvas *c1 = new TCanvas("c1");
   TH1F *h = new TH1F("h","My histogram",100,-3,3);
   h->FillRandom("gaus",6000);
   h->Fit("gaus");
   c1->Update();

   TCanvas *c2 = new TCanvas("c2","contours",10,10,600,800);
   c2->Divide(1,2);
   c2->cd(1);
   /*get first contour for parameter 1 versus parameter 2*/
   TGraph *gr12 = (TGraph*)gMinuit->Contour(40,1,2);
   gr12->Draw("alp");
   c2->cd(2);
   /*Get contour for parameter 0 versus parameter 2  for ERRDEF=2*/
   gMinuit->SetErrorDef(4); //note 4 and not 2!
   TGraph *gr2 = (TGraph*)gMinuit->Contour(80,0,2);
   gr2->SetFillColor(42);
   gr2->Draw("alf");
   /*Get contour for parameter 0 versus parameter 2 for ERRDEF=1*/
   gMinuit->SetErrorDef(1);
   TGraph *gr1 = (TGraph*)gMinuit->Contour(80,0,2);
   gr1->SetFillColor(38);
   gr1->Draw("lf");
}

void plotMuFromToys(std::string inputFile, std::string selectString="fit_status==0"){

	// Some Global preferences
        gROOT->SetStyle("Plain");
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gStyle->SetOptFit(1111);
	gStyle->SetOptStat(1111);
	gStyle->SetPalette(1,0);

	TFile *fi_ = TFile::Open(inputFile.c_str());
        TTree *tree_sb = (TTree*) fi_->Get("tree_fit_sb");

        TH1F *mures = new TH1F("mures","",25,-2.,2.);
        
        mures->SetLineColor(kBlue+3);
        mures->SetMarkerStyle(kOpenCircle);
        mures->SetLineWidth(2);
        
        mures->GetXaxis()->SetTitle("#mu - 1");
        mures->GetYaxis()->SetTitle(Form("no toys (%d total)",int(tree_sb->GetEntries())));
        mures->GetYaxis()->SetTitleOffset(1.05);
        mures->GetXaxis()->SetTitleOffset(0.9);
        mures->GetYaxis()->SetTitleSize(0.05);
        mures->GetXaxis()->SetTitleSize(0.05);

       	TCanvas *c = new TCanvas("c","",960,800);
        tree_sb->Draw("mu>>mures",selectString.c_str());
        mures->Fit("gaus");
        mures->GetFunction("gaus")->SetLineColor(kCyan+3);
        mures->Draw("pe1");
        c->SaveAs("mlfit/mu_residual.pdf");
}

void myfit()
{
   TString dir = gSystem->UnixPathName(__FILE__);
   dir.ReplaceAll("myfit.C","../hsimple.C");
   dir.ReplaceAll("/./","/");
   if (!gInterpreter->IsLoaded(dir.Data())) gInterpreter->LoadMacro(dir.Data());
   TFile *hsimple = (TFile*)gROOT->ProcessLineFast("hsimple(1)");
   if (!hsimple) return;

   TCanvas *c1 = new TCanvas("c1","the fit canvas",500,400);

   TH1F *hpx = (TH1F*)hsimple->Get("hpx");

// Creates a Root function based on function fitf above
   TF1 *func = new TF1("fitf",fitf,-2,2,3);

// Sets initial values and parameter names
   func->SetParameters(100,0,1);
   func->SetParNames("Constant","Mean_value","Sigma");

// Fit histogram in range defined by function
   hpx->Fit(func,"r");

// Gets integral of function between fit limits
   printf("Integral of function = %g\n",func->Integral(-2,2));
}

vector<double> one_fit(bool do_fit)
{
    TH1F* h = new TH1F("h","h", 20, -10, 10);
    h->FillRandom("gaus");
    vector<double> ret;
    if (do_fit) {
        h->Fit("gaus","");
        TF1* g = h->GetFunction("gaus");
        h->Draw();
        canvas->Modified();
        canvas->Update();
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(g->GetParameter(ind));
        }
    }
    else {
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(0.0);
        }
    }
    
    ret.push_back(h->GetMean());
    ret.push_back(h->GetRMS());
    delete h;
    return ret;
}

void ThrowPSEME(int me_i, bool use_old=false){

  // ====================================================
  //                       Create a fit
  // ====================================================
  TF1 *ftot = new TF1("fitFunction",fitFunction,-100,0,2);

  //threw one data PSE  and fill it
  if (use_old){
    cout<<"OLD"<<endl;
    siggnal = MEWH_old[me_i];
    background = MEWBB_old[me_i];
  }else{
    cout<<"NEW"<<endl;
    siggnal = MEWH[me_i];
    background = MEWBB[me_i];
  }

  TH1F * Err = new TH1F("Err","Err",1000,0,1);

  for (int pse=0;pse<100;pse++){

    //Get one PSE
    TH1F * data  =  MakeOnePSE(siggnal,1000,background,2000);
    
    //fit it
    data->Fit("fitFunction","b");   
    Err->Fill(ftot->GetParError(1));
    delete data;

  }
  Err->Draw();

}//ThrowPSEME

void peaks(Int_t np=10) {
   npeaks = TMath::Abs(np);
   TH1F *h = new TH1F("h","test",500,0,1000);
   //generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
   Int_t p;
   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 1;
      par[3*p+3] = 10+gRandom->Rndm()*980;
      par[3*p+4] = 3+2*gRandom->Rndm();
   }
   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
   f->SetNpx(1000);
   f->SetParameters(par);
   TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
   c1->Divide(1,2);
   c1->cd(1);
   h->FillRandom("f",200000);
   h->Draw();
   TH1F *h2 = (TH1F*)h->Clone("h2");
   //Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,2,"",0.10);
   printf("Found %d candidate peaks to fit\n",nfound);
   //Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,20,"same");
   if (hb) c1->Update();
   if (np <0) return;

   //estimate linear background using a fitting method
   c1->cd(2);
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h->Fit("fline","qn");
   //Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Double_t *xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
      par[3*npeaks+4] = 3;
      npeaks++;
   }
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   //we may have more than the default 25 parameters
   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h2->Fit("fit");
}

void Fit::process(const string &output_file)
{
    file = new TFile(output_file.c_str(), "RECREATE");
    TF1 *f1 = new TF1("f1", "gaus", fit_start, fit_end);
    hist->Fit("f1", "R");
    file->cd();
    hist->Write();
    file->Close();
}

void fit(const string &input_file)
{
    TChain *chain = new TChain("SmearEnergy");
    chain->Add(input_file.c_str());
    TH1F *h = new TH1F("energy", "Energy", 500, 0, 70);
    chain->Draw("smearEnergy>>energy", "smearEnergy>0");
    TF1 *g1 = new TF1("g1", "gaus", 4, 12);
    TF1 *g2 = new TF1("g2", "gaus", 14, 32);
    TF1 *g3 = new TF1("g3", "gaus", 36, 54);
    g1->SetLineColor(2);
    g2->SetLineColor(3);
    g3->SetLineColor(4);
    h->Fit(g1, "R");
    cout << "Integral(4,12): " << g1->Integral(4, 12) << endl;
    h->Fit(g2, "R+");
    cout << "Integral(14,32): " << g2->Integral(14, 32) << endl;
    h->Fit(g3, "R+");
    cout << "Integral(36,54): " << g3->Integral(36, 54) << endl;
}

void proton_pion(){
	r = new TRandom();
	TCanvas *c1 = new TCanvas();
	TH1F *Epion = new TH1F("Epion", "Epion", 100, 0, 140);
	TH1F *Npion = new TH1F("Npion", "Npion", 100, 0, 50);

	
	for (int i = 0; i < 100000; ++i)
	{
		double E_init_proton = 1000; // GeV
		double E_proton;
		double E_init_pion;

		twosplit(E_init_proton, E_proton, E_init_pion);
		
		double pion_branch[n_branches];
		
		Nsplit(E_init_pion, pion_branch, n_branches);
		int n_pion = 0;
		for (int i = 0; i < n_branches; ++i)
		{
			double energy_in = pion_branch[i];
			while(energy_in > 0){
				double E_a, E_b;
				twosplit(energy_in, E_a, E_b);
				double pion_E = E_a + M_pion;
				Epion->Fill(pion_E);

				energy_in = E_b - M_pion;
				if (energy_in>0)
				{
					n_pion++;
				}
				// double frac = r->Rndm();
				// double A = E_init_pion*frac;
				// pizza_mass = pizza_mass*(1-frac)-tax;
				// n_pis++;
				// h1->Fill(A);		
					
			}

		}
		Npion->Fill(n_pion);
			
	}
	c1->Divide(2,1);
	
	c1->cd(1);
	Epion->Draw();

	c1->cd(2);
	Npion->Draw();
	Npion->Fit("gaus");
}

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

void ZVXT1HF0()
{

  NCanvas(1,1,"data");
  NCanvas(1,1,"ratio");

  TH1F * HDatHF0;
  TH1F * HDat;

  TFile *fzee = new TFile("ZDiffOutputfile_DataOnly_2010_2011.root");
  // TFile *fzee = new TFile("ZDiffOutputfile.root");
  HDat  =  (TH1F*)fzee->Get("NVTX1_InstLumiPerBx_DATA10");
  
  HDatHF0  =  (TH1F*)fzee->Get("HF0NVTX1_InstLumiPerBx_DATA10");
  data->cd(1);
  HDat->Draw();
  NSetTitle(HDat,"Luminosity [10^{-30} cm^{-2} s^{-1}]", "Entries");
  HDat->SetTitle("Number of Z produced in events without PU");
  HDatHF0->Draw("SAME HIST");
  NHSetMarker(HDat,2,20,0.4);
  TLegend *legend = new TLegend(0.6,0.7,0.9,0.8);
  legend->SetTextFont(72);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);


  legend->AddEntry(HDat,"All Z","p");
  legend->AddEntry(HDatHF0,"Z with minE_HF=0","l");
  legend->Draw();


  TH1F * Hra = (TH1F * ) HDatHF0->Clone();
  //  Hra->GetXaxis->SetRange(0,12);
  NStat(Hra,0);
  Hra->Divide(HDat);
  NLogAxis(0,1);
  ratio->cd(1);

  Hra->Draw();
  Hra->SetTitle("Fraction of Z with (PU=0 && E_HF=0)");
  Hra->Fit("expo","","",0.,1.2);
  NSetTitle(Hra,"Luminosity [10^{-30} cm^{-2} s^{-1}]", "Fraction");
  NLogAxis(0,1);
  NText(.05,0.1, "0.62");
  NText(0.3,0.1, "2.47");
  NText(0.7,0.1, "4.93");
  NText(1.1,0.1, "7.4");
  NText(1.4,0.1, "<Number of Int.>");


}