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

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

TGraphErrors* ExtractCalCurve(vector<calData>& FitPositions, vector<Float_t>& Energies)
//void ExtractCalCurve(vector<Float_t>& FitPositions, vector<Float_t>& Energies)
{
    //TF1* calFit = new TF1("calFit","[0]*x + [1]",0,800E3);
    TGraphErrors* calCurve = new TGraphErrors(Energies.size());
    
    //Get highest energy peak and estimate [0]
    std::vector<calData>::iterator pit;
    std::vector<Float_t>::iterator eit;
    
    pit = max_element(FitPositions.begin(),FitPositions.end(),CompareByadc);
    eit = max_element(Energies.begin(),Energies.end());
    
    //Extrapolate coefficient and place first point
    int currPoint = 1;
    calCurve->SetPoint(currPoint,(*pit).fadc,*eit);
    calCurve->SetPointError(currPoint,(*pit).efadc,0.0);
    TFitResultPtr r = calCurve->Fit("pol1","SQ");   
    Float_t a = r->Parameter(1);

    Float_t CurrentPeak = 0.0;
    Float_t CurrentEnergy = 0.0;
    Float_t CurrentEnergyEst = 0.0;
    // Loop through found peaks and locate closest estimated energy
    // Assume fitted peaks are already ordered from lowest to highest
    for(std::vector<calData>::iterator i = --(FitPositions.end()); i!=FitPositions.begin(); --i)
    {
        currPoint++;
        CurrentPeak = (*i).fadc;
        CurrentEnergyEst = CurrentPeak*a;
        Float_t CurrentDelta = 800E3;
        for(std::vector<Float_t>::iterator j = Energies.begin(); j!=Energies.end(); j++)
	{
		if( abs(*j - CurrentEnergyEst) < CurrentDelta)
		{
			CurrentDelta = abs(*j - CurrentEnergyEst);
                        CurrentEnergy = *j;
		}
	} 
	
	calCurve->SetPoint(currPoint,CurrentPeak,CurrentEnergy);
        calCurve->SetPointError(currPoint,(*i).efadc,CurrentDelta);
        r = calCurve->Fit("pol1","SQ");
        a = r->Parameter(1);
    }
    r->Print("V");
    return calCurve; 
}

void draw() {

  gStyle->SetOptStat(0000);

  TGraphErrors* gr = new TGraphErrors("papel_um.txt", "%lg %lg %lg %lg");
  gr->Draw("A");

  gr->GetXaxis()->SetTitle("Diametro (mm)");
  gr->GetYaxis()->SetTitle("Massa (u.m.)");
  gr->SetTitle("Dimensao Fractal - Papel");
  gr->SetMarkerStyle(20);
  gr->SetMarkerSize(gr->GetMarkerSize()/2.);

  TCanvas* c1 = new TCanvas();
  gr->Draw("AP");
  c1->Print("gr1.png");


  TGraphErrors* grlog = new TGraphErrors();
  for (int i=0; i<gr->GetN(); i++) {
    grlog->SetPoint(i, TMath::Log(gr->GetX()[i]), TMath::Log(gr->GetY()[i]));
    grlog->SetPointError(i, gr->GetEX()[i]/gr->GetX()[i], 0.0);
  }
  grlog->Draw("A");
  grlog->GetXaxis()->SetTitle("Logaritmo do diametro (diametro em mm)");
  grlog->GetYaxis()->SetTitle("Logaritmo da massa (massa em unidade arbitraria)");
  grlog->SetTitle("Dimensao Fractal - Papel");
  grlog->SetMarkerStyle(20);
  //  grlog->SetMarkerSize(gr->GetMarkerSize()/2.);
  grlog->Fit("pol1");
  TCanvas* c2 = new TCanvas();
  grlog->Draw("AP");
  c2->Print("gr2.png");

}

void LSF() {
  float x[100], y[100], ex[100], ey[100];
  char filename[80];
  printf("Enter filename.\n");
  scanf("%s",filename);

  FILE *f = fopen(filename,"r");
  if (f == NULL) {
    sprintf(line,"File %s not found.\n",filename);
    printf(line);
    return;
  }

  // Read data
  int i = 0;
  while (fscanf(f,"%f %f %f",&(x[i]),&(y[i]),&(ey[i])) == 3 ) {
    i++;
  }

  printf("Read %d data points.\n",i);

  TGraphErrors *tge = new TGraphErrors(i,x,y,ex,ey); 
  tge->SetMarkerStyle(20);
  tge->SetMarkerColor(kBlue);
  tge->SetLineColor(kBlue);

  c1 = new TCanvas("c1","c1",600,600);

  tge->Draw("AP");
  tge->Fit("pol1");
}

instrAsy() {
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  
  TGraphErrors *gr =new  TGraphErrors;
  gr->SetName("insyAsy");
  gr->SetTitle("Instrumental Asymmetries");
  //gr->SetMarkerColor(ampCol[iam]);
  gr->SetMarkerSize(0.8);
  gr->SetMarkerStyle(21);
  TGraphErrors *gr1=gr;
  

  TString  wrkDir="final/";
  TString  fname=wrkDir+"defaultB-H/asyVer1.hist.root";
  TFile *inpH=new TFile(fname);
  assert(inpH->IsOpen());
  //inpH->ls();
  const int  nSel=5;
  char *selL[nSel]={ "a0*All","b0*All","a2*All","b2*All","c1*All"};
  int isel;
  for (isel=0;isel<nSel;isel++) {
    TGraphErrors *gr =(TGraphErrors *) inpH->Get(selL[isel]);
    assert(gr);
    //gr->Print();
    //gr->Draw("AP");
    float x1=0;
    if(isel<2) x1=2189;
    gr->Fit("pol0","","",x1,3000);
    TF1 *ff=gr->GetFunction("pol0");
    if(ff==0) continue; // no fit was made
    float val=ff->GetParameter(0);
    float err=ff->GetParError(0);
    printf("pol0= %f +/- %f\n",val,err);
        int n=gr1->GetN();
      float x=isel+1;
      gr1->SetPoint(n,x,val);
      gr1->SetPointError(n,0,err); 
    }


  gr1->Print();

  // return;
  gr1->Draw("AP");

  // save Tgraph
  TString fname=wrkDir+"instAsy.hist.root";
  TFile *outH=new TFile(fname,"RECREATE");
  assert(outH->IsOpen());
  printf("save outH -->%s\n", fname.Data());
  gr1->Write();
  outH->ls();
  outH->Write();


    return;
}

void YZFit(){   
    if (yerr[0] == 0 && yerr[1] == 0 && yerr[2] == 0 && zerr[0] == 0 && zerr[1] == 0 && zerr[2] == 0) {cout << "ERRORE: Incertezze non specificate" << endl;}
    else {
      delete yzfitfunc;
      delete yzgraph;
    yzgraph = new TGraphErrors(3,yv,zv,yerr,zerr);
    yzfitfunc = new TF1("yzfittingfunction","[0]+[1]*x",-1000,1000);
    yzfitfunc->SetParameters(300,-6);
    yzgraph->Fit(yzfitfunc,"0S");
    }
  };

int main(int narg,char **arg)
{
  int i;
  double x,y;
  ifstream file;
  TApplication myapp("App",NULL,NULL);
  TCanvas tela;
  TGraphErrors graf;
  TF1 fun("fun","[0]*exp(-x*x/(2*[1]*[1]))");

  if(narg<2)
    {
      cerr<<"usa: "<<arg[0]<<" file"<<endl;
      exit(0);
    }

  file.open(arg[1]);

  if(file.good()==0)
    {
      cout<<"File: "<<arg[1]<<" inesistente!"<<endl;
      exit(0);
    }

  i=0;
  double dum,N=1;
  while(file>>x>>dum>>y)
    {
      if(i==0) N=y;
      graf.SetPoint(i,x,y);
      graf.SetPointError(i,0,0.0001);
      i++;
    }
  
  fun.SetParameter(0,N);
  fun.SetParameter(1,1);

  fun.SetParName(0,"N");
  fun.SetParName(1,"S");

  graf.Draw("AP");

  graf.Fit("fun");

  //graf.GetXaxis()->SetLimits(0,x);

  tela.Modified();
  tela.Update();

  myapp.Run(true);

  return 0;
}

void gerrors() {
   
   TCanvas *mycanvas = new TCanvas("mycanvs","A Simple Graph with error bars",200,10,700,500);

//   mycanvas->SetFillColor(42);
   mycanvas->SetGrid();
//   mycanvas->GetFrame()->SetFillColor(21);
//  mycanvas->GetFrame()->SetBorderSize(12);

 const int  n_points =10;
 #double x_vals [ n_points]= {1,2,3,4,5,6,7,8,9,10};
 double  y_vals [ n_points]= {6 ,12 ,14 ,20 ,22 ,24 ,35 ,45 ,44 ,53};
 double y_errs [ n_points]= {5 ,5 ,4.7 ,4.5 ,4.2 ,5.1,2.9,4.1,4.8,5.43};

   TGraphErrors *gr = new TGraphErrors(n_points,x_vals,y_vals,NULL,y_errs);
//   TGraphErrors graph(n_points,x_vals,y_vals,NULL,y_errs);
  gROOT -> SetStyle("Plain");
   gr->SetTitle("TGraphErrors Example; lengtht  [cm];Arb. Units");
   gr->SetMarkerColor(kBlue);
   gr->SetMarkerStyle(kOpenCircle);
   gr->SetLineColor ( kBlue ) ;
   gr->Draw("ALP");
 
 //Define a linear function
TF1 f("Linear law" ,"[0]+x*[1]" ,.5 ,10.5) ; 
 
// Let's make the funcion line nicer 
f.SetLineColor(kRed);
 f.SetLineStyle(2);

 // Fit it to the graph and draw it 
  gr->Fit(&f); 
      f.DrawClone("Same"); 

 // Build and Draw a legend 
      TLegend leg(.1 ,.7 ,.3 ,.9 ,"Lab. Lesson 1"); 
     leg.SetFillColor (0) ; 
    gr->SetFillColor (0) ; 
     leg.AddEntry(gr,"Exp. Points"); 
     leg.AddEntry(&f,"Th. Law"); 
     leg.DrawClone("Same"); 
         // Draw an arrow on the canvas 
       TArrow arrow(8,8,6.2,23,0.02,"----|>"); 
       arrow.SetLineWidth (2) ; 
       arrow.DrawClone ( ) ; 
         // Add some text to the plot 
 TLatex text(8.2,7.5,"#splitline{Maximum}{Deviation}"); 
         text.DrawClone ( ) ; 

   mycanvas->Update();
   mycanvas -> Print("example.pdf");
}

std::pair<double,double> plotVelocity(std::vector<double> positions, std::vector<double> means, int strip, std::vector<double> errors, char thresholdLabel)
{
	TCanvas* c = new TCanvas();
	TGraphErrors* vGraph = new TGraphErrors(positions.size(), &means[0], &positions[0], 0, &errors[0] );
	vGraph->SetMarkerStyle(20);
	vGraph->SetMarkerSize(1);
	vGraph->Draw("AP");
	vGraph->Fit("pol1");
	TF1* fit = vGraph->GetFunction("pol1");
	std::stringstream buf;
	buf << strip << "_" << thresholdLabel;
	c->SaveAs( (buf.str() + ".png").c_str() );
	std::pair<double,double> resultOfFit = std::make_pair<double,double> ( (fit->GetParameter(1))*-0.2,  (fit->GetParError(1))*-0.2 );
	return resultOfFit;
}

  void XZFit(){   
    if (xerr[0] == 0 && xerr[1] == 0 && xerr[2] == 0 && zerr[0] == 0 && zerr[1] == 0 && zerr[2] == 0) {cout << "ERRORE: Incertezze non specificate" << endl;}
    else {
      //  TGraphErrors *graph1 = new TGraphErrors(3,x,y,xerr,yerr);
    //  TF1* fitfunc1 = new TF1("fittingfunction","[0]+[1]*x",-100,100);
    //   xyfitfunc = fitfunc1;
    // xygraph = graph1;
      delete xzfitfunc;  //Dealloca prima di riallocarne uno nuovo
      delete xzgraph;
xzfitfunc = new TF1("xzfittingfunction", "[0]+[1]*x",-100,100);
        xzgraph = new TGraphErrors(3,xv,yv,xerr,zerr);
    xzfitfunc->SetParameters(30,-2);
    xzgraph->Fit(xzfitfunc,"0S");
    }
  };

void super(){

    TH1::StatOverflows(kTRUE);
    
    cout << "\n\nStarting process...\n\n";

    TGraphErrors *data = readData("magneticField");

    TF1 *f1 = new TF1("fitRight", "[0]+[1]*TMath::Abs(TMath::Sin( [2]*(x-[3]) )/( [2]*(x-[3]) ) )", -200, 200);
    f1->SetParameter(0, 0.);
    f1->SetParameter(1, 0.36);
    f1->SetParameter(2, 0.05);
    f1->SetParameter(3, -15);
    f1->SetParName(0, "C");
    f1->SetParName(1, "N");
    f1->SetParName(2, "#pi#Phi/#Phi_{0}I");
    f1->SetParName(3, "x_{0}");
    f1->SetLineColor(kBlue);

    TF1 *f2 = new TF1("fitLeft", "[0]+[1]*TMath::Abs(TMath::Sin( [2]*(x-[3]) )/( [2]*(x-[3]) ) )");
    f2->SetParameter(0, 50./2000.);
    f2->SetParameter(1, 700./2000.);
    f2->SetParameter(2, 0.06283);
    f2->SetParameter(3, -15);
    f2->SetParName(0, "C");
    f2->SetParName(1, "N");
    f2->SetParName(2, "A");
    f2->SetParName(3, "A");
    f2->SetLineColor(kBlue);

    // data->Fit(f1, "ME", "", 18, 60);
    // data->Fit(f2, "ME", "", -200, -33);

    data->Fit(f1, "ME", "", -90, 70);


    Double_t pi = 3.1415926538, L = 7.5e-6, l = 1.75e-9, lambda = 39.78e-9, conv = 0.54e-4, t=2*lambda+l;
    Double_t fQauntaRight = pi*L*t*conv/f1->GetParameter(2);//, fQauntaLeft = pi*L*t*conv/f2->GetParameter(2);
    Double_t fQuantaError = f1->GetParError(2)*pi*L*t*conv/f1->GetParameter(2)/f1->GetParameter(2);

    cout << "Flux quanta from right is " << fQauntaRight << " +- " << fQuantaError << endl;//<< " and from left is " << fQauntaLeft << endl;

    TCanvas *can = new TCanvas("Canvas", "Canvas", 1600, 900);
    graph(data, "", can, "Solenoid current (mA)", "Zero-voltage current (mA)", "magField");


}

void TGraphFit(){
   //
   // Draw a graph with error bars and fit a function to it
   //
   gStyle->SetOptFit(111) ; //superimpose fit results
   // make nice Canvas
   TCanvas *c1 = new TCanvas("c1" ,"Daten" ,200 ,10 ,700 ,500) ;
   c1->SetGrid( ) ;
   //define some data points ...
   const Int_t n = 10;
   Float_t x[n] = {-0.22, 0.1, 0.25, 0.35, 0.5, 0.61, 0.7, 0.85, 0.89, 1.1};
   Float_t y[n] = {0.7, 2.9, 5.6, 7.4, 9., 9.6, 8.7, 6.3, 4.5, 1.1};
   Float_t ey[n] = {.8 ,.7 ,.6 ,.5 ,.4 ,.4 ,.5 ,.6 ,.7 ,.8};
   Float_t ex[n] = {.05 ,.1 ,.07 ,.07 ,.04 ,.05 ,.06 ,.07 ,.08 ,.05};
   // and hand over to TGraphErros object
   TGraphErrors *gr = new TGraphErrors(n,x,y,ex,ey);
   gr->SetTitle("TGraphErrors with Fit") ;
   gr->Draw("AP");
   // now perform a fit (with errors in x and y!)
   gr->Fit("gaus");
   c1->Update();
}

int
eicsmear_check_calorimeter( TString filename_mc,
			    TString filename_mc_smeared = "")
{
  /* Uncomment this line when running without compilation */
  gSystem->Load("libeicsmear");

  /* Open input files. */
  TFile *file_mc = new TFile(filename_mc, "OPEN");
  TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");

  /* Get trees from files. */
  TTree *tree = (TTree*)file_mc->Get("EICTree");
  TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");

  /* Add friend to match branches in trees. */
  tree->AddFriend(tree_smeared);

  erhic::EventPythia *event  = NULL;
  Smear::Event       *eventS = NULL;

  tree->SetBranchAddress("event", &event);
  tree->SetBranchAddress("eventS", &eventS);

  /* Define base cut for reconstructed final state particles */
  TCut cut_base("particles.KS==1 && Smeared.particles.p > 0 && (particles.id==11 || particles.id==22)");

  /* colors array */
  unsigned colors[7] = {1,2,3,4,6,7,14};

  /* select how many eta settings to plot */
  unsigned etas_plotmax = 7;

  /* Create vector of theta values to include */
  vector< double > etas;
  etas.push_back(-2.75);
  etas.push_back(-2.25);
  etas.push_back(-1.75);
  etas.push_back(-0.25);
  etas.push_back( 0.25);
  etas.push_back( 1.75);
  etas.push_back( 2.25);

  /* Create fit function */
  TF1* f_momres = new TF1("f_momres", "sqrt( ([0])**2 + ([1]/(x**0.5))**2 )" );

  cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;

  /* Create scratch canvas */
  TCanvas *cscratch = new TCanvas("cscratch");

  /* Create framehistogram */
  TH1F* hframe = new TH1F("hframe","",100,0,40);
  hframe->GetYaxis()->SetRangeUser(0,0.15);
  hframe->GetXaxis()->SetTitle("Energy (GeV/c)");
  hframe->GetYaxis()->SetTitle("#sigma_{E}/E");

  /* Create 2Dhistogram */
  unsigned nbinsp = 40;
  TH2F* h2tmp = new TH2F("h2tmp","",nbinsp,0,(float)nbinsp,50,-1,1);
  h2tmp->GetXaxis()->SetTitle("Energy (GeV)");
  h2tmp->GetYaxis()->SetTitle("(#Delta E)/E");

  /* create combined canvas plot */
  TCanvas *c1 = new TCanvas();
  hframe->Draw();

  /* Create legend */
  TLegend* leg_eta = new TLegend( 0.25, 0.40, 0.45, 0.90);

  /* Create resolution-vs-momentum plot with fits for each selected theta value */
  for ( int i = 0; i < etas.size(); i++ )
    {
      /* Switch to scratch canvas */
      cscratch->cd();

      double eta = etas.at(i);

      cout << "\n***Eta = " << eta << endl;

      /* Define range of theta because float comparison with fixed value doesn't work
         too well for cuts in ROOT trees */
      double eta_min = eta - 0.25;
      double eta_max = eta + 0.25;

      /* Cut for tree */
      TCut cutx( Form("particles.p > 1 && ( particles.eta > %f && particles.eta < %f )", eta_min, eta_max) );
      cout << cutx.GetTitle() << endl;

      /* Fill 2D histogram and fit slices with gaussian to get track resolution */
      h2tmp->Reset();
      tree->Draw("(Smeared.particles.E-particles.E)/particles.E:particles.E >> h2tmp",cut_base && cutx,"");
      h2tmp->FitSlicesY();

      /* Get histogram with sigma from gauss fits */
      TH1D* h2tmp_2 = (TH1D*)gDirectory->Get("h2tmp_2");

      /* Fill fit parameters in TGraphErrors */
      TGraphErrors *gres = new TGraphErrors(nbinsp);
      gres->SetMarkerColor(colors[i]);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  
    vector<float> thetabinsp;
    vector<float> thetabinsm;

    for(int ibin=14;ibin<=30;ibin++)    thetabinsp.push_back(ibin*0.05-1);
    for(int ibin=10;ibin<=26;ibin++)    thetabinsm.push_back(ibin*0.05-1);
    
 
    TLatex *t=new TLatex();
    t->SetNDC();

    TGraphErrors *gduthp[nfiles];
    TGraphErrors *gduthm[nfiles];
        
    TF1 *fduthp[nfiles];
    TF1 *fduthm[nfiles];

    stringstream sduthp1[nfiles];
    stringstream sduthp2[nfiles];
    
    stringstream sduthm1[nfiles];
    stringstream sduthm2[nfiles];

    bool fit = false;

    for(unsigned int ifile = 0 ; ifile < nfiles ; ifile++ ){
      
      can[idx]->cd(1);
      
      hduthp[ifile]->SetName(Form("hduthp_%i",ifile));
      gduthp[ifile] = getTGraphFromTH2(hduthp[ifile],thetabinsp, colors[ifile] );
      fduthp[ifile]=new TF1(Form("fduthp_%i",ifile),"pol1");
      fduthp[ifile]->SetLineColor( colors[ifile] );
      if(fit) gduthp[ifile]->Fit(fduthp[ifile]);
         
      float dwp     = fduthp[ifile]->GetParameter(1);
      float dwperr  = fduthp[ifile]->GetParError(1);
      float bp      = fduthp[ifile]->GetParameter(0);
      float bperr   = fduthp[ifile]->GetParError(0);
      
      sduthp1[ifile] << "#DeltaW = " << fround(dwp,3) << " #pm " << fround(dwperr,3) << " #mum" << endl;
      sduthp2[ifile] << "#Deltatan(LA) = " << fround(bp/(235.-dwp),3) << " #pm " << fround(bperr/(235.-dwp),4) << endl;
      
      gduthp[ifile]->SetTitle("TOB (V+)");
      gduthp[ifile]->Draw("AP");
      gduthp[ifile]->GetXaxis()->SetTitle("<tan(#theta_{trk})-tan(#theta_{LA})>");
      gduthp[ifile]->GetYaxis()->SetTitle("<#Deltau> [#mum]");
      //gduthp[ifile]->GetXaxis()->SetLimits(-1,1);
      gduthp[ifile]->SetMinimum(-25);
      gduthp[ifile]->SetMaximum(20);
      
      
      can[idx]->cd(2);

      hduthm[ifile]->SetName(Form("hduthm_%i",ifile));
      gduthm[ifile] = getTGraphFromTH2(hduthm[ifile],thetabinsm, colors[ifile] );
      fduthm[ifile]=new TF1(Form("fduthm_%i",ifile),"pol1");
      fduthm[ifile]->SetLineColor( colors[ifile] );
      if(fit) gduthm[ifile]->Fit(fduthm[ifile]);
         
      float dwm     = fduthm[ifile]->GetParameter(1);
      float dwmerr  = fduthm[ifile]->GetParError(1);
      float bm      = fduthm[ifile]->GetParameter(0);
      float bmerr   = fduthm[ifile]->GetParError(0);
      
      sduthm1[ifile] << "#DeltaW = " << fround(dwm,3) << " #pm " << fround(dwmerr,3) << " #mum" << endl;

void fit_test(char name[])
{
	TCanvas *c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);
	c1->SetGrid();
	
	std::vector<double> xv;
	std::vector<double> yv;
	
	std::vector<double> xverr;
	std::vector<double> yverr;
		
	Double_t x, y, xerr, yerr;
	FILE *f = fopen(name,"r");

	int counter = 1;
	
	while (!feof(f))
	{ 
		fscanf(f,"%lf %lf lf %lf\n", &x, &y, &xerr, &yerr);
		xv.push_back(x*4); // in ns
		yv.push_back(y);
		
		xverr.push_back(xerr*4);
		yverr.push_back(yerr);
	}	
	

			TGraphErrors * gr = new TGraphErrors(xv.size(), &xv[0], &yv[0], &xverr[0], &yverr[0]);
			TF1 *fitFcn = new TF1("fitFcn", fitFunction, 50, 300,5);
			
			double base_line = 0;
			for(int j = 0; j < 20; j++)
			{
				base_line += yv[j];			
			}
			
			base_line /= 20;
			
			//cout << base_line << endl;
			
			
			fitFcn->SetParLimits(0, 15, 40); 
			fitFcn->SetParLimits(1, 85, 100); 
			fitFcn->SetParLimits(2, 5, 40); 
			fitFcn->SetParLimits(3, 1, 25);
			
			
			fitFcn->SetParameter(4, base_line);
			fitFcn->SetParLimits(4, base_line, base_line);
			
			
			//fitFcn->SetParameter(3, 153.622);
			//fitFcn->SetParameter(4, 502.345);
			//fitFcn->SetParameter(5, 86.9206);
			
			gr->Fit("fitFcn", "R");	
			gr->SetMarkerColor(4);
			gr->SetMarkerStyle(kFullCircle);
			//gr->Draw("APE");
			gr->Draw("APE");				

			//gSystem->Sleep(3000);			
			
			c1->Update();
			
			for(int j =0; j < xv.size(); j++)
			{
				cout << xverr[j] << "\t" << yverr[j] << endl;				
			}
							
			xv.clear();
			yv.clear();

		
		
		
	
}