#include "CbmRichMCbmToTShifter.h"

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


#include "CbmRichDigi.h"
#include "TLatex.h"
#include "CbmDrawHist.h"
#include "CbmRichGeoManager.h"
#include "CbmRichPoint.h"
#include "CbmRichConverter.h"

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

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

using boost::assign::list_of;

#define RichZPos 348.

CbmRichMCbmToTShifter::CbmRichMCbmToTShifter() :
                        FairTask("CbmRichMCbmToTShifter"),
                        fRichDigis(nullptr),
                        fHM(nullptr),
                        fEventNum(0),
                        fOutputDir("result_ToTOffset")
{

}

InitStatus CbmRichMCbmToTShifter::Init()
{
    std::cout << "CbmRichMCbmToTShifter::Init" << std::endl;

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

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

    InitHistograms();

    return kSUCCESS;
}

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

    fHM->Create1<TH1D>("fhNofEvents","fhNofEvents;Entries", 1, 0.5, 1.5);
    fHM->Create1<TH1D>("fhToT_all_begin","fhToT_all_begin;Entries", 601 , 19.975, 50.025);
    fHM->Create1<TH1D>("fhToT_all_end","fhToT_all_end;Entries", 301 , 19.975, 35.025);
    
    for (int i=1; i<33;++i){
        fHM->Create1<TH1D>(Form("fhToT_7321_%d",i),Form("fhToT_7321_%d;Entries",i), 601 , 19.975, 50.025);
        fHM->Create1<TH1D>(Form("fhToT_7211_%d",i),Form("fhToT_7211_%d;Entries",i), 601 , 19.975, 50.025);
        fHM->Create1<TH1D>(Form("fhToT_7220_%d",i),Form("fhToT_7220_%d;Entries",i), 601 , 19.975, 50.025);
        fHM->Create1<TH1D>(Form("fhToT_7320_%d",i),Form("fhToT_7320_%d;Entries",i), 601 , 19.975, 50.025);
        fHM->Create1<TH1D>(Form("fhToT_7221_%d",i),Form("fhToT_7221_%d;Entries",i), 601 , 19.975, 50.025);
        
        fHM->Create1<TH1D>(Form("fhToT_32ch_%d",i),Form("fhToT_32ch_%d;Entries",i), 601 , 19.975, 50.025);
    }
}


void CbmRichMCbmToTShifter::Exec(Option_t* /*option*/)
{
    fEventNum++;
    fHM->H1("fhNofEvents")->Fill(1);
    
    int nofDigis = fRichDigis->GetEntriesFast();

    std::cout << "CbmRichMCbmToTShifter, event No. " <<  fEventNum << std::endl;
    std::cout << "\t Number of Digis in TS: " <<  nofDigis << std::endl;
    
    
    for (int i=0; i < nofDigis; ++i){
        CbmRichDigi* digi = static_cast<CbmRichDigi*>(fRichDigis->At(i));
        std::cout << "Address: 0x"<<std::hex <<getDirichAddress(digi->GetAddress())<<std::dec <<"   " << digi->GetToT() << std::endl;
       // if ( getDirichAddress(digi->GetAddress()) == 0x7321 && (getDirichChannel(digi->GetAddress())== 1))  
        fHM->H1(Form("fhToT_%x_%i",getDirichAddress(digi->GetAddress()),getDirichChannel(digi->GetAddress())))->Fill(digi->GetToT());
        if (getDirichAddress(digi->GetAddress()) != 0x7321) fHM->H1(Form("fhToT_32ch_%d",getDirichChannel(digi->GetAddress())))->Fill(digi->GetToT());
        fHM->H1("fhToT_all_begin")->Fill(digi->GetToT());
    }

}

void CbmRichMCbmToTShifter::DrawHist()
{
    std::cout.precision(4);
 
    SetDefaultDrawStyle();
    double nofEvents = fHM->H1("fhNofEvents")->GetEntries();
    fHM->ScaleByPattern("fh_.*", 1./nofEvents);

    {
        fHM->CreateCanvas("rich_mcbm_fhNofEvents","rich_mcbm_fhNofEvents", 600 , 600);
        DrawH1(fHM->H1("fhNofEvents"));
    }
    
    {
        TCanvas* c = fHM->CreateCanvas("fhToT_7321","fhToT_7321", 1000 , 1000);
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            DrawH1(fHM->H1(Form("fhToT_7321_%d",i)));
        }
    }

    {
        TCanvas* c = fHM->CreateCanvas("fhToT_7211","fhToT_7211", 1000 , 1000);
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            DrawH1(fHM->H1(Form("fhToT_7211_%d",i)));
        }
    }
    
        {
        TCanvas* c = fHM->CreateCanvas("fhToT_7220","fhToT_7220", 1000 , 1000);
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            DrawH1(fHM->H1(Form("fhToT_7220_%d",i)));
        }
    }
    
        {
        TCanvas* c = fHM->CreateCanvas("fhToT_7320","fhToT_7320", 1000 , 1000);
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            DrawH1(fHM->H1(Form("fhToT_7320_%d",i)));
        }
    }
    
    {
        TCanvas* c = fHM->CreateCanvas("fhToT_7221","fhToT_7221", 1000 , 1000);
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            DrawH1(fHM->H1(Form("fhToT_7221_%d",i)));
        }
    }
    
    {
        fHM->CreateCanvas("fhToT_all_begin","fhToT_all_begin", 1000 , 1000);
        DrawH1(fHM->H1("fhToT_all_begin"));

    }
    
    {
        fHM->CreateCanvas("fhToT_all_end","fhToT_all_end", 1000 , 1000);
        fHM->H1("fhToT_all_end")->Draw("HIST");
        
        Double_t par[9];
        TF1 *g1    = new TF1("g1","gaus",22,23.8);
        TF1 *g2    = new TF1("g2","gaus",24.3,26);
        TF1 *total = new TF1("total","gaus(0)+gaus(3)",22,26);
        total->SetLineColor(2);
        fHM->H1("fhToT_all_end")->Fit(g1,"R");
        fHM->H1("fhToT_all_end")->Fit(g2,"R+");
        g1->GetParameters(&par[0]);
        g2->GetParameters(&par[3]);
        total->SetParameters(par);
        fHM->H1("fhToT_all_end")->Fit(total,"R+");
        //g1->Draw("SAME");
        //g2->Draw("SAME");
        total->Draw("SAME");
    }    
    
    {
        TCanvas* c = fHM->CreateCanvas("fhToT_32ch","fhToT_32ch", 1600 , 1600);
        
        c->Divide(6,6);
        for(int i=1; i<33;++i){
            c->cd(i);
            fHM->H1(Form("fhToT_32ch_%d",i))->Draw("HIST");
//             TLine *line = new TLine(1,c->GetUymin(),80,c->GetUymax());
//             line->SetLineColor(kRed);
//             line->Draw("");
        }

    }    
}



void CbmRichMCbmToTShifter::Finish()
{   
    
    std::cout<<"Drawing Hists..."<<std::endl;
    
    for (int i=1;i<33;++i){
        Double_t peak = FindPeak(Form("fhToT_32ch_%d",i));
        std::cout<< "PEAK "<< i  <<": " << peak << std::endl;
        std::cout<< "ToT 7321: "<< i  <<": " << FindPeak(Form("fhToT_7321_%d",i)) << std::endl;
        std::cout<< "ToT 7211: "<< i  <<": " << FindPeak(Form("fhToT_7211_%d",i)) << std::endl;
        std::cout<< "ToT 7220: "<< i  <<": " << FindPeak(Form("fhToT_7220_%d",i)) << std::endl;
        std::cout<< "ToT 7320: "<< i  <<": " << FindPeak(Form("fhToT_7320_%d",i)) << std::endl;
        std::cout<< "ToT 7221: "<< i  <<": " << FindPeak(Form("fhToT_7221_%d",i)) << std::endl;
        
        Int_t nofBins = fHM->H1(Form("fhToT_32ch_%d",i))->GetNbinsX();
   
        for (Int_t j=1;j<=nofBins; ++j){
            Double_t cont = fHM->H1(Form("fhToT_32ch_%d",i))->GetBinContent(j);
            fHM->H1("fhToT_all_end")->Fill(fHM->H1(Form("fhToT_32ch_%d",i))->GetBinCenter(j)-peak+25.0,cont);
        }
        
    }
    
    DrawHist();
    fHM->SaveCanvasToImage(fOutputDir,"png,pdf");
    //fHM->WriteToFile();
}


void CbmRichMCbmToTShifter::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);
}

Double_t CbmRichMCbmToTShifter::FindPeak( const std::string histName){
//    Double_t peak = -999.9;
    std::vector<Int_t> peaks;
//    Double_t height  = 0.0;
    
    const Int_t width = 5; // how many left and right from Center. 2 is 5 Bins
    
    //std::array tmpCont[width];
   
    Int_t nofBins = fHM->H1(histName)->GetNbinsX();
    std::cout <<"NofBins: "<< nofBins << std::endl;
    Double_t slope = 0.;
    Double_t slopeOld = 0.;
    //TODO: Find Peak most right!! ; First Find all Peaks
    
    // divide Histogram in Blocks of 2*width+1 to calculate slopes.
    for (Int_t i=1;i<=nofBins; i= i+2*width+1){
//         Double_t cont = fHM->H1(histName)->GetBinContent(i);
//         if (cont > height) {
//             peak = i;
//             height = cont;
//         }
        
        if (i<= width || i > (nofBins-width)) continue; // skip borders for beginning
        
        // calculate the slope
        slopeOld = slope;
        slope = 0.;
        for (int k=1;k<(width*2+1);++k){
            slope += (fHM->H1(histName)->GetBinContent(i- width+k) - fHM->H1(histName)->GetBinContent(i- width+k-1));
        }
        slope /= (width*2);
        
        Double_t slopeDiff = slopeOld - slope;
        
        //cut on the slope: from + to -: Peak ; Cut on significance of peak
        if (slopeOld > 0.1 && slope < -0.1 && (slopeDiff > 2.8)) {
            Double_t tmpMax = 0.;
            Int_t pos = -99;
            // Find Maximum in Slope Block
            for (int m=0;m<4*width+2;++m){
                Double_t tmplocMax = fHM->H1(histName)->GetBinContent(i- (width*3+1)+m);
                if (tmplocMax > tmpMax ) {
                    tmpMax = tmplocMax;
                    pos    = i- (width*3+1)+m;
                }
            }
            
            peaks.emplace_back(pos); // peak
            //std::cout<<"PEAK_intern: "<< fHM->H1(histName)->GetBinCenter(i) << "  PeakCalc: " << fHM->H1(histName)->GetBinCenter(pos) << "  slope:" << slope << "  slopeOld:" << slopeOld << "  slopeDiff: " << slopeDiff << std::endl;
        }
    }
    if (peaks.size() == 0) return -1;
    if (peaks.back() <  0) return -2;
    return fHM->H1(histName)->GetBinCenter(peaks.back());
    
}

ClassImp(CbmRichMCbmToTShifter)