#include "CbmRichMCbmQaReal.h"

#include "TH1D.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TStyle.h"
#include "TEllipse.h"
#include "TMarker.h"
#include "TGeoNode.h"
#include "TGeoManager.h"
#include "TGeoBBox.h"


#include "CbmRichDigi.h"
#include "CbmRichHit.h"
#include "CbmRichRingFinderHoughImpl.h"
#include "TLatex.h"
#include "CbmDrawHist.h"
#include "CbmTrackMatchNew.h"
#include "CbmRichRing.h"
#include "CbmRichHit.h"
#include "CbmMatchRecoToMC.h"
#include "CbmRichGeoManager.h"
#include "CbmRichPoint.h"
#include "utils/CbmRichDraw.h"
#include "CbmGlobalTrack.h"
#include "CbmTrdTrack.h"
#include "CbmTofHit.h"
#include "CbmTofDigi.h"

#include "CbmRichConverter.h"

#include "CbmUtils.h"
#include "CbmHistManager.h"

#include <iostream>
#include <string>
#include <boost/assign/list_of.hpp>
#include <sstream>

using namespace std;
using boost::assign::list_of;

CbmRichMCbmQaReal::CbmRichMCbmQaReal() :
                        FairTask("CbmRichMCbmQaReal"),
                        fHM(nullptr),
                        fEventNum(0),
                        fOutputDir(""),
                        fRichDigis(nullptr),
                        fRichHits(nullptr),
                        fRichRings(nullptr),
                        fStsDigis(nullptr),
                        fT0Digis(nullptr),
                        fTofDigis(nullptr),
                        fNofDrawnRings(0)

{

}

InitStatus CbmRichMCbmQaReal::Init()
{
    cout << "CbmRichMCbmQaReal::Init" << endl;

    FairRootManager* ioman = FairRootManager::Instance();
    if (nullptr == ioman) { Fatal("CbmRichMCbmQaReal::Init","RootManager not instantised!"); }

    fRichDigis =(TClonesArray*) ioman->GetObject("CbmRichDigi");
    if (nullptr == fRichDigis) { Fatal("CbmRichMCbmQaReal::Init", "No Rich Digis!");}

    fRichHits =(TClonesArray*) ioman->GetObject("RichHit");
    if (nullptr == fRichHits) { Fatal("CbmRichMCbmQaReal::Init", "No Rich Hits!");}

    fRichRings =(TClonesArray*) ioman->GetObject("RichRing");
    if (nullptr == fRichRings) { Fatal("CbmRichMCbmQaReal::Init", "No Rich Rings!");}

    fStsDigis =(TClonesArray*) ioman->GetObject("CbmStsDigi");
    if (nullptr == fStsDigis) { Fatal("CbmRichMCbmQaReal::Init", "No STS Digis!");}

    fTofDigis =(TClonesArray*) ioman->GetObject("CbmTofDigi");
    if (nullptr == fTofDigis) { Fatal("CbmRichMCbmQaReal::Init", "No Tof Digis!");}

    fT0Digis =(TClonesArray*) ioman->GetObject("CbmT0Digi");
    if (nullptr == fT0Digis) { Fatal("CbmRichMCbmQaReal::Init", "No T0 Digis!");}

    InitHistograms();

    return kSUCCESS;
}

void CbmRichMCbmQaReal::InitHistograms()
{
    fHM = new CbmHistManager();

    fHM->Create1<TH1D>("fhNofEvents","fhNofEvents;Entries", 1, 0.5, 1.5);

    // nof objects per timeslice
    fHM->Create1<TH1D>("fhNofRichDigisInTimeslice","fhNofRichDigisInTimeslice;# RICH digis / timeslice;Entries", 150, -0.5, 1499.5);
    fHM->Create1<TH1D>("fhNofRichHitsInTimeslice","fhNofRichHitsInTimeslice;# RICH hits / timeslice;Entries", 150, -0.5, 1499.5);
    fHM->Create1<TH1D>("fhNofRichRingsInTimeslice","fhNofRichRingsInTimeslice;# RICH rings / timeslice;Entries", 15, -0.5, 14.5);

    // RICH hits
    fHM->Create2<TH2D>("fhRichHitXY","fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries", 100, -20, 20,  100, -20, 20);

    // RICH digis, the limits of log histograms are set in Exec method
    fHM->Create1<TH1D>("fhRichDigisTimeLog","fhNofRichDigisTimeLog;Time [ns];Entries", 200, 0., 0.);
    fHM->Create1<TH1D>("fhRichDigisTimeLogZoom","fhNofRichDigisTimeLogZoom;Time [ns];Entries", 200, 0., 0.);

    fHM->Create1<TH1D>("fhTofDigisTimeLog","fhTofDigisTimeLog;Time [ns];Entries", 200, 0., 0.);
    fHM->Create1<TH1D>("fhTofDigisTimeLogZoom","fhTofDigisTimeLogZoom;Time [ns];Entries", 200, 0., 0.);

    fHM->Create1<TH1D>("fhStsDigisTimeLog","fhStsDigisTimeLog;Time [ns];Entries", 200, 0., 0.);
    fHM->Create1<TH1D>("fhStsDigisTimeLogZoom","fhStsDigisTimeLogZoom;Time [ns];Entries", 200, 0., 0.);

    fHM->Create1<TH1D>("fhT0DigisTimeLog","fhT0DigisTimeLog;Time [ns];Entries", 200, 0., 0.);
    fHM->Create1<TH1D>("fhT0DigisTimeLogZoom","fhT0DigisTimeLogZoom;Time [ns];Entries", 200, 0., 0.);

    // RICH rings
    fHM->Create2<TH2D>("fhRichRingXY","fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 100, -20, 20,  100, -20, 20);
    fHM->Create1<TH1D>("fhRichRingRadius","fhRichRingRadius;Ring radius [cm];Entries", 100, 0., 10.);
    fHM->Create1<TH1D>("fhNofHitsInRing","fhNofHitsInRing;# hits in ring;Entries", 50, -0.5, 49.5);
}



void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
{
    fEventNum++;
    fHM->H1("fhNofEvents")->Fill(1);

    cout << "CbmRichMCbmQaReal, event No. " <<  fEventNum << endl;

    if (fEventNum == 1) {
        double minTime = std::numeric_limits<double>::max();
        for (int i = 0; i < fRichDigis->GetEntries(); i++) {
            CbmRichDigi* richDigi = static_cast<CbmRichDigi*>(fRichDigis->At(i));
            if (richDigi->GetTime() < minTime) minTime = richDigi->GetTime();
        }

        double dT = 40e9;
        double dTZoom1 = 0.4e9;
       // double dTZoom2 = 0.004e9;
        fHM->H1("fhRichDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
        fHM->H1("fhRichDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);

        fHM->H1("fhTofDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
        fHM->H1("fhTofDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);

        fHM->H1("fhStsDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
        fHM->H1("fhStsDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);

        fHM->H1("fhT0DigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
        fHM->H1("fhT0DigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
    }

    int nofRichDigis = fRichDigis->GetEntries();
    fHM->H1("fhNofRichDigisInTimeslice")->Fill(nofRichDigis);
    for (int i = 0; i < nofRichDigis; i++) {
        CbmDigi* digi = static_cast<CbmDigi*>(fRichDigis->At(i));
        fHM->H1("fhRichDigisTimeLog")->Fill(digi->GetTime() );
        fHM->H1("fhRichDigisTimeLogZoom")->Fill(digi->GetTime() );
    }

    int nofTofDigis = fTofDigis->GetEntries();
    for (int i = 0; i < nofTofDigis; i++) {
        CbmDigi* digi = static_cast<CbmDigi*>(fTofDigis->At(i));
        fHM->H1("fhTofDigisTimeLog")->Fill(digi->GetTime() );
        fHM->H1("fhTofDigisTimeLogZoom")->Fill(digi->GetTime() );
    }

    int nofStsDigis = fStsDigis->GetEntries();
    for (int i = 0; i < nofStsDigis; i++) {
        CbmDigi* digi = static_cast<CbmDigi*>(fStsDigis->At(i));
        fHM->H1("fhStsDigisTimeLog")->Fill(digi->GetTime() );
        fHM->H1("fhStsDigisTimeLogZoom")->Fill(digi->GetTime() );
    }

    int nofT0Digis = fT0Digis->GetEntries();
    for (int i = 0; i < nofT0Digis; i++) {
        CbmDigi* digi = static_cast<CbmDigi*>(fT0Digis->At(i));
        fHM->H1("fhT0DigisTimeLog")->Fill(digi->GetTime() );
        fHM->H1("fhT0DigisTimeLogZoom")->Fill(digi->GetTime() );
    }


    int nofRichHits = fRichHits->GetEntries();
    fHM->H1("fhNofRichHitsInTimeslice")->Fill(nofRichHits);
    for (int iH = 0; iH < nofRichHits; iH++) {
        CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
        fHM->H2("fhRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
    }

    RichRings();
}

void CbmRichMCbmQaReal::RichRings()
{
    int nofRichRings = fRichRings->GetEntries();
    fHM->H1("fhNofRichRingsInTimeslice")->Fill(nofRichRings);
    for (int i = 0; i < nofRichRings; i++) {
        CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(i));
        if (ring == nullptr) continue;
        DrawRing(ring);
        fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
        fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
        fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
    }
}

void CbmRichMCbmQaReal::DrawHist()
{
    cout.precision(4);

    SetDefaultDrawStyle();
    double nofEvents = fHM->H1("fhNofEvents")->GetEntries();
    fHM->ScaleByPattern("fh_.*", 1./nofEvents);

    {
        TCanvas* c = fHM->CreateCanvas("rich_mcbm_fhNofEvents","rich_mcbm_fhNofEvents", 600 , 600);
        DrawH1(fHM->H1("fhNofEvents"));
    }

    {
        TCanvas* c = fHM->CreateCanvas("rich_mcbm_nofObjectsInTimeslice","rich_mcbm_nofObjectsInTimeslice", 1200 , 1200);
        c->Divide(2,2);
        c->cd(1);
        DrawH1(fHM->H1("fhNofRichDigisInTimeslice"));
        c->cd(2);
        DrawH1(fHM->H1("fhNofRichHitsInTimeslice"));
        c->cd(3);
        DrawH1(fHM->H1("fhNofRichRingsInTimeslice"));
    }

    {
        TCanvas* c = fHM->CreateCanvas("rich_mcbm_XY","rich_mcbm_XY", 1200 , 600);
        c->Divide(2,1);
        c->cd(1);
        DrawH2(fHM->H2("fhRichHitXY"));
        c->cd(2);
        DrawH2(fHM->H2("fhRichRingXY"));
    }


    {
        TCanvas* c = fHM->CreateCanvas("rich_mcbm_richDigisTimeLog","rich_mcbm_richDigisTimeLog", 1200 , 1200);
        c->Divide(1,2);
        c->cd(1);
        DrawH1({fHM->H1("fhStsDigisTimeLog"), fHM->H1("fhTofDigisTimeLog"), fHM->H1("fhT0DigisTimeLog"), fHM->H1("fhRichDigisTimeLog")},
                {"STS", "TOF", "T0", "RICH"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
        gPad->SetLeftMargin(0.1);
        gPad->SetRightMargin(0.10);
        fHM->H1("fhStsDigisTimeLog")->GetYaxis()->SetTitleOffset(0.7);
        fHM->H1("fhStsDigisTimeLog")->SetMinimum(0.9);
        c->cd(2);
        DrawH1({fHM->H1("fhStsDigisTimeLogZoom"), fHM->H1("fhTofDigisTimeLogZoom"),   fHM->H1("fhT0DigisTimeLogZoom"), fHM->H1("fhRichDigisTimeLogZoom")},
                {"STS", "TOF", "T0", "RICH"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
        gPad->SetLeftMargin(0.1);
        gPad->SetRightMargin(0.1);
        fHM->H1("fhStsDigisTimeLogZoom")->GetYaxis()->SetTitleOffset(0.7);
        fHM->H1("fhStsDigisTimeLogZoom")->SetMinimum(0.9);
    }

    {
        TCanvas* c = fHM->CreateCanvas("rich_mcbm_rings","rich_mcbm_rings", 1200 , 600);
        c->Divide(2,1);
        c->cd(1);
        DrawH1(fHM->H1("fhRichRingRadius"));
        c->cd(2);
        DrawH1(fHM->H1("fhNofHitsInRing"));
    }
}

void CbmRichMCbmQaReal::DrawRing(
                              CbmRichRing* ring)
{

    if (fNofDrawnRings > 20) return;
    fNofDrawnRings++;
    stringstream ss;
    ss << "Event" << fNofDrawnRings;
    fNofDrawnRings++;
    TCanvas *c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 500, 500);
    c->SetGrid(true, true);
    TH2D* pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -5., 5., 1, -5., 5);
    pad->SetStats(false);
    pad->Draw();

    // find min and max x and y positions of the hits
    // in order to shift drawing
    double xmin = 99999., xmax = -99999., ymin = 99999., ymax = -99999.;
    for (int i = 0; i < ring->GetNofHits(); i++){
        Int_t hitInd = ring->GetHit(i);
        CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
        if (nullptr == hit) continue;
        if (xmin > hit->GetX()) xmin = hit->GetX();
        if (xmax < hit->GetX()) xmax = hit->GetX();
        if (ymin > hit->GetY()) ymin = hit->GetY();
        if (ymax < hit->GetY()) ymax = hit->GetY();
    }
    double xCur = (xmin + xmax) / 2.;
    double yCur = (ymin + ymax) / 2.;

    //Draw circle and center
    TEllipse* circle = new TEllipse(ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, ring->GetRadius());
    circle->SetFillStyle(0);
    circle->SetLineWidth(3);
    circle->Draw();
    TEllipse* center = new TEllipse(ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, .1);
    center->Draw();

    // Draw hits
    for (int i = 0; i < ring->GetNofHits(); i++){
        Int_t hitInd = ring->GetHit(i);
        CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
        if (nullptr == hit) continue;
        TEllipse* hitDr = new TEllipse(hit->GetX() - xCur, hit->GetY() - yCur, .25);
        hitDr->SetFillColor(kRed);
        hitDr->Draw();
    }

    //Draw information
    stringstream ss2;
    ss2 << "(x,y,r, n)=(" << setprecision(3) << ring->GetCenterX() << ", " << ring->GetCenterY()<< ", " << ring->GetRadius() << ", " << ring->GetNofHits()<<")";
    TLatex* latex = new TLatex(-4., 4., ss2.str().c_str());
    latex->Draw();
}


void CbmRichMCbmQaReal::Finish()
{
    DrawHist();
    //fHM->SaveCanvasToImage(fOutputDir);
    fHM->WriteToFile();
}


void CbmRichMCbmQaReal::DrawFromFile(
        const string& fileName,
        const string& outputDir)
{
    fOutputDir = outputDir;

    if (fHM != nullptr) delete fHM;

    fHM = new CbmHistManager();
    TFile* file = new TFile(fileName.c_str());
    fHM->ReadFromFile(file);

    DrawHist();

    fHM->SaveCanvasToImage(fOutputDir);
}


ClassImp(CbmRichMCbmQaReal)