// -------------------------------------------------------------------------
// -----                PndTofHitProducerIdeal source file             -----
// -------------------------------------------------------------------------

#include <cmath>

#include "TClonesArray.h"
#include "TGeoManager.h"
#include "FairRootManager.h"
#include "PndTofHitProducerIdeal.h"
#include "PndTofHit.h"
#include "TGeoBBox.h"
//#include "PndTofHitInfo.h"
#include "PndTofPoint.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairGeoVector.h"
#include "TVector3.h"


// -----   Default constructor   -------------------------------------------
PndTofHitProducerIdeal::PndTofHitProducerIdeal() :
  FairTask("Ideal PndTof Hit Producer") 
{
	fBranchName 	= "TofPoint";
	fBranchName2 	= "TofSciFPoint";
}
// -------------------------------------------------------------------------

// -----   Default constructor   -------------------------------------------
PndTofHitProducerIdeal::PndTofHitProducerIdeal(Double_t dt, Double_t dt2) :
  FairTask("Ideal PndTof Hit Producer") 
{
	fBranchName 	= "TofPoint";
	fBranchName2 	= "TofSciFPoint";
	fdt= dt;	fdt2= dt2;
	
}
// -------------------------------------------------------------------------

// -----   Destructor   ----------------------------------------------------
PndTofHitProducerIdeal::~PndTofHitProducerIdeal() 
{ 
}

// -----   Public method Init   --------------------------------------------
InitStatus PndTofHitProducerIdeal::Init() 
{
  // Get RootManager
  FairRootManager* ioman = FairRootManager::Instance();
  
  if ( ! ioman ) 
    {
      std::cout << "-E- PndTofHitProducerIdeal::Init: "
	   << "RootManager not instantiated!" << std::endl;
      return kFATAL;
    }

  // Get input array
  fPointArray = (TClonesArray*) ioman->GetObject(fBranchName);

  if ( ! fPointArray ) 
    {
      std::cout << "-W- PndTofHitProducerIdeal::Init: "
	   << "No TofPoint array!" << std::endl;
      return kERROR;
  }

   fsciFPointArray = (TClonesArray*) ioman->GetObject(fBranchName2);

  if ( ! fsciFPointArray ) 
    {
      std::cout << "-W- PndTofHitProducerIdeal::Init: "
	   << "No TofSciFPoint array!" << std::endl;
      return kERROR;
  }

  // Create and register output array
  fHitArray = new TClonesArray("PndTofHit");
  ioman->Register("TofHit", "Tof", fHitArray, kTRUE);

  //std::cout << "-I- PndTofHitProducerIdeal: Intialisation successfull" << std::endl;
  //return kSUCCESS;


 // Create and register output array
  fsciFHitArray = new TClonesArray("PndTofHit");
  ioman->Register("TofSciFHit", "Tof", fsciFHitArray, kTRUE);

  std::cout << "-I- PndTofHitProducerIdeal: Intialisation successfull" << std::endl;
  return kSUCCESS;
}
// -------------------------------------------------------------------------
void PndTofHitProducerIdeal::SetParContainers()
{
  // Get Base Container
  FairRun* ana = FairRun::Instance();
  FairRuntimeDb* rtdb=ana->GetRuntimeDb();
  //fGeoPar = (PndGeoTofPar*)(rtdb->getContainer("PndGeoTofPar"));

}


// -----   Public method Exec   --------------------------------------------
void PndTofHitProducerIdeal::Exec(Option_t* opt) 
{
  // Reset output array
  if ( ! fHitArray ) 
    Fatal("Exec", "No HitArray");
  
  fHitArray->Clear();
  fsciFHitArray->Clear();
  
  // Declare some variables
  PndTofPoint *point = 0;
  PndTofPoint *sciFpoint = 0;
  Int_t 
    detID = 0,       // Detector ID
    trackID = 0;     // Track index

  Int_t 
    scidetID = 0,       // Detector ID
    scitrackID = 0;     // Track index

  Double_t time = 0.;Double_t scitime = 0.;
  Double_t t2 = 0.;
  Double_t t1 = 0.;
  Double_t phdt2 = 0.,tdc12 =0.,tz=0.;
  Double_t phdt1 = 0.,tsig = 0.;

  // Loop over TofPoints
  Int_t 
    nPoints = fPointArray->GetEntriesFast();

  for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) 
    {
      point = (PndTofPoint*) fPointArray->At(iPoint);
//	std::cout << " Ideal Hit Producer -Point-: " << point << std::endl;
      if ( ! point) 
	continue;

      // Detector ID
      detID = point->GetVolumeID();

      // MCTrack ID
      trackID = point->GetTrackID();

      FairGeoVector posCInL, posCOut, meanPos, meanPosL;
      Double_t ZLoc;Double_t Ztdc[3];
      Double_t DistZ[3];

      GetLocalHitPoints(point, posCInL,meanPos);

      TVector3 size = GetSensorDimensions(point->GetDetName().Data());

      if (meanPos.getZ()>0.)
      {
	  tz = size.z()-meanPos.getZ();
	  phdt1 = 1.492*(tz)/30;
	  phdt2 = 1.492*(2*size.z()-tz)/30;
	  }
	  if(meanPos.getZ()<0.){
	    tz = -size.z()-meanPos.getZ();
	  phdt1 = 1.492*(2*size.z()+tz)/30;
	  phdt2 = -1.492*(tz)/30;
	 }
      if(meanPos.getZ()==0.)phdt1=phdt2=1.492*(size.z())/30;
      
      //   std::cout<<" z tdc "<<tz<<" z point "<<meanPos.getZ()<<std::endl;
      //       std::cout<<" t1 phot "<<phdt1<<" t2 phot "<<phdt2<<std::endl;
      
      tdc12 = (phdt1-phdt2);
      tsig = 0.08;
      smear(tdc12,tsig);

      
      if(phdt2>phdt1){
	ZLoc = ((tdc12*30)/(2*1.492))+size.z();
	DistZ[2]=size.z()-ZLoc;
	//std::cout<<" t2>t1"<<std::endl;
	
      }
      
      if(phdt1>phdt2)
	{ZLoc = ((tdc12*30)/(2*1.492))-size.z();
	  DistZ[2]=-size.z()-ZLoc;
	  //std::cout<<" t1>t2"<<std::endl;
	}
      

      if(phdt1==phdt2)DistZ[2]=0.;

      DistZ[0]=DistZ[1]=0.;

      //std::cout<<" z loc "<<ZLoc<<" distz "<<DistZ[2]<<std::endl;
      
      TGeoHMatrix trans = GetTransformation(point->GetDetName().Data());
      trans.LocalToMaster(DistZ, Ztdc);
      
      
      //std::cout<<" z tdc "<<Ztdc[2]<<" z point "<<point->GetZin()<<std::endl;
      //        ____________/__
      // phdt2  |      @   / | phdt1, (phdt1-phdt2)c/2n=x
      //        |_________/__|
      //           L-x   / x 

      //phdt1 phdt 2 time signals produced at each side of the
      //scintillator bar(scintillation process not treated
      // i introduce some aprox. for the TDC signal
      //to calculate the position in the bar by knowing tdc time. 


      //refraction index polipropilrne 1.492 poliviniltoluene 1.58
      
      
      
      TVector3 dpos; 

      dpos.SetXYZ(3.,0.25,fabs(Ztdc[2]-point->GetZin()));

      // coord. of the center of the slab in lab.c.s
      // z coord. determined by measuring t1-t2
      // at each side of the bar
      TVector3 position(posCInL.getX(), 
			posCInL.getY(), 
			Ztdc[2]);
      //point->Position(pos);
   
      
      time = point->GetTime();
      
      t1 = 0.080;//100 ps time resolution
      smear(time,t1);
 
      // Create new hit
      new ((*fHitArray)[iPoint]) PndTofHit(trackID, detID, 
					   point->GetDetName(),time, t1, 
					   position,dpos,iPoint, 
					   point->GetEnergyLoss());
      //std::cout << "Hit created for module: " << point->GetDetName() << std::endl;
      
      
      
    }   // Loop over MCPoints

 
  // Loop over TofsciFPoints
  

  Int_t nFPoints = fsciFPointArray->GetEntriesFast();

  for (Int_t iFPoint = 0; iFPoint < nFPoints; iFPoint++) 
    {
      sciFpoint = (PndTofPoint*) fsciFPointArray->At(iFPoint);
      //std::cout << " Ideal Hit Producer -Point-: " <<  sciFpoint << std::endl;
      if ( ! point) 
	continue;

      // Detector ID
      scidetID = sciFpoint->GetVolumeID();

      // MCTrack ID
      scitrackID = sciFpoint->GetTrackID();
      TVector3 dpos2;
       dpos2.SetXYZ(0.1,0.1,0.1);
      TVector3 pos(sciFpoint->GetXin(), 
		   sciFpoint->GetYin(), 
		   sciFpoint->GetZin());
      
      
      scitime = sciFpoint->GetTime();
      
      t2 = 0.2;//100 ps time resolution
      smear(scitime,t2);
      //std::cout<<" time gauss"<<scitime<<" "<<t2<<std::endl;
      // Create new hit
      new ((*fsciFHitArray)[iFPoint]) PndTofHit(scitrackID, scidetID, 
						sciFpoint->GetDetName(),
						scitime, t2, pos,dpos2,iFPoint, 
						sciFpoint->GetEnergyLoss());

      
      
      
      
    }
  
  // Event summary
  std::cout << "-I- PndTofHitProducerIdeal: " << nPoints << " TofPoints, "
       << nPoints << " Hits created." << std::endl;
   // Event summary
  std::cout << "-I- PndTofHitProducerIdeal: " << nFPoints << " SciFTofPoints, "
	    << nFPoints << " sciF Hits created." << std::endl;

}
// -------------------------------------------------------------------------
void PndTofHitProducerIdeal::smear(Double_t& time, Double_t& dt)
{
/// smear a 3d vector

  Double_t t = time;
  //std::cout<<" time "<<time<<std::endl;
  Double_t sigt;
  
  sigt=gRandom->Gaus(0,dt);
  t += sigt;
  time = t;
  return;
}

void PndTofHitProducerIdeal::GetLocalHitPoints(PndTofPoint* myPoint, 
					       FairGeoVector& myHitIn,FairGeoVector& myInL)
{
 	 
  if (fVerbose > 1)
    std::cout << "GetLocalHitPoints" << std::endl;
  TGeoHMatrix trans = GetTransformation(myPoint->GetDetName().Data());
  
  Double_t posIn[3];
  Double_t posOut[3];
  Double_t posInLocal[3];
  Double_t posOutLocal[3];
  Double_t posCInLocal[3];
  Double_t posCOutLocal[3];
  Double_t posCIn[3];
  Double_t posCOut[3];
  
  posIn[0] = myPoint->GetXin();
  posIn[1] = myPoint->GetYin();
  posIn[2] = myPoint->GetZin();
  
  posOut[0] = myPoint->GetXout();
  posOut[1] = myPoint->GetYout();
  posOut[2] = myPoint->GetZout();
   
  //trans.MasterToLocal(posIn, posInLocal);
  
  if (fVerbose > 1){
    for (Int_t i = 0; i < 3; i++)
      std::cout << "posIn "<< i << ": " << posIn[i] << std::endl;
    //std::cout << "posInLOcal "<< i << ": " << posInLocal[i] << std::endl;
    trans.Print("");
  }
  trans.MasterToLocal(posIn, posInLocal);
  
  posCInLocal[0]=0.;
  posCInLocal[1]=0.;
  posCInLocal[2]=0.;
  
  trans.LocalToMaster(posCInLocal, posCIn);
  
  if (fVerbose > 1) {
    for (Int_t i = 0; i < 3; i++){
      std::cout << "posCInLocal "<< i << ": " << posCInLocal[i] << std::endl;
      std::cout << "posCIn "<< i << ": " << posCIn[i] << std::endl;
      std::cout << "posInLocal "<< i << ": " << posInLocal[i] << std::endl;
    }
  }

  myHitIn.setVector(posCIn);
  myInL.setVector(posInLocal);

}


 	
 	
TGeoHMatrix PndTofHitProducerIdeal::GetTransformation(std::string detName) const
{
  gGeoManager->cd(detName.c_str());
  TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix();
  if (fVerbose > 1)
    transMat->Print("");
  return *transMat;
}

TVector3 PndTofHitProducerIdeal::GetSensorDimensions(std::string detName) const
	{
 	  gGeoManager->cd(detName.c_str());
 	  TGeoVolume* actVolume = gGeoManager->GetCurrentVolume();
 	  TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape());
 	  TVector3 result;
 	  result.SetX(actBox->GetDX());
 	  result.SetY(actBox->GetDY());
 	  result.SetZ(actBox->GetDZ());
 	 
 	  //result.Dump();
 	 
 	  return result;
 	}

ClassImp(PndTofHitProducerIdeal)