#include "TCanvas.h"
#include "TBrowser.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMath.h"
#include <fstream>
#include "TFrame.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TFile.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TBox.h"
#include "TRandom.h"
#include "TObject.h"
#include "TObjString.h"
#include <iostream>
#include <cstdio>
#include <string>
#include <sstream>
#include "TLine.h"
#include "TTree.h"
#include "TBrowser.h"
#include "TF1.h"
#include <TStyle.h>
#include <THStack.h>
#include <TPad.h>
#include "TRandom.h"




void plot(){
    // replace this with name of your singles file
    TFile *f=new TFile("Single_Si29_ded.root","read");
    TString histo, norm, protons;
    
    ///file with a list of the histogram names for one detector
    // format Zn_T_u3_C6D6_runxxxxxx Zn_A_u3_SILI_runxxxxxx Zn_h_info_C6D6_runxxxxxx
    ifstream inni("Si29_C6D64_HistList.dat");
    TCanvas *c=new TCanvas();

    /// replace this with the tof around a resonance, i.e. binlow is the start of the resonance, binhigh is the end of the resonance
    // you can choose any resonance, best is a large one which is usally at big tof values
    /// this will be different for Au and for Zn.

float binlow=2250000;
float binhigh=2400000;

    TLegend *legend = new TLegend(0.65, 0.38, 0.89, 0.89);
    legend->SetBorderSize(1);
    legend->SetFillColor(0);
    legend->SetTextSize(0.03);
    int count2=0;
    int count_protons=0;
    
    
    //// replace this with the number of lines in your file "list"
int lines=59;

            float value[lines];
            float value2[lines];
            float valuesili[lines];
            float valuepkup[lines];
    float staterror[lines];
    while(1){
        
        inni>>histo>>norm>>protons;
        if ( ! inni ) break;

        
        char *s;
        if(count2==0) s="hist";
        if(count2>0)s="hist same";
        count2++;
        cout<<histo<<" "<<norm<<endl;
        
        
        TH1F *h=(TH1F*)f->Get(histo);
        //TH1F *hnorm=(TH1F*)f->Get(norm);
        TH1F *hpro=(TH1F*)f->Get(protons);
        
        
        int bin1=h->FindBin(binlow); //Ge70Bins
        int bin2=h->FindBin(binhigh);
        
        // this gives the total number of count in the resonance
        
        int counter2=h->Integral(bin1,bin2);

      
    
    value[count2-1]=float(counter2)/hpro->GetBinContent(4)*7E12;
        staterror[count2-1]=sqrt(float(counter2))/hpro->GetBinContent(4)*7E12;
             value2[count2-1]=float(counter2);
             valuesili[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(8));
            valuepkup[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(7));
   
        h->Scale(1/hpro->GetBinContent(4)*7E12);
        h->Rebin(50);
        h->SetLineColor(count2);
        h->DrawCopy(s);
        legend->AddEntry(h,histo, "l");
    }
    float sum=0;
    float weight=0;
     float weightsum=0;   
    for(int i=0;i<count2;i++){

        weight=value2[i]/(value[i]*value[i]);

       sum=sum+value[i]*weight;


        weightsum=weightsum+weight;

        cout<<value[i]<<" "<<weight<<" "<<i<<endl;
    }


cout<< sum/weightsum <<endl;
   
    float deviation=0;
double stddev=0;

    for(int i=0;i<count2;i++){
        deviation=(sum/weightsum-value[i])/sum*weightsum*100;

stddev=(sum/weightsum-value[i])*(sum/weightsum-value[i])+stddev;

        cout<<" "<<deviation<<" % "<<sqrt(value2[i])/value2[i]*100<<" % counting "<<deviation/(sqrt(value2[i])/value2[i]*100)<< "   "<<i<<endl;
    }
    deviation=sqrt(1/count2*deviation);
    cout<<count2<<" "<<sum/count2<<endl;

double stddev2=sqrt(1/(float(count2-1))*stddev);

//STOP PROTON SCALING

cout<<"std dev is "<<stddev2/(sum/weightsum)*100<<" %"<<endl;
    legend->Draw();
cout<<"                                "<<endl;
    for(int i=0;i<count2;i++){

        cout<<" "<<valuesili[i]<<"  Prot/SILI    "<<valuepkup[i]<<"  Prot/PKUP   "<<valuesili[i]/valuepkup[i]<<" SILI/PKUP"<<endl;
    }


TH1F *hprosi=new TH1F("","",50,0,50);
TH1F *hprook=new TH1F("","",50,0,50);
TH1F *hsipk=new TH1F("","",50,0,50);
    TH1F *hcounts=new TH1F("","",50,0,50);

for(int i=0;i<count2;i++){
hprosi->SetBinContent(i+1,valuesili[i]);
hprook->SetBinContent(i+1,valuepkup[i]);
hsipk->SetBinContent(i+1,valuesili[i]/valuepkup[i]);
    hcounts->SetBinContent(i+1,value[i]);
    hcounts->SetBinError(i+1,staterror[i]);
    }

TCanvas *c3=new TCanvas;
c3->Divide(2,2);
c3->cd(1);
hprosi->Draw();
c3->cd(2);
hprook->Draw();
c3->cd(3);
hsipk->Draw();
    c3->cd(4);
    hcounts->Draw();
    
    
    TCanvas *c5=new TCanvas;
    hcounts->Draw();
}










void plotau(){
    // replace this with name of your singles file
    TFile *f=new TFile("Single_Au_ded.root","read");
    TString histo, norm, protons;
    
    ///file with a list of the histogram names for one detector
    // format Zn_T_u3_C6D6_runxxxxxx Zn_A_u3_SILI_runxxxxxx Zn_h_info_C6D6_runxxxxxx
    ifstream inni("Au_C6D64_HistList.dat");
    TCanvas *c=new TCanvas();

    /// replace this with the tof around a resonance, i.e. binlow is the start of the resonance, binhigh is the end of the resonance
    // you can choose any resonance, best is a large one which is usally at big tof values
    /// this will be different for Au and for Zn.

float binlow=12000000;
float binhigh=14000000;

    TLegend *legend = new TLegend(0.65, 0.38, 0.89, 0.89);
    legend->SetBorderSize(1);
    legend->SetFillColor(0);
    legend->SetTextSize(0.03);
    int count2=0;
    int count_protons=0;
    
    
    //// replace this with the number of lines in your file "list"
int lines=59;

            float value[lines];
            float value2[lines];
            float valuesili[lines];
            float valuepkup[lines];
    float staterror[lines];
    while(1){
        
        inni>>histo>>norm>>protons;
        if ( ! inni ) break;

        
        char *s;
        if(count2==0) s="hist";
        if(count2>0)s="hist same";
        count2++;
        cout<<histo<<" "<<norm<<endl;
        
        
        TH1F *h=(TH1F*)f->Get(histo);
        //TH1F *hnorm=(TH1F*)f->Get(norm);
        TH1F *hpro=(TH1F*)f->Get(protons);
        
        
        int bin1=h->FindBin(binlow); //Ge70Bins
        int bin2=h->FindBin(binhigh);
        
        // this gives the total number of count in the resonance
        
        int counter2=h->Integral(bin1,bin2);

      
    
    value[count2-1]=float(counter2)/hpro->GetBinContent(4)*7E12;
        staterror[count2-1]=sqrt(float(counter2))/hpro->GetBinContent(4)*7E12;
             value2[count2-1]=float(counter2);
             valuesili[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(8));
            valuepkup[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(7));
   
        h->Scale(1/hpro->GetBinContent(4)*7E12);
        h->Rebin(50);
        h->SetLineColor(count2);
        h->DrawCopy(s);
        legend->AddEntry(h,histo, "l");
    }
    float sum=0;
    float weight=0;
     float weightsum=0;
    for(int i=0;i<count2;i++){

        weight=value2[i]/(value[i]*value[i]);

       sum=sum+value[i]*weight;


        weightsum=weightsum+weight;

        cout<<value[i]<<" "<<weight<<" "<<i<<endl;
    }


cout<< sum/weightsum <<endl;
   
    float deviation=0;
double stddev=0;

    for(int i=0;i<count2;i++){
        deviation=(sum/weightsum-value[i])/sum*weightsum*100;

stddev=(sum/weightsum-value[i])*(sum/weightsum-value[i])+stddev;

        cout<<" "<<deviation<<" % "<<sqrt(value2[i])/value2[i]*100<<" % counting "<<deviation/(sqrt(value2[i])/value2[i]*100)<< "   "<<i<<endl;
    }
    deviation=sqrt(1/count2*deviation);
    cout<<count2<<" "<<sum/count2<<endl;

double stddev2=sqrt(1/(float(count2-1))*stddev);

//STOP PROTON SCALING

cout<<"std dev is "<<stddev2/(sum/weightsum)*100<<" %"<<endl;
    legend->Draw();
cout<<"                                "<<endl;
    for(int i=0;i<count2;i++){

        cout<<" "<<valuesili[i]<<"  Prot/SILI    "<<valuepkup[i]<<"  Prot/PKUP   "<<valuesili[i]/valuepkup[i]<<" SILI/PKUP"<<endl;
    }


TH1F *hprosi=new TH1F("","",50,0,50);
TH1F *hprook=new TH1F("","",50,0,50);
TH1F *hsipk=new TH1F("","",50,0,50);
    TH1F *hcounts=new TH1F("","",50,0,50);

for(int i=0;i<count2;i++){
hprosi->SetBinContent(i+1,valuesili[i]);
hprook->SetBinContent(i+1,valuepkup[i]);
hsipk->SetBinContent(i+1,valuesili[i]/valuepkup[i]);
    hcounts->SetBinContent(i+1,value[i]);
    hcounts->SetBinError(i+1,staterror[i]);
    }

TCanvas *c3=new TCanvas;
c3->Divide(2,2);
c3->cd(1);
hprosi->Draw();
c3->cd(2);
hprook->Draw();
c3->cd(3);
hsipk->Draw();
    c3->cd(4);
    hcounts->Draw();
    
    
    
    
    TCanvas *c5=new TCanvas;
    hcounts->Draw();
}