#include<iostream>

void sigma()
{

  gROOT->SetStyle("Plain");
  //create histogram
  TH1F *hist= new TH1F("hist","Electronenergies",500,0,600000);
  //create a chain with tree Atlfast/51 
  TChain chain("CBNT/t3333");  
  // add the files            
  chain.Add("ntuples/elec_100_eta1/ntuple*.root");
  chain.Add("ntuples/elec_250_eta1/ntuple*.root");
  chain.Add("ntuples/elec_25_eta1/ntuple*.root"); 
  chain.Add("ntuples/elec_10_eta1/ntuple*.root");
  chain.Add("ntuples/elec_150_eta1/ntuple*.root");
  chain.Add("ntuples/elec_200_eta1/ntuple*.root");
  chain.Add("ntuples/elec_500_eta1/ntuple*.root");
  // chain.Add("ntuples/elec_50_eta0/ntuple*.root");
  chain.Add("ntuples/elec_50_eta1/ntuple*.root");
  //chain.Add("ntuples/elec_50_eta1-5/ntuple*.root"); 
  // chain.Add("ntuples/elec_50_eta2/ntuple*.root");
  //chain.Add("ntuples/elec_50_eta2-5/ntuple*.root");
  //define events
  Long64_t nbrEntries = chain.GetEntries();  
  //cout<<"Antall events: "<<nbrEntries<<endl;
  
  float eg_et[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
  
  //make loop for all events  
  for(int i=0;i<nbrEntries;i++)
  {
    chain.GetEntry(i);
    // cout<<" Elektron energi for event nr: "<<i<<" "<<eg_et[0]<<" GeV"<<endl;
    //fill histogram with energi      
    hist->Fill(eg_et[0]);   
  }  
  canvas1 = new TCanvas("canvas1","Energy distribusjon",500,500);
  
  hist->SetXTitle("Energy [MeV]");
  hist->SetYTitle("Events");
  hist->Draw();
  gStyle->SetOptStat(1000011);
  
  //Fitting gauscurves to each energi distribution
 
  TF1 *fit1 = new TF1("fit1","gaus",1000,17000);
  hist->Fit("fit1","R");
  TF1 *fit2 = new TF1("fit2","gaus",17000,35000);
  hist->Fit("fit2","R+");
  TF1 *fit3 = new TF1("fit3","gaus",35000,65000);
  hist->Fit("fit3","R+");
  TF1 *fit4 = new TF1("fit4","gaus",80000,120000);
  hist->Fit("fit4","R+");
  TF1 *fit5 = new TF1("fit5","gaus",125000,175000);
  hist->Fit("fit5","R+");
  TF1 *fit6 = new TF1("fit6","gaus",180000,220000);
  hist->Fit("fit6","R+");
  TF1 *fit7 = new TF1("fit7","gaus",225000,275000);
  hist->Fit("fit7","R+");
  TF1 *fit8 = new TF1("fit8","gaus",450000,550000);
  hist->Fit("fit8","R+");
  // TF1 *fit9 = new TF1("fit9","gaus",420,470);
  //hist->Fit("fit9","R+");
  // TF1 *fit10 = new TF1("fit10","gaus",470,550);
  //hist->Fit("fit10","R+");
   
 //create table of sigmas  
  float p3[8];  
   p3[0] = fit1->GetParameter(2);
   p3[1] = fit2->GetParameter(2);
   p3[2] = fit3->GetParameter(2);
   p3[3] = fit4->GetParameter(2);
   p3[4] = fit5->GetParameter(2);
   p3[5] = fit6->GetParameter(2);
   p3[6] = fit7->GetParameter(2);
   p3[7] = fit8->GetParameter(2);
   // p3[8] = fit9->GetParameter(2);
   // p3[9] = fit10->GetParameter(2);
 
  
  //print sigma on screen 
  cout<<" p3 (sigma) er: "<<p3[0]<<endl;
  cout<<" p3 (sigma) er: "<<p3[1]<<endl;
  cout<<" p3 (sigma) er: "<<p3[2]<<endl;
  cout<<" p3 (sigma) er: "<<p3[3]<<endl;
  cout<<" p3 (sigma) er: "<<p3[4]<<endl;
  cout<<" p3 (sigma) er: "<<p3[5]<<endl;
  cout<<" p3 (sigma) er: "<<p3[6]<<endl;
  cout<<" p3 (sigma) er: "<<p3[7]<<endl;
  /*
  cout<<" p3 (sigma) er: "<<p3[8]<<endl;
  cout<<" p3 (sigma) er: "<<p3[9]<<endl;
  */
  

  //TH2F *hist2= new TH2F("hist2","sigma/energi fordeling" ,1100,0,1100,1000,0,13);

  //create table with energidistribution 
   float Energi[8] = {10000,25000,50000,100000,150000,200000,250000,500000};

   float frac[8] = {(p3[0]*100)/Energi[0],(p3[1]*100)/Energi[1],(p3[2]*100)/Energi[2],(p3[3]*100)/Energi[3],(p3[4]*100)/Energi[4],(p3[5]*100)/Energi[5],(p3[6]*100)/Energi[6],(p3[7]*100)/Energi[7]};
  
   // TGraph *graph = new TGraph(8,Energi,p3);
 
   TGraphErrors *graph = new TGraphErrors(8,Energi,frac);
  
  /*   
  //Printing each point (x,y) on screen just to check  
  
  for(int t=0;t<12;t++)
  {
    cout<<"Energi["<<t<<"] : "<<Energi[t]<<endl;
    cout<<"p3["<<t<<"] : "<<p3[t]<<endl;
      
    //hist2->Fill(Energi[t],p3[t]);
   }
   */


   

  canvas2 = new TCanvas("canvas2","Sigma/Energy",400,400);

   Double_t fitf(Double_t *x, Double_t *par)
   {
     Double_t arg = 0;
     Double_t arg1 = 0;
     if (x[0] != 0) arg = (par[0]*par[0])/x[0];
     if (x[0] != 0) arg1 = (par[2]*par[2])/(x[0]*x[0]);
     Double_t fitval = sqrt( arg + par[1]*par[1] + arg1 );
   }  
    
   TF1 *fitf = new TF1("fitf",fitf,0,550000,3);
   
   fitf->SetParameters(10,0.5,0.4);
   fitf->SetParNames("a","b","c");


   
 /*
  //Histogram istead of graph
  
   hist2->SetXTitle("energi [GEV]");
   hist2->SetYTitle("sigma");
   hist2->SetMarkerStyle(21);
   hist2->SetMarkerSize(1);
   hist2->Draw();
  */
  //gStyle->SetOptFit(1111);
  //gStyle->SetTitleOffset(1.5);
  //create the graph and draw it
  //graph->new TGraphErrors(8,Energi,frac);
    graph->SetTitle("#frac{#sigma_{E}}{E} vs E");
    graph->GetXaxis()->SetTitle("Energy [MeV]");
    graph->GetYaxis()->SetTitle("#frac{#sigma_{E}}{E} [%]");
    graph->SetMarkerStyle(21);
    graph->Draw("ACZ*");
    graph->Fit("fitf","r"); 
   
  





}