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

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 fitBarea(TString infname="",bool doweight = 1)
{
  if (doweight==0) weight="1";
  if (infname=="") infname=inputdata.Data();
  TFile *inf = new TFile(infname.Data());
  TTree *nt = (TTree*) inf->Get("ntKp");

  TFile *infMC = new TFile(inputmc.Data());
  TTree *ntGen = (TTree*)infMC->Get("ntGen");
  TTree *ntGen2 = (TTree*)inf->Get("ntGen");
  TTree *ntMC = (TTree*)infMC->Get("ntKp");
    
  ntGen->AddFriend(ntMC);
  ntGen2->AddFriend(ntMC);
    
  const int nBins = 1;
  double ptBins[nBins+1] = {10,60};
//  const int nBins = 1;
//  double ptBins[nBins+1] = {10,60};
  TH1D *hPt = new TH1D("hPt","",nBins,ptBins);

  for (int i=0;i<nBins;i++)
    {
      TF1 *f = fit(nt,ntMC,ptBins[i],ptBins[i+1]);
      double yield = f->Integral(5,6)/0.02;
      double yieldErr = f->Integral(5,6)/0.02*f->GetParError(0)/f->GetParameter(0);
      hPt->SetBinContent(i+1,yield/(ptBins[i+1]-ptBins[i]));
      hPt->SetBinError(i+1,yieldErr/(ptBins[i+1]-ptBins[i]));
    }
  
}

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{

  Double_t Lsum=0.;
  Double_t Nevt=dataColl.size();
  Double_t fs = par[0];
  Double_t fb = par[1];

  double tmppar[12];
  for(int i=0; i<12; i++){
    tmppar[i] = par[i+2];
  }
  TF1 *funcS = new TF1("funcS", exp_conv, -1., 20., 12);
  funcS->SetParameters(tmppar);
  Double_t Signorm = funcS->Integral(-1., 20.);
  TF1 *funcB = new TF1("funcB", expinv_power, -1., 20., 12);
  funcB->SetParameters(tmppar);
  Double_t Bkgnorm = funcB->Integral(-1., 20.);

  SigPDFnorm = Signorm;
  BkgPDFnorm = Bkgnorm;

  for ( int i=0; i<dataColl.size(); i++ ) {
    //PDF for signal and background
    Double_t x = dataColl[i];
    Double_t Ls = exp_conv(&x,tmppar)/Signorm;
    Double_t Lb = expinv_power(&x,tmppar)/Bkgnorm;
    //Get Log Likelihood		  
    Lsum += TMath::Log( (fs*Ls + fb*Lb) / (fs+fb) );
  }
  //f=2*( -1*Lsum + (fs+fb) - Nevt*TMath::Log(fs+fb) );
  f=2*( -1*Lsum + (fs+fb) - Nevt*TMath::Log(fs+fb) );
}

void plot_diff_xsec()
{
  TF1* fXsec = new TF1("fXsec", "diff_xsec(x)", 1, 0);
  fXsec->SetLineWidth(3);
  fXsec->SetTitle("#nu_{#mu}-e differential cross section (E_{#nu} = 2 GeV)");
  fXsec->Draw();
  fXsec->GetXaxis()->SetTitle("cos#theta");
  fXsec->GetXaxis()->CenterTitle();
  fXsec->GetYaxis()->SetTitle("d#sigma/dcos#theta (cm^{-2})");
  fXsec->GetYaxis()->CenterTitle();
  fXsec->GetYaxis()->SetTitleOffset(1.4);
  gPad->SetLogy(1);
  
  double totXsec = fXsec->Integral(0,1);
  for(int i = 0; i < 100; i++)
  {
    cout << fXsec->Integral(0.01*i,1)/totXsec << endl;
  }
  cout << totXsec << " " << fXsec->Integral(.99,1) << endl;
  
  ROOT::Math::GSLIntegrator ig(1.E-6,1.E-6,1000); 
  ROOT::Math::WrappedTF1 wf(*fXsec);
  ig.SetFunction(wf);
  cout << ig.Integral(.9,1) << endl;
  
  //~ TCanvas* c2 = new TCanvas("c2");
  //~ TF1* fCXsec = new TF1("fCXsec", "cumul_xsec(x)", 0, 90);
  //~ fCXsec->Draw();
  //~ cout << cumul_xsec(45) << endl;
}

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 fitDdecreasewid(TString collsyst="PbPb", Float_t centMin=0, Float_t centMax=100, TString outputfile="outfMasshisto/mass")
{
  collisionsystem = collsyst;
  infname = outputfile;
  centmin = centMin;
  centmax = centMax;

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

  TF1* fit (Float_t ptmin, Float_t ptmax);
  ofstream fout(Form("outYield/decreasewid_%s_cent_%.0f_%.0f.dat",collsyst.Data(),centmin,centmax));
  for(int i=0;i<nBins;i++)
    {
      TF1* f = fit(ptBins[i],ptBins[i+1]);
      Float_t yield = f->Integral(minhisto,maxhisto)/binwidthmass;
      fout<<ptBins[i]<<" "<<ptBins[i+1]<<" "<<yield<<endl;
    }
  fout.close();
}

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

double intXS(double *ecm, double *par)
{
    double res=par[3];
    double hs =par[5];

    TF1 *fXS = new TF1("fXS", rescs, 3.05, 3.12, 8);
    fXS->SetParameter(0, par[0]); //mass of Jpsi
    fXS->SetParameter(1, par[1]); //total width of Jpsi
    fXS->SetParameter(2, par[2]); //Gamma_ee
    fXS->SetParameter(3, par[3]); //Gauss resolution
    fXS->SetParameter(4, ecm[0]); //initial energy
    fXS->SetParameter(5, par[4]); // BfKsX
    fXS->SetParameter(6, par[5]); // h(s) for sigma_mumu 
    fXS->SetParameter(7, par[6]); // phase
    double my_intXS=fXS->Integral(TMath::Min(ecm[0],par[0])-20*res, TMath::Max(ecm[0],par[0])+20*res);


    double beta = 2.0*alpha/pi*(2.0*log(ecm[0]/me)-1.0);
    double delta= 0.75*beta+beta*beta/24.0*(1.0/3*log(ecm[0]/me/me)+2.0*pi*pi-37.0/4);
    double epsilon = 1.0e-8;
    double emin = 1.4;
    double UPL = 1.0-4.0*emin*emin/(ecm[0]*ecm[0]);
    TF1 *fXS_funs = new TF1("fXS_funs", funs, 0.0, 1.0, 1);
    fXS_funs->SetParameter(0, ecm[0]); //initial energy
    double my_intXS_qed1 = fXS_funs->Integral(epsilon, UPL);
    double my_intXS_qed2 = TMath::Power(epsilon,beta)*qedxs(ecm[0]*ecm[0]);//(86.8/(ecm[0]*ecm[0])), born cross section of mu+mu-
    TF1 *fXS_funh = new TF1("fXS_funh", funh, 0.0, 1.0, 1);
    fXS_funh->SetParameter(0, ecm[0]); //initial energy
    double my_intXS_qed3 = fXS_funh->Integral(0.0, UPL);
    double my_intXS_qed = hs*( my_intXS_qed3); //

    return my_intXS+my_intXS_qed;

}

void fitMassHisto(TString inputfile="outfiles/MassHisto", TString outputfile="outfiles/YieldHisto",
                  Float_t centmin=0., Float_t centmax=100.)
{
  infname = inputfile;
  centMin = centmin;
  centMax = centmax;

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

  TF1* fit(TH1D* h, TH1D* hMCSignal, TH1D* hMCSwapped, Float_t ptmin, Float_t ptmax, Int_t j);
  TH1D** hYield = new TH1D*[4];
  for(int i=0;i<4;i++) hYield[i] = new TH1D(Form("hYield_%s",tfend[i].Data()),"",nPtBins,ptBins);
  for(int i=0;i<nPtBins;i++)
    {
      TFile* infile = new TFile(Form("%s_cent_%.0f_%.0f_pt_%.0f_%.0f.root",infname.Data(),centMin,centMax,ptBins[i],ptBins[i+1]));
      TH1D* hMCSignal = (TH1D*)infile->Get("hMCSignal");
      TH1D* hMCSwapped = (TH1D*)infile->Get("hMCSwapped");
      for(int j=0;j<4;j++)
        {
          TH1D* h = (TH1D*)infile->Get(Form("h_%s",tfend[j].Data()));
          TF1* f = fit(h,hMCSignal,hMCSwapped,ptBins[i],ptBins[i+1],j);
          Double_t yield = f->Integral(minhisto,maxhisto)/binwidthmass;
          Double_t yieldErr = f->Integral(minhisto,maxhisto)/binwidthmass*f->GetParError(0)/f->GetParameter(0);
          hYield[j]->SetBinContent(i+1,yield);
          hYield[j]->SetBinError(i+1,yieldErr);
          delete h;
        }
      delete hMCSwapped;
      delete hMCSignal;
      delete infile;
    }

  TFile* outf = new TFile(Form("%s_cent_%.0f_%.0f.root",outputfile.Data(),centmin,centmax),"recreate");
  outf->cd();
  for(int i=0;i<4;i++) hYield[i]->Write();
  outf->Close();
  delete []hYield;
  delete outf;

  cout<<endl;
}

Double_t Likelihood(Double_t *x, Double_t *p)
{
   TF1 *out = new TF1("Outer",Outer,0.,B0 + 5.*sB,1);   
   Double_t low = B0 > 5.*sB ? B0 - 5.*sB : 0.;
   Double_t retval = out->Integral(low,B0+5.*sB,x,epsilon);
   delete out;
   return retval;
}

Double_t NtrueOverNmeas (Double_t ptMin, Double_t ptMax) {

  //Double_t n = 6;
  //Double_t n = 5.78121e+00; // from ref level
  Double_t n = 5.95950e+00; // from reco level
  Double_t c1 = 5.23359e-01;
  Double_t c2 = 1.15226e-01;
  Double_t c3 = -7.96949e-03;

  TF1 *fTrue = new TF1("fMeas","([1]+[2]*log(x)+[3]*log(x)*log(x)+[2]+2*[3]*log(x))/exp([0]*log(x*([1]+[2]*log(x)+[3]*log(x)*log(x))))",0,500); 
  fTrue->SetParameters(n,c1,c2,c3);
  TF1 *fMeas = new TF1("fTrue","1/exp([0]*log(x))",0,500);
  fMeas->SetParameter(0,n);

  return fTrue->Integral(ptMin,ptMax)/fMeas->Integral(ptMin,ptMax);

}

void GetRandomTest(){
	double k0=1.39;
	double k1 = 0.425;
	double theta0 = 3.41;
	double theta1 = 1.30;
	int iNpart=2;
	double k_=k0+k1*(iNpart-2);
        double theta_=theta0+theta1*TMath::Log(iNpart-1);
	TF1 *f = new TF1("f","TMath::GammaDist(x,[0],0,[1])",0,200);
	f->SetParameters(k1,theta_);
	cout<<"Value at 0 = "<<f->Eval(0)<<endl;
	cout<<"Value at 1e-11 = "<<f->Eval(1e-11)<<endl;
	cout<<"Integral 1= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	f->SetRange(1e-12,200);
	cout<<"Integral 2= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	cout<<"fXmin = "<<f->GetXmin()<<"\tfXmax = "<<f->GetXmax()<<"\tfNpx = "<<f->GetNpx()<<endl;
	f->SetNpx(1e5);
	TCanvas *c1 = new TCanvas();
	c1->SetLogy();
	cout<<"f mean = "<<f->Mean(0,200)<<endl;
	cout<<"math mean = "<<f->GetParameter(0)*f->GetParameter(1)<<endl;
	TH1D* h = new TH1D("h","h",1000,0,200);
	for(int i=0;i<1e6;i++){
		double para = f->GetRandom();
		h->Fill(para);
	}
	h->Scale(1.0/h->Integral()*1000/200);
	h->GetYaxis()->SetRangeUser(1e-10,1);
	h->SetMarkerStyle(24);
	h->SetMarkerColor(4);
	h->SetMarkerSize(1.1);
	TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
        leg->SetFillColor(0);
        leg->SetFillStyle(0);
        leg->SetBorderSize(0);
        leg->SetTextFont(42);
        leg->SetTextSize(0.03);
        leg->AddEntry(f,"function","lp");
        leg->AddEntry(h,"filled histogram","lp");
	h->Draw("P");
	f->Draw("same");
        leg->Draw("same");
	cout<<"h mean = "<<h->GetMean(1)<<endl;
	c1->Print("TestGetRandom.png");
	}

double GetChiSquareProb(double ndf, double chi2)
{
  TF1* func = new TF1("func","1/(pow(2,[0]/2.0)*TMath::Gamma([0]/2.0))*pow(x,[0]/2.0-1.0)*TMath::Exp(-x/2.0)",0,1000);//pdf of chi-squared distribution
  func->SetParameter(0,ndf);
  double integral = func->Integral(0,chi2);
  double prob = 1-integral;
  std::cout << "chi2/df = " << chi2/ndf << " with p-value = " << prob << std::endl;
  return prob;
}

void plotV2vsflow(){
gStyle->SetOptFit(kFALSE);
gStyle->SetOptStat(kFALSE);
TCanvas *c1 = new TCanvas("c1","c1",600,600);
TLegend *leg = new TLegend(0.2,0.1,0.4,0.4);
leg->SetTextSize(0.05);
leg->SetBorderSize(0);
leg->SetFillColor(0);
const int ninflow = 6;
int xbin=0;
TString dir;
double V2_Prod[ninflow],V2err_Prod[ninflow];
const int mult = 150;

TF1 *PtDistr  = new TF1("PtDistr","0.03*(exp (-(x/0.594540))+0.00499506*exp (-(x/1.89391)))", 0.3,6.0);       //Real Data
TF1 *V2vsPt = new TF1("V2vsPt","((x/3.31699)^2.35142/(1+(x/3.49188)^3.54429))*(.00005+(1/x)^1.50600)",0.3,6.0);
TF1 *V2 = new TF1("V2","0.03*(exp (-(x/0.594540))+0.00499506*exp (-(x/1.89391)))*((x/3.31699)^2.35142/(1+(x/3.49188)^3.54429))*(.00005+(1/x)^1.50600)",0.3,6.0);
double InV2 = V2->Integral(0.3,6.0)/PtDistr->Integral(0.3,6.0);
double inV2[ninflow]={};
for(int i=0;i<ninflow;i++){
        inV2[i] = InV2 * i / 5 ; 
	dir=Form("pPbDataV20%dm%d",i,mult);
	TFile *mergedV_Prod = TFile::Open(Form("%s/mergedV_Prod.root",dir.Data()));
        if(i==)
	TFile *mergedV_Prod = TFile::Open(Form("%s/mergedV_Sum.root",dir.Data()));
	TVectorD *vecMult = (TVectorD*)mergedV_Prod->Get("totmultall");
	TVectorD *vecNevent = (TVectorD*)mergedV_Prod->Get("Nevent");
	TVectorD *vecV2_Prod=(TVectorD*)mergedV_Prod->Get(Form("Vmean",xbin));
	TVectorD *vecV2err_Prod=(TVectorD*)mergedV_Prod->Get(Form("deltaVmean",xbin));
	V2_Prod[i]=(*vecV2_Prod)[xbin];
	V2err_Prod[i]=(*vecV2err_Prod)[xbin];
        cout<<inV2[i]<<endl;
//	V2sp[i]=InV2/V2_Prod[i];
//	V2spError[i]=InV2/V2_Prod[i]/V2_Prod[i]*(*vecV2err_Prod)[xbin];
//	cout<<V2sp[i]<<"pm"<<V2spError[i]<<"\t";
	}
        V2_Prod[1]=0;
cout<<endl<<InV2<<endl;
TH1D* hFrame = new TH1D("","",2000,-1,1);
hFrame->GetXaxis()->SetRangeUser(0., 0.08);
hFrame->GetYaxis()->SetRangeUser(0., 0.08);
hFrame->Draw();
TGraphErrors *grProd=new TGraphErrors(ninflow,inV2,V2_Prod,0,V2err_Prod);
grProd->SetMarkerColor(1);
grProd->SetMarkerSize(1.6);
grProd->SetLineColor(1);
grProd->SetMarkerStyle(20);
hFrame->GetXaxis()->SetTitle("Input V2");
hFrame->GetYaxis()->SetTitle("Calc V2");
hFrame->Draw();
grProd->Draw("Psame");
TLine *l = new TLine(0,0,0.06,0.06);
l->SetLineStyle(2);
l->Draw("same");
c1->SaveAs("V2vsflow.gif");
}

void MyEff(double n, double N, double& eff, double& errp, double& errn)
{

  TF1* gaus = new TF1("gaus","TMath::Gaus(x,0,1,1)");
  double lowerB = gaus->Integral(-1,0);
  double upperB = 1- lowerB;

  eff = n/N;
  if(n > N)
    {
      errp = errn = -999;
      return;
    }

  double a = n+1;
  double b = N-n+1;
  
  double integralPeak = incompletebeta(a, b, eff);

  if(integralPeak >= lowerB && integralPeak <= upperB)
    {
      double xmin = inverseBeta(integralPeak - lowerB, a, b);
      double xmax = inverseBeta(integralPeak + lowerB, a, b);
      errp = xmax-eff;
      errn = eff-xmin;
      return;

    }
  else if(integralPeak < lowerB)
    {
      double xmax = inverseBeta(2*lowerB, a, b);      
      errp = xmax-eff;
      errn = eff;
      return;

    }

  else if(integralPeak > upperB)
    {
      double xmin = inverseBeta(1-2*lowerB, a, b);
      errp = 1-eff;
      errn = eff-xmin;
      return;      
    }

  else{

    cout << "There is a bug !!!" << endl;
    errp = errn = -999;
    return;

  }


}