#include "PndMvdAllDataEventAna.h"
#include "PndMvdMCPoint.h"
#include "PndMvdPixel.h"
#include "TrackCand.h"
#include "TVector3.h"
#include "TGeoManager.h"
#include "TGeoBBox.h"
#include "TLegend.h"
#include "TGeoMatrix.h"
#include "PndMvdFileNameCreator.h"

#include "TGeoTrack.h"
#include <sstream>

ClassImp(PndMvdAllDataEventAna);

PndMvdAllDataEventAna::PndMvdAllDataEventAna(TString fileName):PndMvdEventAna(fileName)
{
  Init(fileName);
  SetCanvasColumns(6);
  

  geoH = new PndMvdGeoHandling(fileName);
  
  fRecoVolume = gGeoManager->MakeSphere("RecoHit",gGeoManager->GetMedium("vacuum"),0,0.1);
  fRecoVolume->SetLineColor(kRed);
  
  TGeoVolume* mcHit = gGeoManager->MakeSphere("MCHit",gGeoManager->GetMedium("vacuum"),0,0.1);
  mcHit->SetLineColor(kGreen);
  
  TGeoShape* stripShape = gGeoManager->GetVolume("StripSensor")->GetShape();
  TGeoVolume* stripVol = new TGeoVolume("StripSensorRed",stripShape, gGeoManager->GetMedium("vacuum"));
  stripVol->SetLineColor(kRed);
  gGeoManager->AddVolume(stripVol);
  stripVol = gGeoManager->GetVolume("StripSensorRed");
  
  TGeoShape* pixShape1 = gGeoManager->GetVolume("SensorActiveAreao[NrFE=8ooE]")->GetShape();
  TGeoVolume* pixVol1 = new TGeoVolume("PixelSensor8",pixShape1, gGeoManager->GetMedium("vacuum"));
  pixVol1->SetLineColor(kRed);
  gGeoManager->AddVolume(pixVol1);
  pixVol1 = gGeoManager->GetVolume("PixelSensor8");
  
  TGeoShape* pixShape2 = gGeoManager->GetVolume("SensorActiveAreao[NrFE=4]")->GetShape();
  TGeoVolume* pixVol2 = new TGeoVolume("PixelSensor4",pixShape2, gGeoManager->GetMedium("vacuum"));
  pixVol2->SetLineColor(kRed);
  gGeoManager->AddVolume(pixVol2);
  pixVol2 = gGeoManager->GetVolume("PixelSensor4");

  TGeoShape* pixShape3 = gGeoManager->GetVolume("SensorActiveAreao[NrFE=5]")->GetShape();
  TGeoVolume* pixVol3 = new TGeoVolume("PixelSensor5",pixShape3, gGeoManager->GetMedium("vacuum"));
  pixVol3->SetLineColor(kRed);
  gGeoManager->AddVolume(pixVol3);
  pixVol3 = gGeoManager->GetVolume("PixelSensor5");
  
  TGeoVolume* stripVolG = new TGeoVolume("StripSensorGreen",stripShape, gGeoManager->GetMedium("vacuum"));
  stripVolG->SetLineColor(kGreen);
  gGeoManager->AddVolume(stripVolG);
  stripVolG = gGeoManager->GetVolume("StripSensorGreen");
  
  TGeoVolume* pixVolG1 = new TGeoVolume("PixelSensor8Green",pixShape1, gGeoManager->GetMedium("vacuum"));
  pixVolG1->SetLineColor(kGreen);
  gGeoManager->AddVolume(pixVolG1);
  pixVolG1 = gGeoManager->GetVolume("PixelSensor8Green");
  
  TGeoVolume* pixVolG2 = new TGeoVolume("PixelSensor4Green",pixShape2, gGeoManager->GetMedium("vacuum"));
  pixVolG2->SetLineColor(kGreen);
  gGeoManager->AddVolume(pixVolG2);
  pixVolG2 = gGeoManager->GetVolume("PixelSensor4Green");

  TGeoVolume* pixVolG3 = new TGeoVolume("PixelSensor5Green",pixShape3, gGeoManager->GetMedium("vacuum"));
  pixVolG3->SetLineColor(kGreen);
  gGeoManager->AddVolume(pixVolG3);
  pixVolG3 = gGeoManager->GetVolume("PixelSensor5Green");
    
  fMvdTopVolume = gGeoManager->GetVolume("MVDoOption1");
  fGeoList = new PndEventDisplay();
  
  fGeoList->AddNewGroup("RecoHits",new PndGeoHitList("RecoHits","MVDoOption1", fRecoVolume));
  fGeoList->AddNewGroup("MCHits",new PndGeoHitList("MCHits","MVDoOption1", mcHit));
  fGeoList->AddNewGroup("StripHits", new PndGeoHitList("StripHits","MVDoOption1",stripVol));
  fGeoList->AddNewGroup("PixHits8", new PndGeoHitList("PixHits8","MVDoOption1",pixVol1));
  fGeoList->AddNewGroup("PixHits4", new PndGeoHitList("PixHits4","MVDoOption1",pixVol2));
  fGeoList->AddNewGroup("PixHits5", new PndGeoHitList("PixHits5","MVDoOption1",pixVol3));
  fGeoList->AddNewGroup("StripHitsGreen", new PndGeoHitList("StripHitsGreen","MVDoOption1",stripVolG));
  fGeoList->AddNewGroup("PixHits8Green", new PndGeoHitList("PixHits8Green","MVDoOption1",pixVolG1));
  fGeoList->AddNewGroup("PixHits4Green", new PndGeoHitList("PixHits4Green","MVDoOption1",pixVolG2));
  fGeoList->AddNewGroup("PixHits5Green", new PndGeoHitList("PixHits5Green","MVDoOption1",pixVolG3));
    
  fPixelCon = new PndMvdCalcFePixel(76,84,10);
}

void PndMvdAllDataEventAna::Init(TString fileName)
{
  PndMvdFileNameCreator nameCreator(fileName.Data());
  fFile = new TFile(fileName.Data());
  fTree = (TTree*)fFile->Get("cbmsim");
  
  fTree->AddFriend("digi=cbmsim",nameCreator.GetDigiFileName().c_str()); 
  fTree->AddFriend("reco=cbmsim",nameCreator.GetRecoFileName().c_str());
  fTree->AddFriend("trackF=cbmsim",nameCreator.GetTrackFindingFileName().c_str());
  
  fHitArray=new TClonesArray("PndMvdMCPoint");
  fDigiArray = new TClonesArray("PndMvdDigiPixel");
  fClusterArray = new TClonesArray("PndMvdClusterCand");
  fRecoArray = new TClonesArray("PndMvdHit");
  fGeoTrackArray = new TClonesArray("TGeoTrack");
  fTrackFArray = new TClonesArray("TrackCand");

  fTree->SetBranchAddress("MVDPoint",&fHitArray);
  fTree->SetBranchAddress("MVDPixelHit",&fDigiArray);
  fTree->SetBranchAddress("MVDClusterCand",&fClusterArray);
  fTree->SetBranchAddress("MVDClusterHit",&fRecoArray);
  fTree->SetBranchAddress("MVDTrackCand",&fTrackFArray);
  if (fTree->FindBranch("GeoTracks") != 0)
    fTree->SetBranchAddress("GeoTracks", &fGeoTrackArray);

  fAllHitPerClusterHistos = new TH1I("hitPerCluster","hitPerCluster",21,0,20);
  fAllHitResolutionHistos = new TH1D("AllHitResolution","AllHitResolution", 1000,0,0.1);
  f3DMCHisto = new TH3D("h3D","h3D", 50,-15,15, 50, -15,15, 50, -30,30);
  f3DRecoHisto = new TH3D("h3D","h3D", 50,-15,15, 50, -15,15, 50, -30,30);
  
}

void PndMvdAllDataEventAna::InitBranch()
{
}

void PndMvdAllDataEventAna::InitHistos()
{
}

void PndMvdAllDataEventAna::ClearHistos()
{
  ClearAllHMaps();
  ClearAllVectors();

  fGeoList->ClearHits("RecoHits");
  fGeoList->ClearHits("MCHits");
  fGeoList->ClearHits("StripHits");
  fGeoList->ClearHits("PixHits4");
  fGeoList->ClearHits("PixHits5");
  fGeoList->ClearHits("PixHits8");
  fGeoList->ClearHits("StripHitsGreen");
  fGeoList->ClearHits("PixHits4Green");
  fGeoList->ClearHits("PixHits5Green");
  fGeoList->ClearHits("PixHits8Green");
}

void PndMvdAllDataEventAna::AnaHits()
{
  fHitArray->Clear();
  fDigiArray->Clear();
  fClusterArray->Clear();
  fRecoArray->Clear();

  fTree->GetEntry(fActiveEvent);

  std::cout << "hitArray: " << fHitArray->GetEntries() << std::endl;
  std::cout << "digiArray: " << fDigiArray->GetEntries() << std::endl;
  std::cout << "clusterArray: " << fClusterArray->GetEntries() << std::endl;
  std::cout << "recoArray: " << fRecoArray->GetEntries() << std::endl;
    
  PrintHitArray();
  PrintDigiArray();
  PrintClusterArray();
  PrintRecoArray();

  FillHitHistos();
  FillDigiHistos();
  FillClusterHistos();
  FillRecoHistos();
  FillHitPerClusterHistos();
  FillHitResolutionHistos();
  Fill3DHisto();
  FillHitProjHistos();
  
  Create3DGeoHits();
}

void PndMvdAllDataEventAna::PrintHitArray()
{
  for (Int_t i = 0; i < fHitArray->GetEntries(); i++){
    std::cout << i << ": ";
    PndMvdMCPoint* myPoint = (PndMvdMCPoint*)fHitArray->At(i);
    std::cout << *myPoint;
  }
}
void PndMvdAllDataEventAna::PrintDigiArray()
{
  for (Int_t i = 0; i < fDigiArray->GetEntries(); i++){
    PndMvdDigiPixel* myDigi = (PndMvdDigiPixel*)fDigiArray->At(i);
    std::cout << i << ": ";
    //myDigi->Print();
    std::cout << *myDigi;
  }
}
void PndMvdAllDataEventAna::PrintClusterArray()
{
  for (Int_t i = 0; i < fClusterArray->GetEntries(); i++){
    std::cout << "Cluster " << i << ": " << std::endl;
    PndMvdClusterCand* myCluster = (PndMvdClusterCand*)fClusterArray->At(i);
    //myCluster->Print();
    std::vector<Int_t> myClusterList = myCluster->GetClusterList();
    for (Int_t j = 0; j < myClusterList.size(); j++){
      PndMvdDigiPixel* myDigi = (PndMvdDigiPixel*)fDigiArray->At(myClusterList[j]);
      std::cout << myClusterList[j] << ": ";
      std::cout << *myDigi;
    }
  }

}
void PndMvdAllDataEventAna::PrintRecoArray()
{
  for (Int_t i = 0; i < fRecoArray->GetEntries(); i++){
    std::cout << i << ": ";
    PndMvdHit* myHit = (PndMvdHit*)fRecoArray->At(i);
    std::cout << *myHit;
  }
}


TVector3 PndMvdAllDataEventAna::GetLocalHitPoints(TString detName, TVector3 input)
{
  TVector3 result;
  Double_t in[3];
  Double_t local[3];
  
  in[0] = input.X();
  in[1] = input.Y();
  in[2] = input.Z();

  gGeoManager->cd(geoH->GetPath(detName.Data()));
    TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix();

  transMat->MasterToLocal(in, local);
    
  TGeoVolume* actVolume = gGeoManager->GetCurrentVolume();
  TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape());

  result.SetX(local[0] + actBox->GetDX());
  result.SetY(local[1] + actBox->GetDY());
  result.SetZ(local[2] + actBox->GetDZ());

  return result;
}

void PndMvdAllDataEventAna::FillHitHistos()
{
  for (Int_t i = 0; i < fHitArray->GetEntries(); i++){
    PndMvdMCPoint *myPoint = (PndMvdMCPoint*)fHitArray->At(i);
    TString detName = myPoint->GetDetName();
  //if (hit->GetDetName().Contains("119_2")){
    if (fHistos[detName] == 0){
      fHistos[detName] = new TH2I("MCHisto",detName.Data(), 1000,0,1000,201,0,200);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    TH2* tempHisto = (TH2*)(fHistos[detName]);
    gGeoManager->cd(geoH->GetPath(detName.Data()));
      TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix();
    
    TVector3 in(myPoint->GetX(), myPoint->GetY(), myPoint->GetZ());
    TVector3 out(myPoint->GetXOut(), myPoint->GetYOut(), myPoint->GetZOut());
    TVector3 inLocal, outLocal;
    inLocal = GetLocalHitPoints(detName, in);
    outLocal = GetLocalHitPoints(detName, out);
  
    fHistos[detName]->Fill(inLocal.X()*100, inLocal.Y()*100);
    fHistos[detName]->Fill(outLocal.X()*100, outLocal.Y()*100);
    
//    std::cout << i << ": " << std::endl;
//     for (Int_t i = 0; i < 3; i++){
//          std::cout << "posInLocal "<< i << ": " << inLocal(i) << std::endl;
//            std::cout << "posOutLocal "<< i << ": " << outLocal(i) << std::endl;
//        }
  }
} 
void PndMvdAllDataEventAna::FillDigiHistos()
{
  for (Int_t i = 0; i < fDigiArray->GetEntries(); i++){
    PndMvdDigiPixel *hit = (PndMvdDigiPixel*)fDigiArray->At(i);
  //if (hit->GetDetName().Contains("119_2")){
    if (fDigiHistos[hit->GetDetName()] == 0){
      fDigiHistos[hit->GetDetName()] = new TH2I("digiHisto",hit->GetDetName().Data(), 1000,0,1000,201,0,200);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    TH2* tempHisto = (TH2*)(fDigiHistos[hit->GetDetName()]);
    PndMvdPixel myFePixel(hit->GetDetName().Data(), hit->GetFE(), hit->GetPixelColumn(), hit->GetPixelRow(), hit->GetCharge());
    PndMvdPixel mySensorPixel = fPixelCon->CalcSensorHit(myFePixel);
    tempHisto->Fill(mySensorPixel.GetCol(), mySensorPixel.GetRow(),mySensorPixel.GetCharge());//, (Double_t)(hit->GetCharge()));
  }
}
void PndMvdAllDataEventAna::FillClusterHistos()
{
  for (Int_t i = 0; i < fClusterArray->GetEntries(); i++){
    PndMvdClusterCand *cluster = (PndMvdClusterCand*)fClusterArray->At(i);
  //if (hit->GetDetName().Contains("119_2")){
    std::vector<Int_t> clusterList = cluster->GetClusterList();
    PndMvdDigiPixel* hit = (PndMvdDigiPixel*)fDigiArray->At(clusterList[0]);
    TString detName = hit->GetDetName();
    if (fClusterHistos[detName] == 0){
      fClusterHistos[detName] = new TH2I("clusterHisto",detName.Data(), 1000,0,1000,201,0,200);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    TH2* tempHisto = (TH2*)(fClusterHistos[detName]);
    for (Int_t j = 0; j < clusterList.size(); j++){
      hit = (PndMvdDigiPixel*)fDigiArray->At(clusterList[j]);
      PndMvdPixel myFePixel(detName.Data(), hit->GetFE(), hit->GetPixelColumn(), hit->GetPixelRow(), hit->GetCharge());
      PndMvdPixel mySensorPixel = fPixelCon->CalcSensorHit(myFePixel);
      tempHisto->Fill(mySensorPixel.GetCol(), mySensorPixel.GetRow());
    }
  }
}

void PndMvdAllDataEventAna::FillRecoHistos()
{
  for (Int_t i = 0; i < fRecoArray->GetEntries(); i++){
    PndMvdHit *myHit = (PndMvdHit*)fRecoArray->At(i);
    TString detName = myHit->GetDetName();
  //if (hit->GetDetName().Contains("119_2")){
    if (fRecoHistos[detName] == 0){
      fRecoHistos[detName] = new TH2I("RecoHisto",detName.Data(), 1000,0,1000,201,0,200);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    TH2* tempHisto = (TH2*)(fRecoHistos[detName]);
    
    TVector3 in(myHit->GetX(), myHit->GetY(), myHit->GetZ());
    TVector3 inLocal;
    inLocal = GetLocalHitPoints(detName, in);
  
    fRecoHistos[detName]->Fill(inLocal.X()*100, inLocal.Y()*100);
    
//    std::cout << i << ": " << std::endl;
//     for (Int_t i = 0; i < 3; i++){
//          std::cout << "posInLocal "<< i << ": " << inLocal(i) << std::endl;
//        }
  }
}

void PndMvdAllDataEventAna::FillHitPerClusterHistos()
{
  for (Int_t i = 0; i < fClusterArray->GetEntries(); i++){
    PndMvdClusterCand* cand = (PndMvdClusterCand*)fClusterArray->At(i);
    std::vector<Int_t> hits = GetHitPerCluster(cand);
    PndMvdMCPoint* point = (PndMvdMCPoint*)(fHitArray->At(hits[0]));
    TString detName = point->GetDetName();
    if (fHitPerClusterHistos[detName] == 0){
      fHitPerClusterHistos[detName] = new TH1I("HitPerCluster",detName.Data(), 21,0,20);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    fHitPerClusterHistos[detName]->Fill(hits.size());
    fAllHitPerClusterHistos->Fill(hits.size());
    std::cout << "Hits in Cluster: " << hits.size() << " for Det: " << detName.Data() << std::endl;
  }
}

void PndMvdAllDataEventAna::FillHitResolutionHistos()
{
  for (Int_t i = 0; i < fRecoArray->GetEntries(); i++){
    PndMvdHit *myHit = (PndMvdHit*)fRecoArray->At(i);
    TString detName = myHit->GetDetName();
    std::cout << "HitResolution for: " << detName << std::endl;
    TVector3 recoPos = myHit->GetPosition();
    std::cout << "RecoPos: " << recoPos.X() << " " << recoPos.Y() << " " << recoPos.Z() << std::endl;
    std::cout << "Hit RefIndex: " << myHit->GetRefIndex() << std::endl;
    PndMvdClusterCand *myCand = (PndMvdClusterCand*)(fClusterArray->At(myHit->GetRefIndex()));
    std::vector<Int_t> points = GetHitPerCluster(myCand);
    for (Int_t i = 0; i < points.size(); i++)
      std::cout << "ClusterPoints: " << points[i] << std::endl;
    TVector3 hitPos = CalcMeanHitPos(points);
    std::cout << "HitPos: " << hitPos.X() << " " << hitPos.Y() << " " << hitPos.Z() << std::endl;
    if (fHitResolutionHistos[detName] == 0){
      fHitResolutionHistos[detName] = new TH1D("HitRes",detName.Data(), 1000,0,0.01);
      //fDrawOption[hit->GetDetName()] = "colz";
    }
    TVector3 resPos = recoPos - hitPos;
    fHitResolutionHistos[detName]->Fill(resPos.Mag());
    fAllHitResolutionHistos->Fill(resPos.Mag());
    std::cout << "HitRes: " << resPos.Mag() << std::endl;
  }
}

void PndMvdAllDataEventAna::Fill3DHisto()
{
  for (Int_t i = 0; i < fRecoArray->GetEntries(); i++){
    PndMvdHit* myHit = (PndMvdHit*)fRecoArray->At(i);
    f3DRecoHisto->Fill(myHit->GetX(), myHit->GetY(), myHit->GetZ());
  }
  for (Int_t j = 0; j < fHitArray->GetEntries(); j++){
      PndMvdMCPoint* myPoint = (PndMvdMCPoint*)fHitArray->At(j);
      f3DMCHisto->Fill(myPoint->GetX(), myPoint->GetY(), myPoint->GetZ());
  }
    
}

void PndMvdAllDataEventAna::FillHitProjHistos()
{
  ClearHistoVector(&fRecoHisxy);
  ClearHistoVector(&fRecoHisrz);
  unsigned int detID, hitID;
  TVector3 vec;
    
  for (Int_t i = 0; i < fTrackFArray->GetEntries(); i++){
    TrackCand* trackC = (TrackCand*)fTrackFArray->At(i);
    TH2D* myHistoXY = new TH2D("recohisxy","MVD Reco Points, xy view",400,-15.,15.,400,-15.,15.);
    myHistoXY->SetMarkerStyle(7);
    myHistoXY->SetMarkerColor(i+1);
    
    TH2D* myHistoRZ = new TH2D("recohisrz","MVD Reco Points, rz view",400,-20.,20.,400,-20.,20.);
    myHistoRZ->SetMarkerStyle(7);
    myHistoRZ->SetMarkerColor(i+1);
            
    fRecoHisxy.push_back(myHistoXY);
    fRecoHisrz.push_back(myHistoRZ);
    for (Int_t j = 0; j < trackC->getNHits(); j++){
      trackC->getHit(j, detID, hitID);
      PndMvdHit* myHit = (PndMvdHit*)fRecoArray->At(hitID);
      vec.SetXYZ(myHit->GetX(), myHit->GetY(), myHit->GetZ());
      myHistoXY->Fill(vec.x(), vec.y());
      myHistoRZ->Fill(vec.z(), vec.Perp());
    }
  }
}


TVector3 PndMvdAllDataEventAna::CalcMeanHitPos(std::vector<Int_t> points)
{
  TVector3 result;
  Double_t energy = 0;
  for (Int_t i = 0; i < points.size(); i++){
    PndMvdMCPoint *myPoint = (PndMvdMCPoint*)fHitArray->At(points[i]);
    TVector3 posIn = myPoint->GetPosition();
    std::cout << "posIn: " << posIn.X() << " " << posIn.Y() << " " << posIn.Z() << std::endl;
    TVector3 posOut = myPoint->GetPositionOut();
    std::cout << "posOut: " << posOut.X() << " " << posOut.Y() << " " << posOut.Z() << std::endl;

    TVector3 posHit = posIn + posOut;
    posHit *= 0.5;
    
    std::cout << "posHit: " << posHit.X() << " " << posHit.Y() << " " << posHit.Z() << std::endl;
    //posHit *= myPoint->GetEnergyLoss();
    energy += myPoint->GetEnergyLoss();
    result += posHit;
  }
  //Double_t div = 1/energy;
  Double_t div = 1/(Double_t)points.size();
  result *= div;
  return result;
}

void PndMvdAllDataEventAna::DrawHitHisto(TString detName, TCanvas* extCan, Int_t pad)
{
  if (extCan)
    extCan->cd(pad);
  
  if (fHistos[detName] != 0)
    fHistos[detName]->Draw();
}

void PndMvdAllDataEventAna::DrawDigiHisto(TString detName, TCanvas* extCan, Int_t pad)
{
  if (extCan)
    extCan->cd(pad);

  if(fDigiHistos[detName] != 0)
    fDigiHistos[detName]->Draw("colz");
}
void PndMvdAllDataEventAna::DrawClusterHisto(TString detName, TCanvas* extCan, Int_t pad)
{
  if (extCan)
      extCan->cd(pad);
    
  if(fClusterHistos[detName] != 0)
    fClusterHistos[detName]->Draw("colz");
}
void PndMvdAllDataEventAna::DrawRecoHisto(TString detName, TCanvas* extCan, Int_t pad)
{
  if (extCan)
      extCan->cd(pad);
    
  if (fRecoHistos[detName] != 0)
    fRecoHistos[detName]->Draw();
}

void PndMvdAllDataEventAna::DrawAllHistos(TString detName, TCanvas* extCan)
{
  if (extCan)
      fCan1 = extCan;
  if (fCan1 == 0)
    fCan1 = new TCanvas();
  DrawHitHisto(detName, fCan1, 1);
  DrawDigiHisto(detName, fCan1, 2);
  DrawClusterHisto(detName, fCan1, 3);
  DrawRecoHisto(detName, fCan1, 4);
  fCan1->cd(5);
  fHitPerClusterHistos[detName]->Draw();
  fCan1->cd(6);
  fHitResolutionHistos[detName]->Draw();
  fCan1->Modified(kTRUE);
  fCan1->Update();
}

void PndMvdAllDataEventAna::DrawAllHistos(Int_t index, TCanvas* extCan)
{
  std::map<TString,TH1*>::const_iterator ki;
  Int_t histoSize = fHistos.size();
  
  if (index >= histoSize) return;
  ki = fHistos.begin();
  for (Int_t i = 0; i < index; i++)
    ki++;
  TString detName = ki->first;
  std::cout << "DetName: " << detName.Data() << std::endl;
  DrawAllHistos(detName, extCan);
}

void PndMvdAllDataEventAna::Draw3D(TString opt, TCanvas* extCan, Int_t pad)
{
  if (extCan == 0)
      extCan = new TCanvas();
  extCan->cd(pad);
  f3DMCHisto->SetMarkerStyle(8);
  f3DMCHisto->SetMarkerSize(0.5);
  f3DMCHisto->SetMarkerColor(2);
  f3DMCHisto->Draw(opt.Data());
  f3DRecoHisto->SetMarkerStyle(8);
  f3DRecoHisto->SetMarkerSize(0.5);
  f3DRecoHisto->SetMarkerColor(1);
  f3DRecoHisto->Draw("same");

}

void PndMvdAllDataEventAna::DrawAllTracks(TCanvas* extCan, Int_t pad)
{
  if (extCan == 0)
    extCan = new TCanvas();
  extCan->cd(pad);
  
  for (Int_t i = 0; i < fGeoTrackArray->GetEntries(); i++){
    TGeoTrack* myTrack = (TGeoTrack*) fGeoTrackArray->At(i);
    myTrack->Print();
    myTrack->Draw();
  }
}

void PndMvdAllDataEventAna::DrawHitTracks(TCanvas* extCan, Int_t pad)
{
  if (extCan == 0)
    extCan = new TCanvas();
  extCan->cd(pad);
  
  std::map<Int_t, Int_t> trackUsed;
  for (Int_t i = 0; i < fHitArray->GetEntries(); i++){
    PndMvdMCPoint* myPoint = (PndMvdMCPoint*)fHitArray->At(i);
    std::cout << "Hit " << i << " TrackID " << myPoint->GetTrackID() << std::endl;
    if (trackUsed[myPoint->GetTrackID()] == 0){
      for (Int_t j = 0; j < fGeoTrackArray->GetEntries(); j++){
        TGeoTrack* myTrack = (TGeoTrack*)fGeoTrackArray->At(j);
        if (myPoint->GetTrackID() == myTrack->GetId()){
          myTrack->Draw();
          myTrack->Print();
        }
      }
      trackUsed[myPoint->GetTrackID()] = 1;
    }
  }
}

void PndMvdAllDataEventAna::DrawTopVolume(TCanvas* extCan, Int_t pad, const char* opt)
{
  if (extCan)
    extCan->cd(pad);
//  fMvdTopVolume->Draw(""); //ogl for OpenGL Viewer
  fMvdTopVolume->Draw(opt);
}

void PndMvdAllDataEventAna::DrawHisrz(TCanvas* extCan, Int_t pad)
{
  DrawHistoVec(&fRecoHisxy, extCan, pad);
}

void PndMvdAllDataEventAna::DrawHisxy(TCanvas* extCan, Int_t pad)
{
  DrawHistoVec(&fRecoHisrz, extCan, pad);
}

void PndMvdAllDataEventAna::DrawHistoVec(std::vector<TH1*> *vec, TCanvas* extCan, Int_t pad) const
{
  if (extCan)
    extCan->cd(pad);
  if (vec->size() > 0)
      (*vec)[0]->Draw();
    for (Int_t i = 1; i < vec->size(); i++)
      (*vec)[i]->Draw("same");
}

void PndMvdAllDataEventAna::DrawHisto(TH1* histo, TCanvas* extCan, Int_t pad) const
{
  if (extCan)
    extCan->cd(pad);
  histo->Draw();
}

void PndMvdAllDataEventAna::DrawEvent(bool tracks, TCanvas* extCan)
{
  if (extCan)
        fCan2 = extCan;
  if (fCan2 == 0)
    fCan2 = new TCanvas();
  
  fCan2->Divide(2,2);
  
  fCan2->cd(1);
  DrawTopVolume();

  if (tracks)
    DrawHitTracks();
  
  fCan2->cd(2);
  fAllHitResolutionHistos->Draw();
  
  fCan2->cd(3);
  DrawHistoVec(&fRecoHisxy);
  
  fCan2->cd(4);
  DrawHistoVec(&fRecoHisrz);
}

std::vector<Int_t> PndMvdAllDataEventAna::GetHitPerCluster(PndMvdClusterCand* clusterCand)
{
  std::vector<Int_t> result;
  std::vector<Int_t> digiPos = clusterCand->GetClusterList();
  bool isInResult = false;
  for (Int_t i = 0; i < digiPos.size(); i++){
    PndMvdDigiPixel* digiHit = (PndMvdDigiPixel*)(fDigiArray->At(digiPos[i]));
    Int_t mcID = digiHit->GetIndex();
    std::cout << " -I- GetHitPerCluster: mcID: " << mcID << std::endl;
      for (Int_t j = 0; j < result.size() && isInResult == false; j++){
        //std::cout << "Result: " << result[j] << std::endl;
        if (mcID == result[j])
          isInResult = true;
      }
    if (isInResult == false)
      result.push_back(mcID);
    isInResult = false;
  }
  for (Int_t k = 0; k < result.size(); k++)
    std::cout << " Result: " << k << ": " << result[k] << std::endl;
  return result;
}

void PndMvdAllDataEventAna::Create3DGeoHits()
{
  fGeoList->SetHits("RecoHits",fRecoArray);
  
  for (Int_t i = 0; i < fHitArray->GetEntriesFast(); i++){
    PndMvdMCPoint* myPoint = (PndMvdMCPoint*)fHitArray->At(i);
    fGeoList->AddHit("MCHits", myPoint->GetX(), myPoint->GetY(), myPoint->GetZ());
    std::cout << "gGeoManager cd: " << myPoint->GetDetName() << " " << gGeoManager->cd(geoH->GetPath(myPoint->GetDetName())) << std::endl;
    TGeoHMatrix* mat = gGeoManager->GetCurrentMatrix();
    mat->Print();
    if (myPoint->GetDetName().Contains("71")){
      std::cout << "StripHits added!" << std::endl;
      fGeoList->AddHit("StripHits", mat);
      fGeoList->AddHit("StripHitsGreen", mat, kFALSE);
    }
    else if (myPoint->GetDetName().Contains("89")){
      std::cout << "PixHits4 added!" << std::endl;
      fGeoList->AddHit("PixHits4", mat);
      fGeoList->AddHit("PixHits4Green", mat, kFALSE);
    }
    else if (myPoint->GetDetName().Contains("56")){
      std::cout << "PixHits8 added!" << std::endl;
      fGeoList->AddHit("PixHits8", mat);
      fGeoList->AddHit("PixHits8Green", mat, kFALSE);
    }
  }
}


void PndMvdAllDataEventAna::ClearHistoMaps(std::map<TString, TH1*>* myMaps) const 
{
  for (std::map<TString,TH1*>::const_iterator ki = myMaps->begin(); ki != myMaps->end(); ki++){
    ki->second->Delete();
    (*myMaps)[ki->first] = 0;
  }
  myMaps->clear();
  std::cout << " Maps.size: " << myMaps->size() << std::endl;
}

void PndMvdAllDataEventAna::ClearAllHMaps()
{
  ClearHistoMaps(&fHistos);
  ClearHistoMaps(&fRecoHistos);
  ClearHistoMaps(&fDigiHistos);
  ClearHistoMaps(&fClusterHistos);
  ClearHistoMaps(&fHitPerClusterHistos);
  ClearHistoMaps(&fHitResolutionHistos);
}

void PndMvdAllDataEventAna::ClearAllVectors()
{
  ClearHistoVector(&fHitHistoVec);
  ClearHistoVector(&fDigiHistoVec);
  ClearHistoVector(&fClusterHistoVec);
  ClearHistoVector(&fRecoHistoVec); 
}
void PndMvdAllDataEventAna::ClearHistoVector(std::vector<TH1*>* myVectors) const 
{
  for (Int_t i = 0; i < myVectors->size(); i++)
  {
    delete ((*myVectors)[i]);
  }
  myVectors->clear();
}

std::vector<TString> PndMvdAllDataEventAna::GetModulesHit()
{
    std::vector<TString> result;
    for(std::map<TString, TH1*>::const_iterator ki= fRecoHistos.begin(); ki != fRecoHistos.end(); ++ki){
      std::cout << "First: " << ki->first << std::endl;
      std::cout << "Path: " << geoH->GetPath(ki->first) << std::endl;
      result.push_back(geoH->GetPath(ki->first));
    }
    return result;
}