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

void LongiProfile(int E, std::string runPeriod)
{
  char inputFile[200];
  sprintf(inputFile,"%s%d%s","/lyoserv/home/ilc/steen/resultRootFile/tb_data/DHCAL_",Run(E,runPeriod),".root");
  std::cout << "inputFile = " << inputFile << std::endl;
  PionTimeCalibration* timeCalib=new PionTimeCalibration(inputFile);
  timeCalib->ShowerBarycenterCut(runPeriod);
  timeCalib->Loop();
  
  std::cout<< timeCalib->longitudinalVec.size()<<std::endl;
  TF1 *func;
  char outputFile[200];
  sprintf(outputFile,"%s%d%s","/home/steen/timeCalib/LongiProfile/calib_",Run(E,runPeriod),".txt");
  fstream out;
  out.open(outputFile,ios::out);
  for(unsigned int k=0; k<48; k++){
    TH2D* hRadPro=new TH2D("hName","",100,0,10,4000,0,200);
    func=new TF1("func","pol2",0,10);
    for(unsigned int i=0; i<timeCalib->nhit1.size(); i++){
      hRadPro->Fill(timeCalib->time.at(i),timeCalib->longitudinalVec.at(i).layer[k]);
    }
    hRadPro->ProfileX()->Fit(func,"NQ");
    out << k 
	<< " " << func->GetParameter(0) << " " << func->GetParError(0) 
	<< " " << func->GetParameter(1) << " " << func->GetParError(1)
	<< " " << func->GetParameter(2) << " " << func->GetParError(2)
	<< " " << func->GetChisquare()/func->GetNDF() << std::endl;
    delete hRadPro;
  }
  out.close();
  delete func;
  delete timeCalib;
}

/************************************************************
*  Fit ADC spectrum
************************************************************/
int FitAdcHist(TH1F *h, float data[NDATA], const int print ) {
  TF1 *fadc;
  TPaveText *adcPT;

  // Protect for null histograms
  if( h == 0 || h->GetEntries()<1000. ) return -1;
  printf("Fitting %s\n",h->GetName());

  gStyle->SetOptStat("emr"); 
  adcPT = new TPaveText(0.58,0.315,0.98,0.635,"NDC");

//  Refit ADC
  Double_t n    = h->GetEntries();
  Double_t rms  = h->GetRMS();
  Double_t mean = h->GetMean();
  Double_t ymin = mean-2*rms;
  if( ymin < 60. ) ymin=60.;
  Double_t ymax = mean+1.5*rms;
  fadc = new TF1("adc_fun",skewnormal,ymin,ymax,4);
  fadc->SetParameters(n,mean,rms,1.5);
  fadc->SetParLimits(1,mean-rms,mean+0.5*rms);
  fadc->SetParLimits(2,rms*0.5,rms*1.5);
  fadc->SetParLimits(3,1.,4.);
  h->Fit("adc_fun","r"); 
//  If skew is too high refit with fixed skew
  if( fadc->GetParameter(3) > 2.5 ) {
    printf("REFIT %s\n",h->GetName());
    fadc->SetParameter(3,1.6);
    fadc->SetParLimits(3,1.6,1.6);
    h->Fit("adc_fun","r") ;
  }
// if have enough hits use fit results  
  if( n > MINHITS ) {
    data[13] = fadc->GetParameter(1);
    data[14] = fadc->GetParError(1);
    data[15] = fadc->GetParameter(2);
    data[16] = fadc->GetParError(2);
    data[17] = fadc->GetParameter(3);
    data[18] = fadc->GetParError(3);
    double dof = fadc->GetNDF();
    if( dof > 0. ) data[19] = fadc->GetChisquare()/dof;
    adcPT->AddText(Form("Peak %.1lf #pm %.1lf",data[13],data[14]));
    adcPT->AddText(Form("Width %.1lf #pm %.1lf",data[15],data[16]));
    adcPT->AddText(Form("Skew %.2lf #pm %.2lf",	data[17],data[18]));
    if( dof > 0. ) adcPT->AddText(Form("#Chi^{2}/DoF     %.2lf",data[19]));
    adcPT->SetTextColor(2);
    adcPT->Draw();
// if have too few hits use fit stats instead of fit results  
  } else {
    data[12] = mean;
    data[13] = h->GetMeanError();
    data[14] = rms;
    data[15] = h->GetRMSError();
    data[16] = h->GetSkewness();
  }
  if( print ) c1->Print(Form("%s_%s.png",fname,h->GetName()));
  return 0;
}    //end FitAdcHist

	void FillFromFunction(TF1& infunc){
		norm = normscale * infunc.GetParameter(0);
		norm_err = normscale * infunc.GetParError(0);
		gamma = infunc.GetParameter(2);
		gamma_err = infunc.GetParError(2);
		ang = infunc.GetParameter(1);
		ang_err = infunc.GetParError(1);
		chi2 = infunc.GetChisquare();
		ndf = infunc.GetNDF();
	}

TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");

   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}

void FillGoodRun(int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
     for(int icent=0;icent<ncent;icent++){
 //     for(int ihar=0;ihar<nhar;ihar++){
 //         if(ihar!=1) continue;
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      //std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
        //fin = TFile::Open(Form("/phenix/plhf/xuq/taxi/%s%s/%d/data/%s.root",dataset.Data(),pro.Data(),taxi,GetRun(irun).Data()));
         fin = TFile::Open(Form("/gpfs/mnt/gpfs02/phenix/plhf/plhf1/xuq/phenix/flow/Run16dAu/work/62GeV/output/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_0_0_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
   //   }
     }
}

void Chi(int num, double* arr)
{
  //double chi[3] = arr;
  //double* pointer = χ
  //if(num > -1)
    {
      TF1* f = new TF1("f", "5*(TMath::Power(x,-4))",1,100);
      TH1D* h1 = new TH1D("h1","h1", 100, -0.5, 99.5);

      for(int i = 0; i < 1000000; i++)
	{
	  h1->Fill(f->GetRandom(1,100));
	  //cout << f.GetRandom(1,100) << endl;
	}
      
      //cout << __FILE__ << __LINE__<< endl;

      TF1* pwrLw = new TF1("pwrLw", "[0]*TMath::Power(x,[1])",1,1000);//changed to 0 parameter
      pwrLw->SetParameter(0,5);
      pwrLw->SetParameter(1,-4);
      //cout << pwrLw->GetParameter(0) << "\t" << pwrLw->GetParameter(1) << endl;

      //cout << __FILE__ << __LINE__<< endl;
      
      h1->Fit("pwrLw","0","",2,25);
      TF1* func = h1->GetFunction("pwrLw");
      double chi = func->GetChisquare();
      double dof = func->GetNDF();
      arr[0] = chi/dof;
      double par0 = func->GetParameter(0);
      double par1 = func->GetParameter(1);
      //cout << par0 << "\t" << par1 << endl;
      arr[1] = CalcChiSqr(h1,par0,par1,2,25);
      //chi[2] = CalcChiSqr(h1,5,-4,2,25);
      //cout << "Outputing values from Chi()" << endl;
      //cout << arr[0] << "\t" << arr[1] << endl;//"\t" << chi[2] << endl;
      //DoCleanUp(h1, f, pwrLw);
    }
    //arr = chi;
    //return chi;
}

TF1* getLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
{
   FitResults[0]         = -0.5;  //MPV
   FitResults[1]         =  0;    //MPV error
   FitResults[2]         = -0.5;  //Width
   FitResults[3]         =  0;    //Width error
   FitResults[4]         = -0.5;  //Fit Chi2/NDF

   // perform fit with standard landau
   TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
   MyLandau->SetParameter(1,300);
   InputHisto->Fit("MyLandau","0QR WW");

   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
   FitResults[0]         = MyLandau->GetParameter(1);  //MPV
   FitResults[1]         = MyLandau->GetParError(1);   //MPV error
   FitResults[2]         = MyLandau->GetParameter(2);  //Width
   FitResults[3]         = MyLandau->GetParError(2);   //Width error
   FitResults[4]         = MyLandau->GetChisquare() / MyLandau->GetNDF();  //Fit Chi2/NDF

   return MyLandau;
}

int main(int argc, char *argv[])
{
  int binCount = 100;
  if (argc > 2) {
    std::cout << "too many argument. exit." << std::endl;
    std::exit(1);
  } else if (argc == 2) {
    binCount = std::atoi(argv[1]);
    --argc;
  }
  TApplication app("app", &argc, argv);
  gStyle->SetOptStat(0);

  TH1D *hist = new TH1D("hist", "ns vs count", binCount, minNanoSec, maxNanoSec);

  char filename[256];
  GetFileName(filename, "../data/muon_lifetime/muon_lifetime.dat");
  std::ifstream ifs(filename);
  for (int i = 0; !ifs.eof(); ++i) {
    Double_t data;
    ifs >> data;
    if (data) {
      NanoSecWithError ns;
      ConvertTdcChannelToNanoSec(i%8, data, &ns);
      hist->Fill(ns.time);
    }
  }
  hist->Draw();

  TF1 *fit = new TF1("fit", "[0] * exp(- x / [1]) + [2]", 1000, 20000);
  fit->SetParameters(200, 2100, 40);
  hist->Fit(fit, "R+");
  std::cout << "chi^2/ndf: ";
  std::cout << fit->GetChisquare() / fit->GetNDF() << std::endl;

  app.Run();
  return 0;
}

void FillGoodRun(int icent,int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
         fin = TFile::Open(Form("/store/user/qixu/flow/Run16dAu/62GeV/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
}

float GoodRun(int icent, int ihar, int isub, int irun){
    float pi = acos(-1);
    TF1 *fun = new TF1("fun","pol0",-pi,pi);
    TString str;
    TFile *fin;

     ofstream fout;
        if(isub==1){
         str = "FVTX1S";
        }
        else if(isub==2){
         str = "FVTX2S";
        }
        else return -9999;
         fin = TFile::Open(Form("Run15pAu200MinBias/output_fvtxwithcntEP_%d.root",GetRun(irun)));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          hpsitemp = (TH1F*)fin->Get(Form("psi_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
      if(hpsi->GetEntries()>10000){
	hpsi->SetMarkerStyle(20);
	hpsi->SetMarkerSize(0.6);
	hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	float par=fun->GetParameter(0);
	hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
        fin->Close();
	return fun->GetChisquare()/fun->GetNDF();
      }
      else{
        fin->Close();
        return -9999;
      }
    }

void Nhit5by5Calibration(int E, std::string runPeriod)
{
  char inputFile[200];
  sprintf(inputFile,"%s%d%s","/lyoserv/home/ilc/steen/resultRootFile/tb_data/DHCAL_",Run(E,runPeriod),".root");
  std::cout << "inputFile = " << inputFile << std::endl;
  PionTimeCalibration* timeCalib=new PionTimeCalibration(inputFile);
  timeCalib->ShowerBarycenterCut(runPeriod);
  timeCalib->Loop();
  
  std::cout<< timeCalib->nhit5by5.size()<<std::endl;
  TF1 *func;
  char outputFile[200];
  sprintf(outputFile,"%s%d%s","/home/steen/timeCalib/Nhit5by5/calib_",Run(E,runPeriod),".txt");
  fstream out;
  out.open(outputFile,ios::out);
  TH2D* hRadPro=new TH2D("hName","",100,0,10,600,0,600);
  func=new TF1("func","pol1",0,10);
  for(unsigned int i=0; i<timeCalib->nhit5by5.size(); i++){
    hRadPro->Fill(timeCalib->time.at(i),timeCalib->nhit5by5.at(i));
  }
  //TCanvas *cc=new TCanvas();
  //hRadPro->Draw("colz");
  //cc->Update();
  //cc->WaitPrimitive();
  hRadPro->ProfileX()->Fit(func,"NQ");
  //cc->Update();
  //cc->WaitPrimitive();
  out << func->GetParameter(0) << " " << func->GetParError(0) 
      << " " << func->GetParameter(1) << " " << func->GetParError(1) 
      << " " << func->GetChisquare()/func->GetNDF() << std::endl;
  delete hRadPro;
  out.close();
  delete func;
  delete timeCalib;
  //delete cc;
}

void linfit()
{
   gStyle->SetOptFit();

   Double_t x2[1000], y2[1000], ex2[1000], ey2[1000];
   Int_t n2 = 10;

   TRandom3 random;
   Double_t x0, delta;
   Double_t dy = 10;
   Double_t sigma = 10;

   // parameters of line y = k*x + b
   Double_t k;
   Double_t b;
   // errors on parameters
   Double_t dk;
   Double_t db;
   // intersection with x-axis
   Double_t xint;
   Double_t dx;

   TF1 *fitfunction = 0;

   delta = 100./n2;
   x0 = 0 - delta;
   for (int i=0; i<n2; ++i) {
      x0 += delta;
      x2[i] = x0;
      y2[i] = random.Gaus(x2[i], sigma);
      ex2[i] = 0;
      ey2[i] = dy;         // all weights are equal
      ey2[i] = TMath::Sqrt(y2[i]);
   }

   TGraphErrors *gr2 = new TGraphErrors(n2, x2,y2, ex2, ey2);
   gr2->SetNameTitle("gr2", "gr2");
   gr2->SetMarkerStyle(20);
   gr2->SetMarkerColor(4);
   gr2->SetLineColor(4);

   new TCanvas;
   gr2->Draw("ap");
   gr2->Fit("pol1");

   fitfunction = gr2->GetFunction("pol1");
   b = fitfunction->GetParameter(0);
   db = fitfunction->GetParError(0);
   k = fitfunction->GetParameter(1);
   dk = fitfunction->GetParError(1);
   xint = -b / k;
   dx = TMath::Abs(b/k) * TMath::Sqrt( (db/b)*(db/b) + (dk/k)*(dk/k) );
   gr2->SetTitle(Form("intersection with x-axis: %0.3f #pm %0.3f", xint,dx));
   cout<< "--> intersection with x-axis: " << xint << " +- " << dx <<endl;

   cout<< "--> chi2ndf = " << fitfunction->GetChisquare() / fitfunction->GetNDF();
   cout<<endl<<endl;

   // fast straight line fit y = am + bm*x
   // Source: Kalashnikova, Kozodaev, Detektory Elementarnyx chastic.

   //-- Float_t version of the data
   Float_t xf[1000], yf[1000], eyf[1000];
   Int_t nf = n2;
   for (int i=0; i<nf; ++i) {
      xf[i] = x2[i];
      yf[i] = y2[i];
      eyf[i] = ey2[i];
   }

   Double_t wx = 0;
   Double_t wy = 0;
   Double_t wxy = 0;
   Double_t wx2 = 0;
   Double_t W = 0;

   wx = wy = wxy = wx2 = W = 0;
   for (int i=0; i<gr2->GetN(); ++i)
   {
      Double_t x = gr2->GetX()[i];
      Double_t y = gr2->GetY()[i];
      Double_t yerror = gr2->GetEY()[i];
      Double_t w = 1./(yerror*yerror);       // weights are 1/ey

      W += w;
      wx += w*x;
      wx2 += w*x*x;
      wy += w*y;
      wxy += w*x*y;
   }

   Double_t Discriminant = W*wx2 - wx*wx;
   Double_t am = (wy*wx2 - wxy*wx) / Discriminant;    // y = am + bm*x
   Double_t bm = (W*wxy - wx*wy) / Discriminant;

   Double_t chi2 = 0;
   Double_t chi2ndf = 0;
   for (int i=0; i<gr2->GetN(); ++i)
   {
      Double_t x = gr2->GetX()[i];
//.........这里部分代码省略.........

void photonContainmentAnalysis() {

  gStyle->SetOptFit(1111);
  
  TFile *infile = new TFile("contCorrAnalyzer.root");
  infile->ls();
  
  TH1F *theHistos[2];
  theHistos[0] = (TH1F*)infile->Get("contCorrAnalyzer/EB_e25EtrueReference");
  theHistos[1] = (TH1F*)infile->Get("contCorrAnalyzer/EE_e25EtrueReference");
  
  // fitting the distributions to extract the parameter
  TF1 *gausa;
  TF1 *cb_p;
  for(int myH=0; myH<2; myH++) {
  
    // histos parameters
    int peakBin   = theHistos[myH]->GetMaximumBin();
    double h_norm = theHistos[myH]->GetMaximum();
    double h_rms  = theHistos[myH]->GetRMS();
    double h_mean = theHistos[myH]->GetMean();
    double h_peak = theHistos[myH]->GetBinCenter(peakBin);
    
    // gaussian fit to initialize
    gausa = new TF1 ("gausa","[0]*exp(-1*(x-[1])*(x-[1])/2/[2]/[2])",h_peak-10*h_rms,h_peak+10*h_rms);
    gausa ->SetParameters(h_norm,h_peak,h_rms);
    theHistos[myH]->Fit("gausa","","",h_peak-3*h_rms,h_peak+3*h_rms);
    double gausNorm  = gausa->GetParameter(0);
    double gausMean  = gausa->GetParameter(1);
    double gausSigma = fabs(gausa->GetParameter(2));
    double gausChi2  = gausa->GetChisquare()/gausa->GetNDF();
    if (gausChi2>100){ gausMean = h_peak; gausSigma = h_rms; }
    
    // crystalball limits
    double myXmin = gausMean - 7.*gausSigma;
    double myXmax = gausMean + 5.*gausSigma;
    
    // crystalball fit
    cb_p = new TF1 ("cb_p",crystalball,myXmin,myXmax, 5) ;
    cb_p->SetParNames ("Mean","Sigma","alpha","n","Norm","Constant");
    cb_p->SetParameter(0, gausMean);
    cb_p->SetParameter(1, gausSigma);
    cb_p->SetParameter(2, 1.);
    cb_p->SetParLimits(2, 0.1, 5.);
    cb_p->FixParameter(3, 5.);
    cb_p->SetParameter(4, gausNorm);
    theHistos[myH]->Fit("cb_p","lR","",myXmin,myXmax);
    theHistos[myH]->GetXaxis()->SetRangeUser(0.95,1.05); 
    double matrix_gmean      = cb_p->GetParameter(0);
    double matrix_gmean_err  = cb_p->GetParError(0);
    
    c1->Update();
    if(myH == 0) { 
      c1->SaveAs("e25EtrueReference_EB.eps");
      std::cout << "E25 barrel containment: " << matrix_gmean << " +/- " << matrix_gmean_err << std::endl; 
    }
    if(myH == 1) { 
      c1->SaveAs("e25EtrueReference_EE.eps");
      std::cout << "E25 endcap containment: " << matrix_gmean << " +/- " << matrix_gmean_err << std::endl; 
    }
  }
}

/** Make Va1 response
    @param n 
    @param B 
    @param dc 
    @param errors 
    @ingroup simple_script
*/
void 
VA1Response(Int_t n=4, Float_t B=6, Float_t dc=.01, Bool_t errors=kFALSE,
	    Bool_t doFit=kFALSE) 
{
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  gStyle->SetLabelFont(132, "xyz");
  gStyle->SetTitleFont(132, "xyz");
  gStyle->SetTitleSize(0.08, "y");
  gStyle->SetTitleOffset(0.5, "y");
  gStyle->SetTitleSize(0.06, "x");
  gStyle->SetTitleOffset(0.7, "x");

  TCanvas* c = new TCanvas("c", "c", 800, 500);
  c->SetFillColor(0);
  c->SetBorderMode(0);
  c->SetBorderSize(0);
  c->SetTopMargin(0.05);
  c->SetRightMargin(0.05);
  c->SetGridx();
  c->SetGridy();

  TF1* response = new TF1("response", "[0] * (1 - exp(-[1] * x))", 0, 1.4);
  response->SetParameters(1, B);
  response->SetParNames("A", "B");
  response->SetLineColor(2);
  
  TGraph* graph = 0;
  if (n >= 2) {
    if (errors) graph = new TGraphErrors(n);
    else        graph = new TGraph(n);
    for (Int_t i = 0; i < n; i++) {
      Float_t t = Float_t(i + 1) / n;
      Float_t q = gRandom->Gaus(response->Eval(t), dc);
      graph->SetPoint(i, t, q);
      if (errors) ((TGraphErrors*)graph)->SetPointError(i, 0, dc);
    }
  }

  response->Draw();
  response->GetHistogram()->GetYaxis()->SetRangeUser(0, 1.05);
  response->GetHistogram()->GetXaxis()->SetRangeUser(0, 1.1);
  response->GetHistogram()->GetXaxis()->SetNdivisions(6, kTRUE);
  response->GetHistogram()->GetYaxis()->SetNdivisions(10, kTRUE);
  response->GetHistogram()->SetXTitle("t");
  response->GetHistogram()->SetYTitle(Form("1-e^{-%3.1f t}", B));
  
  if (graph) {
    
    
    graph->Draw("P*");
    TString fitOpt("E");
    if (!errors) fitOpt.Append("W");

    if (doFit) {
      TF1* fit = new TF1("fit",  "[0] * (1 - exp(-[1] * x))",  0, 1);
      fit->SetParameters(.5, B/2);
      fit->SetParNames("A", "B");
      fit->SetLineColor(3);
      graph->Fit("fit", fitOpt.Data());
      graph->Fit("fit", fitOpt.Data());
      
      std::cout << "Chi^2/NDF = " << fit->GetChisquare() << "/" << fit->GetNDF()
		<< " = " << std::flush;
      if (fit->GetNDF() == 0) 
	std::cout << " undefined!" << std::endl;
      else
	std::cout << (fit->GetChisquare() / fit->GetNDF()) << std::endl;
      std::cout << "f(t) = " 
		<< fit->GetParameter(0) << "+/-" << fit->GetParError(0) 
		<< " * (1 - exp(" 
		<< fit->GetParameter(1) << "+/-" << fit->GetParError(1) 
		<< " * t))" << std::endl;
    }
  }

  c->Modified();
  c->Update();
  c->cd();
  c->SaveAs("va1_response.png");
}

int main(int argc, char **argv)
{
	TFile *f;
	TFile *newf;
	TNtuple *ntuple;
	TNtuple *newntuple;
	TH1F *h1;
	TF1 *func;
	TString Metal;
	TString dataLoc;
	Int_t dataSL;
	
	if(argc < 5 || argc > 6){
		std::cout << "The number of arguments is incorrect" << std::endl;
		return 0;
	}
	
	dataLoc = argv[1];
	PHI_MIN = (Double_t) std::stod(argv[2]);
	PHI_MAX = (Double_t) std::stod(argv[3]);
	N_PHI = (Int_t) std::stoi(argv[4]);
	if(argc >= 6) Metal = argv[5];
		
	dataSL = dataLoc.Length();
	if(dataLoc(dataSL-1, 1) != "/"){
		dataLoc = dataLoc + "/";
	}	
	
	std::cout << "The following settings are being used" << std::endl;
	std::cout << "Data Directory = " << dataLoc << std::endl;
	std::cout << PHI_MIN << " < PHI < " << PHI_MAX << ", N_PHI = " << N_PHI << std::endl;
	if(argc == 6) std::cout << "Running only for Metal "<< Metal << std::endl;
	else std::cout << "Running all Metals" << std::endl;	
	
	Float_t Q2, Xb, Zh, Pt, Phi, Val, Err;
	Float_t A, Ac, Acc;
	Float_t AErr, AcErr, AccErr;
	Float_t ChiSQ;

	Int_t nentries, empty;

	if(Metal == "") N_METAL = 6;
	else N_METAL = 1;
	
	for(Int_t met = 0; met < N_METAL; met++){
		if(Metal == ""){
			if(met == 0) Metal = "C";
			else if(met == 1) Metal = "Fe";
			else if(met == 2) Metal = "Pb";
			else if(met == 3) Metal = "D_C";
			else if(met == 4) Metal = "D_Fe";
			else if(met == 5) Metal = "D_Pb";
		}
		if(dataLoc == "")
			f = new TFile("../Gen5DimData/" + Metal + "_5_dim_dist.root", "READ");
		else
			f = new TFile(dataLoc + Metal + "_5_dim_dist.root", "READ");
		ntuple = (TNtuple*) f->Get("fit_data");

		nentries = ntuple->GetEntries();
		ntuple->SetBranchAddress("Xb", &Xb);
		ntuple->SetBranchAddress("Q2", &Q2);
		ntuple->SetBranchAddress("Xb", &Xb);
		ntuple->SetBranchAddress("Zh", &Zh);
		ntuple->SetBranchAddress("Pt", &Pt);
		ntuple->SetBranchAddress("Phi", &Phi);
		ntuple->SetBranchAddress("Val", &Val);
		ntuple->SetBranchAddress("Err", &Err);

		newntuple = new TNtuple("AAcAcc_data", "AAcAcc_data", "Q2:Xb:Zh:Pt:A:AErr:Ac:AcErr:Acc:AccErr:ChiSQ");
		func = new TF1("fit", "[0]+TMath::Cos(x*TMath::Pi()/180)*[1]+TMath::Cos(2*x*TMath::Pi()/180)*[2]");

		newf = new TFile(Metal + "newphihist.root", "RECREATE");
		newf->cd();

		for(Int_t i = 0; i < nentries; i = i + N_PHI){
			ntuple->GetEntry(i);
			h1 = new TH1F((const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt), 
							(const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt), N_PHI, PHI_MIN, PHI_MAX);
			empty = 0;
			for(Int_t j = 1; j <= N_PHI; j++){
				ntuple->GetEntry(i+j-1);
				h1->SetBinContent(j, Val);
				if(h1->GetBinContent(j) == 0)
					empty++;
				h1->SetBinError(j, Err*1.04);
			}
			if(empty <= (N_PHI - MIN_BIN)){
				h1->Fit(func, "q");
				h1->Write();
				if(func->GetNDF() != 0){
					ChiSQ = func->GetChisquare();
					A = func->GetParameter(0);
					AErr = func->GetParError(0);
					Ac = func->GetParameter(1);
					AcErr = func->GetParError(1);
					Acc = func->GetParameter(2);
					AccErr = func->GetParError(2);
					newntuple->Fill(Q2, Xb, Zh, Pt, A, AErr, Ac, AcErr, Acc, AccErr, ChiSQ);
				}
//.........这里部分代码省略.........