//dEdx-Plot
//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 "TMath.h"
#include "TVector3.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>
#include <GFTrack.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 ExecutedEdxplot(bool verbose, int runnumber, int datasetnumber, int numevents, bool truncatedmean = true, double cutlow=0., double cuthigh=.4, unsigned int dEdxinputdirectory = 0){


    int lowtrunc = int ((cutlow+0.005) *100);
    int hightrunc = int ((cuthigh+0.005) *100);

    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(s)
    TFile *reco = new TFile(inpathname);

    // create canvas
    TCanvas *cdEdx = new TCanvas("cdEdx","dEdx",1000,700);
    TCanvas *cTracklength = new TCanvas("cTracklength","underlying distribution for cut",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->Divide(2,2);
    cdEdx->SetGrid();


    TString histname = "dEdx";
    if (truncatedmean){
        histname+=" truncated mean ("+doubletostring(cutlow)+", "+doubletostring(cuthigh)+") -";
    }
    else{
        histname+=" simple mean -";
    }
    histname+= " single track events";

    TString histnamePID = histname+ ": e and Pion cut";
    TString histnameE1 = histname + ": e-cher cut(>400)";
    TString histnameE2 = histname + ": e- double cut(>400)";
    TString histnamePi = histname + ": Pionen";

    // create histogram
    TH1F *histall = new TH1F ("histall", histname, 600,0,600);
    TH1F *histPID = new TH1F ("histPID", histnamePID, 600,0,600);
    TH1F *histelEsinglecut = new TH1F ("histEsinglecut", histnameE1, 600,0,600);
    TH1F *histEledoublecut = new TH1F ("histEdoublecut", histnameE2, 600,0,600);
    TH1F *histpion = new TH1F ("histpion", histnamePi, 400,0,400);
    TH1F *histtracklength = new TH1F ("histtracklength", "TrackLength", 160,0,80);
    // create Tree
    TTree *tpcTree = (TTree*)reco->Get("cbmsim");


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

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

    // set branches
    if(dEdxinputdirectory>0){
        //merge_681_13_3_dEdxDATA.reco.root
        TString DEDXindirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdxDATA/";
        TString DEDXinfilename = "merge_" + inttostring(runnumber) + "_"+inttostring(datasetnumber)+"_"+inttostring(dEdxinputdirectory)+"_dEdxDATA.reco.root";
        TString DEDXinpathname = DEDXindirname + DEDXinfilename;
        tpcTree->AddFriend("dedxcbmsim = cbmsim",DEDXinpathname);
        tpcTree->SetBranchStatus("dEdx.*",1);
        tpcTree->SetBranchAddress("dedxcbmsim.dEdx",&dedx);
    }
    else{
        tpcTree->SetBranchStatus("dEdx.*",1);
        tpcTree->SetBranchAddress("dEdx",&dedx);
    }

    tpcTree->SetBranchAddress("TrackPostFit",&tracks);
    tpcTree->SetBranchAddress("GEMEvent",&gemevent);

    // loop over dedx

    //events


    if (numevents == 0){
        numevents = tpcTree->GetEntries();
    }
    std::cout<<std::endl;
    for( Int_t ev = 0; ev <numevents; ev++){

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

        dedx->Clear();
        tracks->Clear();
        gemevent->Clear();

        tpcTree->GetEvent(ev);

        GEMEvent* gemEvent = (GEMEvent*)gemevent->At(0); //there is only one element in gemevent
        Int_t NumTpcdEdx = dedx->GetEntries();
        if (NumTpcdEdx == 0) continue;


        int numberoftracks = tracks->GetEntries();
        if(verbose){
            std::cout << "number of tracks per event:" << numberoftracks << std::endl;
        }
        //select only single track events
        if (numberoftracks!=1){
                continue;
            }

        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);
            }
            histtracklength->Fill(tracklength);

            // select only tracks that go all over the detector
            if (tracklength < 35){ ///
                continue;
            }
            TVector3 position = ((GFTrack*) tracks->At(0))->getCardinalRep()->getPos();
            TVector3 momentum = ((GFTrack*) tracks->At(0))->getCardinalRep()->getMom();

            //select only tracks going in very forward direction
            TVector3 normalXZplane (0,1,0);
            double trackphi = TMath::ASin(momentum.Dot(normalXZplane)/momentum.Mag());

//            if (TMath::Abs(trackphi > 0.005)){
//                continue;
//            }

            //select only tracks with z-position above 2.713 cm

            if (position.Z() < 2.7){
                continue;
            }

            if (position.Z()>10.6){
	      continue;
	    }

            double mean = 0;
            if (truncatedmean){
                mean = tpcdedx->truncMean(cutlow,cuthigh);
            }
            else{
                mean = tpcdedx->simpleMean();
            }
            if(verbose){
                std::cout<< "num of dedx entries:" << NumTpcdEdx <<std::endl;
                std::cout<< "num of dedx entries in single dedx:" << tpcdedx->nEntries() <<std::endl;
                for(unsigned int i = 0; i < tpcdedx->nEntries(); i++){
                    std::cout<<"dedx" << idedex<<"       de: " << tpcdedx->getdx(i)<< ", dx: "<<tpcdedx->getdE(i)<<std::endl;
                }
                std::cout<<mean<<std::endl;
            }
            histall->Fill(mean);
            // is it an electron?
            if (gemEvent->GetCherenkovValue()>400){
                histelEsinglecut->Fill(mean);
                histPID->Fill(mean);
            }
            if (gemEvent->GetCherenkovValue()>400 && gemEvent->GetPbGlassValue()>400){
                histEledoublecut->Fill(mean);
            }

            // is it a pion?
            if (gemEvent->GetCherenkovValue() < 150){
                histpion->Fill(mean);
                histPID->Fill(mean);
            }
        }
    }


    //define output filename
    TString dirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/";
    TString filename = "merge_" + inttostring(runnumber) + "_"+ inttostring(22) + "_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);
    histall->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histall->GetYaxis()->SetTitle("# of entries");
//    histall->Draw();
    histall->Write();
    cdEdx->cd(2);
    histpion->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histpion->GetYaxis()->SetTitle("# of entries");
//    histpion->Draw();
    histpion->Write();
    cdEdx->cd(3);
    histPID->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histPID->GetYaxis()->SetTitle("# of entries");
//    histPID->Draw();
    histPID->Write();

    cdEdx->cd(4);
    histelEsinglecut->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histelEsinglecut->GetYaxis()->SetTitle("# of entries");
//    histelEsinglecut->Draw();
    histelEsinglecut->Write();

    histEledoublecut->GetYaxis()->SetTitle("# of entries");
//    histEledoublecut->Draw();
    histEledoublecut->Write();

//    cdEdx->SaveAs(pathname1.Data());
//    cdEdx->SaveAs(pathname2.Data());
//    cdEdx->SaveAs(pathname3.Data());
//    cdEdx->Write();

    cTracklength->cd();
//    histtracklength->Draw();
    histtracklength->Write();
    std::cout<<std::endl;
//    cTracklength->SaveAs(pathnamelength1.Data());
//    cTracklength->SaveAs(pathnamelength2.Data());
//    cTracklength->SaveAs(pathnamelength3.Data());

    savefile->Close();
  }