// $Id: DemoTouchable.cc,v 1.7 2014/06/28 15:56:59 kelsey Exp $
////////////////////////////////////////////////////////////////////////
//  File:        DemoTouchable.cc                                     //
//                                                                    //
//  Description: Subclass of G4VTouchable to provide additional data  //
//		 about associated volume.                             //
//		 Adapted from Pedro Arce, GAMOS/GmTouchable           //
//                                                                    //
//  Author:      Michael Kelsey (SLAC)                                //
//  Date:        30 September 2011                                    //
//                                                                    //
//  20111003  M. Kelsey -- Fix bug in diagnostic report of transform  //
//  20111006  M. Kelsey -- Add accessor for full transformations and  //
//		access to volumes.                                    //
//  20120402  M. Kelsey -- Use DemoGeometryTools::delimiter for "/"   //
//  20140627  M. Kelsey -- Add missing G4SystemOfUnits.hh.            //
//  20160204  M. Kelsey -- Can't use typedef name as variable name!   //
//  20180530  M. Kelsey -- Check for null PV in BuildTouchable().     //
//  20190614  Replace ::exit() with DemoAbortJob.                     //
////////////////////////////////////////////////////////////////////////

#include "Demog4base/CDMSTouchable.hh"
#include "Demog4base/CDMSAbortJob.hh"
#include "Demog4base/CDMSGeometryTools.hh"
#include "G4Cons.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4PVParameterised.hh"
#include "G4PVPlacement.hh"
#include "G4Point3D.hh"
#include "G4SystemOfUnits.hh"
#include "G4Tubs.hh"
#include "G4VPVParameterisation.hh"
#include "G4VPhysicalVolume.hh"
#include "G4VSolid.hh"
#include "G4VTouchable.hh"
#include "G4Vector3D.hh"
#include "G4VisAttributes.hh"
#include <vector>


// Constructors

DemoTouchable::CDMSTouchable(std::vector<G4VPhysicalVolume*> pvlist,
			     std::vector<G4int>& ancestorsCopyNo,
			     G4int verbose)
  : verboseLevel(verbose), theCopyNo(0), theVolume(0) {
  BuildTouchable(pvlist, ancestorsCopyNo);
}

DemoTouchable::CDMSTouchable(const G4VTouchable& g4touch, G4int verbose)
  : verboseLevel(verbose), theCopyNo(0), theVolume(0) {
  BuildTouchable(&g4touch);
}

DemoTouchable::CDMSTouchable(const G4VTouchable* g4touch, G4int verbose)
  : verboseLevel(verbose), theCopyNo(0), theVolume(0) {
  BuildTouchable(g4touch);
}

DemoTouchable::~CDMSTouchable() {}


// Return shape associated with this volume

G4VSolid* DemoTouchable::GetSolid() const {
  return theVolume ? theVolume->GetLogicalVolume()->GetSolid() : 0;
}


// Fill with data collected by DemoGeometryTools::GetTouchable()

void DemoTouchable::BuildTouchable(std::vector<G4VPhysicalVolume*> pvlist,
				   std::vector<G4int>& ancestorsCopyNo) {
  if (pvlist.size() == 0) {
    G4cerr << "DemoTouchable::BuildTouchable: empty volume list" << G4endl;
    return;
  }

  // Physical volume is the first one in vector, others are ancestors of it
  theVolume = pvlist[0];
  if (!theVolume) {
    G4cerr << "DemoTouchable::BuildTouchable got null volume" << G4endl;
    return;
  }

  bool isReplica = theVolume->IsReplicated();
  theName = theVolume->GetName();

  if (verboseLevel) 
    G4cout << "DemoTouchable::BuildTouchable " << theName << G4endl;

  theCopyNo = isReplica ? ancestorsCopyNo[0] : theVolume->GetCopyNo();
  theMaterialName = theVolume->GetLogicalVolume()->GetMaterial()->GetName();

  //--- One Exception to the rule: if it is a division of a tubs or cons along kRho, each copy has different solid parameters
  if (isReplica) {
    EAxis axis;
    G4int nReplicas;
    G4double width;
    G4double offset;
    G4bool consuming;
    theVolume->GetReplicationData(axis, nReplicas, width, offset, consuming);
    if (axis == kRho) {
      DemoAbortJob("CDMSTouchable::BuildTouchable",
		   "Division along radius not supported.");
    }
  }

  //------ Get local and global traslation and rotation
  //----- Loop over ancestors to get the global traslation and rotation
  //----- Build the long name
  G4int nAnces = pvlist.size();

  if (verboseLevel>1) G4cout << "nAnces Build " << nAnces << G4endl;

  DemoGeometryTools* geomTool = CDMSGeometryTools::GetInstance();

  G4Transform3D globalTrans = G4Transform3D::Identity;
  theLongName = "";
  G4String ancestorName;

  G4int ii;
  for(ii = 0; ii < nAnces; ii++) {
    G4VPhysicalVolume* pv = pvlist[ii];
    //    G4cout << ii << " trans pv " << pv << G4endl;
    G4int copyNoAncestor =
      pv->IsReplicated() ? ancestorsCopyNo[ii] : pv->GetCopyNo();

    ancestorName = geomTool->delimiter + geomTool->AppendCopyNo(pv->GetName(), copyNoAncestor);

    theLongName.prepend(ancestorName);

    G4Transform3D trans = CalculateTransformation(pv, ancestorsCopyNo[ii]);
    globalTrans = trans * globalTrans;
    //    G4cout << " trans " << trans.getTranslation() << G4endl;

    if (verboseLevel>2) {
      G4cout << " GLOBAL TRANSFORMATION FOR " << ancestorName
	     << "\n TRANS " << globalTrans.getTranslation()
	     << "\n ROT " << globalTrans.getRotation()
	     << "\n CENTRE " << G4Point3D(0.,0.,0.).transform(globalTrans)
	     << "\n X " << G4Vector3D(1.,0.,0.).transform(globalTrans)
	     << "\n Y " << G4Vector3D(0.,1.,0.).transform(globalTrans)
	     << "\n Z " << G4Vector3D(0.,0.,1.).transform(globalTrans)
	     << G4endl;
    }

    if (ii == 0) {
      theLocalPos = trans.getTranslation();
      theLocalRotMat = trans.getRotation();
    }
    
    if (verboseLevel>2) {
      G4cout << "transform rot " << theLocalRotMat
	     << " transl " << theLocalPos << G4endl;
    }
  }
  
  theGlobalPos = globalTrans.getTranslation();
  theGlobalRotMat = globalTrans.getRotation();
  
  if (verboseLevel>1) {
    G4cout << " GLOBAL TRANSFORMATION\n TRANS " << theGlobalPos
	   << "\n ROT " << theGlobalRotMat << G4endl;
  }
}

// Fill with data from GEANT4-created object

void DemoTouchable::BuildTouchable(const G4VTouchable* g4touch) {
  if (!g4touch) return;

  std::vector<G4VPhysicalVolume*> pvlist;
  std::vector<G4int> ancestorsCopyNo;

  G4int nAncest = g4touch->GetHistoryDepth();
  for(G4int ii = 0; ii <= nAncest; ii++) {
    pvlist.push_back(g4touch->GetVolume(ii));
    ancestorsCopyNo.push_back(g4touch->GetCopyNumber(ii));
  }

  BuildTouchable(pvlist, ancestorsCopyNo);
  return;
}


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G4Transform3D DemoTouchable::CalculateTransformation(G4VPhysicalVolume* pv,
						     G4int index) {
  G4Transform3D trans;

  if (verboseLevel) 
    G4cout << " DemoTouchable::CalculateTransformation " << pv->GetName()
	   << " index " << index << G4endl;

  //----- If it is simple placement, just get the translation and 
  if (! pv->IsReplicated()) {
    if (pv->GetMotherLogical() == 0) { //World volume

      if (verboseLevel>2) {
	G4cout << " DemoTouchable::CalculateTransformation "
	       << " No replica: WORLD volume " << G4endl;
      }

      trans = G4Transform3D();
    } else  if (pv->GetRotation() == 0) { 
      trans = G4Transform3D(G4RotationMatrix(), pv->GetObjectTranslation());

      if (verboseLevel>2) {
	G4cout << " DemoTouchable::CalculateTransformation "
	       << " No replica: rotm 0 " << G4endl;
      }
   } else {
      trans = G4Transform3D(pv->GetObjectRotationValue(),
			    pv->GetObjectTranslation());

      if (verboseLevel>2) {
	G4cout << " DemoTouchable::CalculateTransformation "
	       << " No replica: " << G4endl;
      }
    }
    //----- If it is replica 
  } else if (pv->GetParameterisation() == 0){
    EAxis axis;
    G4int nReplicas;
    G4double width;
    G4double offset;
    G4bool consuming;
    pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
    
    //---- Translation is calculated from index of replica
    //---- Rotation is 0 always, except if axis is kPhi     
    G4RotationMatrix rotm;
    G4ThreeVector pos(0.,1.,0.);  //is the one used by kRho, for example
    G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
    G4String sType = solid->GetEntityType();
    //-    G4cout << pv->GetName() << " mother pv " << pv->GetMother() << " sType " << sType << " axis " << axis << G4endl;
    switch (axis){
    case kXAxis:
      pos.setY(0.);
      pos.setX(offset + width * (index+0.5));
      break;
    case kYAxis:
      pos.setY(offset + width * (index+0.5));
      break;
    case kZAxis:
      pos.setY(0.);
      pos.setZ(offset + width * (index+0.5));
      break;
    case kRho:
      //--- If it is a full cone of tube, position is at center, if it is section of cone or tube, position is at radius 
      if (sType == "G4Tubs") { 
	G4Tubs* solTubs = static_cast<G4Tubs*>(solid);
	if ( solTubs->GetDeltaPhiAngle() != 360.*deg) { 
	  //??depends on how you implement this parameterisation	  pos = G4ThreeVector(0.,0.,0.);
	  pos.setMag(offset + width * (index+0.5));
        }
      } else if (sType == "G4Cons") { 
	G4Cons* solCons = static_cast<G4Cons*>(solid);
	if (solCons->GetDeltaPhiAngle() != 360.*deg) { 
	  //??depends on how you implement this parameterisation	  pos = G4ThreeVector(0.,0.,0.);
	  pos.setMag(offset + width * (index+0.5));
        } 
      } else {
	DemoAbortJob("CDMSTouchable::CalculateTransformation",
		     "Replica along kRho for not G4Tubs,G4Cons = "+sType);
      }
      break;
    case kPhi:
      // If it is a cone or tube section, centre is at half between minimum and maximum radii, rotation is that with which the section points to the center
      // If it is the full cone or tube (pathological case where number of divisions is 1) then centre is 0 and rotation identity
      if (sType == "G4Tubs") { 
	G4Tubs* solTubs = static_cast<G4Tubs*>(solid);
	if (solTubs->GetDeltaPhiAngle() != 360.*deg) { 
	  pos = G4ThreeVector(0.,0.,0.);
	  G4double phi =  offset + width * (index+0.5);
	  rotm.rotateZ(phi); 
        }
      } else if (sType == "G4Cons") { 
	G4Cons* solCons = static_cast<G4Cons*>(solid);
	if (solCons->GetDeltaPhiAngle() != 360.*deg) { 
	  pos = G4ThreeVector(0.,0.,0.);
	  G4double phi =  offset + width * (index+0.5);
	  rotm.rotateZ(phi); 
        } 
      } else {
	DemoAbortJob("CDMSTouchable::CalculateTransformation",
		     "Replica along kPhi for not G4Tubs,G4Cons = "+sType);
      }
      break; 
    default:
      break;
    } 

    trans = G4Transform3D(rotm, pos);

    if (verboseLevel>2) {
      G4cout << " DemoTouchable::CalculateTransformation "
	     << " Replica: " << G4endl;
    }

    //----- If it is parameterisation
  } else {
   //--- ComputeTransformation without modifying pv
    G4PVParameterised* pvparam = (G4PVParameterised*)(pv);
    G4VPVParameterisation* param = pvparam->GetParameterisation();

    //create dummy PV just to extract the transformation
    G4RotationMatrix rm0;
    G4PVPlacement* pvDummy =
      new G4PVPlacement(&rm0, G4ThreeVector(), "DUMMY", pv->GetLogicalVolume(),
			0, false, 0);

    param->ComputeTransformation(index, pvDummy);
    trans = G4Transform3D(*(pvDummy->GetRotation()), pvDummy->GetTranslation());
    delete pvDummy;


    if (verboseLevel>2) {
      G4cout << " DemoTouchable::CalculateTransformation "
	     << " parameterisation:  " << G4endl;
    }
  }

  if (verboseLevel>1) {
    G4cout << " DemoTouchable::CalculateTransformation " << index
	   << "\n rotation " << trans.getRotation()
	   << "\n translation " << trans.getTranslation() << G4endl;
  }

  return trans;
}


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
std::ostream& operator<<(std::ostream& os, const DemoTouchable& touch) {
  os << "DemoTouchable: " << touch.GetLongName()
     << " solid type " << touch.GetSolid()->GetEntityType()
     << " material " << touch.GetMaterialName()
     << " Global position " << touch.GetGlobalPosition()
     << " rotation " << touch.GetGlobalRotation()
     << " Local position " << touch.GetLocalPosition();

  return os;
}