C++ TGraph::GetFunction方法代码示例

TGraph* autogain152(TH1 *hist) {

   hist->GetXaxis()->SetRangeUser(200.,16000.);
   TSpectrum *s = new TSpectrum();
   Int_t nfound = s->Search(hist,6,"",0.08); //This will be dependent on the source used.
   printf("Found %d candidate peaks to fit\n",nfound);
   if(nfound > 6)
      nfound = 6;

   std::vector<float> vec;
   for(int x=0;x<nfound;x++)
      vec.push_back(s->GetPositionX()[x]);

   std::sort(vec.begin(),vec.end());

   Float_t energies[] = {121.7830, 244.6920, 344.276, 778.903, 964.131, 1408.011};
   TGraph* slopefit = new TGraph(nfound, &(vec[0]), energies);

   printf("Now fitting: Be patient\n");
   slopefit->Fit("pol1");
   if(slopefit->GetFunction("pol1")->GetChisquare() > 10.) {
      slopefit->RemovePoint(slopefit->GetN()-1);
      slopefit->Fit("pol1");
   }
   TChannel *chan = 0;
   slopefit->Draw("AC*");

   return slopefit;
}

Double_t beta2_for_5sigma(const Int_t fNPts, const Double_t fRangeMin, const Double_t fRangeMax, const Double_t fN_sig_100, const Double_t fN_bkg_100, const Double_t fSigma_N_bkg, const string& fTitle) 
{
 
 Double_t x[fNPts], y[fNPts];

 Double_t step = (fRangeMax - fRangeMin)/(fNPts-1);
   
 for (Int_t i = 0; i < fNPts; i++) {
 
  x[i] = fRangeMin + step*i;
  y[i] = ScP(x[i], fN_sig_100, fN_bkg_100, fSigma_N_bkg);
 }
 
 TCanvas *c_temp = new TCanvas("c_temp","",1120,800);
 c_temp->cd();

 string title = fTitle + ";#beta^{2};S_{cP}";

 TH2F *bg_temp = new TH2F("bg_temp",title.c_str(), 100, fRangeMin, fRangeMax, 100, 0.8*y[0], 1.2*y[fNPts-1]);
 bg_temp->SetStats(kFALSE);
 bg_temp->SetTitleOffset(1.,"X");
 bg_temp->SetTitleOffset(1.,"Y");
 bg_temp->Draw();
 
 TGraph *scP = new TGraph(fNPts, x, y);
 scP->SetMarkerSize(1.);
 scP->SetMarkerStyle(24);
 scP->SetMarkerColor(kRed);
 scP->Draw("P"); 
 scP->Fit("pol2");

 TF1 *fit = (TF1*)scP->GetFunction("pol2");
 Double_t beta2 = fit->GetX(5);
 Double_t N_sig = beta2*fN_sig_100;
 Double_t N_s_b = N_sig + fN_bkg_100;
 
 cout<<">> beta2 for 5 sigma discovery for "<<fTitle<<" = "<<beta2<<"\n";
 cout<<">> ** N_sig = "<<N_sig<<"\n";
 cout<<">> ** N_bkg = "<<fN_bkg_100<<"\n";
 cout<<">> ** N_s_b = "<<N_s_b<<"\n";
 
 if(fSigma_N_bkg != 0) c_temp->SaveAs((fTitle + "_significance_beta_sys.png").c_str());
 else c_temp->SaveAs((fTitle + "_significance_beta.png").c_str());

 delete scP;
 delete bg_temp;
 delete c_temp;
 
 return beta2;
}

/// used in display_beamprofile() : not working !
void getEllipseParameters(const float * x_data, const float * y_data, const unsigned int N, float& x_width, float& y_width, float& angle) {
	// In order to fit a good ellipse on the scattered plot :
	// 1) The TH2 is copied into a TGraph, to fit it with y(x) = ax => to retrieve the angle
	// 2) Rotation of the Graph to get the RMS in X and Y
	// 3) Creation of the final ellipse, with the good widths and angle

	TCanvas * ca0 = new TCanvas;
	ca0->Divide(2,1);
	ca0->cd(1);

	TGraph * draft = new TGraph(N,x_data,y_data);
	draft->Draw("AP");
	draft->Fit("pol1","Q");
	TF1 * pol1 = draft->GetFunction("pol1");
	pol1->Draw("same");

	// gets the angle [rad]
	angle = asin(1.) - atan(pol1->GetParameter(1));
	
	double x_datarot[N], y_datarot[N];
	for (unsigned int i=0; i<N; i++) {
		x_datarot[i]=  x_data[i]*cos(angle) - y_data[i]*sin(angle);
		y_datarot[i]=  x_data[i]*sin(angle) + y_data[i]*cos(angle);
	}

	ca0->cd(2);
	TGraph * draft2 = new TGraph(N,x_datarot,y_datarot);
	draft2->Draw("AP");
	x_width =  draft2->GetRMS(1); 
	y_width =  draft2->GetRMS(2);
	angle = 180-90*angle/asin(1.);
//	draft->Draw("AP");
	ca0->cd(1);
	TEllipse * ell = new TEllipse(draft->GetMean(1),draft2->GetMean(2),x_width*3,y_width*3);
	ell->SetTheta(angle);
	ell->Draw("same");

	//cout << "x = " << x_width << "\t y = " << y_width << "\t angle = " << angle << endl;

//	delete draft2;
//	delete draft;
//	delete ca0;
	return;
}

void DoFit(TH2 *hist) {
   TChannel::DeleteAllChannels();
   TChannel::ReadCalFile("GrifCal.cal");
   printf("made %i channels.\n",TChannel::GetNumberOfChannels());
   TNucleus nuc("152eu");
   for(int x = 1; x <= 64; ++x) {
      printf(" x = %i\n",x);
      TH1 *p = GetProjectionY(hist,x);
      if(p->GetEntries() < 100)
         continue;
      TGraph* graph = autogain(p,&nuc);
      TChannel *chan = TChannel::GetChannelByNumber(x);
      if(!chan)
         continue;
      chan->DestroyCalibrations();
      chan->AddENGCoefficient(graph->GetFunction("pol1")->GetParameter(0));
      chan->AddENGCoefficient(graph->GetFunction("pol1")->GetParameter(1));
      chan->SetIntegration(125);
   }

   TChannel::WriteCalFile("NewGrifCal.cal");

}

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

int AngOptimize()
{
  ifstream infile("angSNR.dat");
  ofstream outfile("OptAng.dat",ios::app);

  istringstream iss;
  char line[1000];

  cout<<"aaaaaaaa"<<endl;
  const int nene=21;
  const int nep = 28;
  double epn[nene][nep];
  double sigN[nene][nep];
  double bckN[nene][nep];
  double snr[nene][nep];
  int iene=-1;
  int iep=-1;
  double enes[nene];
  double opteps[nene];
  char title[nene][1000];

  TF1 *fit = new TF1("fit","[0]*(x-[2])+[1]/(x-[2])+[3]",178,179.6);
  fit->SetParameter(0,1.5);
  fit->SetParameter(1,1);
  fit->SetParameter(2,180);
  fit->SetParLimits(2,180,180)
  fit->SetParameter(3,40);

  bool start=0;
  if (!infile.is_open()) return -1;
  while (!infile.eof()){
    //cout<<"b"<<endl;
    infile.getline(line,1000);
	if (line[0]!='K') 
	{
	  if (!start) continue;
	  if (line[0]=='#' || line[0]==' ') continue;
	}
	else {
	  iep=-1; start=1;
	  iene++;
	  if (iene>=21) break;
	  string tit;
	  iss.clear();
	  iss.str(line);
          iss >> tit;
	  std::cout<<tit<<endl;
	  sprintf(title[iene],"%s",tit.c_str());
	  continue;
	}
	iep++;
	//cout<<line<<endl;
	iss.clear();
	iss.str(line);
	double ratio,signal,background;
	iss >> ratio >> signal >> background; // 3 sigma
	//iss >> signal >> background; // 5 sigma

	epn[iene][iep] = ratio;
	sigN[iene][iep] = signal;
	bckN[iene][iep] = background;
	
	if (sigN[iene][iep]==0) snr[iene][iep] = 0;
	else snr[iene][iep] = sigN[iene][iep]/sqrt(sigN[iene][iep]+bckN[iene][iep]);
	//cout<< epn[iene][i]<<sigN[iene][i]<<bckN[iene][i]<<endl;
  }

  for (iene=0;iene<nene;iene++){
	TGraph *graph = new TGraph(nep,epn[iene],snr[iene]);
	graph->SetTitle(title[iene]);
	graph->SetMarkerStyle(25);
	graph->GetYaxis()->SetTitle("S/#sqrt{S+N}");
	graph->GetXaxis()->SetTitle("#theta");
	TCanvas *ca1 = new TCanvas(title[iene], title[iene]);
	graph->Draw("AP");
	//graph->Fit("pol4","","",178,179.9);
////////TF1 *fpol2 = graph->GetFunction("pol2");
////////fit->SetParameter(0,fpol2->GetParameter(0));
////////fit->SetParameter(1,fpol2->GetParameter(1));
////////fit->SetParameter(2,fpol2->GetParameter(2));
	//graph->Fit(fit,"","",178,179.6);
	graph->Fit(fit,"","",172.0+iene*6.0/21.0,179.9);
	double optep = fit->GetMaximumX(176,179.6);
	outfile<<"Optmized theta at "<< title[iene] << " is "<<optep<<std::endl;
	enes[iene] = atof(&title[iene][3])/10000.;
	opteps[iene] = optep;
	char ofname[1000];
	sprintf(ofname,"output/Angopt_%s.pdf",title[iene]);
	ca1->Print(ofname);
  } 
  TGraph* gep = new TGraph(iene,enes,opteps);
  gep->SetTitle("Optimized #theta");
  gep->SetMarkerStyle(25);
  gep->GetYaxis()->SetTitle("S/#sqrt{S+N}");
  gep->GetXaxis()->SetTitle("#theta");
  TCanvas *ca1 = new TCanvas("gep","Optimized #theta");
  gep->Draw("AP");
  gep->Fit("pol1");
  ca1->Print("output/Angopt.pdf");
  double c0 = gep->GetFunction("pol1")->GetParameter(0);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  float muMiGlobalChi;
  int muMiNMuonhits;
  bool muMiGoodMu;
  bool muMiTrkArb;
  bool muMiTMOneStaTight;
  bool muMiHighPurity;
  
  
  int nx = 14, ny = 15000;
  //double xBins[22] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,30.0}; //19
  double xBins[15] = {3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,15.0,18.0,22.0,30.0}; //19
  
  //double xBins[4] = {3.0,6.5,12.0,30.0}; //19
  //double xBins[6] = {6.5,9.0,12.0,15.0,20.0,30.0}; //19
  
  //double xBins[13] = {6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,15.0,18.0,21.0,24.0,30.0}; //13
  //double xBins[8] = {6.0,9.0,12.0,15.0,18.0,21.0,24.0,30.0}; //19
  //TH2D *hPR2D = new TH2D("hPR2D","hPR2D",nx,xBins,ny,-5.0,10.0);
  
  //TH1F *hPR = new TH1F("hPR","hPR;Inv_mass",100,-2.5,3.2);
  TH1F *hPR = new TH1F("hPR","hPR;p_{T}",nx,xBins);
  //TH1F *hPR = new TH1F("hPR","hPR;l_{J/#psi} (mm);counts",1500,-5,10);
  hPR->Sumw2();
  hPR->SetLineColor(1);
  //TH1F *hPRfunc = new TH1F("hPRfunc","hPRfunc;Inv_mass",100,-2.8,3.2);
  TH1F *hPRfunc = new TH1F("hPRfunc","hPRfunc;p_{T}",nx,xBins);
  //TH1F *hPRfunc = new TH1F("hPRfunc","hPRfunc;l_{J/#psi} (mm);counts",1500,-5,10);
  hPRfunc->Sumw2();
  hPRfunc->SetLineColor(2);
  //TH1F* hPRCut = (TH1F*) hPR->Clone("hPRCut");
  
  TFile *Fitfile = new TFile("Plots/2015/3D/ctauFit/Jpsi_pp_eff_0.9_Rap_1.6-2.4_Pt_3.0-30.0.root");
  TGraph *FitGr = (TGraph*)Fitfile->Get("Graph");
  TF1 *func1 = (TF1*)FitGr->GetFunction("FitFn");
  
  new TCanvas;
  TAxis *XaxctauFW = FitGr->GetXaxis();
  XaxctauFW->SetLimits(0.0,30.0);
  TAxis *YaxctauFW = FitGr->GetYaxis();
  YaxctauFW->SetLimits(0.0,0.1);
  FitGr->GetXaxis()->SetRangeUser(0.0,30.0);
  FitGr->GetYaxis()->SetRangeUser(0.0,0.1);
  //FitGr->Draw("AP");
  //FitGr->Draw();
  func1->Draw();

  initOniaTree(iTree);
  Long64_t nentries = iTree->GetEntries();
  nentries = 500000;          
  cout<<"Total Entries in data  Tree  "<<"  "<<nentries<< "====="<<endl;
  Long64_t nbytes = 0, nb = 0;
 
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    nb = iTree->GetEntry(jentry);   
    
    for (int iQQ=0; iQQ<Reco_QQ_size; iQQ++) {
      TLorentzVector *qq4mom = (TLorentzVector*) Reco_QQ_4mom->At(iQQ);
      TLorentzVector *qq4mupl = (TLorentzVector*) Reco_QQ_mupl_4mom->At(iQQ);
      TLorentzVector *qq4mumi = (TLorentzVector*) Reco_QQ_mumi_4mom->At(iQQ);
      
      invariantMass = qq4mom->M();
      JpsiPt = qq4mom->Pt();
      JpsiRap = qq4mom->Rapidity();

      muPlPt = qq4mupl->Pt();
      muMiPt = qq4mumi->Pt();

void dumpProfile() {
  
 std::string fileName_ = "Profile_SM10.root";
 TFile *shapeFile_ = TFile::Open(fileName_.c_str(),"old");
 TProfile* PROF_704 = (TProfile*) shapeFile_->Get("SHAPE_XTAL_704");

 ofstream out;
 out.open("dat.txt");

 int nBinsHisto_ = 250;
 std::vector<double> shapeArray(nBinsHisto_,0.0);

 double max = -999;
 int imax = 0;
 for(int ibin=0; ibin < nBinsHisto_; ++ibin) 
   {
     out << "shapeArray[" << ibin << "] = " << PROF_704->GetBinContent(ibin+1) << " ; \n";
     shapeArray[ibin] = PROF_704->GetBinContent(ibin);
     std::cout << "Original shape, ns = " << ibin << " shape = " << shapeArray[ibin] << std::endl;
     if ( shapeArray[ibin] > max ) {
       max = shapeArray[ibin];
       imax = ibin;
     }

   }//loop

 out.close();

 double xMinHisto_ = -1.;
 double xMaxHisto_ = 9.;
 double binw = (xMaxHisto_ - xMinHisto_)/(shapeArray.size());
 int nbins = shapeArray.size()/10;

 float low =  xMinHisto_+(double)(imax-nbins/2+0.5)*binw;
 float up = xMinHisto_+(double)(imax+nbins/2+0.5)*binw;
 
 double* x = new double[nbins];
 double* y = new double[nbins];
 for (int i = 0; i < nbins; i++) {
   x[i] = xMinHisto_ + (double)(imax - nbins/2 + i + 0.5)*binw;
   y[i] = shapeArray[imax - nbins/2 + i];
   std::cout << " x,y = " << x[i] << " " << y[i] << " " << (double)(imax - nbins/2 + i + 0.5) << std::endl;
 }
 TGraph* graph = new TGraph(nbins, x, y);
 graph->Fit("pol3", "V");//"Q 0");
 TF1* fFit = graph->GetFunction("pol3");
 double tMax = fFit->GetMaximumX();

 std:;cout << "Maxiumum = " << tMax << std::endl;

 gStyle->SetOptFit(1111);

 TCanvas *MyC = new TCanvas("MyC","Test canvas",1); 
 MyC->Divide(2,1); 
 MyC->cd(1); 
 PROF_704->Draw(); 
 MyC->cd(2); 
 fFit->Draw(); 
 MyC->SaveAs("PROF_704.jpg");
 
}

/*
root -l 'sigmapeak.C+(0.025)'
// 
Processing sigmapeak.C+(0.025)...
background only: bkg_nsigma_0_99 = -0.352803
sig_nsigma_0_59  = -0.608621
sig_nsigma_0_99  = 2.1472
sig_nsigma_60_99 = 4.14042
sig_nsigma_70_90 = 5.57689
sig_nsigma_75_85 = 8.056
*/
void sigmapeak(Double_t signal_area=0.025)
{
   Double_t bkg_mean = 0;
   Double_t bkg_sigma = 0.001;

   Double_t x[100];
   Double_t y[100];
   Int_t np = 100;
   
   TRandom3 rand;
   
   for (int i=0; i<np; ++i) {
      x[i] = i;
      y[i] = rand.Gaus(bkg_mean,bkg_sigma);
   }
   
   TGraph* gr = new TGraph(np,x,y);
   gr->SetNameTitle("gr",Form("#sigma = %0.1f",bkg_sigma));
   gr->SetMarkerStyle(7);
   gr->SetMarkerColor(46);
   new TCanvas();
   gr->Draw("ap");

   Double_t bkg_nsigma_0_99 = Nsigma(gr->GetY(), 0,99, bkg_sigma);
   cout<< "background only: bkg_nsigma_0_99 = " << bkg_nsigma_0_99 <<endl;

   // add signal

   Double_t signal_mean = 80;
   Double_t signal_sigma = 3;

   for (int i=0; i<gr->GetN(); ++i) {
      Double_t xx = (gr->GetX()[i] - signal_mean)/signal_sigma;
      Double_t arg = 0.5*xx*xx;
      Double_t exp = arg < 50? TMath::Exp(-arg): 0;
      Double_t signal = signal_area/(TMath::Sqrt(TMath::TwoPi())*signal_sigma) * exp;
      gr->SetPoint(i, gr->GetX()[i], gr->GetY()[i] + signal);
   }

   gr->SetTitle(Form("#sigma_{bkg} = %0.3f signal: area = %0.3f mean = %0.0f sigma = %0.1f", bkg_sigma,signal_area,signal_mean,signal_sigma));

   gr->Draw("apl");
   gr->Fit("gaus", "R", "", signal_mean - 5*signal_sigma, signal_mean+5*signal_sigma);
   Double_t fit_area = 2.5 * gr->GetFunction("gaus")->GetParameter("Constant") * gr->GetFunction("gaus")->GetParameter("Sigma");
   cout<< "Area under fitted gaussian = " << fit_area <<endl;

   // titmax();
   gPad->Modified();    // to create box (NB: the pad was not drawn yet at this point!)
   gPad->Update();
   TPaveText* tit = (TPaveText*)gPad->GetPrimitive("title");
   tit->SetX1NDC(0.);
   tit->SetX2NDC(1.);
   tit->SetY1NDC(0.9);
   tit->SetY2NDC(1.);
   gPad->Modified();    // to update the pad
   gPad->Update();

   Double_t sig_nsigma_0_59 = Nsigma(gr->GetY(), 0,59, bkg_sigma);
   cout<< "sig_nsigma_0_59  = " << sig_nsigma_0_59 <<endl;

   Double_t sig_nsigma_0_99 = Nsigma(gr->GetY(), 0,99, bkg_sigma);
   cout<< "sig_nsigma_0_99  = " << sig_nsigma_0_99 <<endl;

   Double_t sig_nsigma_60_99 = Nsigma(gr->GetY(), 60,99, bkg_sigma);
   cout<< "sig_nsigma_60_99 = " << sig_nsigma_60_99 <<endl;

   Double_t sig_nsigma_70_90 = Nsigma(gr->GetY(), 70,90, bkg_sigma);
   cout<< "sig_nsigma_70_90 = " << sig_nsigma_70_90 <<endl;

   Double_t sig_nsigma_75_85 = Nsigma(gr->GetY(), 75,85, bkg_sigma);
   cout<< "sig_nsigma_75_85 = " << sig_nsigma_75_85 <<endl;

   Double_t ys5[100];
   smooth5(np, gr->GetY(), ys5);
   TGraph* gr5 = new TGraph(np, x, ys5);
   gr5->SetNameTitle("gr5","smoothed on 5 points");
   gr5->SetMarkerStyle(7);
   gr5->SetMarkerColor(2);
   gr5->SetLineColor(2);
   new TCanvas;
   gr5->Draw("apl");

   Double_t ys7[100];
   smooth7(np, gr->GetY(), ys7);
   TGraph* gr7 = new TGraph(np, x, ys7);
   gr7->SetNameTitle("gr7","smoothed on 7 points");
   gr7->SetMarkerStyle(7);
   gr7->SetMarkerColor(8);
   gr7->SetLineColor(8);
//.........这里部分代码省略.........