//dEdx-Plot of Straggling Function (assuming the tracks are all similar)
//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 ExecutedEdxLandau(bool verbose, int runnumber, int datasetnumber, unsigned int nevents = 0){

    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 *cdEdxel = new TCanvas("cdEdxel","dEdx electrons",1000,700);
    TCanvas *cdEdxpi = new TCanvas("cdEdxpi","dEdx pions",1000,700);
    TString histname = "dE/dx ";

    // create histogram
    TH1F *histelectron = new TH1F ("histelectron", histname+"electrons", 1000,0,2000);
    TH1F *histpion = new TH1F ("histpion", histname+"pions", 1000,0,2000);


    // 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);



    //events

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


    for( Int_t ev = 0; ev < nevents; 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(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 (tracklength> 50 || tracklength<35)
                continue;

            if(verbose){
                std::cout << "tracklength( 35<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){
                for(unsigned int i = 0; i < tpcdedx->nEntries(); i++){
                    double dEdx = tpcdedx->getdE(i)/tpcdedx->getdx(i);
                    histelectron->Fill(dEdx);
                    if(verbose){
                        std::cout<<"ELECTRON dedx " << i<<"       de: " << tpcdedx->getdx(i)<< ", dx: "<<tpcdedx->getdE(i) << "dE/dx: " << dEdx<<std::endl;
                    }
                }
            }

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



    cdEdxel->cd();
    histelectron->GetXaxis()->SetTitle("dE/dx");
    histelectron->GetYaxis()->SetTitle("# of entries");
    histelectron->Draw();
    cdEdxpi->cd();
    histpion->GetXaxis()->SetTitle("dE/dx");
    histpion->GetYaxis()->SetTitle("# of entries");
    histpion->Draw();
    cdEdxel->SaveAs("LandauE.pdf");
    cdEdxpi->SaveAs("LandauP.pdf");
}