// -------------------------------------------------------------------------
// -----                        PndDrc source file                     -----
// -----               Created 11/10/06 by Annalisa Cecchi             -----
// -----                                                               -----
// -------------------------------------------------------------------------



#include <iostream>
#include <valarray>
#include <vector>
#include <string>
#include <list>
using std::list;
#include <limits>
#include <fstream>
#include <utility>
using std::pair;
#include <map>


#include "TClonesArray.h"
#include "TGeoMCGeometry.h"
#include "TGeoManager.h"
#include "TLorentzVector.h"
#include "TParticle.h"
#include "TVirtualMC.h"

#include "Math/Vector3D.h"
using ROOT::Math::XYZVector;

#include "Math/Point3D.h"
using ROOT::Math::XYZPoint;

#include "Math/Transform3D.h"
using ROOT::Math::Transform3D;

#include "Math/RotationX.h"
using ROOT::Math::RotationX;
#include "Math/RotationY.h"
using ROOT::Math::RotationY;
#include "Math/RotationZ.h"
using ROOT::Math::RotationZ;
#include "Math/Rotation3D.h"
using ROOT::Math::Rotation3D;

#include "CbmGeoInterface.h"
#include "CbmGeoLoader.h"
#include "CbmGeoNode.h"
#include "PndGeoDrc.h"
#include "CbmGeoRootBuilder.h"
#include "CbmStack.h"
#include "PndDrc.h"
#include "PndDrcPoint.h"
#include "CbmRootManager.h"
#include "CbmVolume.h"

// add on for debug
#include "CbmGeoG3Builder.h"
#include "CbmRuntimeDb.h"
#include "PndGeoDrcPar.h"
#include "TObjArray.h"
#include "CbmRun.h"


#include "PndDrcPhoton.h"
#include "PndDrcOptReflAbs.h"
#include "PndDrcSurfAbs.h"
#include "PndDrcSurfPolyFlat.h"
#include "PndDrcOptReflSilver.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatLithotecQ0.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptVol.h"
#include "PndDrcOptDevManager.h"
#include "PndDrcUtil.h"


// -----   Default constructor   -------------------------------------------
PndDrc::PndDrc() {
  fDrcCollection = new TClonesArray("PndDrcPoint");
  fPosIndex   = 0;
  fflag = 0;
  volDetector = 0;
  
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
PndDrc::PndDrc(const char* name, Bool_t active)
  : CbmDetector(name, active) {
    fDrcCollection = new TClonesArray("PndDrcPoint");
    fPosIndex   = 0;
    fflag   = 0;
    volDetector = 0;
    fEventID    = -1;
    fSenId1 =0;
    fSenId2 =0;
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
PndDrc::~PndDrc() {
  if (fDrcCollection) {
    fDrcCollection->Delete();
    delete fDrcCollection;
  }
  
}
// -------------------------------------------------------------------------



// -----   Public method Intialize   ---------------------------------------
void PndDrc::Initialize() {
  

  CbmDetector::Initialize();
  CbmRun       *sim  = CbmRun::Instance();
  CbmRuntimeDb *rtdb = sim->GetRuntimeDb();
  PndGeoDrcPar *par  = (PndGeoDrcPar*)(rtdb->getContainer("PndGeoDrcPar"));


  // create a clone of each bar, which will handle the photon propagation.
  //
  //         p4----------p7
  //        /|          /|                 y
  //       / |         / |                 |  z
  //      /  p5-------/--p6                | /
  //     /  /        /  /                  |/
  //    /  /        /  /             x ----0
  //   /  /        /  /
  //  /  /        /  /
  // p0---------p3  /
  // | /         | /
  // |/          |/
  // p1---------p2 






  int verbosity = 0; // 0=quiet ... 5=talk too much

  XYZPoint p[8];   // 8 space points of radiator bar in sequence of above

  PndDrcOptDevManager* manager = new PndDrcOptDevManager();
  manager->setVerbosity(verbosity);


  TObjArray    *sensNodes = par->GetGeoSensitiveNodes();

  CbmGeoNode *fm1= (CbmGeoNode *) sensNodes->FindObject("drc05"); //screen at the end of the bar
  //CbmGeoNode *fm2= (CbmGeoNode *) sensNodes->FindObject("drc06"); //screen far from the bar
  fSenId1=fm1->getMCid();
  //fSenId2=fm2->getMCid();
  

  for( int inode=0; inode<sensNodes->GetEntries(); inode++) 
    {// for inode=
      CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (sensNodes->At(inode));	
      if ( !node ) continue;
      
      TString name = node->getName();
      TString shapeName = node->getShapePointer()->GetName();

      TString name_clone = "none";

      cout<<" node,name "<<inode<<" "<<name.Data()<<endl;//###

      if (shapeName.CompareTo("BOX"))
	{//if box
	  name_clone  = name;
	  name_clone += "_clone";
	    

	  CbmGeoTransform trans = node->getTransform();

	  for (int ii=0; ii<node->getNumPoints(); ii++)
	    {
	      CbmGeoVector vec = *(node->getPoint(ii));
	      CbmGeoVector vec1 = trans.transFrom(vec);
	      p[ii].SetXYZ(vec1.X(),vec1.Y(),vec1.Z());
	      //cout<<" points: "<<vec1.X()<<" "<<vec1.Y()<<" "<<vec1.Z()<<endl;
	    }

	  // assemble bar

	  // Declare 6 flat surfaces with 4 edge points.
	  PndDrcSurfPolyFlat a1,a2,a3,a4,a5,a6;

	  a1.setVerbosity(verbosity);
	  a1.addPoint(p[0]);
	  a1.addPoint(p[1]);
	  a1.addPoint(p[2]);
	  a1.addPoint(p[3]);
	  a1.setPixel();                // screen
	  a1.setName("a_upstream");

	  a2.setVerbosity(verbosity);
	  a2.addPoint(p[1]);
	  a2.addPoint(p[5]);
	  a2.addPoint(p[6]);
	  a2.addPoint(p[2]);
	  a2.setName("a_side1");

	  a3.setVerbosity(verbosity);
	  a3.addPoint(p[0]);
	  a3.addPoint(p[4]);
	  a3.addPoint(p[5]);
	  a3.addPoint(p[1]);
	  a3.setName("a_side2");

	  a4.setVerbosity(verbosity);
	  a4.addPoint(p[3]);
	  a4.addPoint(p[2]);
	  a4.addPoint(p[6]);
	  a4.addPoint(p[7]);
	  a4.setName("a_side3");

	  a5.setVerbosity(verbosity);
	  a5.addPoint(p[0]);
	  a5.addPoint(p[3]);
	  a5.addPoint(p[7]);
	  a5.addPoint(p[4]);
	  a5.setName("a_side4");
	  
	  PndDrcOptReflSilver refl;

	  a6.setVerbosity(verbosity);
	  a6.setReflectivity(refl);   // mirror
	  a6.addPoint(p[7]);
	  a6.addPoint(p[6]);
	  a6.addPoint(p[5]);
	  a6.addPoint(p[4]);
	  a6.setName("a_downstream");

	  // create a volume consiting of surfaces
	  // create a material the bar will consist of
	  PndDrcOptVol bar;
	  PndDrcOptMatLithotecQ0 quartz;
	  bar.setVerbosity(verbosity);
	  bar.setOptMaterial(quartz);
	  bar.addSurface(a1);
	  bar.addSurface(a2);
	  bar.addSurface(a3);
	  bar.addSurface(a4);
	  bar.addSurface(a5);
	  bar.addSurface(a6);
	  bar.setName("bar");
	  //bar.setCopyNumber(inode);

	  // create a screen (simplified photon detector) where photons end and 
	  // get the status Drc::Measured (no further propagation). For this you need 
	  // one single plane.
	  

	  // Build a optical system consisting out of several 
	  // volumes, mirrors and screens. This layer has the advantage, that 
	  // a device consisting out of many equal subsystems
	  // like a bar box, easily can be reproduced. See one of the 
	  // test examples in ./drcprop.

	  PndDrcOptDevSys opt_system;
	  opt_system.setVerbosity(verbosity);
	  opt_system.addDevice(bar);
	  opt_system.setNameCopyNumber(name_clone.Data(),inode);
	  manager->addDeviceSystem(opt_system);

	  
	}



    }
  // print out for debugging all names and pointers stored in manager
  // manager->print();

}

// -------------------------------------------------------------------------
void PndDrc::BeginEvent(){
  cout<<" \n\n>>>>>>>>>>>>>>>>>>>>new event "<<endl;
  fEventID++;
}



// -----   Public method ProcessHits  --------------------------------------
Bool_t PndDrc::ProcessHits(CbmVolume* vol) {

  //cout << "PndDrc::ProcessHits " << vol->GetName() << endl;


  if (vol->getMCid()==fSenId2 && gMC->IsTrackEntering()) 
    {


      TString nam = gMC->CurrentVolName();
      Int_t  fpdgCode = gMC->TrackPid(); 
 
      //Int_t ntr = gMC->GetStack()->GetNtrack();
      Int_t  fEventID = gMC->CurrentEvent();
      Int_t  fCopyNo = vol->getCopyNo();
      //gMC->CurrentVolID(fCopyNo);
      fflag = 1;
  
      fmass   = gMC->TrackMass();   // mass (GeV)
      fcharge = gMC->TrackCharge(); // charge?
    
   
      if (fpdgCode == 50000050) 
	{ //Cherenkov photons
        
	  Int_t  fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); //track ID     
	  Double_t fTime    = gMC->TrackTime() * 1.0e09;
	  Double_t fLength = gMC->TrackLength();
	  Double_t fELoss = gMC->Edep();
	  TLorentzVector fPos, fMom;
	  gMC->TrackPosition(fPos);
	  gMC->TrackMomentum(fMom); // GeV/c
      
	  AddHit(fTrackID,
		 fCopyNo,
		 TVector3(fPos.X(),   fPos.Y(),   fPos.Z()),
		 TVector3(fMom.Px(),  fMom.Py(),  fMom.Pz()),
		 fTime,
		 fLength,
		 fELoss,
		 fcharge,
		 fmass,
		 fflag,
		 nam,
		 fpdgCode,
		 fEventID);
	}  
    }
   

  ResetParameters();

  return kTRUE; 
  


}

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



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



// -----   Public method Register   -------------------------------------------
void PndDrc::Register() {
  CbmRootManager::Instance()->Register("PndDrcPoint","drc", fDrcCollection, kTRUE);
}
// ----------------------------------------------------------------------------



// -----   Public method GetCollection   --------------------------------------
TClonesArray* PndDrc::GetCollection(Int_t iColl) const {
   if (iColl == 0) return fDrcCollection;
  return NULL;
}
// ----------------------------------------------------------------------------




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

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



// -----   Public method Reset   ----------------------------------------------
void PndDrc::Reset() {
   fDrcCollection->Clear();
 
  fPosIndex = 0;
}
// ----------------------------------------------------------------------------


// guarda in CbmRootManager::CopyClones
// -----   Public method CopyClones   -----------------------------------------
void PndDrc::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) {
  Int_t nEntries = cl1->GetEntriesFast();
  cout << "-I- PndDrc: " << nEntries << " entries to add." << endl;
  TClonesArray& clref = *cl2;
  PndDrcPoint* oldpoint = NULL;
  for (Int_t i=0; i<nEntries; i++) {
    oldpoint = (PndDrcPoint*) cl1->At(i);
    Int_t index = oldpoint->GetTrackID() + offset;
    oldpoint->SetTrackID(index);
    new (clref[fPosIndex]) PndDrcPoint(*oldpoint);
    fPosIndex++;
  }
  cout << " -I- PndDrc: " << cl2->GetEntriesFast() << " merged entries."
       << endl;
}
// ----------------------------------------------------------------------------



// -----   Public method ConstructGeometry   ----------------------------------
void PndDrc::ConstructGeometry() {

  CbmGeoLoader*    drcgeoLoad = CbmGeoLoader::Instance();
  CbmGeoInterface* drcgeoFace = drcgeoLoad->getGeoInterface();
  PndGeoDrc*      drcGeo = new PndGeoDrc();
  drcGeo->setGeomFile(GetGeometryFileName());
  drcgeoFace->addGeoModule(drcGeo);

  Bool_t rc = drcgeoFace->readSet(drcGeo);
  if (rc) drcGeo->create(drcgeoLoad->getGeoBuilder());
  TList* volList = drcGeo->getListOfVolumes();

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

  TListIter iter(volList);
  CbmGeoNode* node   = NULL;
  CbmGeoVolume *aVol=NULL;

  while( (node = (CbmGeoNode*)iter.Next()) ) {
    aVol = dynamic_cast<CbmGeoVolume*> ( node );
    if ( node->isSensitive()  ) {
      fSensNodes->AddLast( aVol );
    }else{
      fPassNodes->AddLast( aVol );
    }
  }

  ProcessNodes ( volList );
  par->setChanged();
  par->setInputVersion(fRun->GetRunId(),1);

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


// -----   Private method AddHit   --------------------------------------------
PndDrcPoint* PndDrc::AddHit(Int_t trackID, Int_t copyNo, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eLoss, Double_t charge, Double_t mass, Bool_t flag, TString nam, Int_t pdgCode, Int_t eventID) {
 
  TClonesArray& clref = *fDrcCollection;
  Int_t size = clref.GetEntriesFast();
  if (fVerboseLevel>1) 
    cout << "-I- PndDrc: Adding Point at (" << pos.X() << ", " << pos.Y() 
	 << ", " << pos.Z() << ") cm, detector " << copyNo << ", flag "<<flag<<", track "
	 << trackID <<" event "<<eventID << ", energy loss " << eLoss*1e06 << " keV" << endl;
  return new(clref[size]) PndDrcPoint(trackID, 
				      copyNo, 
				      pos, 
				      mom, 
				      time, 
				      length, 
				      eLoss, 
				      charge,
                                      mass, 
				      flag, 
				      nam,
				      pdgCode,
				      eventID);


}


// ----


ClassImp(PndDrc)