//Creates Illustration to show working process of CFD
//author: P.Gadow
//status: done

#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 <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <sstream>

#include <TpcSample.h>
#include <TpcEvent.h>
#include <TpcCluster.h>
#include <TpcDigi.h>
#include <GFTrack.h>

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



void ExecuteCFDIllu(bool verbose, int runnumber, int nevents, bool singletrackevents){

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

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

    // create canvas
    TCanvas *c1 = new TCanvas("c1","samples: amp vs t",200,10,700,500);
    c1->SetFillColor(00);
    c1->SetGrid();

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

    // create arrays
    TClonesArray *clusters = new TClonesArray("TpcCluster");
    TClonesArray *digis = new TClonesArray("TpcDigi");
    TClonesArray *samples = new TClonesArray("TpcSample");
    TClonesArray *tracks = new TClonesArray("GFTrack");

    // set branches
    tpcTree->SetBranchAddress("TpcSample",&samples);
    tpcTree->SetBranchAddress("TpcDigi",&digis);
    tpcTree->SetBranchAddress("TpcCluster",&clusters);
    tpcTree->SetBranchAddress("TpcTrackPreFit",&tracks);

    // loop over events, clusters, digis and samples
    //events
    if (nevents == 0){
        nevents = tpcTree->GetEntries();
    }

    for( Int_t ev = 0; ev <nevents; ev++){

        if(ev%1==0){
            cout<<"#################### Event done: "<<ev<<"##############"<<endl;
        }
        tpcTree->GetEvent(ev);

        int numberoftracks = tracks->GetEntries();
        std::cout << "number of tracks per event:" << numberoftracks << std::endl;

        //select only single track events
        if (singletrackevents){
            if (numberoftracks!=1)
                continue;
        }

        Int_t NumTpcClusters = clusters->GetEntries();
        if (NumTpcClusters == 0) continue;

        //clusters - gehe array von clustern durch, die zu einem event gehoeren
        for (Int_t clu = 0; clu < NumTpcClusters; ++clu){
            TpcCluster* tpccluster = (TpcCluster*) clusters->At(clu);
            Int_t tpcdigispercluster =  tpccluster->size();

            if (verbose){
                cout<<"Cluster: "<<clu<<endl;
            }
            if (tpcdigispercluster == 0) continue;


            //finde die zugehoerigen digi IDs
            std::vector<unsigned int> tpcdigiIDs = tpccluster->getDigiIDs();

            for(unsigned int iDigi = 0; iDigi < tpcdigiIDs.size();++iDigi){
                TpcDigi* tpcdigi = (TpcDigi*) digis->At(tpcdigiIDs[iDigi]);

                double tpcdigiT = tpcdigi->t();
                Double_t digiTx[]  = {tpcdigiT,tpcdigiT};


                // finde die zugehoerigen sample IDs
                const std::map<unsigned int,double>* tpcsampleIDs = tpcdigi->getSampleMap();
                std::map<unsigned int,double>::const_iterator it;
                if(tpcsampleIDs->size()==0) {
                    continue;
                }

                if (verbose){
                    cout<<"Digi: "<<iDigi<<" t: "<<tpcdigiT<<" nSample: "<<tpcsampleIDs->size()<<endl;
                }

                std::vector<int> sampleT;
                std::vector<int> sampleAmp;

                std::vector<int> sampleTdelayed;
                std::vector<int> sampleAmpdelayed;

                std::vector<int> sampleAmpdiff;

                int size = 0;
                double maxamp = 0;
                double minT = 10000000; //just a very large number
                double maxT = 0.;
                for (it = tpcsampleIDs->begin(); it!=tpcsampleIDs->end(); ++it){

                    if(verbose){
                        cout<<"sampleId "<<it->first<<" weight: "<<it->second<<endl;
                    }

                    if (it->first>4294967294) {
                        cerr << "sampleID is too great!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
                        continue;
                    }
                    TpcSample* tpcsample = (TpcSample*) samples->At(it->first);
                    ++size;
                    double ampl = tpcsample->amp();
                    double time = tpcsample->t();
                    sampleT.push_back(time);
                    sampleAmp.push_back(ampl-0);
                    if (ampl>maxamp){
                        maxamp = ampl;
                    }
                    if (time>maxT){
                        maxT = time;
                    }
                    if (time<minT){
                        minT = time;
                    }

                    if(verbose){
                        cout<<"Sample: pad: "<< tpcsample->padId()<<" t: "<<tpcsample->t()
                                <<" A: "<<tpcsample->amp()<<endl;
                    }
                    // delayed sample and difference
                    if(it->first>2){
                        TpcSample* tpcsampledelayed = (TpcSample*) samples->At(it->first-2);
                        sampleAmpdelayed.push_back((tpcsampledelayed->amp()-0)*2);
                        sampleTdelayed.push_back(time);
                        sampleAmpdiff.push_back((tpcsampledelayed->amp()-0)*2-ampl+0);

                    }
                    else{
                        sampleAmpdiff.push_back(-ampl);
                    }
                }



                //original graph
                Double_t digiTy[]  = {0.,maxamp};
                const int* t = sampleT.data();
                const int* amp = sampleAmp.data();
                TGraph* gr = new TGraph(size,t,amp);
                gr->GetXaxis()->SetRangeUser(minT-40,maxT+15);
                gr->GetYaxis()->SetRangeUser(-1.5*maxamp,2.2*maxamp);
                TGraph* gr2 = new TGraph(2,digiTx,digiTy);
                //delayed graph
                const int* ampd = sampleAmpdelayed.data();
                const int* td = sampleTdelayed.data();
                TGraph* grdelay = new TGraph(size,td,ampd);

                //difference graph
                const int* ampdiff = sampleAmpdiff.data();
                TGraph* grdiff = new TGraph(size,t,ampdiff);


                string outdirname = "plots/shapeofdigi/CFDIllu/";
                string outfilename = "CFDillustrationmerge_" + inttostring(runnumber) + "_digiplot_" + inttostring(ev) + "_" + inttostring(clu) + "_"+ inttostring(iDigi);
                string outpathname1 = outdirname + outfilename + ".C";
                string outpathname2 = outdirname + outfilename + ".pdf";
                // Draw Plot
                gr->SetLineColor(2);
                gr->SetLineWidth(4);
                gr->SetMarkerColor(4);
                gr->SetMarkerStyle(21);
                gr->SetTitle("Constant Fraction Discriminator");
                gr->GetXaxis()->SetTitle("time t");
                gr->GetYaxis()->SetTitle("amplitude amp");
                gr->Draw("AP");
                grdelay->SetLineColor(0);
                grdelay->SetLineWidth(0);
                grdelay->SetMarkerStyle();
                grdelay->Draw("CP");
                grdiff->SetLineColor(5);
                grdiff->SetLineWidth(2);
                grdiff->Draw("CP");
                gr2->SetLineColor(2);
                gr2->SetLineWidth(4);
                gr2->Draw("C");


                c1->Update();
                c1->GetFrame()->SetFillColor(00);
                c1->GetFrame()->SetBorderSize(12);
                c1->Modified();
                c1->SaveAs(outpathname1.data());
                c1->SaveAs(outpathname2.data());
            }

        }

    }

}