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

/*============================================================================*/
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");
}

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

// Loop through detector layers to report 
// various parameters
residu(TString baseName, double min=-1.5, double max=1.5)
{
  double xMin, xMax;
  TString name;
  name = baseName;
  name += "_a";     TH1D * amp  = new TH1D(name,name,55,0,54.);
  name = baseName;
  name += "_mean";  TH1D * mean = new TH1D(name,name,55,0,54.);
  name = baseName;
  name += "_rms";   TH1D * rms = new TH1D(name,name,55,0,54.);
  TString thisL;
  TH1D *h=0;
  TF1 * f;
  int ii;
  for(int i=0; i<51;i++)
    {
      thisL=baseName+i;
      if (i<6)
	{
	  xMin = -0.2;
	  xMax = 0.2;
	}
      else
	{
	  xMin = min;
	  xMax = max;
	}
      h=(TH1D*)gDirectory->Get(thisL);
      h->Fit("gaus","","",xMin,xMax);
      double entries = h->GetEntries();
      if (entries>0)
	{
	  ii = i+1;
	  f = h->GetFunction("gaus");
	  double a = f->GetParameter("Constant"); double ea = f->GetParError(0);
	  double m = f->GetParameter("Mean");double em = f->GetParError(1);
	  double r = f->GetParameter("Sigma");double er = f->GetParError(2);
	  amp->SetBinContent(ii,a);amp->SetBinError(ii,ea);
	  mean->SetBinContent(ii,m);mean->SetBinError(ii,em);
	  rms->SetBinContent(ii,r);rms->SetBinError(ii,er);
	}
    }
  cout << "Fitting completed" << endl;

  c1->Clear();
  c1->Divide(1,3);
  c1->cd(1); amp->Draw();
  c1->cd(2); mean->Draw();
  c1->cd(3); rms->Draw();
  
}

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

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

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

void ThrowPSEME(int me_i, bool use_old=false){

  // ====================================================
  //                       Create a fit
  // ====================================================
  TF1 *ftot = new TF1("fitFunction",fitFunction,-100,0,2);

  //threw one data PSE  and fill it
  if (use_old){
    cout<<"OLD"<<endl;
    siggnal = MEWH_old[me_i];
    background = MEWBB_old[me_i];
  }else{
    cout<<"NEW"<<endl;
    siggnal = MEWH[me_i];
    background = MEWBB[me_i];
  }

  TH1F * Err = new TH1F("Err","Err",1000,0,1);

  for (int pse=0;pse<100;pse++){

    //Get one PSE
    TH1F * data  =  MakeOnePSE(siggnal,1000,background,2000);
    
    //fit it
    data->Fit("fitFunction","b");   
    Err->Fill(ftot->GetParError(1));
    delete data;

  }
  Err->Draw();

}//ThrowPSEME

TF1* fit_xy_res(double &x, double &xe, double &max, bool do_fit, int idel, int ipad, TCanvas *tc,
		TH1F* h, const char *meth, const char *var_tag, const char *det_tag, const char *var, const char *det_desc) {

  tc->cd(1+ipad);

  //cout << "var_tag = " << var_tag << endl;
  TLatex *tt = new TLatex();
  tt->SetTextSize(0.08);
  tt->SetTextColor(4);
  tt->SetNDC(true);

  double fmx = 2.2;
  double fit_min= h->GetMean()-fmx*h->GetRMS();
  double fit_max= h->GetMean()+fmx*h->GetRMS();
  TF1* f = new TF1(Form("fxres_%s_%d_%s",det_tag,idel,meth), "gaus",fit_min,fit_max);

  if (do_fit) {

    h->Fit(f,"R+Q");
    x = f->GetParameter(2);
    xe = f->GetParError(2);

    if (f->GetParameter(1) < f->GetParError(1))
      cout << "   OK: ";
    else
      cout << " OOPS: ";
    cout << var_tag << " " << meth << " " << det_tag << idel;
    cout << " mean= " << f->GetParameter(1) << " pm " << f->GetParError(1);
    cout << " sigma= " << f->GetParameter(2) << " pm " << f->GetParError(2) << endl;

  } else {

    max = TMath::Max(max, x);
    h->Draw();
    //tt->DrawLatex(0.12,0.15, h->GetTitle());
  }

  tt->DrawLatex(0.65,0.8, meth);
  tt->DrawLatex(0.65,0.7, var);
  tt->DrawLatex(0.65,0.6, det_desc);
  tt->DrawLatex(0.65,0.5, Form("#delta = %d%%",idel-5 ));

  max = TMath::Max(max, x);
  gPad->SetGridx();
  return f;
}

void histspec(TH1* hist, Double_t &mean, Double_t &ermean, Double_t &rms, 
		Double_t range=4, Double_t rangef=2.5, int mode=0) {

  Double_t xmin,xmax,rmsmin,meanf,ermeanf,rmsf;
  char ctit[145];

  mean=hist->GetMean();
  rms=hist->GetRMS();
  xmin=hist->GetXaxis()->GetXmin();
  xmax=hist->GetXaxis()->GetXmax();
  int nb=hist->GetNbinsX();
  rmsmin=(xmax-xmin)/(nb+1);
  if (rms==0) rms=rmsmin;

  for (int i=0;i<2;i++) {
    hist->SetAxisRange(mean-range*rms,mean+range*rms);
    mean=hist->GetMean();
    rms=hist->GetRMS();
    if (rms==0) rms=rmsmin;
  }
  
  if (rangef>0) {
    for (int i=0;i<2;i++) {
      hist->SetAxisRange(mean-rangef*rms,mean+rangef*rms);
      mean=hist->GetMean();
      rms=hist->GetRMS();
      if (rms==0) rms=rmsmin;
    }
  }
  
  ermean=hist->GetMeanError();
  if (abs(mode)>0) {
    if (mode<0) {
      //sprintf(ctit,"cHistSpecs%4d",nHistSpecs++);
      //TCanvas *ccxx = new TCanvas("ccxx","ccxx",900,0,450,450);
      //hist->SetMinimum(0);
      //hist->Draw();
      hist->Fit("gaus","LI","same",mean-rangef*rms,mean+rangef*rms);
      //ccxx->Update();
      for (int i=0;i<1e3;i++) for (int j=0;j<1e3;j++) {} 
    }
    else hist->Fit("gaus","LI0","",mean-rangef*rms,mean+rangef*rms);
    TF1 *f = hist->GetFunction("gaus");
    meanf=f->GetParameter(1);
    ermeanf=f->GetParError(1);
    rmsf=f->GetParameter(2);
    if (fabs(mean-meanf)/(fabs(mean+meanf)+0.0001)<0.8 && 
	fabs(ermean-ermeanf)/(fabs(ermean+ermeanf)+0.0001)<0.8) {
      mean=meanf;
      ermean=ermeanf;
      rms=rmsf;
    }
  }
  hist->SetAxisRange(xmin,xmax);
  return;
}

double Fit511Photopeak(TH1* h, double* error=0)
{
  TSpectrum spec(1);
  spec.Search(h);
  h->GetXaxis()->SetTitle("Energy [photoelectrons]");
  h->Draw("e");
  TH1* bg = spec.Background(h);
  TH1* sig = (TH1*)(h->Clone());
  sig->SetLineColor(kGreen);
  
  sig->Add(bg,-1);
  sig->Draw("same e");
  sig->Fit("gaus","m","same e");
  TF1* gaus = sig->GetFunction("gaus");
  if(gaus)
    gaus->SetLineColor(kGreen);
  
  bg->SetLineColor(kRed);
  bg->Draw("same e");
  
  TLine* line = new TLine(gaus->GetParameter(1),0,
			  gaus->GetParameter(1),h->GetMaximum());
  line->SetLineColor(kBlue);
  line->SetLineWidth(2);
  line->Draw();
  

  double yield = spec.GetPositionX()[0]/epeak;
  double err = 0;
  
  cout<<"Results from TSpectrum: \n\t"
      <<"Peak = "<<spec.GetPositionX()[0]<<" p.e.; Light Yield = "
      <<yield<<" p.e./keV"<<endl;
  if(gaus){
    yield = gaus->GetParameter(1)/epeak;
    err = gaus->GetParError(1)/epeak;
    cout<<"Results from BG Subtracted Gaus Fit: \n\t"
	<<"Peak = "<<gaus->GetParameter(1)<<" p.e.; Light Yield = "
	<<yield<<" +/- "<<err<<" p.e./keV"<<endl;
    err = max(err, TMath::Abs(yield-spec.GetPositionX()[0]/epeak));
    
  }
  TLegend* leg = new TLegend(.6,.6,.9,.9);
  leg->AddEntry(h,"Raw Spectrum","lpe");
  leg->AddEntry(bg,"Background","lpe");
  leg->AddEntry(sig,"Signal","lpe");
  char title[20];
  sprintf(title,"Yield = %.2f pe/keV",yield);
  leg->AddEntry(line, title ,"l");
  leg->Draw();

  if(error) *error = err;
  return yield;
}

drawGraph( TGraph* gr, const Char_t* ytitle = 0 ) 
{
	gr->Draw( "AP*" );
	TH1* h = gr->GetHistogram( );

	std::cout << gr->GetMean( ) << std::endl;

	gr->Fit( "pol0", "Q", "SAME" );
	TF1* fit = gr->GetFunction( "pol0" );	
	Double_t chi2 	= fit->GetChisquare( );
	Double_t par0 	= fit->GetParameter( 0 );
	Double_t err 	= fit->GetParError( 0 );	

	TLatex* l = new TLatex;
	l->SetNDC( );
	l->DrawLatex( 0.85, 0.80, Form( "#splitline{Entries = %i}{RMS = %2.3f}",
		gr->GetN(), gr->GetMean(2)) );

	TLatex* l2 = new TLatex;
	l2->SetNDC( );
	l2->DrawLatex( 0.85, 0.7, Form("#splitline{fit = %2.3f#pm%2.3f}{#chi^{2} = %2.3f}", 
		par0, chi2, err) );
		
//	TLegend* leg = new TLegend( 0.8, 0.9, 0.95, 0.95 );
//	leg->SetFillColor( kWhite );
//	leg->SetBorderSize( 0 );
//	leg->AddEntry( "gr", "CPU+AFAR+MBED+ATWD+AMPS", "P" );
//	leg->Draw( );

	if( h != 0 ) 
	{
		h->GetXaxis()->SetTimeFormat("%m/%d %H:%M");
		h->GetXaxis()->SetNdivisions( 10, 10, 0 );
		h->GetXaxis()->SetTimeDisplay( 1 );
		h->GetXaxis()->SetTitle("Date / Time (UTC)");
		h->GetXaxis()->SetLabelSize( 0.06 );
		h->GetYaxis()->SetLabelSize(0.07);
		h->GetYaxis()->SetTitleSize(0.06);
		h->GetXaxis()->SetTitleSize(0.055);
		h->GetYaxis()->SetTitleOffset( 0.9 );
		h->GetXaxis()->SetTitleOffset( 0.9 );
		h->SetMarkerSize( 5 );

		if( ytitle != 0 ) 
		{
			h->GetYaxis()->SetTitle(ytitle);
      	}
	}
   
	if( gPad != 0 ) 
	{
		gPad->Update();
	}
}

scalePair energyCorrectionDiff(TCut centCut1, float lowPt, float highPt, float lowEta, float highEta,TCut addCut) { 
   
   TString fname1 = "forest/barrelHiForestPhoton_MCphoton50_51k.root";
   if ( lowPt > 90 ) 
      fname1 = "forest/barrelHiForestPhoton_MCphoton80_28k.root";
   TFile *f1  =new TFile(fname1.Data());
   TTree *photon1 = (TTree*)f1->Get("yongsunPhotonTree");
   photon1->AddFriend("yEvt=yongsunHiEvt"    ,fname1.Data());
   photon1->AddFriend("ySkim=yongsunSkimTree"   ,fname1.Data());
   photon1->AddFriend("yHlt=yongsunHltTree"     ,fname1.Data());
   
   TCut collisionCut = "ySkim.pcollisionEventSelection==1";
   
   TCut hoeCut = "hadronicOverEm<0.2";
   TCut isoCut = "cc4 + cr4 + ct4j20 < 5 && sigmaIetaIeta<0.011";
   TCut ptCut  = Form("genMatchedPt>%.f && genMatchedPt <%.f",lowPt, highPt);
   TCut etaCut = Form("abs(eta)>%f && abs(eta)<%f",lowEta,highEta);
      
   TCut finalCut1 = genMatchCut1 && collisionCut && centCut1 && hoeCut && isoCut && ptCut && etaCut && addCut ; 
   TString variable1 = "pt/genMatchedPt";
   
   TH1D* hScale = new TH1D("hScale","",100,.5,1.5);
   TH1D* hdpt = new TH1D("hdpt","",100,-20,20);
   photon1->Draw(Form("%s>>%s",variable1.Data(),hScale->GetName()), finalCut1);
   photon1->Draw(Form("pt-genMatchedPt>>%s",hdpt->GetName()),       finalCut1);
   cout << "cut = " << finalCut1.GetTitle() <<endl;
   hScale->Draw();
   hdpt->Draw();
   TF1* ff =  cleverGaus(hScale);
   
   scalePair ret;
   double *ps = ff->GetParameters();
   ret.val = ps[1];
   ret.err = ff->GetParError(1);
   // resErr = ff->GetParError(2);
   ret.res  = ps[2];
   ret.resErr = ff->GetParError(2);  
   ret.absVal = hdpt->GetMean();
   cout <<"scale = " << ret.val << " +-" << ret.err << endl;
   return ret;
}

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

double findDelPhiAngle(TH1D *input, int parorErr){
  
  TF1 *fitter = new TF1("fitter","gaus",-2,2);
  //input->Scale(1./(input->GetEntries()));
  fitter->SetParameters(1./input->GetMaximum(),-0.1,1.0);
  input->Fit("fitter","qN0","",-1.5,1.5);
  if(parorErr==1){
    return fitter->GetParameter(1);
  }
  else if(parorErr==2){
    return fitter->GetParError(1);
  }
  else return 0;
}

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