#include<iostream>

void eta_energi()
{

  gROOT->SetStyle("Plain");
  //create histogram
  TH1F *hist1= new TH1F("hist1","Electronenergies",100,30000,53000);
  
  //create a chain with tree Atlfast/51 
  TChain chain("CBNT/t3333");  
  // add the files            
  
  chain.Add("ntuples/elec_50_eta0/ntuple*.root");
  //define events
  Long64_t nbrEntries = chain.GetEntries();  
  //cout<<"Antall events: "<<nbrEntries<<endl;
  
  float eg_et[2000];
  float eg_eta[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
   chain.SetBranchAddress("eg_eta",&eg_eta);
  //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      
    hist1->Fill(eg_et[0]);   
  }  
  

  TH1F *hist2= new TH1F("hist2","Electronenergies",100,35000,60000);
  TChain chain("CBNT/t3333");  
  chain.Add("ntuples/elec_50_eta1/ntuple*.root");
  //define events
  Long64_t nbrEntries = chain.GetEntries();  
  //cout<<"Antall events: "<<nbrEntries<<endl;
  
  float eg_et[2000];
  float eg_eta[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
   chain.SetBranchAddress("eg_eta",&eg_eta);
  //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      
    hist2->Fill(eg_et[0]);   
  }  


  TH1F *hist3= new TH1F("hist3","Electronenergies",3,45000,110000);
  TChain chain("CBNT/t3333");  
  chain.Add("ntuples/elec_50_eta1-5/ntuple*.root"); 
  //define events
  Long64_t nbrEntries = chain.GetEntries();  
  //cout<<"Antall events: "<<nbrEntries<<endl;
  
  float eg_et[2000];
  float eg_eta[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
   chain.SetBranchAddress("eg_eta",&eg_eta);
  //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      
    hist3->Fill(eg_et[0]);   
  }  



  TH1F *hist4= new TH1F("hist4","Electronenergies",50,40000,55000);
  TChain chain("CBNT/t3333");  
  chain.Add("ntuples/elec_50_eta2/ntuple*.root");
  //define events
  Long64_t nbrEntries = chain.GetEntries();  
  //cout<<"Antall events: "<<nbrEntries<<endl;
  
  float eg_et[2000];
  float eg_eta[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
  chain.SetBranchAddress("eg_eta",&eg_eta);
  //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      
    hist4->Fill(eg_et[0]);   
  }  


  TH1F *hist5= new TH1F("hist5","Electronenergies", 50,30000,60000);
  TChain chain("CBNT/t3333");  
  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];
  float eg_eta[2000];
  //Specify the address where to read 
  chain.SetBranchAddress("eg_et",&eg_et);
   chain.SetBranchAddress("eg_eta",&eg_eta);
  //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      
    hist5->Fill(eg_et[0]);   
  }  
   



  TCanvas *c = new TCanvas ("c", "Energifordeling", 200, 220, 800, 600);
  c->Divide(2,3);
  c->cd(1); 
  hist1->Draw();
  TF1 *fit1 = new TF1("fit1","gaus");
  hist1->Fit("fit1");
  
  c->cd(2);
  hist2->Draw();
  TF1 *fit2 = new TF1("fit2","gaus");
  hist2->Fit("fit2");

  c->cd(3);
  hist3->Draw();
  TF1 *fit3 = new TF1("fit3","gaus");
  hist3->Fit("fit3");

  c->cd(4);  
  hist4->Draw();
  TF1 *fit4 = new TF1("fit4","gaus");
  hist4->Fit("fit4");
  
  c->cd(5);  
  hist5->Draw();
  TF1 *fit5 = new TF1("fit5","gaus");
  hist5->Fit("fit5");
//hist->SetXTitle("Eta");
//hist->SetYTitle("Events");
//hist->Draw();
//gStyle->SetOptStat(1000011);
  
  //Fitting gauscurves to each energi distribution
 
     
 //create table of sigmas  
  float p3[5];  
  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);
  
 
  
  //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;
  
 

  float frac[5] =  {p3[0]/sqrt(50000),p3[1]/sqrt(50000),p3[2]/sqrt(50000),p3[3]/sqrt(50000),p3[4]/sqrt(50000)};


  //create table with energidistribution 
  float Eta[5] = {0,1,1.5,2,2.5};

     
   // TGraph *graph = new TGraph(8,Energi,p3);
 
   TGraphErrors *graph = new TGraphErrors(5,Eta,frac);
  
     
  //Printing each point (x,y) on screen just to check  
   /*  
  for(int t=0;t<5;t++)
  {
    cout<<"frac["<<t<<"] : "<<frac[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 );
   }  
  
// Sets initial values and parameter names
   
   TF1 *fitf = new TF1("fitf",fitf,5,.1,.1);
   
   fitf->SetParameters(10,0.5,0.4);
   fitf->SetParNames("a","b","c");

   

  gStyle->SetOptStat(111111); 
  gStyle->SetOptFit(1011);
  //gStyle->GetYaxisOffset(1.1);
  //create the graph and draw it
  //graph->new TGraphErrors(8,Energi,frac);
    graph->SetTitle("#frac{#sigma_{E}}{#sqrt{E}} vs Eta");
    graph->GetXaxis()->SetTitle("Eta [Pseudorapidity]");
    graph->GetYaxis()->SetTitle("#frac{#sigma_{E}}{#sqrt{E}} [%*#sqrt{GeV}]");
    graph->SetMarkerStyle(21);
    //graph->SetPointError(bin);
    
    graph->Draw("ACZ*");
    graph->Fit("fitf","r");  

   
 
  
 
}