//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//   Primitive Eventdisplay of tracks
// load macro with .L <macroname>.C+ to compile it and then enter ExecuteTrackdisplay and press <tab> to see the arguments
// runnumber: which run? datanumber: which reconstruction of this run? nevents: process how many events? singletrackevents: only process single track events? dirname: output directory
//
// -> there is a simple python GUI, which could be improved
//
// Environment:
//      GEM ALICE IROC prototype data analysis in fopiroot
//
// Author List:
//      Philipp Gadow          (philipp.gadow@mytum.de)
//
//
//-----------------------------------------------------------

#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 "TH2F.h"
#include "TH2I.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 <TpcPadPlane.h>
#include <TpcPadShapePool.h>
#include <GFTrack.h>
#include <GFTrackCand.h>

//---------------------------------------------
//Helper functions
//---------------------------------------------

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

//---------------------------------------------

void ExecuteTrackdisplay(bool verbose, int runnumber, string infilename, int evnumber,
                         string dirname = "plots/trkdisplay/",
                         TString poolfile = "/nfs/hicran/project/panda/SIM/pgadow/fopiROOT_trunk/tpc/ALICE/par/padshape_iroc.dat",
                         TString padplanefilename = "/nfs/hicran/project/panda/SIM/pgadow/fopiROOT_trunk/tpc/ALICE/par/padplane_iroc.dat",
                         ){

    TString indirname = "/nfs/mds/data/tpc/alice/ps_test_beam/RecoOut/";
    TString inpathname = indirname + infilename;

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


    // create canvas
    TCanvas *c1 = new TCanvas("c1","Simple Event Display",400,0,1600,1100);
    c1->Divide(3,2);
    c1->SetFillColor(00);
    c1->SetGrid();

    // histogram names
    TString histtrkdisplayname = "Simple Track Display ev #" + inttostring(evnumber);
    TString histtrkdisplaydiginame = "Track Display (pad hit) ev #" + inttostring(evnumber);
    TString histtrkdisplayclusname = "Simple Track Display (cluster) ev #" + inttostring(evnumber);
    TString histevdisplayname = "Simple Event Display ev #" + inttostring(evnumber);
    TString histevdisplaydiginame = "Simple Event Display (pad hit) ev #" + inttostring(evnumber);
    TString histevdisplayclusname = "Simple Event Display (cluster) ev #" + inttostring(evnumber);

    // create histogram
    TH2F *histtrkdisplay = new TH2F ("histtrkdisplay",histtrkdisplayname, 63,85.225,131.725,108,-22,22);
    TH2F *histtrkdisplaydigi = new TH2F ("histtrkdisplaydigi",histtrkdisplaydiginame, 63,85.225,131.725,108,-22,22);
    TH2F *histtrkdisplayclus = new TH2F ("histtrkdisplayclus",histtrkdisplayclusname, 63,85.225,131.725,108,-22,22);
    TH2F *histevdisplay = new TH2F ("histevdisplay",histevdisplayname, 63,85.225,131.725,108,-22,22);
    TH2F *histevdisplaydigi = new TH2F ("histevdisplaydigi",histevdisplaydiginame, 63,85.225,131.725,108,-22,22);
    TH2F *histevdisplayclus = new TH2F ("histevdisplayclus",histevdisplayclusname, 63,85.225,131.725,108,-22,22);


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

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


    // 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("TrackPostFit",&tracks);

    // create padplane


    TpcPadShapePool *pool = new TpcPadShapePool(poolfile);
    TpcPadPlane* padplane = new TpcPadPlane(padplanefilename,pool);



    tpcTree->GetEvent(evnumber);

    //###################################################
    //    event display
    //###################################################

    Int_t NumTpcClusters = clusters->GetEntries();


    for (Int_t clu = 0; clu < NumTpcClusters; ++clu){
        TpcCluster* tpccluster = (TpcCluster*) clusters->At(clu);
        double clusx = tpccluster->pos().x();
        double clusy = tpccluster->pos().y();
        histevdisplayclus->Fill(clusx,clusy, tpccluster->amp());
    }   

    Int_t NumTpcDigis = digis->GetEntries();
    for (Int_t digi = 0; digi < NumTpcDigis; ++digi){
        TpcDigi* tpcdigi = (TpcDigi*) digis->At(digi);
        double digix = 0.;
        double digiy = 0.;
        padplane->GetPadXY(tpcdigi->padId(),digix, digiy);
        histevdisplaydigi->Fill(digix,digiy, tpcdigi->amp());
    }

    Int_t NumTpcSamples = samples->GetEntries();
    for (Int_t sample = 0; sample < NumTpcSamples; ++sample){
        TpcSample* tpcsample = (TpcSample*) samples->At(sample);
        double padx = 0.;
        double pady = 0.;
        padplane->GetPadXY(tpcsample->padId(),padx, pady);
        histevdisplay->Fill(padx,pady, tpcsample->amp());
    }


    //###################################################
    //    track display
    //###################################################

    int numberoftracks = tracks->GetEntries();

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

    //loop over tracks
    for(Int_t itr=0;itr<numberoftracks;++itr){
        GFTrack* trk=(GFTrack*)tracks->At(itr);
        GFTrackCand cand = trk->getCand();
        std::vector<unsigned int> hitIDs = cand.getHitIDs(-1); //gets IDs of clusters in track

        // get all digis in the track
        std::vector<const TpcDigi*> tpcdigis;           //vector with digis of the track
        for(int i=0; i<hitIDs.size(); i++){
            TpcCluster* cl;
            cl = (TpcCluster*) clusters->At(hitIDs[i]);
            double trkclusx = cl->pos().x();
            double trkclusy = cl->pos().y();
            histtrkdisplayclus->Fill(trkclusx,trkclusy, cl->amp());

            const std::map<unsigned int,double>* digiMap = cl->getDigiMap(); //all the digis in the cluster

            for(std::map<unsigned int, double>::const_iterator it=digiMap->begin();it!=digiMap->end();++it){
                TpcDigi* aDigi=(TpcDigi*)digis->At((*it).first);
                aDigi->amp(aDigi->amp()*(*it).second); //if the digi is shared assign it a fractional amplitude (fraction is in the 2nd entry in map, usually .5)
                tpcdigis.push_back(aDigi);
            }
        }


        // loop over digis in track
        for (unsigned int idigi = 0; idigi<tpcdigis.size(); ++idigi){

            //finde die zugehoerigen digi IDs
            TpcDigi* trktpcdigi = (TpcDigi*) tpcdigis[idigi];
            double trkdigix = 0.;
            double trkdigiy = 0.;
            padplane->GetPadXY(trktpcdigi->padId(),trkdigix, trkdigiy);
            histtrkdisplaydigi->Fill(trkdigix,trkdigiy, trktpcdigi->amp());

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

            for (it = tpcsampleIDs->begin(); it!=tpcsampleIDs->end(); ++it){

                TpcSample* trktpcsample = (TpcSample*) samples->At(it->first);
                double trkpadx = 0.;
                double trkpady = 0.;
                padplane->GetPadXY(trktpcsample->padId(),trkpadx, trkpady);
                histtrkdisplay->Fill(trkpadx,trkpady, trktpcsample->amp());
            }
        }
        //define output filename
        string filename = "merge_" + inttostring(runnumber) + "_ev_"+inttostring(evnumber)+"_" + inttostring(itr) + "_trackdisplay";

        string pathname1 = dirname + filename+".png";
        string pathname2 = dirname + filename+".C";

        //  Draw the histograms and save plot
        c1->cd(1);
        histevdisplay->GetXaxis()->SetTitle("x (cm)");
        histevdisplay->GetYaxis()->SetTitle("y (cm)");
        histevdisplay->Draw("colz");
        c1->cd(2);
        histevdisplaydigi->GetXaxis()->SetTitle("x (cm)");
        histevdisplaydigi->GetYaxis()->SetTitle("y (cm)");
        histevdisplaydigi->Draw("colz");
        c1->cd(3);
        histevdisplayclus->GetXaxis()->SetTitle("x (cm)");
        histevdisplayclus->GetYaxis()->SetTitle("y (cm)");
        histevdisplayclus->Draw("colz");
        c1->cd(4);
        histtrkdisplay->GetXaxis()->SetTitle("x (cm)");
        histtrkdisplay->GetYaxis()->SetTitle("y (cm)");
        histtrkdisplay->Draw("colz");
        c1->cd(5);
        histtrkdisplaydigi->GetXaxis()->SetTitle("x (cm)");
        histtrkdisplaydigi->GetYaxis()->SetTitle("y (cm)");
        histtrkdisplaydigi->Draw("colz");
        c1->cd(6);
        histtrkdisplayclus->GetXaxis()->SetTitle("x (cm)");
        histtrkdisplayclus->GetYaxis()->SetTitle("y (cm)");
        histtrkdisplayclus->Draw("colz");

        c1->SaveAs(pathname1.data());
        c1->SaveAs(pathname2.data());


    }
}