//dEdx-NumberofSamples
//Author: P. Gadow
//status: version 1.0

#include <TClonesArray.h>
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TFrame.h"
#include "TF1.h"
#include "TLegend.h"
#include "TArrow.h"
#include "TLatex.h"
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TH1F.h"

#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <sstream>
#include <cmath>

#include <TpcSample.h>
#include <TpcEvent.h>
#include <TpcCluster.h>
#include <TpcDigi.h>
#include <TpcdEdx.h>
#include <GEMEvent.h>

std::string inttostring(Int_t input)
{
    std::string s;
    std::stringstream out;
    out << input;
    s = out.str();
    return s;
}

string doubletostring(double input)
{
    string s;
    stringstream out;
    out << input;
    s = out.str();
    return s;
}


void ExecutedEdxNumberofdEdxSamples(bool verbose, int runnumber, int datasetnumber, unsigned int nev){

    TString indirname = "/nfs/mds/data/tpc/alice/ps_test_beam/RecoOut/";
    TString infilename = "merge_" + inttostring(runnumber) + "_"+inttostring(datasetnumber)+"_corrected_hough_smoothed.reco.root";

    TString inpathname = indirname + infilename;

    // open file
    TFile *reco = new TFile(inpathname);



    // create canvas
    TCanvas *cdEdx = new TCanvas("cdEdx","dEdx",1000,700);
    cdEdx->Divide(2,1);

    TString histname = "# of dEdx samples ";

    // create histogram
    TH1F *histdistributionel = new TH1F ("histdistributionel", histname + "e", 100,0,100);
    TH1F *histdistributionpi = new TH1F ("histdistributionpi", histname + "pion", 100,0,100);

    // create Tree
    TTree *tpcTree = (TTree*)reco->Get("cbmsim");

    //save time and load only necessary branches
    tpcTree->SetBranchStatus("*",0);
    tpcTree->SetBranchStatus("dEdx.*",1);
    tpcTree->SetBranchStatus("TrackPostFit.*",1);
    tpcTree->SetBranchStatus("GEMEvent.*",1);
    tpcTree->SetBranchStatus("TpcCluster.*",1);

    // create arrays
    TClonesArray *dedx = new TClonesArray("TpcdEdx");
    TClonesArray *tracks = new TClonesArray("GFTrack");
    TClonesArray *gemevent = new TClonesArray("GEMEvent");
    TClonesArray *clusters = new TClonesArray("TpcCluster");

    // set branches
    tpcTree->SetBranchAddress("dEdx",&dedx);
    tpcTree->SetBranchAddress("TrackPostFit",&tracks);
    tpcTree->SetBranchAddress("GEMEvent",&gemevent);
    tpcTree->SetBranchAddress("TpcCluster",&clusters);

    if (nev == 0){
        nev = tpcTree->GetEntries();
    }

    //event

    for (unsigned int ev = 0; ev < nev; ev++){

        if(ev%1000==0){
            std::cout<<"\r#################### Event done: "<<ev<<"##############";
        }

        tpcTree->GetEvent(ev);
        GEMEvent* gemEvent = (GEMEvent*)gemevent->At(0); //there is only one element in gemevent


        Int_t NumTpcdEdx = dedx->GetEntries();

        int numberofclusters = clusters->GetEntries();
        if(verbose){
            std::cout << numberofclusters << std::endl;
        }

        int numberoftracks = tracks->GetEntries();

        if (numberoftracks!=1) continue;

        if(verbose){
            std::cout << "number of tracks per event:" << numberoftracks << std::endl;
        }



        for (Int_t idedex = 0; idedex < NumTpcdEdx; ++idedex){
            TpcdEdx* tpcdedx = (TpcdEdx*) dedx->At(idedex);

            //tracklength
            //detector length in x = 46.5 cm
            double tracklength = 0.;
            for(unsigned int idx=0;idx<tpcdedx->nEntries();++idx){
                tracklength += tpcdedx->getdx(idx);
            }
            if(verbose){
                std::cout << "tracklength( 40<l<50 ? ): " << tracklength << std::endl;
                std::cout <<"truncated mean(0,25): " << tpcdedx->truncMean(0,.25) << ", simple mean: " << tpcdedx->simpleMean() << std::endl;
                std::cout<< "num of dedx entries:" << NumTpcdEdx <<std::endl;
                std::cout<< "num of dedx entries in single dedx:" << tpcdedx->nEntries() <<std::endl;
            }


            if (gemEvent->GetCherenkovValue()>400){
                histdistributionel->Fill(tpcdedx->nEntries());

            }

            if (gemEvent->GetCherenkovValue()<150){
                histdistributionpi->Fill(tpcdedx->nEntries());

            }


            for(unsigned int i = 0; i < tpcdedx->nEntries(); i++){
                double dEdx = tpcdedx->getdE(i)/tpcdedx->getdx(i);
                if(verbose){
                    std::cout<<"dedx" << i<<"       de: " << tpcdedx->getdx(i)<< ", dx: "<<tpcdedx->getdE(i) << "dE/dx: " << dEdx<<std::endl;
                }
            }

        }
    }
    //    //define output filename
    //    //TString dirname = "plots/dEdx/truncations/";
    //    TString dirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/";
    //
    //    TString filename = "merge_" + inttostring(runnumber) + "_"+ inttostring(datasetnumber) + "_dEdx";
    //
    //    if(truncatedmean){
    //        filename+= "_truncmean_"+inttostring(lowtrunc)+"_"+inttostring(hightrunc);
    //    }
    //    else{
    //        filename+= "_simplemean";
    //    }
    //    TString savefilename = dirname + filename + "histograms.root";
    //    TFile *savefile = new TFile(savefilename.Data(),"RECREATE");

    //    TString pathname1 = dirname + filename + ".C";
    //    TString pathname2 = dirname + filename + ".png";
    //    TString pathname3 = dirname + filename + ".pdf";
    //
    //    TString pathnamelength1 = dirname + filename + "_tracklength.C";
    //    TString pathnamelength2 = dirname + filename + "_tracklength.png";
    //    TString pathnamelength3 = dirname + filename + "_tracklength.pdf";

    //    savefile->cd();
    //  Draw the histogram
    cdEdx->cd(1);
    histdistributionel->GetXaxis()->SetTitle("Number of dE/dx samples/track e");
    histdistributionel->GetYaxis()->SetTitle("# of entries");
    histdistributionel->Draw();

    cdEdx->cd(2);
    histdistributionpi->GetXaxis()->SetTitle("Number of dE/dx samples/track pion");
    histdistributionpi->GetYaxis()->SetTitle("# of entries");
    histdistributionpi->Draw();

}