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

TF1 *f_pol4(char *funcn, double min, double max, bool skippeak)
{
  TF1 *retf = new TF1(funcn,&d_pol4,min,max,7);
  if (skippeak) {
    retf->FixParameter(5,0.74);
    retf->FixParameter(6,0.85);
  } else {
    retf->FixParameter(5,0);
    retf->FixParameter(6,0);
  }
  return retf;
}

void Fit(TH1D* hAll,TH1D* hSignal){
  Double_t par[3];
  hSignal->Fit("gaus","Q","");
  hSignal->GetFunction("gaus")->GetParameters(par);

  TF1* fitSB = new TF1("fitSB","gaus(0)+pol2(3)",1.10,1.13);
  fitSB->FixParameter(0,par[0]);
  fitSB->FixParameter(1,par[1]);
  fitSB->FixParameter(2,par[2]);
  hAll->Fit("fitSB","","",1.10,1.13);

  hSignal->GetFunction("gaus") ->SetLineWidth(2);
  hAll   ->GetFunction("fitSB")->SetLineWidth(2);
  hSignal->GetFunction("gaus") ->SetLineColor(kBlue);
  hAll   ->GetFunction("fitSB")->SetLineColor(kBlue);

}

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 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* GetFitFunc_ZYAM_poor(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])+[3]*(x-[2])*(x-[2])*(x-[2])",0.5,1.5); //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->FixParameter(3,0);

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

TH1D* GetJetCorrFunc1D_ZYA1(int itrg, int jass)
{
                    TH1D* hcorrphi = (TH1D*)GetRawCorrFunc1D_ratio(itrg,jass);
                    TH1D* hcorrphi_clone = (TH1D*)hcorrphi->Clone(Form("corrphi_clone_itrg%d_jass%d",itrg,jass));
                    hcorrphi->SetAxisRange(0.5,2.0,"X");
                    double histminY = hcorrphi->GetMinimum();
		    TF1* fitfunc = new TF1("fitfunc","[0]+[1]*x",0.8,1.2);		    
                    fitfunc->SetParameters(histminY,0);
                    fitfunc->FixParameter(1,0);
		    for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit("fitfunc","RNO");                     
		    float level = fitfunc->GetParameter(0);
		    for(int ibin=1;ibin<=hcorrphi_clone->GetNbinsX();ibin++) hcorrphi_clone->SetBinContent(ibin,hcorrphi_clone->GetBinContent(ibin)-level);
                    float max = hcorrphi_clone->GetBinContent(hcorrphi_clone->GetMaximumBin());
                    hcorrphi_clone->SetAxisRange(ymin,max*1.3,"Y");
		    hcorrphi_clone->SetAxisRange(-PI/2.,3.*PI/2.,"X");
                    delete fitfunc;
                    return hcorrphi_clone;
}

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

void vertexAssociationAnalysis(TH1* fraction, unsigned int index, double& slope,double& slopeError,double& constNorm,double& constError,TGraphErrors*& graph,TFitResultPtr fitRes) {
  // do the fit
  if(constrainSlope) {
    fraction->Fit("expo","LS","",xminfit[index],xmaxfit[index]);
    TF1* expo = (TF1*)gROOT->FindObjectAny("expo");
    expo->FixParameter(1,contrainedSlopeParameters[0]+contrainedSlopeParameters[1]*(index+1));
    fitRes = fraction->Fit("expo","LSB","",xminfit[index],xmaxfit[index]);
    expo->SetParLimits(1,-10000,10000);
  } else {
    // fitRes = fraction->Fit("expo","LS","",xminfit[index],xmaxfit[index]);  
    fitRes = fraction->Fit("expo","LS","",getConstrainedMinimum(fraction,0.1,0.3),getConstrainedMaximum(fraction,0.7,0.9));
  }
  // estimate efficiency and purity as a function of the cut
  double *efficiency, *purity, *efferror, *purerror;
  efficiency = new double[nSteps];
  purity     = new double[nSteps];
  efferror   = new double[nSteps];
  purerror   = new double[nSteps];
  unsigned int step = 0;
  for(float cut=cutmin; cut<cutmax; cut += (cutmax-cutmin)/float(nSteps), ++step) {
    estimateEfficiencyAndPurity(fraction,cut,efficiency[step],purity[step],efferror[step],purerror[step],fitRes);
  }
  // output: the graph with efficiency and purity vs cut and the expo slope
  slope      = fraction->GetFunction("expo")->GetParameter(1);
  slopeError = fraction->GetFunction("expo")->GetParError(1);
  
  constNorm  = fraction->GetFunction("expo")->GetParameter(0);
  constError = fraction->GetFunction("expo")->GetParError(0);

  //slope = result->Parameter(1);
  graph = new TGraphErrors(nSteps,purity,efficiency,purerror,efferror);
  //graph = new TGraph(nSteps,purity,efficiency);   
  delete[] efficiency;
  delete[] purity;
  delete[] efferror;
  delete[] purerror;
}

//.........这里部分代码省略.........
        c[i]->cd(1);
        
        /*p definitions
         [0] CB1 yield;
         [1] Common mean of CB and Gaus;
         [2] CB1 sigma;
         [3] CB n;
         [4] CB alpha;
         [5] CB2 yield;
         [6] CB2 sigma;
         [7-10] poly 3;
         [11] v2 signal;
         [12-13] v2 bkg;
         */
        TF1* f = new TF1(Form("f_%d",i), crystalball_function_total, fit_range_low, fit_range_high, 11);
        f->SetLineColor(2);
        f->SetLineWidth(1);
        f->SetParNames("CB1_Yield","common_mean","CB1_sigma","CB_N","CB_Alpha","CB2_Yield","CB2_Sigma","Pol0","Pol1","Pol2","Pol3");

        //first fit data mass signal + bkg
        
        f->SetParameter(0,10000.);
        f->SetParameter(1,JPsi_mass);
        f->SetParameter(2,0.03);
        f->SetParameter(3,1.0);
        f->SetParameter(4,1.0);
        f->SetParameter(5,10000);
        f->SetParameter(6,0.03);
        
        f->SetParLimits(2,0.01,0.1);
        f->SetParLimits(6,0.01,0.1);
        
        //fix alpha & n from MC
        f->FixParameter(4,alpha_fit[i]);
        f->FixParameter(3,n_fit[i]);
        
        f->FixParameter(1,JPsi_mass); //for first few attempt fix mean of gaussian to get reasonable estimation of other pars; later open it up
        h_data->Fit(Form("f_%d",i),"q","",fit_range_low,fit_range_high);
        h_data->Fit(Form("f_%d",i),"q","",fit_range_low,fit_range_high);
        f->ReleaseParameter(1); //now let gaussian mean float
        h_data->Fit(Form("f_%d",i),"L q","",fit_range_low,fit_range_high);
        h_data->Fit(Form("f_%d",i),"L q","",fit_range_low,fit_range_high);
        h_data->Fit(Form("f_%d",i),"L m","",fit_range_low,fit_range_high);
        
        
        //draw D0 signal separately
        TF1* f1 = new TF1(Form("f_sig_%d",i), crystalball_function_signal, fit_range_low, fit_range_high, 7);
        f1->SetLineColor(kOrange-3);
        f1->SetLineWidth(1);
        f1->SetLineStyle(2);
        f1->SetFillColorAlpha(kOrange-3,0.3);
        f1->SetFillStyle(1001);
        f1->FixParameter(0,f->GetParameter(0));
        f1->FixParameter(1,f->GetParameter(1));
        f1->FixParameter(2,f->GetParameter(2));
        f1->FixParameter(3,f->GetParameter(3));
        f1->FixParameter(4,f->GetParameter(4));
        f1->FixParameter(5,f->GetParameter(5));
        f1->FixParameter(6,f->GetParameter(6));
        
        fmasssig[i] = (TF1*)f1->Clone();
        fmasssig[i]->SetName(Form("masssigfcn_pt%d",i));
        fmasssig[i]->Write();
        
        f1->Draw("LSAME");
        

//________________________________________________________________________________
void peaks(TH1* h, Int_t np=10, Int_t ij=0) {
  if (! h) return;
  npeaks = TMath::Abs(np);
  if (! c1)  c1 = new TCanvas();
  else       c1->Clear();
  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
  c1->Divide(1,2);
  c1->cd(1);
  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,5,"",0.05);
  printf("Found %d candidate peaks to fit\n",nfound);
  if (! nfound) return;
  //Estimate background using TSpectrum::Background
  TH1 *hb = s->Background(h,20,"same");
  hb->Draw("same");
  c1->Update();
  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
  if (np <0) return;

  //estimate linear background using a fitting method
  c1->cd(2);
  TF1 *fline = new TF1("fline","pol1",0,1000);
  fline->FixParameter(1,0.);
  h->Fit("fline","qnlw");
  if (c2 && ij > 0) {c2->cd(ij+1); h->Draw(); c2->Update(); c1->cd(2);}
  //Loop on all found peaks. Eliminate peaks at the background level
  Double_t par[3000];
  par[0] = fline->GetParameter(0);
  par[1] = fline->GetParameter(1);
  npeaks = 0;
  Float_t *xpeaks = s->GetPositionX();
  Float_t ymax = 0;
  for (Int_t p=0;p<nfound;p++) {
    Float_t xp = xpeaks[p];
    Int_t bin = h->GetXaxis()->FindBin(xp);
    Float_t yp = h->GetBinContent(bin);
    if (yp-3*TMath::Sqrt(yp) < fline->Eval(xp)) continue;
    par[3*npeaks+2] = yp;
    if (yp > ymax) ymax = yp;
    par[3*npeaks+3] = xp;
    par[3*npeaks+4] = 3;
    npeaks++;
  }
  cout << "Found " << npeaks << " useful peaks to fit" << endl;
  Int_t NP = 0;
  Int_t nbad = 0;
  TString LineH("  {\"");  LineH +=  h->GetName(); LineH += "\"";
  TString Line("");
  struct ParErr_t {Double_t par, err;};
  ParErr_t parErr[10];
  TF1 *fit = 0;
  if (ymax > 2*par[0]) {
    cout << "Now fitting: Be patient" << endl;
    fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
    TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
    fit->SetParameter(0,par[0]);
    fit->FixParameter(1,0.);
    for (Int_t p = 0; p < npeaks; p++) {
      fit->SetParName(3*p+2,Form("A%i",p));
      fit->SetParLimits(3*p+2,0,1e6);
      fit->SetParameter(3*p+2,par[3*p+2]);
      fit->SetParName(3*p+3,Form("#mu%i",p));
      fit->SetParLimits(3*p+3,TMath::Max(0.,par[3*p+3]-2), TMath::Min(240.,par[3*p+3]+2));
      fit->SetParameter(3*p+3,par[3*p+3]);
      fit->SetParName(3*p+4,Form("#sigma%i",p));
      fit->SetParLimits(3*p+4,0.01,20);
      fit->SetParameter(3*p+4,par[3*p+4]);
    }
    fit->SetNpx(1000);
    h2->SetStats(0);
    h2->Fit("fit","l");
    if (c2 && ij > 0) {c2->cd(ij); h2->Draw("same"); c2->Update(); c1->cd(2);}
    fit->GetParameters(par);
    for (Int_t p = 0; p<np;p++) {
      if (p < npeaks && par[3*p+2] > 5*fit->GetParError(3*p+2) && 
	  par[3*p+2] > par[0]) {
	if (TMath::Abs(par[3*p+4]) > 5.0) nbad++;
	//      Line += Form(",%f,%f,%7.2f,%5.2f",par[3*p+2],fit->GetParError(3*p+2),par[3*p+3],TMath::Abs(par[3*p+4]));
	parErr[NP].par = par[3*p+3];
	parErr[NP].err = TMath::Abs(par[3*p+4]);
	for (Int_t i = 0; i < NP; i++) {
	  if (parErr[i].par > parErr[NP].par) {
	    ParErr_t temp = parErr[i];
	    parErr[i] = parErr[NP];
	    parErr[NP] = temp;
	  }
	}
	NP++;
      } 
    }
  }
  for (Int_t p = 0; p < np; p++) {
    if (p < NP) Line += Form(",%7.2f,%5.2f",parErr[p].par,parErr[p].err);
    else  Line += ",0,0";
  }
  Line += "},";
//.........这里部分代码省略.........

TF1* fitDstar5prongs(TTree* nt, Double_t ptmin, Double_t ptmax)
{
  static int count5p=0;
  count5p++;
  
  TCanvas* c = new TCanvas(Form("c_5p_%d",count5p),"",600,600);
  TH1D* h = new TH1D(Form("h_5p_%d",count5p),"",60,0.140,0.160);
  TF1* f = new TF1(Form("f_5p_%d",count5p),"[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*((1-[8])*TMath::Gaus(x,[6],[7])/(sqrt(2*3.14159)*[7])+[8]*TMath::Gaus(x,[6],[9])/(sqrt(2*3.14159)*[9]))",minmass3prong,maxmass3prong);
  nt->Project(Form("h_5p_%d",count5p),"Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f&&Dpt<%f)",weight.Data(),seldata5p.Data(),triggerselection[isData].Data(),ptmin,ptmax));   
    
  f->SetLineColor(4);
  f->SetParameters(0,0,0,0,0,2e2,1.45491e-1,9e-4,0.1,8e-4);
  f->FixParameter(9,15e-4);
  f->FixParameter(6,0.145491);
  f->FixParameter(7,8e-4);
  f->SetParLimits(8,0,1);
  f->SetParLimits(5,0,100000);
  h->Fit(Form("f_5p_%d",count5p),"LL");
  h->Fit(Form("f_5p_%d",count5p),"LL");
  h->Fit(Form("f_5p_%d",count5p),"LL","",minmass3prong,maxmass3prong);
  f->ReleaseParameter(6);
  f->ReleaseParameter(7);
  f->ReleaseParameter(9);
  f->SetParLimits(6,0.144,0.147);
  f->SetParLimits(7,1e-4,9e-4);
  f->SetParLimits(9,1e-4,9e-4);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.142,0.148);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.142,0.16);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.142,0.16);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.141,0.16);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.141,0.16);
  h->Fit(Form("f_5p_%d",count5p),"LL","",0.141,0.16);

  TF1* background = new TF1(Form("background_5p_%d",count5p),"[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x");
  background->SetParameter(0,f->GetParameter(0));
  background->SetParameter(1,f->GetParameter(1));
  background->SetParameter(2,f->GetParameter(2));
  background->SetParameter(3,f->GetParameter(3));
  background->SetParameter(4,f->GetParameter(4));
  background->SetLineColor(4);
  background->SetRange(minmass3prong,maxmass3prong);
  background->SetLineStyle(2);
  
  TF1* mass = new TF1(Form("fmass_5p_%d",count5p),"[0]*((1-[3])*TMath::Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+[3]*TMath::Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
  mass->SetParameters(f->GetParameter(5),f->GetParameter(6),f->GetParameter(7),f->GetParameter(8),f->GetParameter(9));
  mass->SetParError(0,f->GetParError(5));
  mass->SetParError(1,f->GetParError(6));
  mass->SetParError(2,f->GetParError(7));
  mass->SetParError(3,f->GetParError(8));
  mass->SetParError(4,f->GetParError(9));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);
    
  h->SetXTitle("M_{K#pi#pi#pi#pi}-M_{K#pi#pi#pi} (GeV/c^{2})");
  h->SetYTitle("Entries / (1/3 MeV/c^{2})");
  h->SetStats(0);
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();
  h->SetAxisRange(0,h->GetMaximum()*1.4*1.2,"Y");
  h->GetXaxis()->SetTitleOffset(1.3);
  h->GetYaxis()->SetTitleOffset(1.8);
  h->GetXaxis()->SetLabelOffset(0.007);
  h->GetYaxis()->SetLabelOffset(0.007);
  h->GetXaxis()->SetTitleSize(0.045);
  h->GetYaxis()->SetTitleSize(0.045);
  h->GetXaxis()->SetTitleFont(42);
  h->GetYaxis()->SetTitleFont(42);
  h->GetXaxis()->SetLabelFont(42);
  h->GetYaxis()->SetLabelFont(42);
  h->GetXaxis()->SetLabelSize(0.04);
  h->GetYaxis()->SetLabelSize(0.04);
  h->SetMarkerSize(0.8);
  h->SetMarkerStyle(20);
  h->SetStats(0);
  h->Draw("e");

  background->Draw("same");   
  mass->SetRange(0.142,0.152);	
  mass->Draw("same");
  f->Draw("same");
  
  Double_t yield = mass->Integral(minmass5prong,maxmass5prong)/binwidth5prong;
  Double_t yieldErr = mass->Integral(minmass5prong,maxmass5prong)/binwidth5prong*mass->GetParError(0)/mass->GetParameter(0);

  TLatex* tex;
  TLegend* leg = new TLegend(0.60,0.62,0.85,0.88,NULL,"brNDC");
  leg->SetBorderSize(0);
  leg->SetTextSize(0.04);
  leg->SetTextFont(42);
  leg->SetFillStyle(0);
  leg->AddEntry((TObject*)0,"D* D^{0}(K#pi#pi#pi)#pi",NULL);
  leg->AddEntry(h,"Data","pl");
  leg->AddEntry(f,"Fit","l");
  leg->AddEntry(mass,"D*^{+}+D*^{-} Signal","f");
  leg->AddEntry(background,"Combinatorial","l");
  leg->Draw("same");

//.........这里部分代码省略.........

//.........这里部分代码省略.........
            for(int ibin=1; ibin < histos[i]->GetNbinsX(); ibin++){
                histos[i]->SetBinError(ibin, histos[i]->GetBinContent((int)(0.892*(250.0/2.0)))*0.05);
                newHistos[i]->SetBinContent(ibin, histos[i]->GetBinContent(ibin));
                newHistos[i]->SetBinError(ibin, histos[i]->GetBinError(ibin));
            }
            
            pad = (TPad*)canvas[i]->cd(nfile+1);
            histos[i]->SetLineColor(1);
            histos[i]->SetLineWidth(1);
            histos[i]->GetXaxis()->SetRangeUser(0.7, 1.2);
            histos[i]->GetYaxis()->SetRangeUser(0, 1.5*histos[i]->GetBinContent(histos[i]->GetMaximumBin()));
            //histos[i]->SetStats(kFALSE);
            
            //histos[i]->Draw("HIST");

            printf("mean PT: %f\n", meanPT);

            TF1 *fit = new TF1(Form("fitPTbin%d00particle%d", nfile*2+1, i), FitFunRelBW, 0.68, 1.05, 5);
            //TF1 *fit = new TF1(Form("fitPTbin%d00particle%d", nfile*2+1, i), "gaus(0)", 0.86, 0.92);


            fit->SetParNames("BW Area", "Mass", "Width", "PT", "Temp");
            fit->SetParameters(TMath::Power(10.0, (float)(nfile)/1.7), 0.89, 0.1, 0.5, 0.130);
            //fit->SetParNames("BW Area", "Mass", "Width");
            //fit->SetParameters(100, 0.89, 0.0474);
            //fit->SetParLimits(0, -10, 1.5e9);
            Float_t max = histos[i]->GetXaxis()->GetBinCenter(histos[i]->GetMaximumBin());
            //if(max < 0.91 && max > 0.892){
            //    fit->SetParLimits(1, max-0.001, max+0.001);
            //}else{
                fit->SetParLimits(1, 0.82, 0.98);
            //}
            //fit->SetParLimits(2, 0.005, 0.15);
            fit->FixParameter(2, 0.042);
            fit->FixParameter(3, meanPT);
            //fit->SetParLimits(4, 0.05, 0.2);
            fit->FixParameter(4, 0.100001);
            fit->SetLineColor(2);

            printf("%s\n", fit->GetName());

            histos[i]->Fit(Form("fitPTbin%d00particle%d", nfile*2+1, i), "BRIM", "SAME");
            TVirtualFitter *fitter = TVirtualFitter::GetFitter();
           
            histos[i]->SetStats(1);
            histos[i]->Draw();
            gPad->Update();
            pad->Update();
            st = (TPaveStats*)histos[i]->FindObject("stats");
            st->SetX1NDC(0.524);
            st->SetY1NDC(0.680);
            st->SetX2NDC(0.884);
            st->SetY2NDC(0.876);
            //fit->Draw("SAME");
            //histos[i]->Draw();
            gPad->Update();
            pad->Update();
            printf("\n");
    
            diffHistos[i] = (TH1D*)histos[i]->Clone(Form("diffPTbin%d00particl%d", nfile*2+1, i));
            diffHistos[i]->Add(fit, -1);
            diffCanvas[i]->cd(nfile+1);
            diffHistos[i]->Draw("HIST E");
            diffHistos[i]->Write();

            //counting bins

//.........这里部分代码省略.........
  double errEff_cat_S1_1[100]={0},errEff_cat_S1_2[100]={0};
  double errEff_cat_S2_1[100]={0},errEff_cat_S2_2[100]={0};
  char PlotName[500],PlotName1[500], PlotName2[500]; 
  char GPlotName[500],GPlotName1[500], GPlotName2[500];

  for (Int_t ih = 0; ih < Nptbin; ih++) {
   
    cout<<" no of gen dimuons from diMuons Pt histo    "<<diMuonsPt_GenA1[ih]->Integral(1,100)<<endl;
    cout<<" no of gen dimuons from diMuons Mass histo  "<<diMuonsInvMass_GenA1[ih]->Integral(1,100)<<endl;
    

    //from pT histogram
    //gen_pt[ih] =diMuonsPt_GenA1[ih]->IntegralAndError(1,100,genError);
   
    gen_pt[ih] = diMuonsInvMass_GenA1[ih]->IntegralAndError(1,100,genError);
    gen_ptError[ih]= genError;
    
    if(iSpec==1) sprintf(PlotName,"plots/DiMuonMass_PtBin_%d.png",ih);
    if(iSpec==2) sprintf(PlotName,"plots/DiMuonMass_RapBin_%d.png",ih);
    if(iSpec==3) sprintf(PlotName,"plots/DiMuonMass_CentBin_%d.png",ih);
    
    if(iSpec==1) sprintf(PlotName1,"plots/DiMuonMass_PtBin_%d.pdf",ih);
    if(iSpec==2) sprintf(PlotName1,"plots/DiMuonMass_RapBin_%d.pdf",ih);
    if(iSpec==3) sprintf(PlotName1,"plots/DiMuonMass_CentBin_%d.pdf",ih);
    
    if(iSpec==1) sprintf(PlotName2,"plots/DiMuonMass_PtBin_%d.eps",ih);
    if(iSpec==2) sprintf(PlotName2,"plots/DiMuonMass_RapBin_%d.eps",ih);
    if(iSpec==3) sprintf(PlotName2,"plots/DiMuonMass_CentBin_%d.eps",ih);
 
    //giving inetial value for crystall ball fourth parameter 
    diMuonsInvMass_RecA1[ih]->Rebin(2);
    //GAUSPOL->SetParameter(0, diMuonsInvMass_RecA1[ih]->GetMaximum());
    GAUSPOL->SetParameter(0, diMuonsInvMass_RecA1[ih]->Integral(0,50));
    GAUSPOL->FixParameter(1, diMuonsInvMass_RecA1[ih]->GetBinWidth(1));
    
    new TCanvas; 
    diMuonsInvMass_RecA1[ih]->Fit("GAUSPOL","LLMERQ", "", mass_low, mass_high);
   
    double UpsilonMass = GAUSPOL->GetParameter(2);
    double UpsilonWidth = GAUSPOL->GetParameter(3);
   
    double UpsilonYield = GAUSPOL->GetParameter(0);
    double UpsilonYieldError = GAUSPOL->GetParError(0);

    double par[20];
    GAUSPOL->GetParameters(par);
    sprintf(namePt_1B,"pt_1B_%d",ih);
    backfun_1 = new TF1(namePt_1B, Pol2, mass_low, mass_high, 3);
    backfun_1->SetParameters(&par[4]);
   

    double MassLow=(UpsilonMass-3*UpsilonWidth);
    double MassHigh=(UpsilonMass+3*UpsilonWidth);
    
    int binlow =diMuonsInvMass_RecA1[ih]->GetXaxis()->FindBin(MassLow);
    int binhi =diMuonsInvMass_RecA1[ih]->GetXaxis()->FindBin(MassHigh);
    double binwidth=diMuonsInvMass_RecA1[ih]->GetBinWidth(1);
        
    //yield by function 
    //rec_pt[ih] = UpsilonYield;
    //rec_ptError[ih]=UpsilonYieldError;

    cout<<"Rec diMuons from Pt histo "<<diMuonsPt_RecA1[ih]->Integral(1,100)<<endl;  
    cout<<"Rec dimuons from mass "<<diMuonsInvMass_RecA1[ih]->Integral(1,100)<<endl; 
    

//.........这里部分代码省略.........
  hModulationAv->Fit(fmod); 

  printf("Estimates from time delay: Distance = %f +/- %f m -- Angle = %f +/- %f deg\n",fmod->GetParameter(1),fmod->GetParError(1),fmod->GetParameter(2),fmod->GetParError(2));

  h->SetStats(0);

  hDeltaThetaBack->Sumw2();
  hDeltaPhiBack->Sumw2();
  hThetaRelBack->Sumw2();
  hDeltaThetaBack->Scale(0.1);
  hDeltaPhiBack->Scale(0.1);
  hThetaRelBack->Scale(0.1);
  hAngleBack->Scale(0.1);
  hAngle->Add(hAngleBack,-1);

  printf("bin counting: SIGNAL = %f +/- %f\n",hDeltaPhi->Integral()-hDeltaPhiBack->Integral(),sqrt(hDeltaPhi->Integral()));
  rate = (hDeltaPhi->Integral()-hDeltaPhiBack->Integral())/nsecGR*86400;
  rateErr = sqrt(hDeltaPhi->Integral())/nsecGR*86400;


  Float_t val,eval;
  TCanvas *c1=new TCanvas();
  TF1 *ff = new TF1("ff","[0]*[4]/[2]/sqrt(2*TMath::Pi())*TMath::Exp(-(x-[1])*(x-[1])*0.5/[2]/[2]) + [3]*[4]/6/[2]");
  ff->SetParName(0,"signal");
  ff->SetParName(1,"mean");
  ff->SetParName(2,"sigma");
  ff->SetParName(3,"background");
  ff->SetParName(4,"bin width");
  ff->SetParameter(0,42369);
  ff->SetParameter(1,0);
  ff->SetParLimits(2,10,maxwidth);
  ff->SetParameter(2,350); // fix witdh if needed
  ff->SetParameter(3,319);
  ff->FixParameter(4,(tmax-tmin)/nbint); // bin width

  ff->SetNpx(1000);
  
  if(cout) cout->cd();
  h->Fit(ff,"EI","",-10000,10000);
  
  val = ff->GetParameter(2);
  eval = ff->GetParError(2);
  
  printf("significance = %f\n",ff->GetParameter(0)/sqrt(ff->GetParameter(0) + ff->GetParameter(3)));

  h->Draw();
  
  new TCanvas;

  TF1 *func1 = (TF1 *)  h->GetListOfFunctions()->At(0);
  
  func1->SetLineColor(2);
  h->SetLineColor(4);
  
  TPaveText *text = new TPaveText(1500,(h->GetMinimum()+(h->GetMaximum()-h->GetMinimum())*0.6),9500,h->GetMaximum());
  text->SetFillColor(0);
  sprintf(title,"width = %5.1f #pm %5.1f",func1->GetParameter(2),func1->GetParError(2));
  text->AddText(title);
  sprintf(title,"signal (S) = %5.1f #pm %5.1f",func1->GetParameter(0),func1->GetParError(0));
  text->AddText(title);
  sprintf(title,"background (B) (3#sigma) = %5.1f #pm %5.1f",func1->GetParameter(3),func1->GetParError(3));
  text->AddText(title);
  sprintf(title,"significance (S/#sqrt{S+B}) = %5.1f",func1->GetParameter(0)/sqrt(func1->GetParameter(0)+func1->GetParameter(3)));
  text->AddText(title);
  
  text->SetFillStyle(0);

TF1 *GausBF::Bfit(TGraph *gr)
{
  TF1 *func = new TF1("func", FB(), -200, 200);
  gr->Fit("gaus", "q0");
  TF1 *fg0 = gr->GetFunction("gaus");

  Int_t ip = 3;
  for (Int_t i = 0; i < 3; i++) func->SetParameter(i, fg0->GetParameter(i));

  TGraph gd;
  TGraph gm;

  Double_t min = 0, max = 0;
  Double_t xmn = gr->GetY()[0];
  Double_t xmx = gr->GetY()[0];
  Double_t xb  = gr->GetX()[0];
  Double_t db  = 0;

  Double_t mmin = 0.002;

  for (Int_t i = 0; i < gr->GetN(); i++) {
    Double_t x = gr->GetX()[i];
    Double_t y = gr->GetY()[i];
    Double_t d = y-func->Eval(x);

    gd.SetPoint(i, x,  d);
    gm.SetPoint(i, x, -d);

    if (db*d < 0) {
      if ((d > 0 && min < -mmin) || (d < 0 && max > mmin)) {
	Double_t xm = (d > 0) ? xmn : xmx;
	Double_t w1 = x-xm;
	Double_t w2 = xm-xb;
	Double_t w  = (w1 > w2) ? w2 : w1;
	Double_t par[3];
	if (d < 0) {
	  gd.Fit("gaus", "q0", "", xm-w, xm+w);
	  TF1 *fg = gd.GetFunction("gaus");
	  for (Int_t j = 0; j < 3; j++) par[j] = fg->GetParameter(j);
	}
	else {
	  gm.Fit("gaus", "q0", "", xm-w, xm+w);
	  TF1 *fg = gm.GetFunction("gaus");
	  for (Int_t j = 0; j < 3; j++) par[j] = fg->GetParameter(j);
	  par[0] = -par[0];
	}

	for (Int_t j = 0; j < 3; j++) func->SetParameter(ip+j, par[j]);
	ip += 3;
      }
      xb = x;
      if (d < 0) min = 0;
      if (d > 0) max = 0;
    }
    db = d;

    if (d < min) { min = d; xmn = x; }
    if (d > max) { max = d; xmx = x; }
  }
  for (Int_t i = ip; i < Np; i++) func->FixParameter(i, (i%3 == 2) ? 1 : 0);

  gr->Fit(func, "q0");

  return func;
}