// -------------------------------------------------------------------------
// -----                   PndDrcRecoLookupMap source file                -----
// -----               Created 10/11/10  by Maria Patsyuk              -----
// -----                                                               -----
// -----                                                               -----
// -------------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include "stdio.h"

#include "PndGeoDrc.h"
#include "PndDrcRecoLookupMap.h"
#include "FairRootManager.h"
#include "PndMCTrack.h"
#include "PndDrcBarPoint.h"
#include "PndDrcPDPoint.h"
#include "PndDrcHit.h"
#include "PndDrcPDHit.h"
#include "TVector3.h"
#include "TRandom.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairBaseParSet.h"
#include "FairGeoVolume.h"
#include "TString.h"
#include "FairGeoTransform.h"
#include "FairGeoVector.h"
#include "FairGeoMedium.h"
#include "FairGeoNode.h"

#include "PndGeoDrcPar.h"
#include "TFormula.h"
#include "TMath.h"
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "TPDGCode.h"
#include "TGeoManager.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TF1.h"
#include "TExec.h"
#include "TLine.h"
#include "TPolyLine.h"
//#include "PndChPho.h"
#include "PndDrcLutInfo.h"
#include "PndDrcDigiPar.h"
#include "TVectorD.h"

using std::endl;
using std::cout;

// -----   Default constructor   -------------------------------------------
PndDrcRecoLookupMap::PndDrcRecoLookupMap() 
:FairTask("PndDrcRecoLookupMap")
{
  fGeo = new PndGeoDrc();
  fGeoH=NULL;
  fDigiPar = NULL;  
}
// -----   Standard constructor with verbosity level  -------------------------------------------

PndDrcRecoLookupMap::PndDrcRecoLookupMap(Int_t verbose) 
  :FairTask("PndDrcRecoLookupMap")
{
  fVerbose = verbose;  
  fGeo = new PndGeoDrc();
  fGeoH=NULL;
  fDigiPar = NULL;
}
// -----   Destructor   ----------------------------------------------------
PndDrcRecoLookupMap::~PndDrcRecoLookupMap()
{
  if (fGeo) delete fGeo;
  if (fGeoH) delete fGeoH;
  fHistoList->Delete();  
  delete fHistoList;      
}

// -----   Initialization   -----------------------------------------------
InitStatus PndDrcRecoLookupMap::Init()
{
  cout << " ---------- INITIALIZATION ------------" << endl;
  nevents = 0;
  // Get RootManager
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) {
    cout << "-E- PndDrcRecoLookupMap::Init: "
         << "RootManager not instantiated!" << endl;
    return kFATAL;
  }

  // Get input array
  fMCArray = (TClonesArray*) ioman->GetObject("MCTrack");
  if ( ! fMCArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: "
         << "No MCTrack array!" << endl;
    return kERROR;
  }
  // Get input array
  fBarPointArray = (TClonesArray*) ioman->GetObject("DrcBarPoint");
  if ( ! fBarPointArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: "
         << "No DrcBarPoint array!" << endl;
    return kERROR;
  }
  // Get Photon point array
  fPDPointArray = (TClonesArray*) ioman->GetObject("DrcPDPoint");
  if ( ! fPDPointArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: "
         << "No DrcPDPoint array!" << endl;
    return kERROR;
  }
/*  // Get input array
  fHitArray = (TClonesArray*) ioman->GetObject("DrcHit");
  if ( ! fHitArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: "
         << "No DrcHit array!" << endl;
	 return kERROR;
	 }
*/

  // Get digi array
  fDigiArray = (TClonesArray*) ioman->GetObject("DrcDigi");
  if ( ! fDigiArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: " << "No DrcDigi array!" << endl;
    return kERROR;
  }

  // Get input array
  fPDHitArray = (TClonesArray*) ioman->GetObject("DrcPDHit");
  if ( ! fPDHitArray ) {
    cout << "-W- PndDrcRecoLookupMap::Init: "
         << "No DrcPDHit array!" << endl;
    return kERROR;
  }
    
   // Create and register output array
  // fChPhoArray = new TClonesArray("PndChPho");
  // ioman->Register("ChPho","Drc",fChPhoArray, kTRUE);  
  fDrcLutInfoArray = new TClonesArray("PndDrcLutInfo");
  ioman->Register("DrcLutInfo","Drc",fDrcLutInfoArray, kTRUE); 
  
   // Get parameters:
   fpi = TMath::Pi();
   
   fR          = fGeo->radius();
   fzup        = fGeo->barBoxZUp();
   fzdown      = fGeo->barBoxZDown();
   fHThick     = fGeo->barHalfThick();
   fBboxNum    = fGeo->BBoxNum();
   fBarNum     = fGeo->barNum();
   fPipehAngle = fGeo->PipehAngle();
   fBarBoxGap  = fGeo->BBoxGap();
   //fLength     = (180. - 2.*fPipehAngle - fBarBoxGap/fR*(fBboxNum/2. - 1.)/fpi*180.)/(fBboxNum/2.) * fR/ 180.*fpi;
   fDphi       = 2.*(180. - 2*fPipehAngle)/ fBboxNum; //[degrees]
   fLSide      = fGeo->Lside();//(180. - 2.*fPipehAngle - fBarBoxGap/fR*(fBboxNum/2. - 1.)/fpi*180.)/(fBboxNum/2.) * fR/ 180.*fpi;
   fBarWidth   = fGeo->BarWidth();//fLSide/fBarNum;
   
   // vectors of bar surfaces:
   fnX1.SetXYZ(-1., 0.,0.);   
   fnY1.SetXYZ( 0., 1.,0.);   
   
   // Time resolution as a function of photon path
   timeres = new TF1("timeres","[0] + [1]*x", 0., 1000.);
   timeres->SetParameters(23.2773, 1.11008); 

   cout << "-I- PndDrcRecoLookupMap: Intialization successfull" << endl;
   CreateHisto();
   return kSUCCESS;
   
}

// -----   Private method SetParContainers   -------------------------------
void PndDrcRecoLookupMap::SetParContainers() {

  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get DIRC digitisation parameter container
  fDigiPar = (PndDrcDigiPar*)(db->getContainer("DIRCLookupTable"));
  cout << "-I- PndDrcRecoLookupMap::SetParContainers(). read parameters" << endl;
  cout<<"read them!"<<endl;
  //cout << "-I- PndDrcRecoLookupMap: Number of hit pixels "<< fDigiPar->GetNHitPixels()<< endl;
  cout << "-I- PndDrcRecoLookupMap: Number of hit pixels "<< fDigiPar->GetNAmbiguities()<< endl;
  
  if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance();
  fGeoH->SetParContainers(); 
  if(fVerbose>1) Info("SetParContainers","done.");
  return;
}
// -------------------------------------------------------------------------

// -----   Execution of Task   ---------------------------------------------
void PndDrcRecoLookupMap::Exec(Option_t* option)
{
  //if ( ! fChPhoArray ) Fatal("Exec", "No fChPhoArray");
  //fChPhoArray->Clear();
  fGeoH->SetVerbose(fVerbose);
  
  nevents++;
  fDetectorID = 0;
   
  cout<<"EVENT # "<<nevents<<endl;
  // fNHits = 0;
  //ProcessBarHit();
  ProcessPhotonHit();
  //ProcessPhotonMC();
}

//--------------Process Photon MC Points----------------------------------------------------
void PndDrcRecoLookupMap::ProcessPhotonMC()
{
  if (fVerbose > 0) {
    cout <<"PndDrcRecoLookupMap: Number of Detected MC Tracks : "<<fMCArray->GetEntries()<<endl;
    PndMCTrack* ptr = NULL;
    PndMCTrack* trMr=NULL;
    
  }
}
//--------------Process Bar Hits----------------------------------------------------
void PndDrcRecoLookupMap::ProcessBarHit()
{
 PndDrcBarPoint* pt=NULL;
  PndDrcHit* hit=NULL;
  PndMCTrack* tr = NULL;
  
  for(Int_t j=0; j<fHitArray->GetEntriesFast(); j++) {
    hit = (PndDrcHit*)fHitArray->At(j);
    Int_t mcRef= hit->GetRefIndex();
    pt= (PndDrcBarPoint*)fBarPointArray->At(mcRef);
    
    Int_t chtrID= pt->GetTrackID();
    tr = (PndMCTrack*)fMCArray->At(chtrID);
    
    cout<<"nother track px = "<<pt->GetPx()<<", py = "<<pt->GetPy()<<endl;
    
  }
}

//--------------Process Photon Hits----------------------------------------------------
void PndDrcRecoLookupMap::ProcessPhotonHit()
{
  //Int_t nofChPho = 0;
  PndDrcLutInfo lutinfo;
  
  PndDrcPDPoint* Ppt=NULL;
  PndDrcPDHit* pdhit=NULL;
  PndMCTrack* tr = NULL;
  PndMCTrack* trMr;
  	
  if (fVerbose > 0) {
    cout <<"PndDrcRecoLookupMap: Number of Detected Photon MCPDPoints : "<<fPDPointArray->GetEntries()<<endl;
    cout <<"PndDrcRecoLookupMap: Number of Detected Photon PDHits : "<<fPDHitArray->GetEntries()<<endl;
  }
  
 // Loop over PndDrcPDHits
 for(Int_t k=0; k<fPDHitArray->GetEntriesFast(); k++) {
      
    pdhit = (PndDrcPDHit*)fPDHitArray->At(k);
        
    Int_t mcPDRef= pdhit->GetRefIndex();
    
    Ppt = (PndDrcPDPoint*)fPDPointArray->At(mcPDRef);
    
    PndDrcBarPoint *fBarPoint= (PndDrcBarPoint*)fBarPointArray->At(Ppt->GetBarPointID());
    fBarId = fBarPoint->GetDetectorID();
    //cout<<"bar name - "<<fGeoH->GetPath(fBarId)<<endl;
            
    Int_t trID= Ppt->GetTrackID();
    tr = (PndMCTrack*)fMCArray->At(trID);
    
    Int_t trMID=tr->GetMotherID();
    
    //if(trMID > -1){
    trMr = (PndMCTrack*)fMCArray->At(trMID);
    trMr->GetMomentum().Print();
    
    Int_t trMpdg=trMr->GetPdgCode();
      
    if(trMr->GetMotherID()==-1){ // charged track is a primary
    
      lutinfo.SetChPartPdg(trMpdg);
    
      // expected cherenkov angle
      Double_t Mrmass;
      Int_t MoSign;
      if(fabs(trMpdg) == 11){Mrmass = 0.000511; MoSign = 1;}
      if(fabs(trMpdg) == 13){Mrmass = 0.105658;  MoSign = 1;}
      if(fabs(trMpdg) == 211){Mrmass = 0.139570; MoSign = -1;}
      if(fabs(trMpdg) == 321){Mrmass = 0.49368; MoSign = -1;}
      if(fabs(trMpdg) == 2212){Mrmass = 0.938; MoSign = -1;}
      if(fabs(trMpdg) == 50000050){Mrmass = 0.0; MoSign = -1;}
      Double_t Mrmom = trMr->GetMomentum().Mag();      
      CHexp = acos(sqrt(pow(Mrmom,2) + pow(Mrmass,2))/Mrmom/fGeo->nQuartz());
      lutinfo.SetCherenkovMC(CHexp);
      cout<<"+++++++++++++++++++++++++++"<<endl;
      cout<<"CH expected = "<<CHexp<<endl;
    
      fxPHit= pdhit->GetX();
      fyPHit= pdhit->GetY();
      fzPHit= pdhit->GetZ();
      ftime = pdhit->GetTime();      
           
      //get the pixel id
      fpixID = pdhit->GetDetectorID(); 
      cout<<"hit id "<<fpixID<<", x = "<<pdhit->GetX()<<", y = "<<pdhit->GetY()<<", time = "<<ftime<<endl;      
      
      Double_t xPPoi= Ppt->GetX();
      Double_t yPPoi= Ppt->GetY();
      
      lutinfo.SetPath(Ppt->GetLength());     
	   
      // momentum of a photon on the PD Plane     
      fPphoPD.SetXYZ(Ppt->GetPx(), Ppt->GetPy(), Ppt->GetPz());
      //cout<<"Photon momentum on the PDplane: "<<endl;
      //fPphoPD.Print();
      
      {     
	// mother momentum/direction      
	fPMo.SetXYZ(trMr->GetMomentum().X(), trMr->GetMomentum().Y(), trMr->GetMomentum().Z());	
	//fPMo.Print();
	lutinfo.SetChPartDir(fPMo);
	cout<<"mother phi = "<<fPMo.Phi()/3.1415*180.<<endl;            
	if(fB > 0.){       
	  Double_t PtMo = sqrt(pow(trMr->GetMomentum().X(),2)+pow(trMr->GetMomentum().Y(),2));
	  Double_t Rratio = fR*fR/2./pow((PtMo/0.29979/fB)*100.,2);
	  //Double_t phi_extra = fPMo.Phi() + trMpdg/fabs(trMpdg) * acos(1. - Rratio); // muon - ; proton +
	  Double_t phi_extra = fPMo.Phi() + MoSign * trMpdg/fabs(trMpdg) * acos(1. - Rratio); // muon - ; proton +
	  //fHAngleInBDeg = 0.5 * trMpdg/fabs(trMpdg) * (acos(1. - Rratio) /TMath::Pi())*180.;   
	  fHAngleInBDeg = 0.5 * MoSign * trMpdg/fabs(trMpdg) * (acos(1. - Rratio) /TMath::Pi())*180.;
	     // cout<<"B = "<<fB<<", R = "<<fR<<", trMpdg = "<<trMpdg<<", half angle in B = "<<fHAngleInBDeg<<endl;
	  fPMo.SetPhi(phi_extra);
	}	
	lutinfo.SetChPartDirInBar2(fPMo);
      }

      {
	fBarPoint->Momentum(fPMo);
	lutinfo.SetChPartDirInBar(fPMo);	
      }


      cout<<"Mother vector (global cs) : "<<endl;
      TVector3 motherMom = fPMo.Unit();
      motherMom.Print();
         
      // initial momentum of the photon      
      fPphoInit.SetXYZ(tr->GetMomentum().X(), tr->GetMomentum().Y(), tr->GetMomentum().Z());
      cout<<"Initial momentum of the photon :"<<endl;
      fPphoInit.Print();
                              
      // real = generated cherenkov angle
      CHreal = fPphoInit.Angle(fPMo);
      lutinfo.SetCherenkovReal(CHreal);
      //cout<<"CH real (generated) = "<<CHreal<<endl;
    
      Double_t etot = sqrt(pow(fPphoPD.X(),2) + pow(fPphoPD.Y(),2) +pow(fPphoPD.Z(),2));      
      flambdah=197.0*2.0*TMath::Pi()/(etot*1.0E9);//wavelength of photon in nm
      lutinfo.SetLambda(flambdah);
      
      // production point of the photon
      fStartVertex = tr->GetStartVertex();
      
      //cout<<"start: phi = "<<fStartVertex.Phi()<<endl;
      //fStartVertex.Print();
      ftime0 = tr->GetStartTime()/1.0E9;
      lutinfo.SetTime(ftime-ftime0);
      
      //cout<<"fstart time = "<<ftime0<<endl;
      //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      // proceed only if photon was born in the bar for which the LUT is used!!!
      //if(fStartVertex.Phi() >0./180.*fpi && fStartVertex.Phi() < 6./180.*fpi){
      if(fPphoInit.Z() > 0. || fPphoInit.Z() < 0.){ // reflected photons
       
      // Transformation of the photon initial direction to the local bar coodr system:      
      // to get to the bar coord system.
      fPhiRot = InBarCoordSyst(fStartVertex, &fBBver1, &fBBver2, &fBBver3, &fBBver4);
      cout<<"PhiRot = "<<fPhiRot/fpi*180.<<", Mother phi = "<<fPMo.Phi()/fpi*180.<<endl;
      
      // Photon initial momentum in the bar coord. syst.
      fPphoB = fPphoInit.Unit();
      fPphoB.RotateZ(-fPhiRot);
      fkxBar = -fPphoB.Y();
      fkyBar = fPphoB.X();
      fkzBar = fPphoB.Z();
      fPphoB.SetXYZ(fkxBar, fkyBar, fkzBar); 
      cout<<"photon momentum in the bar = "<<endl;
      fPphoB.Print();
      TVector3 fPpB = (fGeoH->MasterToLocalShortId(fPphoInit, fBarId) - fGeoH->MasterToLocalShortId((0.,0.,0.),fBarId)).Unit();
      fPpB.Print();
          
      // mother momentum in the bar' (and bar) coord syst:
      fPMoB = fPMo.Unit();
      fPMoB.RotateZ(-fPhiRot); 
      cout<<"Mother vector (bar' cs): "<<endl;      
      fPMoB.Print();
      fPxMoBar = -fPMoB.Y();
      fPyMoBar = fPMoB.X();
      fPzMoBar = fPMoB.Z();
      TVector3 PMoBar;
      PMoBar.SetXYZ(fPxMoBar, fPyMoBar, fPzMoBar);
      TVector3 PMBar = (fGeoH->MasterToLocalShortId(fPMo, fBarId) - fGeoH->MasterToLocalShortId((0.,0.,0.),fBarId)).Unit();       
      PMBar.Print();     
	   
      // to find Cherenkov Phi (NOT THETA!!!):
      TVector3 kB;
      kB.SetXYZ(fkxBar, fkyBar, fkzBar);
            
      // use lookups to get the kXbar and kYbar of photons:
      NHitPix = fDigiPar->GetNHitPixels();
      //cout<<"N pixels = "<<NHitPix<<endl;
      NAmb = fDigiPar->GetNAmbiguities();
      cout<<"N amb = "<<NAmb<<endl;
      NPixPar = fDigiPar->GetNPixelParam();
      Double_t par[NPixPar]; 
           
      if (fDigiPar->GetParamsForPixel(fpixID, par) == kFALSE) {
        cout<<"NO SUCH HIT IN LUT!!!"<<endl;
        continue;
      }
      
      cout<<"LOOKUP!!!!"<<endl;      
      //cout<<"N Hit Pix = "<<NHitPix<<endl;
      //cout<<"N Amb = "<<NAmb<<endl;
      cout<<"NPixPar = "<<NPixPar<<endl;
/*      cout<<"par 0 = "<<par[0]<<endl; // kx 
      cout<<"par 1 = "<<par[1]<<endl; // ky 
      cout<<"par 2 = "<<par[2]<<endl; // kxB 
      cout<<"par 3 = "<<par[3]<<endl; // kyB 
      cout<<"par 4 = "<<par[4]<<endl; // kxU1
      cout<<"par 5 = "<<par[5]<<endl; // kyU1
      cout<<"par 6 = "<<par[6]<<endl; // kxU2
      cout<<"par 7 = "<<par[7]<<endl; // kyU2
      cout<<"par 8 = "<<par[8]<<endl; // kxU3
      cout<<"par 9 = "<<par[9]<<endl; // kyU3
      cout<<"par 10 = "<<par[10]<<endl; // kxBU1
      cout<<"par 11 = "<<par[11]<<endl; // kyBU1
      cout<<"par 12 = "<<par[12]<<endl; // kxBU2
      cout<<"par 13 = "<<par[13]<<endl; // kyBU2
      cout<<"par 14 = "<<par[14]<<endl; // kxBU3
      cout<<"par 15 = "<<par[15]<<endl; // kyBU3
      cout<<"par 16 = "<<par[16]<<endl; // kxUB 
      cout<<"par 17 = "<<par[17]<<endl; // kyUB 
      cout<<"par 18= "<<par[18]<<endl; // kxUU 
      cout<<"par 19= "<<par[19]<<endl; // kyUU 
      cout<<"par 20= "<<par[20]<<endl; // kxUUU
      cout<<"par 21= "<<par[21]<<endl; // kyUUU
      cout<<"par 22= "<<par[22]<<endl; // kxBUU
      cout<<"par 23= "<<par[23]<<endl; // kyBUU
      cout<<"par 24= "<<par[24]<<endl; // kxBUB 
      cout<<"par 25= "<<par[25]<<endl; // kyBUB 
      cout<<"par 26= "<<par[26]<<endl; // kxUBU
      cout<<"par 27= "<<par[27]<<endl; // kyUBU         
*/     
      // fill in NPixPar/2 * 1 ambiguities:
      Double_t kX, kY, kZ;                                  
      std::vector<TVector3> ambig;
      std::vector<Double_t> CHreco;     
      std::vector<Double_t> Tamb; 
      std::vector<Double_t> Path; 
                
      // length of the vector should be the same as NumAmb*8   !!!!!
      Int_t NumAmb = NAmb;  //=14        
      TVectorD CHdiff(8*NumAmb); 
      TVectorD Adiff(8*NumAmb); // vector to store angle diff btw real and reco photons
            
      for(Int_t i=0; i<NumAmb; i++){
        // in bar' coordinate system	
        kX = par[i*2]; // kXD	
	kY = par[i*2 + 1]; // kYD	
	kZ = -sqrt(1. - pow(kX,2.) - pow(kY,2.));	
	//cout<<"mom:  x = "<<fPMoB.X()<<",  y = "<<fPMoB.Y()<<", z = "<<fPMoB.Z()<<endl;	
	fkBar.SetXYZ(kX, kY, kZ);
	//cout<<"photon from LUT in bar' coordinate system: "<<endl;
	fkBar.Print();
		
	for(Int_t jamb=0; jamb <8; jamb++){
  		if(jamb == 1){fkBar.SetXYZ(-kX, kY, kZ);}
		if(jamb == 2){fkBar.SetXYZ(kX, -kY, kZ);}
		if(jamb == 3){fkBar.SetXYZ(kX, kY, -kZ);}
		if(jamb == 4){fkBar.SetXYZ(-kX, -kY, kZ);}
		if(jamb == 5){fkBar.SetXYZ(kX, -kY, -kZ);}
		if(jamb == 6){fkBar.SetXYZ(-kX, kY, -kZ);}
		if(jamb == 7){fkBar.SetXYZ(-kX, -kY, -kZ);}
							
		// check if there is such enrty in the LUT and if the ambiguity fulfills the total internal reflection requitrement
		if((kX != 0. && kY != 0.) && ((fkBar.Cross(fnX1)).Mag() > 1.00028/fGeo->nQuartz() || 
	   	   (fkBar.Cross(fnY1)).Mag() > 1.00028/fGeo->nQuartz())){
			ambig.push_back(fkBar);
			CHreco.push_back(fkBar.Angle(PMoBar));
			CHdiff(8*i+jamb) = fabs(CHreco[8*i+jamb] - CHexp);
			Tamb.push_back(RecoAmbigTime(fkBar, fStartVertex, &fPath, 0));			
			Path.push_back(fPath);	 
			Adiff(8*i+jamb) = fkBar.Angle(fPphoB);
			//cout<<"reco ch angle = "<<fkBar.Angle(PMoBar)<<endl;
			//cout<<"number of bounces = "<<NumberOfBounces(fStartVertex, fPphoB, fBarId)<<endl;	  
			//cout<<"CHdiff = "<<CHdiff(8*i+jamb)<<", CHreco"<<8*i+jamb<<" = "<<CHreco[8*i+jamb]<<endl;
			
			//fill lutinfo only with credible information:
			lutinfo.AddAngle(fkBar.Angle(PMoBar));
			lutinfo.AddTime(RecoAmbigTime(fkBar, fStartVertex, &fPath, 0));
			lutinfo.AddPath(fPath);
			lutinfo.AddChDiff(CHreco[8*i+jamb] - CHexp);
			//cout<<"fkBar = "<<endl;
			//fkBar.Print();
			lutinfo.AddNOfBounces(NumberOfBounces(fStartVertex, fkBar/*fPphoB*/, fBarId));
			
			fkBarXHist->Fill(kX, fPphoB.X()); 
			fkBarYHist->Fill(kY, fPphoB.Y());				  
		}else{
			ambig.push_back((0.,0.,0.));
			CHreco.push_back(-999.);
			CHdiff(8*i+jamb) = -999.;
			Tamb.push_back(-999.);			
			Path.push_back(-999.);	 
			Adiff(8*i+jamb) = -999.;	  				
		}
	}// for				
      }
      		
	// F I L L     T H E      H I S T O S			
	fhPDHits->Fill(xPPoi, yPPoi);//(fxPHit, fyPHit);	
	fhPDTime->Fill(ftime-ftime0); // substract production time of a photon 	
	
	fhRecoT1->Fill(flambdah, Ppt->GetLength()/(ftime-ftime0));
			
	fhCHreal->Fill(CHreal);
	fhLam->Fill(flambdah);
	  
	if(fB > 0.){
	  fMapHist->Fill(fPMo.Theta()/fpi*180.,fPMo.Phi()/fpi*180. - fHAngleInBDeg);  // assuming that R = 50 cm
	}	     
	if(fB == 0.){
	  fMapHist->Fill(fPMo.Theta()/fpi*180.,fPMo.Phi()/fpi*180.); 
	}
	cout<<"fMAPHIST is filled in theta="<<fPMo.Theta()/fpi*180.<<", phi ="<<fPMo.Phi()/fpi*180. - fHAngleInBDeg<<endl;
	cout<<"mother phi angle = "<<fPMo.Phi()/fpi*180.<<endl;
	cout<<"mother phi in the name of histo = "<<(Int_t)(fPMo.Phi()/fpi*180.*100.)<<endl;
	//#############################################
        TString PhiT = Form("t%g_phi%d", fPMo.Theta()/fpi*180., (Int_t)(fPMo.Phi()/fpi*180.*100.));
	TString PhiTCut = Form("t%g_phi%d_cut", fPMo.Theta()/fpi*180., (Int_t)(fPMo.Phi()/fpi*180.*100.));
	TString PhiTWeight = Form("t%g_phi%d_weight", fPMo.Theta()/fpi*180., (Int_t)(fPMo.Phi()/fpi*180.*100.));
	TString Nbo = Form("t%g_phi%d_nbo", fPMo.Theta()/fpi*180., (Int_t)(fPMo.Phi()/fpi*180.*100.));
		
        Int_t currPhiTh = PhiThetaPoints.size();
        for ( Int_t ihi = 0 ; ihi < PhiThetaPoints.size() ; ihi++ ) {
	  if ( PhiT == PhiThetaPoints[ihi]->GetName() ) {
	    currPhiTh = ihi;
	    break;
	  }
        }
        if ( currPhiTh == PhiThetaPoints.size() ) {
 	  Int_t nbins = (Int_t)(2.*fWidth/0.005);		
	  TH1F* newThetaPhiPoint = new TH1F(PhiT.Data(),"Cherenkov angle difference for theta", nbins, -fWidth,fWidth);
	  PhiThetaPoints.push_back(newThetaPhiPoint);
	  TH1F* newThetaPhiPointCut = new TH1F(PhiTCut.Data(),"Cherenkov angle difference for theta with time cut", nbins, -fWidth,fWidth);
	  PhiThetaPointsCut.push_back(newThetaPhiPointCut);
	  TH1F* newThetaPhiPointWeight = new TH1F(PhiTWeight.Data(),"Cherenkov angle difference for theta with time weights", nbins, -fWidth,fWidth);
	  PhiThetaPointsWeight.push_back(newThetaPhiPointWeight);
	  TH1F* newNbo = new TH1F(Nbo.Data(), "Number of bounces", 1000,0.,500.);
	  NboPoints.push_back(newNbo);	  
        }
        //PhiThetaPoints[currPhiTh]->Fill();      
        //#############################################	
	
	// fill N bounces for all photons:
	NboPoints[currPhiTh]->Fill(NumberOfBounces(fStartVertex, fPphoB, fBarId));
	fhNboLam->Fill(flambdah, NumberOfBounces(fStartVertex, fPphoB, fBarId));
	// fill all the ambiguities with weights:
	// WEIGHTS ARE BASED ON TIMING!!!!!
	for(Int_t j=0; j< CHreco.size(); j++){	  	  
	  fWeight = exp(- pow(Tamb[j] - ftime + ftime0, 2.)/2./pow(timeres->Eval(Path[j])/1000.,2.));
	  //cout<<"CHreco = "<<CHreco[j]<<", d time = "<<Tamb[j]-ftime + ftime0<<", weight = "<<fWeight<<endl;
	  
	  // using time cuts
	  if(fabs(CHdiff(j)) < fWidth){
	    PhiThetaPoints[currPhiTh]->Fill(CHreco[j]-CHexp,1.);
	    //fhDiff->Fill(CHreco(j)-CHexp,1./fInvSum/CHdiff(j));
	    fhDiff->Fill(CHreco[j]-CHexp,1.);
	    //fhDTime->Fill(Tamb[j] - ftime + ftime0);
	    //cout<<"CHreco = "<<CHreco[j]<<", time = "<<Tamb[j]<<", path = "<<Path[j]<<endl;	    
	    //cout<<"time resolution at this path = "<<timeres->Eval(Path[j])/1000.<<", d t = "<<Tamb[j] - ftime + ftime0<<endl;
	    //cout<<"histo "<<PhiThetaPoints[currPhiTh]->GetName()<<" filled with"<<CHreco[j]-CHexp<<endl;
	    // using time cuts
	    if(fabs(Tamb[j] - ftime + ftime0) < fNSigma* timeres->Eval(Path[j])/1000.){
	      PhiThetaPointsCut[currPhiTh]->Fill(CHreco[j]-CHexp);
	      //fhDTimeCut->Fill(Tamb[j] - ftime + ftime0);
	    }
	    //using weights based on time
	    PhiThetaPointsWeight[currPhiTh]->Fill(CHreco[j]-CHexp, fWeight);
	    //fhDTimeWeight->Fill(Tamb[j] - ftime + ftime0, fWeight);
	    
	    //cout<<"diff was filled with "<<CHreco(j)-CHexp<<endl;
	  }	  
	}
	
     // new ((*fChPhoArray)[nofChPho++]) PndChPho(fxPHit,fyPHit,fzPHit,ftime,flambdah,     	//					fPphoB,fPMo,fPphoInit,fPphoPD);	
    //} // if kxBar	
    } // if start vertex lies within the right bar!
    }// photon from primary particle 
             
 }// photon hits 
 new ((*fDrcLutInfoArray)[fDrcLutInfoArray->GetEntriesFast()]) PndDrcLutInfo(lutinfo);
}
//----------------------------------------------------------------------------------------------
Double_t PndDrcRecoLookupMap::InBarCoordSyst(TVector3 start, TVector3 *v1, TVector3 *v2, TVector3 *v3, TVector3 *v4){
    Double_t startPhi = start.Phi()/fpi*180.; // [degrees]
    //cout<<"-I- InBarCoordinateSystem: start phi = "<<start.Phi()/fpi*180.<<endl;
    if(startPhi < 0.){startPhi = 360. + startPhi;}
    //if(startPhi > 360.){startPhi = startPhi - 360.;}
    //cout<<"-I- InBarCoordinateSystem: start phi = "<<startPhi<<endl;
    //Double_t dphi = 2.*(180. - 2*fPipehAngle)/ fBboxNum; //[degrees]
    //cout<<"-I- InBarCoordinateSystem: dphi = "<<dphi<<endl;
    Double_t PhiRot = 0.; //[degrees]
    if(startPhi >= 0. && startPhi < 90.){
      PhiRot = TMath::Floor(startPhi/fDphi) *fDphi + fDphi/2.;
    }
    if(startPhi >= 90. && startPhi < 270.){
      PhiRot = 90.  + fPipehAngle + TMath::Floor((startPhi-90.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
    } 
    if(startPhi >= 270. && startPhi < 360.){
      PhiRot = 270. + fPipehAngle + TMath::Floor((startPhi-270.-fPipehAngle)/fDphi) *fDphi + fDphi/2.;
    }
    //cout<<"-I- InBarCoordinateSystem: PhiRot = "<<PhiRot<<endl;
      
    TVector3 ver1, ver2, ver3, ver4;
    ver1.SetXYZ(fR-(fHThick), fLSide/2. ,0.);
    ver2.SetXYZ(fR+(fHThick), fLSide/2. ,0.);
    ver3.SetXYZ(fR+(fHThick), -fLSide/2. ,0.);
    ver4.SetXYZ(fR-(fHThick), -fLSide/2. ,0.);
    
    ver1.RotateZ(PhiRot/180.*fpi);
    ver2.RotateZ(PhiRot/180.*fpi);
    ver3.RotateZ(PhiRot/180.*fpi);
    ver4.RotateZ(PhiRot/180.*fpi);
    
    *v1 = ver1;
    *v2 = ver2;
    *v3 = ver3;
    *v4 = ver4;   
    
    return PhiRot/180.*fpi;
}
//--------------------------------------------------------------------------------------------
Double_t PndDrcRecoLookupMap::SectorNum(Double_t xhit, Double_t yhit){

  // this function is used for photon hits to create LUTs, so that hits in front of target pipe area 
  // are treated as additional small LUTs
 
  TVector3 hit;
  hit.SetXYZ(xhit, yhit, 0.);  
  Double_t phi = hit.Phi()*180./fpi;
  if(phi < 0.){phi = 360. + hit.Phi()*180./fpi;}
  if(phi >= 0. && phi < 90.-fPipehAngle){
    //return 13.+TMath::Floor(phi/fDphi);
    return 1.+TMath::Floor((90.-fPipehAngle-phi)/fDphi);
  }
  if(phi >= 90.+fPipehAngle && phi<270.-fPipehAngle){
    //return TMath::Floor((phi-2.*fPipehAngle)/fDphi) - 3.; 
    return 9.+TMath::Floor((270.-fPipehAngle-phi)/fDphi);
  }
  if(phi >= 270.+fPipehAngle){
    //return TMath::Floor((phi-4.*fPipehAngle)/fDphi) - 3.;
    return 5.+TMath::Floor((360.-phi)/fDphi);
  }
  // small areas in front of the target pipe:
  if(phi >= 90.-fPipehAngle && phi < 90.+fPipehAngle){
   return 17.;
  }
  if(phi >= 270.-fPipehAngle && phi < 270.-fPipehAngle){
   return 18.;
  }
}

//------   Find Nubmer of Bounces     --------------------------------------
Int_t PndDrcRecoLookupMap::NumberOfBounces(TVector3 start, TVector3 dir, Int_t barId){
    // start - photon production point in global coord system
    // dir - photon direction in bar coord system
    
    //cout<<"-I- NumberOfBounces: dir"<<endl;
    //dir.Print();
    //cout<<"-I- NumberOfBounces: start"<<endl;
    //start.Print();   
       
    // Find coordinates of X0, Y0:
    Double_t Z0, X0, Y0;
    if(dir.Theta() < 3.1415/2.){
      Z0 = -(fabs(fzup) + 2.*fzdown - start.Z());
    }
    if(dir.Theta() >= 3.1415/2.){
      Z0 = -(start.Z() -  fzup);
    }    
    X0 = Z0*TMath::Tan(dir.Theta())*TMath::Cos(dir.Phi());
    Y0 = Z0*TMath::Tan(dir.Theta())*TMath::Sin(dir.Phi());
    //cout<<"-I- NumberOfBounces: tan th = "<<TMath::Tan(dir.Theta())<<", sin ph = "<<TMath::Sin(dir.Phi())<<", cos ph = "<<TMath::Cos(dir.Phi())<<endl;
    //cout<<"-I- NumberOfBounces: X0 = "<<X0<<", Y0 = "<<Y0<<", Z0 = "<<Z0<<endl;
    
    // Find the start position of the photon with respect to the middle of the bar:
    TVector3 startLocal = fGeoH->MasterToLocalShortId(start, barId);
    //cout<<"-I- NumberOfBounces: Xen = "<<startLocal.X() + fBarWidth/2.<<", Yen = "<<startLocal.Y() + fHThick<<endl;
        
    // Find the number of bounces in each direction       
    Double_t N1, N2;   
    FindOutPoint(X0, startLocal.X() + fBarWidth/2., fBarWidth,  &N1, 0);
    FindOutPoint(Y0, startLocal.Y() + fHThick,      2.*fHThick, &N2, 0);
    //cout<<"-I- NumberOfBounces: N1 = "<<N1<<", N2 = "<<N2<<endl;
   
    return (Int_t)N1+N2;
}

//----------------------------------------------------------------------------------------------------------
Double_t PndDrcRecoLookupMap::FindOutPoint(Double_t x0, Double_t xEn, Double_t a, Double_t *NN, Bool_t printt){	
	Double_t m=99.;
	Double_t n=TMath::Floor(x0/a);
	m = n;		
	if(printt){std::cout<<"n = "<<n<<", NN = "<<*NN<<", x0 = "<<x0<<", a = "<<a<<std::endl;}
	Double_t x1 = x0 - n*a;
	if(x0 < 0.){x1 = x0 - (n+1)*a;}
	if(printt){std::cout<<"xy = "<< x1<<std::endl;}
	Double_t xK = 0.;
	if((m/2. - TMath::Floor(m/2.)) == 0.) { // 4etnoe
		if(print){std::cout<<"odd==0"<<std::endl;}
		if(x0 >= 0. && x1 + xEn <= a){xK = x1 + xEn;} 
		if(x0 >= 0. && x1 + xEn >  a){xK = 2*a - x1 - xEn; n = 1. + n;}
		if(x0 < 0. && x1 + xEn >= 0.){xK = a - (x1 + xEn); n = -1. -n;}
		if(x0 < 0. && x1 + xEn < 0.) {xK = a + x1 + xEn; n = -n;}
	if(printt){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}
		
	}

	if((m/2. - TMath::Floor(m/2.)) != 0.) { // ne4etnoe
		if(print){std::cout<<"even!=0"<<std::endl;} 
		if(x0 >= 0. && x1 + xEn <= a){xK = a - (x1 + xEn);} 
		if(x0 >= 0. && x1 + xEn >  a){xK = x1 + xEn - a; n = 1. + n;}
		if(x0 < 0. && x1 + xEn >= 0.){xK = x1 + xEn; n = -1. -n;} 
		if(x0 < 0. && x1 + xEn < 0.) {xK = - (x1 + xEn); n = -n;}
	if(printt){std::cout<<"xK = "<< xK<<", n = "<<n<<std::endl;}

	}
 	
	*NN = n;		
	return xK;
}

//-------------------------------------------------------------------------------------------
void PndDrcRecoLookupMap::DrawBarBox(TVector3 v1, TVector3 v2, TVector3 v3, TVector3 v4){
  Double_t xx[5];
  Double_t yy[5];
  xx[0] = v1.X(); xx[1] = v2.X(); xx[2] = v3.X(); xx[3] = v4.X(); xx[4] = v1.X();	
  yy[0] = v1.Y(); yy[1] = v2.Y(); yy[2] = v3.Y(); yy[3] = v4.Y(); yy[4] = v1.Y();
  TPolyLine *p6 =new TPolyLine(5,xx,yy);
  p6->SetLineWidth(2);
  p6->Draw("same");
}

//--------------------------------------------------------------------------------------------
Double_t PndDrcRecoLookupMap::RecoAmbigTime(TVector3 kb, TVector3 start, Double_t *l, Bool_t printout){  // [ns]  

  // group velocities for oil and quartz for Noil and Nquartz assuming that mean lambda = 410 nm:
  //Double_t nOil = fGeo->nEV();
  Double_t u_oil = 19.8;//30./nOil;
  //Double_t u_quartz = 30./1.46907; // for lambda = 410 nm 
  //Double_t u_quartz = 30./1.47012; // for lambda = 400 nm 
  //Double_t u_quartz = 30./1.46838; // for lambda = 417 nm 
  //Double_t u_quartz = 30./1.47125; // for lambda = 390 nm 
  //Double_t u_quartz = 30./1.47248; // for lambda = 380 nm
  Double_t u_quartz = 19.8;//30./fNquartz;//1.47805; // for lambda = 370 nm
  
  // kb - photon momentum at the production point, first calculate path in QUARTZ BAR:
  // assume that there are only DIRECT photons!!!   
  
  // path in quartz:
  Double_t Z0 = 0.;
  Double_t kz = 0.;
  if(kb.Z() <  0.){Z0 = start.Z() - fGeo->barBoxZUp(); kz = kb.Z();}
  if(kb.Z() >= 0.){Z0 = 2.*fGeo->barBoxZDown() - fGeo->barBoxZUp() - start.Z(); kz = -kb.Z();}
  
  if(printout){
    cout<<"-I- Reco ambig Time: Zstart = "<<start.Z()<<", Z0 = "<<Z0<<endl;
  }
  
  Double_t L0 = Z0/fabs(cos(kb.Theta()));
  kb.SetZ(kz);
  Double_t thetaOil = asin(fGeo->nQuartz()/fGeo->nEV()*sin(kb.Theta()));
  if(printout){
    cout<<"RECO ambig TIME: thetas are: quartz = "<<kb.Theta()/fpi*180.<<", oil = "<<thetaOil/fpi*180.<<endl;
  }
  kb.SetTheta(thetaOil);
  Double_t Z1 = fGeo->EVlen();
  Double_t L1 = Z1/fabs(cos(kb.Theta())); // path in oil
    
  if(printout){
    cout<<"-I- Reco ambig Time: Zev = "<<fGeo->EVlen()<<endl;     
    cout<<"-I- Reco ambig Time: L0 = "<<L0<<", L1 = "<<L1<<endl;
    cout<<"-I- Reco ambig Time: time is "<< L0/u_quartz + L1/u_oil <<endl; // 19.5 cm/ns - average (lambda) group speed
  } 
  
  *l = L0 + L1;
   
  return L0/u_quartz + L1/u_oil; //[ns]  
}

void PndDrcRecoLookupMap::CreateHisto()
{
  // Histogram list
  fHistoList = new TList(); 
  
  fhPDTime = new TH1D("fhPDTime","Time of photons in ns",1000,0,100);  
  fhRecoT1 = new TH2F("fhRecoT1","Reconstructed group velocity using MC track path [ns]", 500, 200.,700., 1000, 0., 40.);  
  fhPDHits = new TH2D("fhPDHits","Photon hits on the PD plane", 2500, -250, 250, 2500, -250, 250);  

  fhDiff = new TH1F("fhDiff","Difference btw reco and expected Cherenkov Angle", 50, -fWidth, fWidth);

  fhCHreal = new TH1F("fhCHreal","Generated cherenkov angle", 100, 0.7, 0.9);
  fhLam = new TH1F("fhLam","Lambda of detected Cherenkov photons", 480, 200., 800.);
  
  fhNboLam = new TH2F("fhNboLam", "Nbo as a function of lambda", 600, 200.,800., 100, 0.,500.);
 
  fHistoList->Add(fhPDTime);  
  fHistoList->Add(fhPDHits); 
  fHistoList->Add(fhDiff); 
  fHistoList->Add(fhCHreal);
  fHistoList->Add(fhLam);
  // also add temp histos
  fHistoList->Add(fhRecoT1);
  fHistoList->Add(fhNboLam);
  
  fPixelSize = 0.65;
   
  fMapHist = new TH2F("fMapHist","Photon yield",30,2.5,152.5,20,0.35,20.35);// 3 bars 41,0.35,20.85); //25,0.5,25.5);  
  fSigHist = new TH2F("fSigHist","fSigHist",30,2.5,152.5, 20,0.35,20.35);// 3 bars 41,0.35,20.85);//7 bars("fSigHist","fSigHist",30,2.5,152.5, 25,0.5,25.5);// 5 bars in bar box
  fCheHist = new TH2F("fCheHist","fCheHist",30,2.5,152.5, 20,0.35,20.35);// 3 bars 41,0.35,20.85); //("fCheHist","fCheHist",30,2.5,152.5, 25,0.5,25.5); 
  fkBarXHist = new TH2F("fkBarXHist","fkBarXHist",200,-1.,1.,200,-1.,1.);
  fkBarYHist = new TH2F("fkBarYHist","fkBarYHist",200,-1.,1.,200,-1.,1.);
  fHistoList->Add(fMapHist);  
  fHistoList->Add(fSigHist);  
  fHistoList->Add(fCheHist);   
  fHistoList->Add(fkBarXHist);
  fHistoList->Add(fkBarYHist);
}

//------------------Write to File----------------------------------------------
void PndDrcRecoLookupMap::WriteToFile()
{	
	fOutputName.Append(":/");
	gDirectory->cd(fOutputName);
        //gDirectory->cd("13_t90_500ev_l0_oil_5_HIT.root:/");
  	
	gDirectory->mkdir("diff");
	gDirectory->cd("diff");
	for(int i=0; i<PhiThetaPoints.size();i++){	      
	      	PhiThetaPoints[i]->Write();	      
	}
		
	gDirectory->cd("../");
	gDirectory->mkdir("diffCut");
	gDirectory->cd("diffCut");
	for(int i=0; i<PhiThetaPointsCut.size();i++){	     
	      	PhiThetaPointsCut[i]->Write();	      		      
	}	
	
	gDirectory->cd("../");
	gDirectory->mkdir("diffWeight");
	gDirectory->cd("diffWeight");
	for(int i=0; i<PhiThetaPointsWeight.size();i++){	     
	      	PhiThetaPointsWeight[i]->Write();	      		      
	}
	
	gDirectory->cd("../");
	gDirectory->mkdir("Nbounces");
	gDirectory->cd("Nbounces");
	for(int i=0; i<NboPoints.size();i++){	     
	      	NboPoints[i]->Write();	      		      
	}
	gDirectory->cd("../");
	TIter next(fHistoList);
		while ( TH1* histo = ((TH1*)next()) ) histo->Write();
	
}

// -----   Finish Task   ---------------------------------------------------
void PndDrcRecoLookupMap::Finish()
{
  cout << "-I- PndDrcRecoLookupMap: Finish" << endl;
      
  // creating a map in every pixel   
  Double_t sig = 0., chr = 0.;    
  TString hnam, num1, num2, num3;
  Float_t xx, yy;
  for(int i=0; i<PhiThetaPoints.size(); i++){
       
       TF1* fit = new TF1("fit","gaus(0)+pol1(3)");
       fit->SetParameters(10.,-0.001,0.017,1.,1.);
       PhiThetaPoints[i]->Fit("fit","V","");
       sig = fabs(fit->GetParameter(2));
       chr = fit->GetParameter(1);       
       hnam = PhiThetaPoints[i]->GetName();       
       //cout<<"hnam - "<<hnam<<endl;       		       				
       num1 = hnam;
       num1.Remove(0,1);
       num1.Remove(num1.Index("_"), (num1.Length()-num1.Index("_")));
       xx = num1.Atof();
       //cout<<"t = "<<xx<<endl;
       num2 = hnam;
       num2.Remove(0,num2.Index("i")+1);
       yy = num2.Atof()/100.-fHAngleInBDeg;
       //cout<<"phi = "<<yy<<endl;       			
	cout<<"i = "<<i<<endl;		
	fSigHist->Fill(xx,yy,sig);
	cout<<"SigHIST was filled with "<<sig<<" in phi = "<<yy <<", theta - "<<xx <<endl;
	fCheHist->Fill(xx, yy, chr);
  } 
/*  
  // fill diff map with low value to avoid the background color
  for(Int_t i=0; i<fMapHist->GetXaxis()->GetNbins(); i++){
    for(Int_t j=0; j<fMapHist->GetYaxis()->GetNbins(); j++){
      if(fMapHist->GetBinContent(i+1,j+1) == 0.){
        fMapHist->SetBinContent(i+1, j+1, -10.);
      }
    }
  }
*/    
  WriteToFile();
  DrawHisto();
}

//----------------------------------------------------------------
void PndDrcRecoLookupMap::DrawHisto()
{
gStyle->SetOptFit(0111);

TCanvas *C1 = new TCanvas("C1","Time", 500,500);
C1->Divide(1,2);
C1->cd(1);
fhPDTime->Draw();
C1->cd(2);
fhRecoT1->Draw();

TCanvas *C2 = new TCanvas("C2","PD Plane", 500, 500);
fhPDHits->SetMarkerStyle(20);
fhPDHits->SetMarkerSize(0.25);
//fhPDHits->SetMarkerColor(2);
fhPDHits->Draw();

TCanvas *C3 = new TCanvas("C3","Diff", 500,500);
fhDiff->GetXaxis()->SetTitle("[rad]");
fhDiff->SetLineWidth(3);
fhDiff->Draw();

// C A N V A S     W I T H      M A P S (lambda and NEntries of pixels)  

TCanvas *C8= new TCanvas("C8","Nentries map",500,500);
SetPlotStyle();
gStyle->SetOptStat(11);
fMapHist->GetXaxis()->SetTitle("#Theta_{track} [degrees]");
fMapHist->GetYaxis()->SetTitle("#phi_{track} [degrees]");
fMapHist->GetZaxis()->SetTitle("nEntries");
//fMapHist->Draw("LEGO2Z");
fMapHist->Scale(0.05);
fMapHist->Draw("COL2Z");
//DrawDetectorLayout();
//DrawBarBox(fBBver1,fBBver2,fBBver3,fBBver4);


TCanvas *C9= new TCanvas("C9","Mean CHreco map",500,500);
SetPlotStyle();
gStyle->SetOptStat(11);
fCheHist->GetXaxis()->SetTitle("#theta track angle [degrees]");
fCheHist->GetYaxis()->SetTitle("#phi track angle [degrees]");
fCheHist->GetZaxis()->SetTitle("<CHreco - CHexp> [rad]");
//fCheHist->Draw("LEGO2Z");
fCheHist->Draw("COL2Z");
//DrawDetectorLayout();
//DrawBarBox(fBBver1,fBBver2,fBBver3,fBBver4);

TCanvas *C10= new TCanvas("C10","Sigma map",500,500);
SetPlotStyle();
gStyle->SetOptStat(11);
//fSigHist->SetMinimum(-0.2);
//fSigHist->SetMaximum(0.2);
fSigHist->GetXaxis()->SetTitle("#theta track angle [degrees]");
fSigHist->GetYaxis()->SetTitle("#phi track angle [degrees]");
fSigHist->GetZaxis()->SetTitle("#sigma_{CHreco - CHexp}");
fSigHist->Draw("COL2Z");
//DrawDetectorLayout();
//DrawBarBox(fBBver1,fBBver2,fBBver3,fBBver4);

TCanvas *C11= new TCanvas("C11","kBar LUT & kBar gen",500,500);
C11->Divide(2,1);
C11->cd(1);
fkBarXHist->GetXaxis()->SetTitle("kxBar from LUT");
fkBarXHist->GetYaxis()->SetTitle("kxBar generated");
fkBarXHist->Draw();
TLine *l1 =new TLine(1.,1.,-1.,-1.);
TLine *l2 = new TLine(-1.,1.,1.,-1.);
l1->SetLineColor(8);
l2->SetLineColor(8);
l1->Draw("same");
l2->Draw("same");
C11->cd(2);
fkBarYHist->GetXaxis()->SetTitle("kyBar from LUT");
fkBarYHist->GetYaxis()->SetTitle("kyBar generated");
fkBarYHist->Draw();
l1->Draw("same");
l2->Draw("same");

TCanvas *C12= new TCanvas("C12","Lam & CHreal",500,500);
C12->Divide(1,2);
C12->cd(1);
fhCHreal->Draw();
C12->cd(2);
fhLam->Draw();

}

// -------------------------------------------------------------------------
ClassImp(PndDrcRecoLookupMap)