#include<stdio.h>
#include<math.h>
#include<string.h>
#include<TFile.h>
#include<TTree.h>
#include<TH1F.h>
#include<TH2F.h>
#include<TCanvas.h>
#include<TPad.h>
#include<TApplication.h>
#include<TBranch.h>
#include<TStyle.h>
#include<iostream>
#include<fstream>
#include<string>

#define SIZE 4000
#define SIG_SIZE 5
#define SIG_TMIN 70
#define SIG_TMAX 120
#define BKG_TMIN 200
#define BKG_TMAX 250

int main(int argc, char *argv[])
{
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(111);
    gStyle->SetPalette(55);
    
    int i,j;
    char in_file[80];
    strcpy (in_file,argv[1]);
    
    std::string input=in_file;
    // std::string rpl="root";
    // input.replace(input.find("txt"), sizeof(rpl)-1, rpl);
    // const char *output = input.c_str();

    TFile *f=new TFile(in_file,"update");
    TTree *htree = (TTree*)f->Get("htree");
    int nentries=htree->GetEntries();      
    std::cout<<nentries<<std::endl;
    
    int count=0;
    int evnt;
    int run;
    int sub;
    int nhits;
    int nclusters;
    int h_cluster[SIZE];
    int h_xtalk[SIZE];
    int h_pmt[SIZE];
    int h_time[SIZE];
    float c_time[SIZE];
    int c_size[SIZE];
    int h_duration[SIZE];
    int h_glx[SIZE];
    int h_gly[SIZE];
    
    int hits1=0;
    int hits2=0;
    int hits3=0;
	
    TH1F *ht1=new TH1F("ht1","",100,0.5,400.5);
    TH1F *ht2=new TH1F("ht2","",100,0.5,400.5);
    TH1F *ht3=new TH1F("ht3","",100,0.5,400.5);
    TH1F *hhits1=new TH1F("hhits1","",100,0.5,100.5);
    TH1F *hhits2=new TH1F("hhits2","",100,0.5,100.5);
    TH1F *hhits3=new TH1F("hhits3","",100,0.5,100.5);
    TH1F *hdur1=new TH1F("hdur1","",100,0.5,100.5);
    TH1F *hdur2=new TH1F("hdur2","",100,0.5,100.5);
    TH1F *hdur3=new TH1F("hdur3","",100,0.5,100.5);
    TH1F *hduro1=new TH1F("hduro1","",100,0.5,100.5);
    TH1F *hduro2=new TH1F("hduro2","",100,0.5,100.5);
    TH1F *hduro3=new TH1F("hduro3","",100,0.5,100.5);
    TH2F *hcor=new TH2F("hcor","",100,0,100,200,0,200);
    TH1F *hpmt1=new TH1F("hpmt1","",400,0,400);
    TH1F *hpmt2=new TH1F("hpmt2","",400,0,400);
    TH1F *hpmt3=new TH1F("hpmt3","",400,0,400);
    TH1F *hpmt1_1=new TH1F("hpmt1_1","",400,0,400);
    TH1F *hpmt2_1=new TH1F("hpmt2_1","",400,0,400);
    TH1F *hpmt3_1=new TH1F("hpmt3_1","",400,0,400);
    TH1F *hdif=new TH1F("hdif","",100,-50,50);
    TH2F *hmap=new TH2F("hmap","",260,0,260,220,0,220);
    TH2F *htime1=new TH2F("htime1","",210,0.,210,100,50.5,150.5);
    TH2F *htime2=new TH2F("htime2","",401,-0.5,400.5,100,50.5,150.5);
    TH1F *hsig=new TH1F("hsig","",250,0.,250.);
    TH1F *hbck=new TH1F("hbck","",250,0.,250.);
    TH1F *hxtk=new TH1F("hxtk","",250,0.,250.);
    TH1F *hasig=new TH1F("hasig","",2800,0.,2800.);
    TH1F *habck=new TH1F("habck","",2800,0.,2800.);
    TH1F *haxtk=new TH1F("haxtk","",2800,0.,2800.);
    
    htree->SetBranchAddress("evnt",&evnt);
    htree->SetBranchAddress("run",&run);
    htree->SetBranchAddress("sub",&sub);
    htree->SetBranchAddress("nhits",&nhits);
    htree->SetBranchAddress("nclusters",&nclusters);
    htree->SetBranchAddress("c_time",&c_time);
    htree->SetBranchAddress("c_size",&c_size);
    htree->SetBranchAddress("h_cluster",&h_cluster);
    htree->SetBranchAddress("h_xtalk",&h_xtalk);
    htree->SetBranchAddress("h_time",&h_time);
    htree->SetBranchAddress("h_duration",&h_duration);
    htree->SetBranchAddress("h_pmt",&h_pmt);
    htree->SetBranchAddress("h_glx",&h_glx);
    htree->SetBranchAddress("h_gly",&h_gly);
	
    // -------------------------------------
    // first loop to define the normalization (number of good events)
    // -------------------------------------
    for(int i=0;i<nentries;i++){
        htree->GetEntry(i);
        if(nhits!=0){ 

            // count good events (>1 RICH cluster in time)
            if(nclusters>=1){
                int found=0;
                for(int j=0;j<nclusters;j++){
                    if(c_time[j]>SIG_TMIN && c_time[j]<SIG_TMAX && c_size[j]>SIG_SIZE){
                        found=1;
                        //std::cout << " --> " << j <<" "<< c_time[j]<< std::endl;
                    }
                }
                if(found>0)count=count+1;
            }

            // -------------------------------------
            // beam conditions (without selection)
            // -------------------------------------
            for(int j=0;j<nhits;j++){
                if(h_time[j]>SIG_TMIN && h_time[j]<SIG_TMAX) {
                  if(h_xtalk[j]==0){
                    hsig->Fill(evnt/1.e3,1./nentries);
                    hasig->Fill(sub*250+evnt/1.e3,1./nentries);
                  }else{
                    hxtk->Fill(evnt/1.e3,1./nentries);
                    haxtk->Fill(sub*250+evnt/1.e3,1./nentries);
                  }
                }
                if(h_time[j]>BKG_TMIN&&h_time[j]<BKG_TMAX) {
                    hbck->Fill(evnt/1.e3,1./nentries);
                    habck->Fill(sub*250+evnt/1.e3,1./nentries);
                }
            }

        }
    }

    std::cout<<"Normalization Tot "<<nentries<<" Good "<<count<<std::endl;

    // -------------------------------------
    // second loop to fill the histos with proper normalization
    // -------------------------------------
    int eve_size=0;
    int eve_time=0;
    float wc=1./count;
    float we=1./nentries;
    for(int i=0;i<nentries;i++){
        htree->GetEntry(i);
        hits1=0;
        hits2=0;
        hits3=0;
        if(nhits!=0){

            int found=0;
            // -------------------------------------
            // select good events (>1 RICH cluster in time)
            // -------------------------------------
            if(nclusters>=1){
                for(int j=0;j<nclusters;j++){
                    if(c_time[j]>SIG_TMIN && c_time[j]<SIG_TMAX && c_size[j]>SIG_SIZE){
                        found=1;
                        eve_size=c_size[j];
                        eve_time=c_time[j];
                    }
                }
            }

            if(found==1){
            for(int j=0;j<nhits;j++){

                ht1->Fill(h_time[j],wc);
                if(h_cluster[j]==0){
                    ht2->Fill(h_time[j],wc);
                    hcor->Fill(h_duration[j],h_time[j]);
                    htime1->Fill(evnt/1.e3,h_time[j]);
                    htime2->Fill(h_pmt[j],h_time[j]);
                }else{
                    ht3->Fill(h_time[j],wc);
                }

                if(h_time[j]>SIG_TMIN && h_time[j]<SIG_TMAX){
                    hdur1->Fill(h_duration[j],wc);
                    hits1=hits1+1;

                    if(h_cluster[j]==0){
                        hdur2->Fill(h_duration[j],wc);

                        hduro1->Fill(h_duration[j],wc);
                        if(h_xtalk[j]==0)hduro2->Fill(h_duration[j],wc);
                        if(h_xtalk[j]>0)hduro3->Fill(h_duration[j],wc);

                        hpmt1->Fill(h_pmt[j]+1,wc);
                        if(h_xtalk[j]==0)hpmt2->Fill(h_pmt[j]+1,wc);
                        if(h_xtalk[j]>0)hpmt3->Fill(h_pmt[j]+1,wc);

                        hdif->Fill(h_time[j]-eve_time,wc);

                        hmap->Fill(h_glx[j],h_gly[j],1.);
                        hits2=hits2+1;
                    }else{
                        hdur3->Fill(h_duration[j],wc);
                        hits3=hits3+1;
                    }
                }

                if(h_time[j]>BKG_TMIN && h_time[j]<BKG_TMAX){
                    if(h_cluster[j]==0){
                        hpmt1_1->Fill(h_pmt[j]+1,wc);
                        if(h_xtalk[j]==0)hpmt2_1->Fill(h_pmt[j]+1,wc);
                        if(h_xtalk[j]>0)hpmt3_1->Fill(h_pmt[j]+1,wc);

                    }
                }
            }

            hhits1->Fill(hits1,wc);
            hhits2->Fill(hits2,wc);
            hhits3->Fill(hits3,wc);
            }
        }
    }
	
    f->Write("",TObject::kOverwrite);

    //TApplication *myapp=new TApplication("myapp",0,0);
   
    TCanvas *c1=new TCanvas("c1","",1400,700);
    c1->Divide(4,2);

    c1->cd(1);
    ht1->SetLineColor(7);
    ht2->SetLineColor(6);
    ht3->SetLineColor(1);
    ht1->Draw("hist");
    ht2->Draw("histsame");
    ht3->Draw("histsame");
    
    c1->cd(2);
    hdur1->SetLineColor(7);
    hdur2->SetLineColor(6);
    hdur3->SetLineColor(1);
    hdur1->Draw("hist");
    hdur2->Draw("histsame");
    hdur3->Draw("histsame");

    c1->cd(3);
    hduro1->SetLineColor(6);
    hduro2->SetLineColor(4);
    hduro3->SetLineColor(3);
    hduro1->Draw("hist");
    hduro2->Draw("histsame");
    hduro3->Draw("histsame");

    c1->cd(4);
    hhits1->SetLineColor(7);
    hhits1->GetYaxis()->SetRangeUser(0,0.1);
    hhits2->SetLineColor(6);
    hhits3->SetLineColor(1);
    hhits1->DrawCopy("hist");
    hhits2->DrawCopy("histsame");
    hhits3->DrawCopy("histsame");
          
    c1->cd(5);
    hcor->DrawCopy("colz");

    c1->cd(6);
    hpmt1->Add(hpmt1_1,-1);
    hpmt2->Add(hpmt2_1,-1);
    hpmt3->Add(hpmt3_1,-1);
    hpmt1->SetLineColor(6);
    hpmt2->SetLineColor(4);
    hpmt3->SetLineColor(3);
    hpmt1->DrawCopy("hist");
    hpmt2->DrawCopy("histsame");
    hpmt3->DrawCopy("histsame");

    c1->cd(7);
    gPad->SetLogz();
    hmap->DrawCopy("colz");

    /*c1->cd(7);
    hsig->SetLineColor(2);
    hbck->SetLineColor(4);
    hsig->DrawCopy("hist");
    hbck->DrawCopy("histsame");
    hasig->DrawCopy("histsame");
    habck->DrawCopy("histsame");*/
 
    c1->cd(8);
    gPad->SetLogz(0);
    htime2->DrawCopy("colz");
 
    //c1->cd(8);
    //hdif->DrawCopy("");
 
    c1->Write();
    c1->Print("c1.png");
    
    //myapp->Run(true);
    f->Close();
    return 0;
}