/////////////////////////////////////////////////////////////
//
//  PndHyp
//
//  Filler of PndHypPoint
//
//  created by A. Sanchez
//
///////////////////////////////////////////////////////////////

#include <iostream>
//#include "ostringstream.h"

#include "TClonesArray.h"
#include "TGeoManager.h"
#include "TLorentzVector.h"
#include "TParticle.h"
#include "THParticle.h"
#include "TVirtualMC.h"
#include "TString.h"
#include "FairGeoTransform.h"

#include "TGeoBBox.h"
#include "FairGeoInterface.h"
#include "FairGeoLoader.h"
#include "TGeoMCGeometry.h"
#include "FairGeoNode.h"
#include "FairGeoMedium.h"
#include "PndGeoHyp.h"
#include "FairGeoRootBuilder.h"
#include "PndStack.h"
#include "PndHyp.h"
#include "PndHypPoint.h"

#include "FairRootManager.h"
#include "FairVolume.h"
// add on for debug

#include "FairRuntimeDb.h"
#include "TObjArray.h"
#include "FairRun.h"
#include "FairRunSim.h"

#include "PndHypGeoHandling.h"

//#include "PndHypDecayer.h"
//#include "HypStatDecay.h"
#include "TArrayI.h"
#include "TMCProcess.h"

#include "TList.h"
#include "TFile.h"

#include <string>
#include <sstream>
using std::cout;
using std::endl;
using std::ostringstream;
class FairVolume;

// -----   Default constructor   -------------------------------------------
PndHyp::PndHyp() {
  fHypCollection        = new TClonesArray("PndHypPoint");
  fHypSecTarCollection  = new TClonesArray("PndHypPoint");
  //fHypSTpipeCollection  = new TClonesArray("PndHypPoint");
  SiId = 0;
  CId = 0;
  CpipeId = 0;
  alId = 0;
  beId = 0;
  fPosIndex   = 0; 
  // fpreflag = 0;  
  //fpostflag = 0;
  fEventID=-1; 

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

// -----   Standard constructor   ------------------------------------------
PndHyp::PndHyp(const char* name, Bool_t active)
  : FairDetector(name, active) {
    fHypCollection        = new TClonesArray("PndHypPoint");
    fHypSecTarCollection  = new TClonesArray("PndHypPoint");
    //fHypSTpipeCollection  = new TClonesArray("PndHypPoint");
    SiId = 0;
     CId = 0;
     alId = 0;
     beId = 0;
    fPosIndex   = 0;
    
    fEventID=-1;
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
PndHyp::~PndHyp() {
  if (fHypCollection) {
    fHypCollection->Delete();
    delete fHypCollection;
  }
  if (fHypSecTarCollection) {
    fHypSecTarCollection->Delete();
    delete fHypSecTarCollection;
  }


 delete fGeoH;
 //delete fread;
  
}
// -------------------------------------------------------------------------



// -----   Public method Intialize   ---------------------------------------
void PndHyp::Initialize() {
  // Init function
  
  FairDetector::Initialize();

 

  if(0==gGeoManager) {
    std::cout<<" -E- No gGeoManager in PndHyp Detector::Initialize()!"<<std::endl;
    abort();
  }
  fGeoH = new PndHypGeoHandling(gGeoManager);

  //fread = new PndHypDecayer("hypBupDecay2.root");
  //fread = new HypStatDecay("statDecay");
  
  TGeoMedium *Si= gGeoManager->GetMedium("HYPsilicon");
  SiId=  Si->GetId();
  
  //----disactivated when geo file is block
  TGeoMedium *C= gGeoManager->GetMedium("HYPdiamond");
  //TGeoMedium *C= gGeoManager->GetMedium("carbon");
  CId=  C->GetId();
  
  TGeoMedium *Cpipe= gGeoManager->GetMedium("HYPcarbon");
  CpipeId=  Cpipe->GetId();
  
  //    TGeoMedium *al= gGeoManager->GetMedium("HYPaluminium");
  //    alId=  al->GetId();
  //    TGeoMedium *be= gGeoManager->GetMedium("HYPberyllium");
  //    beId=  be->GetId();
  

}
// -------------------------------------------------------------------------
void PndHyp::BeginEvent(){
  // Begin of the event
  
}


// -------------------------------------------------------------------------
void PndHyp::PreTrack(){
  // Begin of the event
  
  fTrackStopNxtStep=kFALSE;
  // cout << " PndHyp::PreTrack() "<< endl;
  
}


void PndHyp::SetSpecialPhysicsCuts(){
  // FairRun* fRun = FairRun::Instance();
  
  //Int_t mat = gGeoManager->GetMaterialIndex("HYPdiamond");
  //gMC->Gstpar(mat,"HADR",1.0e-9);
}


// -----   Public method ProcessHits  --------------------------------------

Bool_t PndHyp::ProcessHits(FairVolume* vol) 
  
  
{ 
  
  //fpdgCode = gMC->TrackPid(); 
  Double_t beta, gamma;	TString nam;Int_t nSiL = -1;
  ostringstream FullName,matName;
  
  Int_t medId =  gMC->CurrentMedium();
  TVector3 radt;
  
  
   if(fTrackStopNxtStep) {
   //if( gMC->TrackPid()==3312)//&&(medId==SiId||medId==CId))
   //cout<<"particle  "<< gMC->TrackPid() << "  " <<fTrackStopNxtStep <<endl;
      gMC->StopTrack();
      return kTRUE;
  }
  
   TString nam2 = gMC->CurrentVolName();   
  
 
   /*gMC->TrackMomentum(PiL);
     if (gMC->TrackPid()>1010000000 ||gMC->TrackPid()>1020000000 ||gMC->TrackPid()==-211) cout<<"ProcessHits :  Energy Loss  hyp "<< gMC->TrackPid() << "  "<<vol->getName() <<" "<<PiL.P()<<" "<<gMC->Edep()<<endl;
     if(PiL.P()>3.)cout<<"ProcessHits :  Energy Loss  "<< gMC->TrackPid() << "  " 
     <<vol->getName() <<" "<<PiL.P()<<" "<<gMC->Edep()<<endl;
   */
   
   if (medId==SiId) 
     { //hola
       
       if ( gMC->IsTrackEntering() ) 
	 {
	   fELoss  = 0.;
	   fTime   = gMC->TrackTime() * 1.0e09;
	   fLength = gMC->TrackLength();
	   fmass   = gMC->TrackMass();   // mass (GeV)
	   fcharge = gMC->TrackCharge(); // charge?
	   fpdgCode = gMC->TrackPid(); 
	   fEventID = gMC->CurrentEvent();
	   gMC->TrackPosition(fPosIn);
	   gMC->TrackMomentum(fMomIn);
	   
	   
	 }
       
       // Sum energy loss for all steps in the active volume
       
       fELoss += gMC->Edep();
       
       // Set additional parameters at exit of active volume. Create CbmStsPoint.
       
       TLorentzVector PL;
       gMC->TrackMomentum(PL);
       
       if ( (gMC->IsTrackExiting()    ||
	     gMC->IsTrackStop()       ||
	     gMC->IsTrackDisappeared() )&& gMC->TrackCharge()  ) 
	 {
	   fTrackID  = gMC->GetStack()->GetCurrentTrackNumber();
	   
	   nam = gMC->CurrentVolName();   
	   if ((nam.Contains("Si"))) {
	     sscanf(nam,"stglSi%d#01", &nSiL);
	     // cout << "hyp::ProcessHits> : " << nam <<" # "
	     //    <<nSiL<<" "<<"Hit in "<< gGeoManager->GetPath()<<endl;
	   }
	   
	   //fVolumeID = vol->getMCid();//before it was on
	   //*** now the volume is through the layer number characterised.(X-Z,Z-Y)
	       fVolumeID = nSiL;
	   
	   
	   //**************///
	     
	     //FullName << gGeoManager->GetPath();
	     
	     
	     
	     //cout << "*******  Info from gMC *************" << endl;
	     Int_t cp=-1;
	     Int_t fVolid = gMC->CurrentVolID(cp);
	     
	     
	     // cout << " Vol Name: " << gMC->CurrentVolName() << endl;
	     
	     /*  TString nam2 = gMC->CurrentVolName();   
		 if ((nam2.Contains("Si"))) {
		 sscanf(nam2,"stglSi%d#01", &nSiL);
		 cout << "hyp::ProcessHits> : " << nam2 <<" # "
		 <<nSiL<<" "<<"Hit in "<< gGeoManager->GetPath()<<endl;
		 } */
	     FullName <<gMC->CurrentVolPath();
	     
	     if(0==fGeoH) {
	       std::cout<<" -E- No PndHypGeoHandling loaded."<<std::endl;
	       abort();
	     }
	     
	     
	     gMC->TrackPosition(fPosOut);
	     gMC->TrackMomentum(fMomOut);
	     
	     if (fELoss == 0. ) return kFALSE;
	     
	     radt= fPosOut.Vect();
	     fdist=radt.Perp();
	     beta = fMomOut.Beta();
	     //gamma = fMomOut.Gamma();
	     fPLin = beta;
	     //fPLin = fMomIn.P();
	     fPLout = fMomOut.P();
	     
	     AddHit(fTrackID, fEventID,fVolumeID, 
		    fGeoH->GetID(gMC->CurrentVolPath()),
		    TVector3(fPosIn.X(),   fPosIn.Y(),   fPosIn.Z()),
		    TVector3(fMomIn.Px(),  fMomIn.Py(),  fMomIn.Pz()),
		    TVector3(fPosOut.X(),  fPosOut.Y(),  fPosOut.Z()),
		    TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()),
		    fTime, fLength,fELoss,fcharge,fmass,fpdgCode,
		    fdist,fPLin,fPLout);
	     
	     // Increment number of PndMvd points for TParticle
	     // PndStack* stack = (PndStack*) gMC->GetStack();
	     // stack->AddPoint(kHYP);
	     
	     ResetParameters();
	 }
       
       //return kTRUE;
       
       
       
     }//volSi
   else if ((fpdgCode==3312)&&(nam2.Contains("Ab")))//||nam2.Contains("Si")))
     {  // absorver
       //if ((fpdgCode==3312)&&(medId==CId || medId==SiId))
       
       if ( gMC->IsTrackEntering() ) 
	 {
	   fELoss  = 0.;
	   fTime   = gMC->TrackTime() * 1.0e09;
	   fLength = gMC->TrackLength();
	   fmass   = gMC->TrackMass();   // mass (GeV)
	   fcharge = gMC->TrackCharge(); // charge?
	   fpdgCode = gMC->TrackPid(); 
	   fEventID = gMC->CurrentEvent();
	   gMC->TrackPosition(fPosIn);
	   gMC->TrackMomentum(fMomIn);
	   
	 }
       
       // Sum energy loss for all steps in the active volume
       
       fELoss += gMC->Edep();
     // Gamma, Beta, tau(proper time) of ximnus
       TLorentzVector PL;
       gMC->TrackMomentum(PL);
       beta = PL.Beta();
       
       
       if (beta==0.0) 
	 {    
	   fTrackID  = gMC->GetStack()->GetCurrentTrackNumber();
	   fVolumeID = vol->getMCid();
	   gMC->TrackPosition(fPosOut);
	   gMC->TrackMomentum(fMomOut);
	   
	   //ostringstream matName;
	   TString mat[4]={"CAbs","Si","Be","Al"};
	   if (medId==CId) matName<<"CAbs";
	   if (medId==SiId) matName<<"Si";
	   // cout << "Hit in fullname " << matName.str() <<endl;
	   if (fELoss == 0. ) return kFALSE;
	   
	   //---Kinetic energy:###(PL.P())^2 + Mass^2 -Mass##
	   radt= fPosOut.Vect();
	   fdist=radt.Perp();
	   //beta = fMomOut.Beta();
	   //gamma = fMomOut.Gamma();
	   fPLin = beta;
	   //fPLin = fMomIn.P();
	   fPLout = fMomOut.P();
	   
	   
	   
	   AddSecTarHit(fTrackID, fEventID,fVolumeID,matName.str(),
			TVector3(fPosIn.X(),   fPosIn.Y(),   fPosIn.Z()),
			TVector3(fMomIn.Px(),  fMomIn.Py(),  fMomIn.Pz()),
			TVector3(fPosOut.X(),  fPosOut.Y(),  fPosOut.Z()),
			TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()),
			fTime, fLength,fELoss,fcharge,fmass,fpdgCode,
			fdist,fPLin,fPLout);
	   
	   
	   fTrackStopNxtStep=kTRUE;
	   
	   // Increment number of PndMvd points for TParticle
	   //PndStack* stack = (PndStack*) gMC->GetStack();
	   //stack->AddPoint(kHYP);
	   
	   ResetParameters();
	 }
       
       //for the pipe return kTRUE;
       
       
       
     }//no volSi
   
   
  
   return kTRUE;
  
  
} //ProcessHits

// ----------------------------------------------------------------------------

// -----   Public method EndOfEvent   -----------------------------------------
void PndHyp::EndOfEvent() {
  if (fVerboseLevel)  Print();
  Reset();
}
// ----------------------------------------------------------------------------

// -----   Public method Register   -------------------------------------------
void PndHyp::Register() {
  FairRootManager::Instance()->Register("HypPoint","Hyp", fHypCollection, kTRUE);
  FairRootManager::Instance()->Register("HypSegTarPoint","HypSecTarg", 
				       fHypSecTarCollection, kTRUE);
  
}
// ----------------------------------------------------------------------------

// -----   Public method GetCollection   --------------------------------------
TClonesArray* PndHyp::GetCollection(Int_t iColl) const {
   if (iColl == 0) return fHypCollection;
   if (iColl == 1) return fHypSecTarCollection;
   //if (iColl == 2) return fHypSTpipeCollection;

  return NULL;
}
// ----------------------------------------------------------------------------

// -----   Public method Print   ----------------------------------------------
void PndHyp::Print() const {
    Int_t nHits = fHypCollection->GetEntriesFast();
    cout << "-I- PndHyp: " << nHits << " points registered in this event."
 	<< endl;

    if (fVerboseLevel>1)
      for (Int_t i=0; i<nHits; i++) (*fHypCollection)[i]->Print();
}
// ----------------------------------------------------------------------------



// -----   Public method Reset   ----------------------------------------------
void PndHyp::Reset() {
   fHypCollection->Clear();
   fHypSecTarCollection->Clear();
   // fHypSTpipeCollection->Clear();
 
  fPosIndex = 0;
}
// ----------------------------------------------------------------------------


// guarda in FairRootManager::CopyClones
// -----   Public method CopyClones   -----------------------------------------
void PndHyp::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) {
  Int_t nEntries = cl1->GetEntriesFast();
  //cout << "-I- PndHyp: " << nEntries << " entries to add." << endl;
  TClonesArray& clref = *cl2;
  PndHypPoint* oldpoint = NULL;
  for (Int_t i=0; i<nEntries; i++) {
    oldpoint = (PndHypPoint*) cl1->At(i);
    Int_t index = oldpoint->GetTrackID() + offset;
    oldpoint->SetTrackID(index);
    new (clref[fPosIndex]) PndHypPoint(*oldpoint);
    fPosIndex++;
  }
  cout << " -I- PndHyp: " << cl2->GetEntriesFast() << " merged entries."
       << endl;
}
// ----------------------------------------------------------------------------
 // -----   Public method ConstructGeometry   ----------------------------------
void PndHyp::ConstructGeometry() {
 FairGeoLoader*    geoLoad = FairGeoLoader::Instance();
  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  PndGeoHyp*      hypGeo = new PndGeoHyp();
  hypGeo->setGeomFile(GetGeometryFileName());
  geoFace->addGeoModule(hypGeo);

  Bool_t rc = geoFace->readSet(hypGeo);
  if (rc) hypGeo->create(geoLoad->getGeoBuilder());
  TList* volList = hypGeo->getListOfVolumes();

  // store geo parameter
  FairRun *fRun = FairRun::Instance();
  FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb();
  PndGeoHypPar* par=(PndGeoHypPar*)(rtdb->getContainer("PndGeoHypPar"));
  TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  TObjArray *fPassNodes = par->GetGeoPassiveNodes();

  TListIter iter(volList);
  FairGeoNode* node   = NULL;
  FairGeoVolume *aVol=NULL;

  while( (node = (FairGeoNode*)iter.Next()) ) {
    aVol = dynamic_cast<FairGeoVolume*> ( node );


    if ( node->isSensitive()  ) {
      fSensNodes->AddLast( aVol );
    }else{
      fPassNodes->AddLast( aVol );
    }
  }
    par->setChanged();
    par->setInputVersion(fRun->GetRunId(),1);

  ProcessNodes ( volList );



}
  
 

// -----   Private method AddHit   --------------------------------------------

PndHypPoint* PndHyp::AddHit(Int_t trackID, Int_t evtID,
			    Int_t detID, TString detName,
			    TVector3 pos, TVector3 mom,
			    TVector3 posout, 
			    TVector3 momout,
			    Double_t time,
			    Double_t length, 
			    Double_t eLoss,
			    Double_t charge, Double_t mass,
			    Int_t pdgCode,
			    Double_t dist,
			    Double_t PLin,Double_t PLout) {
  TClonesArray& clref = *fHypCollection;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) PndHypPoint(trackID,  evtID,detID, detName,pos, mom, 
				      posout, momout,
				      time, length, eLoss,charge, 
				      mass,pdgCode,
				      dist,PLin,PLout);
 }



// ----

// -----   Private method AddSecTarHit   --------------------------------------------

PndHypPoint* PndHyp::AddSecTarHit(Int_t trackID, Int_t evtID,
				  Int_t detID, TString detName,
				  TVector3 pos, TVector3 mom,  
				  TVector3 posout, 
				  TVector3 momout, 
				  Double_t time, 
				  Double_t length, 
				  Double_t eLoss,
				  Double_t charge, Double_t mass,
				  Int_t pdgCode,Double_t dist,
				  Double_t PLin,Double_t PLout) {
  TClonesArray& clref = *fHypSecTarCollection;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) PndHypPoint(trackID,  evtID,detID, detName, pos, 
				      mom,  posout, momout,
				      time, length, eLoss,charge, mass, pdgCode,
				      dist,PLin,PLout);
 }



// ----


// ----


ClassImp(PndHyp)