// -------------------------------------------------------------------------
// -----                FairFtofProducerIdeal source file             -----
// -----                  Created by A. Sanchez              -----
// -------------------------------------------------------------------------

#include <cmath>

#include "TClonesArray.h"
#include "TGeoManager.h"
#include "FairRootManager.h"
#include "PndFtofHitProducerIdeal.h"
#include "PndFtofHit.h"
#include "TGeoBBox.h"
//#include "PndFtofHitInfo.h"
#include "PndFtofPoint.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairGeoVector.h"
#include "TVector3.h"


// -----   Default constructor   -------------------------------------------
PndFtofHitProducerIdeal::PndFtofHitProducerIdeal() :
  FairTask("Ideal PndFtof Hit Producer") 
{
	fBranchName 	= "FtofPoint";

}
// -------------------------------------------------------------------------

// -----   Default constructor   -------------------------------------------
PndFtofHitProducerIdeal::PndFtofHitProducerIdeal(Double_t dt, Double_t dt2) :
  FairTask("Ideal PndFtof Hit Producer") 
{
	fBranchName 	= "FtofPoint";

	fdt= dt;	fdt2= dt2;
	
}
// -------------------------------------------------------------------------

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

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

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

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

 

  // Create and register output array
  fHitArray = new TClonesArray("PndFtofHit");
  ioman->Register("FtofHit", "Ftof", fHitArray, kTRUE);

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




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

}


// -----   Public method Exec   --------------------------------------------
void PndFtofHitProducerIdeal::Exec(Option_t* opt) 
{
  // Reset output array
  if ( ! fHitArray ) 
    Fatal("Exec", "No HitArray");
  
  fHitArray->Clear();
 
  
  // Declare some variables
  PndFtofPoint *point = 0;

  Int_t 
    detID = 0,       // Detector ID
    trackID = 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 FtofPoints
  Int_t 
    nPoints = fPointArray->GetEntriesFast();

  for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) 
    {
      point = (PndFtofPoint*) 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]) PndFtofHit(trackID, detID, 
					   point->GetDetName(),time, t1, 
					   position,dpos,iPoint, 
					   point->GetEnergyLoss());
      //std::cout << "Hit created for module: " << point->GetDetName() << std::endl;
      
      
      
    }   // Loop over MCPoints

 
 
  
  // Event summary
  std::cout << "-I- PndFtofHitProducerIdeal: " << nPoints << " FtofPoints, "
       << nPoints << " Hits created." << std::endl;
  

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

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

void PndFtofHitProducerIdeal::GetLocalHitPoints(PndFtofPoint* 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 PndFtofHitProducerIdeal::GetTransformation(std::string detName) const
{
  gGeoManager->cd(detName.c_str());
  TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix();
  if (fVerbose > 1)
    transMat->Print("");
  return *transMat;
}

TVector3 PndFtofHitProducerIdeal::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(PndFtofHitProducerIdeal)