#include "CbmTrdMCQa.h"

#include "CbmHistManager.h"
#include "CbmTrdAddress.h"
#include "CbmTrdPoint.h"
#include "CbmSimulationReport.h"
//#include "CbmTrdMCQaReport.h"

#include "FairRootManager.h"
#include "FairLogger.h"

#include "TClonesArray.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TProfile.h"
#include "TProfile2D.h"

using std::vector;
using std::map;

CbmTrdMCQa::CbmTrdMCQa():
  FairTask()
  , fHM(new CbmHistManager())
  , fTrdPoints(NULL)
  , fMCTracks(NULL)
  , fNofStation(10)
{
}

CbmTrdMCQa::~CbmTrdMCQa()
{
  if ( fHM ) delete fHM;
}

InitStatus CbmTrdMCQa::Init()
{
  LOG(info) << "Trd Setup consist of " << fNofStation << " stations.";
    
  ReadDataBranches();
  CreateHistograms();
  return kSUCCESS;
}

void CbmTrdMCQa::ReadDataBranches()
{  
  FairRootManager* ioman = FairRootManager::Instance();
  if ( NULL == ioman )
    LOG(fatal) << "No FairRootManager!";
  
  fTrdPoints = dynamic_cast<TClonesArray*>(ioman -> GetObject("TrdPoint"));
  if ( NULL == fTrdPoints )
    LOG(error) << "No TrdPoint array!";

  fMCTracks = dynamic_cast<TClonesArray*>(ioman -> GetObject("MCTrack"));
  if ( NULL == fMCTracks )
    LOG(error) << "No MCTrack array!";
}

void CbmTrdMCQa::CreateHistograms()
{
  CreateNofObjectsHistograms();
  CreatePointHistograms();
  fHM -> Create1<TH1F>("h_trd_EventNo_MCQa", "h_trd_EventNo_MCQa", 1, 0, 1.);
}

void CbmTrdMCQa::CreateNofObjectsHistograms()
{

  Int_t nofBins = 100;
  Double_t minX = -0.5;
  Double_t maxX = 99.5;
  string name = "h_trd_NofObjects_";
  fHM -> Create1<TH1F>(name + "Points", name + "Points;Objects per event;Entries", nofBins, minX, maxX);
  
  nofBins = fNofStation;
  minX = -0.5;
  maxX = static_cast<Float_t>(fNofStation)-0.5;
  fHM -> Create1<TH1F>(name + "Points_Station", name + "Points_Station;Station number;Objects per event", nofBins, minX, maxX);
}

void CbmTrdMCQa::CreatePointHistograms()
{
  for (Int_t stationId = 0; stationId < fNofStation; stationId++){
    Int_t nofBins = 100;
    Double_t minX = -0.5;
    Double_t maxX = 99.5;
    fHM -> Create1<TH1F>(Form("h_trd_MultPoints_Station%i",stationId),
			 Form("Mult, Station %i;Objects per event;Entries", stationId),
			 nofBins, minX, maxX);

    fHM -> Create2<TH2F>(Form("h_trd_PointsMap_Station%i",stationId),
       Form("TrdPoint, Station %i;x, cm;y, cm", stationId),
       120, -60., 60., 120, -60., 60.);
    fHM -> Create2<TH2F>(Form("h_trd_PointsMapEvent_Station%i",stationId),
       Form("TrdPoint/cm^{2}, Station %i;x, cm;y, cm", stationId),
       120, -60., 60., 120, -60., 60.);
    fHM -> Create2<TH2F>(Form("h_trd_PointsMapRate_Station%i",stationId),
       Form("TrdPoint/cm^{2}/s, Station %i;x, cm;y, cm", stationId),
       120, -60., 60., 120, -60., 60.);
    fHM -> Create1<TH1F>(Form("h_trd_XPos_Station%i",stationId),"X position;x, cm; Entries", 1200, -60., 60.);
    fHM -> Create1<TH1F>(Form("h_trd_YPos_Station%i",stationId),"Y position;y, cm; Entries", 1200, -60., 60.);
   }
  fHM -> Create1<TH1F>("h_trd_XPos","X position;x, cm; Entries", 1200, -60., 60.);
  fHM -> Create1<TH1F>("h_trd_YPos","Y position;y, cm; Entries", 1200, -60., 60.);
}

void CbmTrdMCQa::Exec(Option_t*){
  ProcessPoints(fTrdPoints);
  fHM -> H1("h_trd_EventNo_MCQa") -> Fill(0.5);
}

void CbmTrdMCQa::ProcessPoints(const TClonesArray * points)
{
  
  fHM -> H1("h_trd_NofObjects_Points") -> Fill(points -> GetEntriesFast());
  
  Double_t pointX=0.;
  Double_t pointY=0.;
//  Double_t pointZ=0.;

  std::map<Int_t, vector<Int_t>> used_map;

  for(Int_t iPoint = 0; iPoint < points -> GetEntriesFast(); iPoint++) {
    const CbmTrdPoint* trdPoint = static_cast<const CbmTrdPoint*>(points -> At(iPoint));
    Int_t stationId = CbmTrdAddress::GetLayerId(trdPoint->GetModuleAddress());;
    fHM -> H1("h_trd_NofObjects_Points_Station") -> Fill(stationId);
    
    pointX = trdPoint -> GetXIn();
    pointY = trdPoint -> GetYIn();
//    pointZ = trdPoint -> GetZIn();

    fHM -> H1(Form("h_trd_XPos_Station%i",stationId)) -> Fill(pointX);
    fHM -> H1(Form("h_trd_YPos_Station%i",stationId)) -> Fill(pointY);

    fHM -> H2(Form("h_trd_PointsMap_Station%i", stationId)) -> Fill(pointX, pointY);
    fHM -> H2(Form("h_trd_PointsMapEvent_Station%i", stationId)) -> Fill(pointX, pointY);
    fHM -> H2(Form("h_trd_PointsMapRate_Station%i", stationId)) -> Fill(pointX, pointY);
    fHM -> H1("h_trd_XPos") -> Fill(pointX);
    fHM -> H1("h_trd_YPos") -> Fill(pointY);
    
    Int_t mcTrackID = trdPoint -> GetTrackID();
 
    if (std::find(used_map[stationId].begin(), used_map[stationId].end(), mcTrackID) == used_map[stationId].end()) {
      used_map[stationId].push_back(mcTrackID);
    }
  }
  fHM -> H1(Form("h_trd_MultPoints_Station%i",0)) -> Fill(used_map[0].size());
  fHM -> H1(Form("h_trd_MultPoints_Station%i",1)) -> Fill(used_map[1].size());
  fHM -> H1(Form("h_trd_MultPoints_Station%i",2)) -> Fill(used_map[2].size());
  fHM -> H1(Form("h_trd_MultPoints_Station%i",3)) -> Fill(used_map[3].size());

}

void CbmTrdMCQa::Finish(){

  Int_t nofEvents = fHM -> H1("h_trd_EventNo_MCQa") -> GetEntries();
  // Do here some scaling of the histograms to have MCPoint per cm^2

  Int_t xbins = (fHM->H2("h_trd_PointsMap_Station0"))->GetXaxis()->GetNbins();
  Float_t xmax = fHM -> H2("h_trd_PointsMap_Station0") -> GetXaxis()->GetXmax();
  Float_t xmin = fHM -> H2("h_trd_PointsMap_Station0") -> GetXaxis()->GetXmin();
  Float_t scaleX = static_cast<Float_t>(xbins)/(xmax - xmin);

  Int_t ybins = fHM -> H2("h_trd_PointsMap_Station0") -> GetYaxis()->GetNbins();
  Int_t ymax = fHM -> H2("h_trd_PointsMap_Station0") -> GetYaxis()->GetXmax();
  Int_t ymin = fHM -> H2("h_trd_PointsMap_Station0") -> GetYaxis()->GetXmin();
  Float_t scaleY = static_cast<Float_t>(ybins)/(ymax - ymin);

  Float_t scale = scaleX * scaleY;
  scale=1.;

  for(Int_t i=0; i<fNofStation; ++i) {
    fHM -> Scale(Form("h_trd_PointsMapEvent_Station%i", i),scale/nofEvents);
    fHM -> Scale(Form("h_trd_PointsMapRate_Station%i", i),10000000.*scale/nofEvents);
   }

  gDirectory -> mkdir("QA/TrdMCQa");
  gDirectory -> cd("QA/TrdMCQa");
  fHM -> WriteToFile();
  gDirectory -> cd("../..");
  //    CbmSimulationReport* report = new CbmTrdMCQaReport(fSetup, fDigitizer);
  //    report -> Create(fHM, fOutputDir);
  //    delete report;
  
  // Compare results with defined benchmark results and raise an ERROR if there are differences
  // Either get default histogramms from a benchmark qa file or as parameters from the parameter container
}

ClassImp(CbmTrdMCQa);