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

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

void fitSys(char *infname = "background_PbPb.dat")
{
    TGraph *g = new TGraph(infname);
    TF1 *f = new TF1("f","[0]+[1]/(x)+[2]/x/x+[3]*x");
    g->Draw("ap");
    g->Fit("f");
    g->Fit("f");
    g->Fit("f");
}

double get_correctionFactorbb(){

 double xx[5] = {39, 62.4, 200, 2760, 5000};
 double yy[5] = {0.00944, 0.0709, 1.81, 94.92, 180};

 TGraph *g = new TGraph(5);
 for(int i=0; i<5; i++){
        g->SetPoint(i,xx[i],yy[i]);
 }

 TCanvas *c = new TCanvas("c","c",600,450);
 g->SetMarkerStyle(20);
 g->Draw("AP");
 c->SetLogy();
 c->SetLogx();

 TF1 *f = new TF1("f","[0]+[1]*x+[2]*x*x",0,300);
 //TF1 *f = new TF1("f","[0]*pow(x,[1])",0,5500);
 //f->SetParameters(0,-5);
 g->Fit(f,"R");

 double corr = (f->Eval(193))/(f->Eval(200));

 return corr;

}

/*============================================================================*/
void gaus1peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[5],epar[5],x[4],y[4];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  setcanvas(1);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x2;
  x[2] = x3;
  x[3] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x2);
  y[1] = hist->GetBinContent(bin2);
  Int_t bin3 = hist->FindBin(x3);
  y[2] = hist->GetBinContent(bin3);
  Int_t bin4 = hist->FindBin(x4);
  y[3] = hist->GetBinContent(bin4);
  TGraph *g = new TGraph(4,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[3]=fpol1->GetParameter(0);
  par[4]=fpol1->GetParameter(1);

  //--**-- Gaussian Peak estimation without background --**--//
  TF1 *fgaus = new TF1("GAUS","gaus",x2,x3);
  hist->Fit(fgaus,"RQN");
  fgaus->GetParameters(&par[0]);

  //--**-- Final Peak Fit with Background --**--//
  TF1 *func = new TF1("FGAUS","gaus(0)+pol1(3)",x1,x4);
  func->SetParameters(par);
  hist->Fit(func,"R+QN");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  Double_t fwhm = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N0 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r0 = epar[0]/par[0];
  Double_t r2 = epar[2]/par[2];
  Double_t eN0= N0*TMath::Sqrt(r0*r0+r2*r2);
  printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
          par[1],epar[1],fwhm,efwhm,N0,eN0);
  //printf("%11.4f %11.4f %11.0f %11.0f\n",
  //        par[1],epar[1],N0,eN0);
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
}

TGraph* autogain(TH1 *hist,TNucleus *nuc) {    //Display The fits on a TPad  

   if(!hist || !nuc)
      return 0;

   nuc->SetSourceData();

   if(nuc->GetA() == 152) {
      return autogain152(hist);
   }

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

   nuc->TransitionList.Sort();

   std::vector<float> engvec;
   TIter iter(&(nuc->TransitionList));
   TObject* obj;
   while(obj = iter.Next()) {
      if(!obj->InheritsFrom("TGRSITransition"))
         continue;
      TGRSITransition *tran = (TGRSITransition*)obj;

      engvec.push_back(static_cast<float>(tran->energy));
      if(engvec.size() == nfound)
         break;
   }

   if(nfound != engvec.size())
      return 0;

   Float_t *posPeaks = s->GetPositionX();
   Float_t *energies = &(engvec[0]);

   for(int x=0;x<nfound;x++) {
      printf("posPeaks[%i] = %f\t\tenrgies[%i] = %f\n",x,posPeaks[x],x,energies[x]);
   }

   TGraph *slopefit = new TGraph(nfound,posPeaks,energies ); 

   printf("Now fitting: Be patient\n");
   slopefit->Fit("pol1");
   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;
}

Double_t beta2_for_exclusion(const Int_t fNPts, const Double_t fxsTh, const Double_t fRangeMin, const Double_t fRangeMax, const Double_t fSigma_L, const Double_t fS_eff, const Double_t fSigma_S_eff, 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] = CLA(100, 100*fSigma_L, fS_eff*x[i], fS_eff*x[i]*fSigma_S_eff, fN_bkg_100, fN_bkg_100*fSigma_N_bkg);
 }
 
 TCanvas *c_temp = new TCanvas("c_temp","",1120,800);
 c_temp->cd();

 string title = fTitle + ";#beta^{2};95% C.L. upper limit on #sigma [pb]";

 TH2F *bg_temp = new TH2F("bg_temp",title.c_str(), 100, fRangeMin, fRangeMax, 100, 0.8*y[fNPts-1], 1.2*y[0]);
 bg_temp->SetStats(kFALSE);
 bg_temp->SetTitleOffset(1.,"X");
 bg_temp->SetTitleOffset(1.,"Y");
 bg_temp->Draw();
 
 TF1 *f1 = new TF1("f1","[0]+[1]/pow(x,[2])",fRangeMin,fRangeMax);
 f1->SetParameters(0.,1.,1.);
 f1->SetParLimits(2, 0.45, 1.1);
 
 TGraph *xsection = new TGraph(fNPts, x, y);
 xsection->SetMarkerSize(1.);
 xsection->SetMarkerStyle(24);
 xsection->SetMarkerColor(kRed);
 xsection->Draw("P");
 xsection->Fit("f1");

 Double_t beta2 = f1->GetX(fxsTh);

 cout<<">> beta^2 for 95% CL exclusion of "<<fTitle<<" = "<<beta2<<"\n";
 
 if(fSigma_N_bkg != 0) c_temp->SaveAs((fTitle + "_exclusion_sys.png").c_str());
 else c_temp->SaveAs((fTitle + "_exclusion.png").c_str());
 
 delete xsection;
 delete f1;
 delete bg_temp;
 delete c_temp;
 
 return beta2;
}

    TF1* computeFunctionFit(TrialDataSet& eSimData)
    {
         TF1 * func = new TF1("fittingFunction", this->fittingFunction, -39, 200, 3);
         //func->SetParameters(eSimData.startingTime - subtractedTime, 1.0);

         // Convert digital readout values to floats
         float floatReadoutValues[7];
         for (int i = 0; i < 7; ++i){floatReadoutValues[i] = eSimData.digitalReadoutValues[i];}


         TGraph *gr = new TGraph(7, eSimData.measurementTimes, floatReadoutValues);

         gr->Fit(func, "QN");
         return func;
    }

void TrPdf::FitLogLog(double min, double max) {
  TGraph* FluxLogLogTmp = new TGraph(GetGraph()->GetN());
  FluxLogLogTmp->SetName("spectrumloglogtmp");
  FluxLogLogTmp->SetTitle("spectrumloglogtmp");
  for (int ii=0; ii<Graph->GetN(); ii++) {
    double a,b;
    GetGraph()->GetPoint(ii,a,b);
    FluxLogLogTmp->SetPoint(ii,log10(a),log10(b));
  }
  TF1* LinFitTmp = new TF1("linfittmp","[0]+[1]*x+[2]*pow(x,2.)+[3]*pow(x,3.)+[4]*pow(x,4.)+[5]*pow(x,5.)",-2.,5.); 
  FluxLogLogTmp->Fit(LinFitTmp,"EQR","",log10(min),log10(max));
  LogLog = new TF1(Form("LogLog_%s",GetName().Data()),
    "pow(10.,[0]+[1]*log10(x)+[2]*pow(log10(x),2.)+[3]*pow(log10(x),3.)+[4]*pow(log10(x),4.)+[5]*pow(log10(x),5.))",1.e-2,1.e5);
  for (int i=0; i<6; i++) LogLog->SetParameter(i,LinFitTmp->GetParameter(i));
  delete LinFitTmp;
  delete FluxLogLogTmp;
}

void GraphCrdcPads(int num=0,int CRDCNum=0){

  //This script will look at the Raw root trees
  //and plot a single crdc pad distribution 

  S800Event * event= new S800Event();

  rawtree->SetBranchAddress("s800event",&event);

  rawtree->GetEntry(num);

  // cout<<"Size of samples "<<event->GetS800()->GetCrdc(0)->GetSample().size()<<endl;
  // cout<<"Size of data "<<event->GetS800()->GetCrdc(0)->GetData().size()<<endl;
  // cout<<"Size of channels "<<event->GetS800()->GetCrdc(0)->GetChannels().size()<<endl;

  int size = event->GetS800()->GetCrdc(CRDCNum)->GetChannels().size();



  TString install =gSystem->Getenv("R00TLeInstall");
  
  TString calfile = install+"/prm/crdccalNone.dat";
  TString pedfile = install+"/prm/crdcpedestals.dat";
  
  S800Calibration calibration;//new S800Settings(install+"/prm/Raw2Cal.dat"));
  calibration.ReadCrdcCalibration(calfile,pedfile);


  //calibration.SetCrdc(event->GetS800()->GetCrdc(0)->GetChannels(),event->GetS800()->GetCrdc(0)->GetData(),0,0,0);

  calibration.CrdcCal(event->GetS800()->GetCrdc(CRDCNum)->GetChannels(),event->GetS800()->GetCrdc(CRDCNum)->GetData(),0);

  int size2=calibration.GetCRDCCal().size();
  TGraph * graph = new TGraph();
  for (int i=0;i<size2;i++){
    if (TMath::IsNaN(calibration.GetCRDCCal()[i])){
      graph->SetPoint(i,i,0);
    }else {
      graph->SetPoint(i,i,calibration.GetCRDCCal()[i]);
    }
  }
  graph->Fit("gaus");
  graph->Draw("A*");
  return;
}

/// 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 rootFit(){

  TCanvas *c1 = new TCanvas("c1","",600,600);
  c1->Divide(1,2);
  TNtuple *T = new TNtuple("T","","x:y:z");
  T->ReadFile("histogramData.dat");
  T->Draw("y:x");

  Double_t *X = T->GetV1();
  Double_t *Y = T->GetV2();

  for(int loop = 0; loop < 120; loop++){
    cout << X[loop] << "  " << Y[loop] << endl;
  }

  c1->Clear();
  TGraph *graphT = new TGraph(120,Y,X);
  TGraph *graphB = new TGraph(120,Y,X);
  c1->cd(1);
  graphT->Draw("APL");
  graphT->Fit("gaus+pol2","RME");

}

//.........这里部分代码省略.........
  tgy1->SetLineColor(kBlue);
  tgy2->SetLineColor(kRed);
  tgy1->SetLineWidth(2);
  tgy2->SetLineWidth(2);
  tgy1->SetMinimum(0);
  tgy2->SetMinimum(0);

  tgy1->GetXaxis()->SetTitle("Distance (cm)");
  tgy2->GetXaxis()->SetTitle("Distance (cm)");
  tgy1->GetYaxis()->SetTitle("Number of photoelectrons");
  tgy2->GetYaxis()->SetTitle("Number of photoelectrons");

  double maxx = 0;
  if(max1>=max2)
    {
      // case 1
      maxx = max1;
      h1->Draw();
      h2->Draw("same");
    }
  else
    {
      // case 2
      maxx = max2;
      h2->Draw();
      h1->Draw("same");
    }

  TF1 *funx1 = new TF1("funx1","[0]+[1]*TMath::Exp(-x/[2])",0,25);
  funx1->SetLineColor(kBlue);
  funx1->SetParameter(0,5.0); // number of photoelectrons in core
  funx1->SetParameter(1,5.0); // number of photoelectrons in clad
  funx1->FixParameter(2,cld); // clad decay constant
  tgy1->Fit(funx1,"","",1,24); // need to use TGraph here, fitting h1 produces "Warning in <Fit>: Fit data is empty"
                               // I don't know why this happening, probably some dumb mistake I'm making
  TF1 *funx2 = new TF1("funx2","[0]+[1]*TMath::Exp((x-25)/[2])",0,25);
  funx2->SetLineColor(kRed);
  funx2->SetParameter(0,5.0); // number of photoelectrons in core
  funx2->SetParameter(1,5.0); // number of photoelectrons in clad
  funx2->FixParameter(2,cld); // clad decay constant
  tgy2->Fit(funx2,"","",1,24); // need to use TGraph here, fitting h2 produces "Warning in <Fit>: Fit data is empty"
                               // I don't know why this happening, probably some dumb mistake I'm making

  // tgy1->Draw("same"); // looks cool but distracting
  // tgy2->Draw("same"); // looks cool but distracting
  funx1->Draw("same");
  funx2->Draw("same");

  double numcore1 = funx1->GetParameter(0);
  double numclad1 = funx1->GetParameter(1);
  double numcore2 = funx2->GetParameter(0);
  double numclad2 = funx2->GetParameter(1);
  double Enumcore1 = funx1->GetParError(0);
  double Enumclad1 = funx1->GetParError(1);
  double Enumcore2 = funx2->GetParError(0);
  double Enumclad2 = funx2->GetParError(1);

  double fracore1 = (numcore1)/(numcore1+numclad1);
  double fracore2 = (numcore2)/(numcore2+numclad2);
  // ---
  double partB;
  partB = sqrt(Enumcore1**2+Enumclad2**2);
  double Efracore1 = fracore1*sqrt((Enumcore1/numcore1)**2+(partB/(numcore1+numclad1))**2);
  partB = sqrt(Enumcore2**2+Enumclad2**2);
  double Efracore2 = fracore2*sqrt((Enumcore2/numcore2)**2+(partB/(numcore2+numclad2))**2);

void test()
{
//Illustrates TVirtualFitter::GetConfidenceIntervals
//This method computes confidence intervals for the fitted function
//Author: Anna Kreshuk

   TCanvas *myc = new TCanvas("myc",
      "Confidence intervals on the fitted function",1200, 500);
   myc->Divide(3,1);

/////1. A graph
   //Create and fill a graph
   Int_t ngr = 100;
   TGraph *gr = new TGraph(ngr);
   gr->SetName("GraphNoError");
   Double_t x, y;
   Int_t i;
   for (i=0; i<ngr; i++){
      x = gRandom->Uniform(-1, 1);
      y = -1 + 2*x + gRandom->Gaus(0, 1);
      gr->SetPoint(i, x, y);
   }
   //Create the fitting function
   TF1 *fpol = new TF1("fpol", "pol1", -1, 1);
   fpol->SetLineWidth(2);
   gr->Fit(fpol, "Q");

   //Create a TGraphErrors to hold the confidence intervals
   TGraphErrors *grint = new TGraphErrors(ngr);
   grint->SetTitle("Fitted line with .95 conf. band");
   for (i=0; i<ngr; i++)
      grint->SetPoint(i, gr->GetX()[i], 0);
   //Compute the confidence intervals at the x points of the created graph
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint);
   //Now the "grint" graph contains function values as its y-coordinates
   //and confidence intervals as the errors on these coordinates
   //Draw the graph, the function and the confidence intervals
   myc->cd(1);
   grint->SetLineColor(kRed);
   grint->Draw("ap");
   gr->SetMarkerStyle(5);
   gr->SetMarkerSize(0.7);
   gr->Draw("psame");

/////2. A histogram
   myc->cd(2);
   //Create, fill and fit a histogram
   Int_t nh=5000;
   TH1D *h = new TH1D("h",
      "Fitted gaussian with .95 conf.band", 100, -3, 3);
   h->FillRandom("gaus", nh);
   TF1 *f = new TF1("fgaus", "gaus", -3, 3);
   f->SetLineWidth(2);
   h->Fit(f, "Q");
   h->Draw();

   //Create a histogram to hold the confidence intervals
   TH1D *hint = new TH1D("hint",
      "Fitted gaussian with .95 conf.band", 100, -3, 3);
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(hint);
   //Now the "hint" histogram has the fitted function values as the
   //bin contents and the confidence intervals as bin errors
   hint->SetStats(kFALSE);
   hint->SetFillColor(2);
   hint->Draw("e3 same");

/////3. A 2d graph
   //Create and fill the graph
   Int_t ngr2 = 100;
   Double_t z, rnd, e=0.3;
   TGraph2D *gr2 = new TGraph2D(ngr2);
   gr2->SetName("Graph2DNoError");
   TF2  *f2 = new TF2("f2",
      "1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+250",-6,6,-6,6);
   f2->SetParameters(1,1);
   for (i=0; i<ngr2; i++){
      f2->GetRandom2(x,y);
      // Generate a random number in [-e,e]
      rnd = 2*gRandom->Rndm()*e-e;
      z = f2->Eval(x,y)*(1+rnd);
      gr2->SetPoint(i,x,y,z);
   }
   //Create a graph with errors to store the intervals
   TGraph2DErrors *grint2 = new TGraph2DErrors(ngr2);
   for (i=0; i<ngr2; i++)
      grint2->SetPoint(i, gr2->GetX()[i], gr2->GetY()[i], 0);

   //Fit the graph
   f2->SetParameters(0.5,1.5);
   gr2->Fit(f2, "Q");
   //Compute the confidence intervals
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint2);
   //Now the "grint2" graph contains function values as z-coordinates
   //and confidence intervals as their errors
   //draw
   myc->cd(3);
   f2->SetNpx(30);
   f2->SetNpy(30);
   f2->SetFillColor(kBlue);
   f2->Draw("surf4");
//.........这里部分代码省略.........

PMTCalibration(){

   // Variable declarations
   // NOTE - Make sure the number of elements in Voltage is correct. If its not, this will
   // confuse the forloop.

   float  Voltage[19] = {1.432,1.434,1.436,1.438,1.440,1.442,1.444,1.446,1.448,1.450,1.452,1.454, 1.456, 1.458, 1.460, 1.462, 1.464, 1.466, 1.468}; // LED Voltages
   float  PedestalFitMin = -6e-12;                // Minimum for pedestal fitting range
   float  PedestalFitMax = 2e-8;                  // Maximum for pedestal fitting range
   float  SignalNorm = 100;                       // Normalization for signal fitting
   float  SignalMean = 5e-10;                     // Mean of signal distribution
   float  SignalWidth = 1e-10;                    // Width of signal distribution
   float  MaxADCforPlot = 7e-10;                     // Maximum of x-axis for linearity plot
   float  MaxNPEforPlot = 60;                   // Maximum of y-axis for linearity plot
   string IsPedestalFit;                          // String for user interface
   string IsFullADCFit;                           // String for user interface
   string IsPlotPretty;                           // String for user interface
   int    dummy;                                  // Dummy variable used in sprintf stuff
   char   filename[50];                           // Array of filenames
   char   plotname[50];                           // Array of ADC plot names
   // Ntuple for storing the derived quantities at each voltage
   TNtuple *ntuple = new TNtuple("ntuple","data for each voltage","V:mean:sigma:meanerr:sigmaerr");

   // Begin loop over all files/voltages

   for (unsigned int i=0; i<sizeof(Voltage)/sizeof(Voltage[0]); i++) {

      // To learn about sprintf:
      // http://www.cplusplus.com/reference/cstdio/sprintf/
      // To learn about "%.3f":
      // http://en.wikibooks.org/wiki/C++_Programming/Code/Standard_C_Library/Functions/printf

      dummy=sprintf(filename,"F3squarepmt_1300v_%.3fv00000.txt",Voltage[i]);
      dummy=sprintf(plotname,"F3squarepmt_1300v_%.3fv00000.pdf",Voltage[i]);
      cout<<"filename: "<<filename<<endl;

      // Put contents of text file into TGraph

      TGraph *OrignalADC = new TGraph(filename);
      OrignalADC->Draw("A*");

      // While-loop to continue iterating over the fit until the user approves it
      // To learn about cin and cout:
      // http://www.cplusplus.com/doc/tutorial/basic_io/

      IsPedestalFit = "N";
      while (IsPedestalFit != "Y") {

	 // Fit a gaussian over the pedestal only. Use fit output to shift ADC plot.
         // The 'gpad' lines are necessary so that the plot continues to display while waiting
	 // for user input (https://root.cern.ch/phpBB3/viewtopic.php?f=3&t=18852)

         PedestalFit = new TF1("h1","gaus",PedestalFitMin,PedestalFitMax);
         OrignalADC->Fit(PedestalFit,"R");
         OrignalADC->Draw("A*");
         gPad->Modified();
         gPad->Update();
         gSystem->ProcessEvents();
         cout << "Is the pedestal fit correctly?"<< endl;
         cout << "Return Y for yes and N for no."<< endl;
         cin >> IsPedestalFit;

	 // If pedestal fit is incorrect, change pedestal fit range and try again

         if (IsPedestalFit != "Y") {

            cout << "Enter the minimum fit range for the pedestal (ie -2e10)"<<endl;
            cin >> PedestalFitMin;
            cout << "Enter the minimum fit range for the pedestal (ie 2e10)"<<endl;
            cin >> PedestalFitMax;

         }

      }


      // Define a New TGraph with x-axis rescaled using the function below

      TGraph *ModifiedADC = rescaleaxis(OrignalADC,-1.0,-1.0*PedestalFit->GetParameter(1));

      // Use another while loop to assure that the pedestal plus signal fit is good

      IsFullADCFit = "N";
      while (IsFullADCFit != "Y") {

         // Fit with two gaussians. 0-2 are signal. 3-5 are pedestal. The 'FixParameter' keeps the pedestal gaussian fixed at 0.

         TF1 *FullADCFit = new TF1("FullADCFit","([0]*exp(-0.5*((x-[1])/[2])^2))+([3]*exp(-0.5*((x-[4])/[5])^2))",-1,1);
         FullADCFit->SetParameters(SignalNorm,SignalMean,SignalWidth,PedestalFit->GetParameter(0), 0.0, PedestalFit->GetParameter(2));
         FullADCFit->FixParameter(4,0);
         ModifiedADC->Fit(FullADCFit,"R");
         ModifiedADC->Draw("A*");
         gPad->Modified();
         gPad->Update();
         cout << "Is the full ADC distribution fit correctly?"<< endl;
         cout << "Return Y for yes and N for no."<< endl;
         cin >> IsFullADCFit;

         // If ADC isn't fit correctly, give the user a chance to change the inital values of the signal fit.

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