From 1abe3fa4454549bcca9d966c063ed65ad7d025f1 Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Thu, 28 Jun 2012 11:37:14 +0000
Subject: [PATCH] Bug fix for detector element transformations

---
 DDCore/include/DD4hep/GeoHandler.h        |  15 ++
 DDCore/src/Detector.cpp                   | 212 ++++++++++------------
 DDCore/src/Geant4Converter.cpp            |  36 ++--
 DDCore/src/LCDDImp.cpp                    |   6 +-
 DDCore/src/LCDDImp.h                      |  10 +-
 DDCore/src/compact/Compact2Objects.cpp    |   2 +-
 DDExamples/CLICSiD/src/EcalBarrel_geo.cpp |   2 +-
 DDExamples/ILDExDet/src/TPCModule.cpp     |  67 +++++--
 8 files changed, 192 insertions(+), 158 deletions(-)

diff --git a/DDCore/include/DD4hep/GeoHandler.h b/DDCore/include/DD4hep/GeoHandler.h
index dc8170360..d50f2b12f 100644
--- a/DDCore/include/DD4hep/GeoHandler.h
+++ b/DDCore/include/DD4hep/GeoHandler.h
@@ -73,6 +73,21 @@ namespace DD4hep {
       GeoHandler& collect(DetElement top, GeometryInfo& info);
       /// Access to collected node list
       Data* release();
+#if  0
+      template <typename C, typename F> static void _handle(const O* ptr, const C& c, F pmf)  {
+	for(typename C::const_iterator i=c.begin(); i != c.end(); ++i)   {
+	  (ptr->*pmf)((*i).first, (*i).second);
+	}
+      }
+#endif
+      template <typename O, typename C, typename F> static void handle(const O* o, const C& c, F pmf)    {
+	for(typename C::const_iterator i=c.begin(); i != c.end(); ++i)
+	  (o->*pmf)((*i).first, (*i).second);
+      }
+      template <typename O, typename C, typename F> static void handle(O* o, const C& c, F pmf)     {
+	for(typename C::const_iterator i=c.begin(); i != c.end(); ++i)
+	  (o->*pmf)((*i).first, (*i).second);
+      }
     };
 
     struct GeoScan {
diff --git a/DDCore/src/Detector.cpp b/DDCore/src/Detector.cpp
index e77e86ad1..96c2ff28e 100644
--- a/DDCore/src/Detector.cpp
+++ b/DDCore/src/Detector.cpp
@@ -12,137 +12,95 @@
 #include "TGeoVolume.h"
 #include "TGeoMatrix.h"
 #include "TGeoManager.h"
+#include <iostream>
 
 using namespace std;
 using namespace DD4hep::Geometry;
 
-static bool traverse_find(TGeoNode* parent, TGeoNode* child, vector<TGeoMatrix*>& trafos) {
-  TIter next(parent->GetVolume()->GetNodes());
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    if ( daughter == child )   {
-      trafos.push_back(daughter->GetMatrix());
+static bool find_child(TGeoNode* parent, TGeoNode* child, vector<TGeoNode*>& path) {
+  if ( parent && child ) {
+    if ( parent == child ) {
+      path.push_back(child);
       return true;
     }
-  }
-  next.Reset();
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    if ( traverse_find(daughter, child, trafos) ) {
-      trafos.push_back(daughter->GetMatrix());
-      return true;
+    TIter next(parent->GetVolume()->GetNodes());
+    for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
+      if ( daughter == child )   {
+	path.push_back(daughter);
+	return true;
+      }
+    }
+    next.Reset();
+    for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
+      bool res = find_child(daughter, child, path);
+      if ( res ) {
+	path.push_back(daughter);
+	return res;
+      }
     }
   }
   return false;
 }
 
-static bool traverse_up(const DetElement& parent, const DetElement& child, vector<TGeoMatrix*>& trafos) {
-  if ( child.ptr() != parent.ptr() ) {
-    for( DetElement par=parent, cld=child; par.isValid(); cld=par, par=par._data().parent ) {
-      PlacedVolume cld_place = cld._data().placement;
-      PlacedVolume par_place = par._data().placement;
-      if ( !traverse_find(par_place.ptr(),cld_place.ptr(),trafos) ) {
-	return false;
-      }
-    }
+static bool collect_detector_nodes(const vector<TGeoNode*>& det_nodes, vector<TGeoNode*>& nodes) {
+  if ( det_nodes.size() < 1 ) {
+    return false;
   }
-  return true;
-}
-#include <iostream>
-static string find_path(TGeoNode* top, TGeoNode* child) {
-  if ( top == child ) {
-    return child->GetName();
-  }
-  TIter next(top->GetVolume()->GetNodes());
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    if ( daughter == child )   {
-      return child->GetName();
-    }
+  if ( det_nodes.size() < 2 ) {
+    return true;
   }
-  next.Reset();
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    string res = find_path(daughter, child);
-    if ( !res.empty() ) {
-      return top->GetName() + ("/"+res);
+  for(size_t i=0, n=det_nodes.size(); i<n-1; ++i) {
+    if ( !find_child(det_nodes[i+1],det_nodes[i],nodes) )  {
+      return false;
     }
   }
-  return "";
+  return true;
 }
 
-static string find_child(TGeoNode* top, TGeoNode* child, vector<TGeoNode*>& path) {
-  if ( top == child ) {
-    path.push_back(child);
-    return child->GetName();
+/// Create cached matrix to transform to positions to an upper level DetElement
+static TGeoHMatrix* create_trafo(const vector<TGeoNode*>& det_nodes)   {
+  if ( det_nodes.size() < 2 ) {
+    return new TGeoHMatrix(*gGeoIdentity);
   }
-  TIter next(top->GetVolume()->GetNodes());
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    if ( daughter == child )   {
-      path.push_back(daughter);
-      return child->GetName();
-    }
+  vector<TGeoNode*> nodes;
+  if ( !collect_detector_nodes(det_nodes,nodes) ) {
+    throw runtime_error("DetElement cannot connect "+string(det_nodes[0]->GetName())+
+			" to child "+string(det_nodes[1]->GetName()));
   }
-  next.Reset();
-  for (TGeoNode *daughter=(TGeoNode*)next(); daughter; daughter=(TGeoNode*)next() ) {
-    string res = find_child(daughter, child, path);
-    if ( !res.empty() ) {
-      path.push_back(daughter);
-      return top->GetName() + ("/"+res);
-    }
+  TGeoHMatrix* mat = new TGeoHMatrix(*gGeoIdentity);
+  for(vector<TGeoNode*>::const_iterator i=nodes.begin(); i!=nodes.end(); ++i) {
+    TGeoMatrix* m = (*i)->GetMatrix();
+    mat->MultiplyLeft(m);
   }
-  //cout << "FAILED child:" << (unsigned long)child << endl;
-  return "";
+  return mat;
 }
 
-/// Create cached matrix to transform to positions to an upper level DetElement
-static TGeoHMatrix* create_trafo(const Ref_t& parent, const Ref_t& child)   {
-  if ( parent.isValid() && child.isValid() )   {
-    TGeoHMatrix* mat = 0;    // Now collect all transformations
-    vector<TGeoMatrix*> trafos;
-    if ( !traverse_up(parent,child,trafos) )   {
-      trafos.clear();
-      if ( !traverse_up(child,parent,trafos) )   {  // Some error, non-existing path!
-	throw runtime_error("DetElement "+string(parent.name())+" is not connected to child "+string(child.name())+" [Geo-Error]");
-      }
-      else {
-	/// The two detector elements were only found by reversing the hierarchy.
-	/// Hence the transformations must be applied in reverse order.
-	for(vector<TGeoMatrix*>::const_reverse_iterator i=trafos.rbegin(); i!=trafos.rend(); ++i) {
-	  if ( !mat ) mat = new TGeoHMatrix(*(*i));
-	  else mat->MultiplyLeft(*i);
-	}
-      }
-    }
-    else {
-      /// Combine the transformations to a single transformation.
-      for(vector<TGeoMatrix*>::const_iterator i=trafos.begin(); i!=trafos.end(); ++i) {
-	if ( !mat ) mat = new TGeoHMatrix(*(*i));
-	else mat->MultiplyLeft(*i);
-      }
+/// Top detector element
+static DetElement _top(DetElement o, vector<TGeoNode*>& det_nodes)   {
+  DetElement par = o, pp = o;
+  while( par.isValid() ) {
+    if ( par.placement().isValid() )   {
+      det_nodes.push_back(par.placement().ptr());
+      pp = par;
     }
-#if 0
-    // Test: find arbitrary child by traversing TGeoNode tree
-    vector<TGeoNode*> path;
-    TGeoNode* top  = gGeoManager->GetTopNode();
-    TGeoNode* node = DetElement(child).placement().ptr();
-    string res = find_child(top,node,path);
-    cout << "Path:" << res << endl;
-#endif
-    return mat;
+    par = par.parent();
   }
-  if ( parent.isValid() )
-    throw runtime_error("DetElement cannot connect "+string(parent.name())+" to not-existing child!");
-  else if ( child.isValid() )
-    throw runtime_error("DetElement cannot connect "+string(child.name())+" to not-existing parent!");
-  else
-    throw runtime_error("DetElement cannot connect nonexisting parent to not-existing child!");
+  return pp;
 }
 
 /// Top detector element
-static DetElement _top(DetElement o)   {
+static DetElement _par(DetElement o, DetElement top, vector<TGeoNode*>& det_nodes)   {
   DetElement par = o, pp = o;
-  while(par.isValid()) {
-    if ( par.placement().isValid() ) pp = par;
+  while( par.isValid() ) {
+    if ( par.placement().isValid() )   {
+      det_nodes.push_back(par.placement().ptr());
+      pp = par;
+    }
+    if ( par.ptr() == top.ptr() ) break;
     par = par.parent();
   }
-  return pp;
+  if ( pp.ptr() == top.ptr() ) return pp;
+  return DetElement();
 }
 
 /// Default constructor
@@ -200,19 +158,9 @@ DetElement::Object::operator Ref_t() {
 /// Create cached matrix to transform to world coordinates
 TGeoMatrix* DetElement::Object::worldTransformation() {
   if ( !worldTrafo ) {
-    DetElement top_det = _top(DetElement(asRef()));
-    vector<TGeoMatrix*> trafos;
-    TGeoHMatrix* mat = create_trafo(top_det,asRef());
-    // Now we got the point in the top-most detector element. We now have
-    // to translate this to the "world volume", which has no transformation matrix anymore
-    TGeoNode* top_node = gGeoManager->GetTopNode();
-    PlacedVolume place = top_det.placement();
-    if ( !traverse_find(top_node, place.ptr(), trafos) ) {
-      // Some error, non-existing path!
-      throw runtime_error("DetElement "+string(parent.name())+" is not connected to top geo node [Geo-Error]");
-    }    
-    for(vector<TGeoMatrix*>::const_iterator i=trafos.begin(); i!=trafos.end(); ++i)
-      mat->MultiplyLeft(*i);
+    vector<TGeoNode*> det_nodes;
+    _top(DetElement(asRef()),det_nodes);
+    TGeoHMatrix* mat = create_trafo(det_nodes);
     worldTrafo = mat;
   }    
   return worldTrafo;
@@ -221,7 +169,10 @@ TGeoMatrix* DetElement::Object::worldTransformation() {
 /// Create cached matrix to transform to parent coordinates
 TGeoMatrix* DetElement::Object::parentTransformation() {
   if ( !parentTrafo )   {
-    parentTrafo = create_trafo(this->parent,asRef());
+    vector<TGeoNode*> det_nodes;
+    det_nodes.push_back(placement.ptr());
+    det_nodes.push_back(DetElement(parent).placement().ptr());
+    parentTrafo = create_trafo(det_nodes);
   }
   return parentTrafo;
 }
@@ -229,7 +180,24 @@ TGeoMatrix* DetElement::Object::parentTransformation() {
 /// Create cached matrix to transform to reference coordinates
 TGeoMatrix* DetElement::Object::referenceTransformation() {
   if ( !referenceTrafo )   {
-    referenceTrafo = create_trafo(this->reference,asRef());
+    vector<TGeoNode*> nodes;
+    DetElement ref(reference);
+    DetElement self(asRef());
+    DetElement elt = _par(self,ref,nodes);
+    if ( elt.isValid() )   {
+      referenceTrafo = create_trafo(nodes);
+    }
+    else   {
+      nodes.clear();
+      DetElement elt = _par(ref,self,nodes);
+      if ( !elt.isValid() )   {
+	throw runtime_error("referenceTransformation: No path from "+string(self.name())+
+			    " to reference element "+string(ref.name())+" present!");
+      }
+      TGeoMatrix* m = create_trafo(nodes);
+      referenceTrafo = new TGeoHMatrix(m->Inverse());
+      delete m;
+    }
   }
   return referenceTrafo;
 }
@@ -258,8 +226,16 @@ std::string DetElement::placementPath() const {
   if ( isValid() ) {
     Object& o = _data();
     if ( o.placementPath.empty() ) {
-      DetElement top =_top(*this);
-      o.placementPath = find_path(top.placement().ptr(),placement().ptr());
+      string res = "";
+      vector<TGeoNode*> nodes, path;
+      _top(*this,nodes);
+      if ( !collect_detector_nodes(nodes,path) ) {
+	throw runtime_error("DetElement cannot determine placement path of "+string(name()));
+      }
+      path.push_back(gGeoManager->GetTopNode());
+      for(vector<TGeoNode*>::const_reverse_iterator i=path.rbegin(); i!=path.rend(); ++i)
+	res += (string("/")+(*i)->GetName());
+      o.placementPath = res;
     }
     return o.placementPath;
   }
diff --git a/DDCore/src/Geant4Converter.cpp b/DDCore/src/Geant4Converter.cpp
index f57928b91..ada85eb5d 100644
--- a/DDCore/src/Geant4Converter.cpp
+++ b/DDCore/src/Geant4Converter.cpp
@@ -38,11 +38,10 @@ using namespace std;
 namespace {
   static string indent = "";
   struct MyTransform3D : public G4Transform3D {
-    MyTransform3D() : G4Transform3D() {}
-    // Set transformation matrix
-    void set(double XX, double XY, double XZ, double DX,
+    MyTransform3D(double XX, double XY, double XZ, double DX,
 	     double YX, double YY, double YZ, double DY,
-	     double ZX, double ZY, double ZZ, double DZ) {}
+	     double ZX, double ZY, double ZZ, double DZ) : G4Transform3D() 
+    {}
     //  { G4Transform3D::setTransform(XX,XY,XZ,DX,YX,YY,YZ,DY,ZX,ZY,ZZ,DZ);  }
   };
 }
@@ -166,12 +165,12 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape)
       const TGeoCompositeShape* s = (const TGeoCompositeShape*)shape;
       const TGeoBoolNode* boolean = s->GetBoolNode();
       TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
-      TGeoMatrix* m   = boolean->GetRightMatrix();
-      G4VSolid* left  = (G4VSolid*)handleSolid(name+"_left", boolean->GetLeftShape());
-      G4VSolid* right = (G4VSolid*)handleSolid(name+"_right",boolean->GetRightShape());
+      TGeoMatrix* m     = boolean->GetRightMatrix();
+      G4VSolid* left    = (G4VSolid*)handleSolid(name+"_left", boolean->GetLeftShape());
+      G4VSolid* right   = (G4VSolid*)handleSolid(name+"_right",boolean->GetRightShape());
       const Double_t *t = m->GetTranslation();
       const Double_t *r = m->GetRotationMatrix();
-
+      
       if ( !left )   {
 	throw runtime_error("G4Converter: No left Geant4 Solid present for composite shape:"+name);
       }
@@ -180,10 +179,9 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape)
       }
 
       if ( m->IsRotation()    )   {
-	MyTransform3D transform;
-	transform.set(r[0],r[1],r[2],t[0],
-		      r[3],r[4],r[5],t[1],
-		      r[6],r[7],r[8],t[3]);
+	MyTransform3D transform(r[0],r[1],r[2],t[0],
+				r[3],r[4],r[5],t[1],
+				r[6],r[7],r[8],t[3]);
 	if (      oper == TGeoBoolNode::kGeoSubtraction )
 	  solid = new G4SubtractionSolid(name,left,right,transform);
 	else if ( oper == TGeoBoolNode::kGeoUnion )
@@ -247,11 +245,10 @@ void* Geant4Converter::handlePlacement(const std::string& name, const TGeoNode*
       G4LogicalVolume*   g4vol  = data().g4Volumes[node->GetVolume()];
       G4LogicalVolume*   g4mot  = data().g4Volumes[node->GetMotherVolume()];
 
-      if ( trafo->IsRotation() )  {
-	MyTransform3D transform;
-	transform.set(rot[0],rot[1],rot[2],trans[0],
-		      rot[3],rot[4],rot[5],trans[1],
-		      rot[6],rot[7],rot[8],trans[3]);
+      if ( trafo->IsRotation() )    {
+	MyTransform3D transform(rot[0],rot[1],rot[2],trans[0],
+				rot[3],rot[4],rot[5],trans[1],
+				rot[6],rot[7],rot[8],trans[3]);
 	g4pv = new G4PVPlacement(transform, // no rotation
 				 g4vol,     // its logical volume
 				 name,      // its name
@@ -299,10 +296,13 @@ void Geant4Converter::create(DetElement top) {
   m_data->clear();
   collect(top,geo);
 
+  // We do not have to handle defines etc.
+  // All positions and the like are not really named.
+  // Hence, start creating the G4 objects for materials, solids and log volumes.
   handle(this, geo.materials, &Geant4Converter::handleMaterial);
   handle(this, geo.solids,    &Geant4Converter::handleSolid);
   handle(this, geo.volumes,   &Geant4Converter::handleVolume);
-
+  // Now place all this stuff appropriately
   for(Data::const_reverse_iterator i=m_data->rbegin(); i != m_data->rend(); ++i)   {
     const Data::mapped_type& v = (*i).second;
     for(Data::mapped_type::const_iterator j=v.begin(); j != v.end(); ++j) {
diff --git a/DDCore/src/LCDDImp.cpp b/DDCore/src/LCDDImp.cpp
index 3e0796bd5..6a29241c1 100644
--- a/DDCore/src/LCDDImp.cpp
+++ b/DDCore/src/LCDDImp.cpp
@@ -52,7 +52,7 @@ Volume LCDDImp::pickMotherVolume(const DetElement&) const  {     // throw if not
 }
 
 LCDD& LCDDImp::addDetector(const Ref_t& x)    { 
-  m_detectors.append_noCheck(x);
+  m_detectors.append(x);
   m_world.add(DetElement(x));
   return *this;
 }
@@ -147,6 +147,7 @@ void LCDDImp::endDocument()  {
   /// Since we allow now for anonymous shapes,
   /// we will rename them to use the name of the volume they are assigned to
   TGeoManager* mgr = gGeoManager;
+  gGeoManager->SetTopVolume(m_worldVol);
   mgr->CloseGeometry();
   m_world.setPlacement(PlacedVolume(mgr->GetTopNode()));
   ShapePatcher(m_world)();
@@ -172,9 +173,8 @@ void LCDDImp::init()  {
   m_trackingVol    = tracking;
   m_materialAir    = material("Air");
   m_materialVacuum = material("Vacuum");
-  m_detectors.append_noCheck(m_world);
+  m_detectors.append(m_world);
   m_world.add(m_trackers);
-  gGeoManager->SetTopVolume(m_worldVol);
 }
 
 void LCDDImp::fromCompact(const std::string& xmlfile) {
diff --git a/DDCore/src/LCDDImp.h b/DDCore/src/LCDDImp.h
index 58775a4e4..72f7f5efc 100644
--- a/DDCore/src/LCDDImp.h
+++ b/DDCore/src/LCDDImp.h
@@ -9,12 +9,14 @@
 
 #ifndef DD4hep_LCDDGEOIMP_H
 #define DD4hep_LCDDGEOIMP_H
+
+// Framework include files
 #include "DD4hep/LCDD.h"
-//#include "XML/XMLElements.h"
 
+// Forward declarations
 class TGeoManager;
 
-// C++ include files
+// C/C++ include files
 #include <map>
 
 /*
@@ -35,6 +37,7 @@ namespace DD4hep {
       
       struct ObjectHandleMap : public HandleMap  {
         ObjectHandleMap() {}
+#if 0
         void append_noCheck(const Ref_t& e) { 
           if ( e.isValid() )  {
             std::string n = e.name();
@@ -44,6 +47,7 @@ namespace DD4hep {
             this->insert(std::make_pair(n,e.ptr()));
           }
         }
+#endif
         void append(const Ref_t& e, bool throw_on_doubles=true) { 
           if ( e.isValid() )  {
             std::string n = e.name();
@@ -163,7 +167,7 @@ namespace DD4hep {
       virtual LCDD& addReadout(const Ref_t& x)          { m_readouts.append(x);   __R;}
       virtual LCDD& addVisAttribute(const Ref_t& x)     { m_display.append(x);    __R;}
       virtual LCDD& addSensitiveDetector(const Ref_t& x){ m_sensitive.append(x);  __R;}
-      virtual LCDD& addDetector(const Ref_t& x);    //  { m_detectors.append_noCheck(x);  __R;}
+      virtual LCDD& addDetector(const Ref_t& x);
 
       virtual LCDD& addAlignment(const Ref_t& x)        { m_alignments.append(x); __R;}
 #undef __R
diff --git a/DDCore/src/compact/Compact2Objects.cpp b/DDCore/src/compact/Compact2Objects.cpp
index 3687af7ff..783fd0568 100644
--- a/DDCore/src/compact/Compact2Objects.cpp
+++ b/DDCore/src/compact/Compact2Objects.cpp
@@ -366,7 +366,7 @@ namespace DD4hep { namespace Geometry {
       e.second->SetTitle(parent.isValid() ? parent.type().c_str() : e.first.c_str());
     }
     else {
-      cout << "Title present: " << e.second->GetTitle() << endl;
+      //cout << "Title present: " << e.second->GetTitle() << endl;
     }
     for_each(children.begin(),children.end(),setChildTitles);
   }
diff --git a/DDExamples/CLICSiD/src/EcalBarrel_geo.cpp b/DDExamples/CLICSiD/src/EcalBarrel_geo.cpp
index 0987eb3ff..8c206e061 100644
--- a/DDExamples/CLICSiD/src/EcalBarrel_geo.cpp
+++ b/DDExamples/CLICSiD/src/EcalBarrel_geo.cpp
@@ -42,7 +42,7 @@ static Ref_t create_detector(LCDD& lcdd, const xml_h& e, SensitiveDetector& sens
   sdet.setPlacement(env_phv);
 
   DetElement    stave_det("stave0",det_id);
-  double dx        = mod_z / std::sin(dphi); // dx per layer
+  double dx = mod_z / std::sin(dphi); // dx per layer
   dx = 0;
     
   // Compute the top and bottom face measurements.
diff --git a/DDExamples/ILDExDet/src/TPCModule.cpp b/DDExamples/ILDExDet/src/TPCModule.cpp
index f484351a5..7e13e90c0 100644
--- a/DDExamples/ILDExDet/src/TPCModule.cpp
+++ b/DDExamples/ILDExDet/src/TPCModule.cpp
@@ -10,6 +10,7 @@
 //  to inherit the special implementations of different PadLayouts.
 //
 //====================================================================
+#include "DD4hep/LCDD.h"
 
 #include "TPCModule.h"
 #include "TPCModuleData.h"
@@ -113,6 +114,23 @@ namespace DD4hep {
     return pad - 1;
   }
 
+  void printPlace(const PlacedVolume& pl) {
+    TGeoMatrix* m = pl->GetMatrix();
+    if ( m ) {
+      const Double_t* tr = m->GetTranslation();
+      if ( tr ) {
+	//cout << "\t\t Tr:" << tr[0] << " " << tr[1] << " " << tr[2] << endl;
+      }
+    }
+  }
+  void check_points(const char* prefix, const Position& global, const double* point) {
+    cout << prefix << global.x   << " " << global.y   << " " << global.z   << endl;
+    if ( global.x != point[0] || global.y != point[1] || global.z != point[2] ) {
+      cout << "Points are different -- wrong transformation!!!!" << endl;
+      cout << "Points are " << point[0]  << " " << point[1] << " " << point[2] << endl;
+    }
+  }
+
   std::vector<double>  TPCModule::getPadCenter (int pad) const {
     if(pad>getNPads())
       throw OutsideGeometryException("getPadCenter: Requested pad not on module querried!");
@@ -129,9 +147,12 @@ namespace DD4hep {
     Double_t point_local[3];
     Double_t point_global[3];
     Double_t point_global_m[3];
+    Double_t point_global_t[3];
     point_local[0]=pad_x;
     point_local[1]=pad_y;
-    point_local[2]=getModuleZPosition();
+
+    // MF: In local coords this should be 0 no ?
+    point_local[2]=0;//getModuleZPosition();
 
     TGeoManager *geom=volume()->GetGeoManager();
     DetElement parent   = _data().parent;
@@ -146,26 +167,44 @@ namespace DD4hep {
     placement()->LocalToMaster(point_local, point_global);
     //one further up the tree. framework should provilde local to world
     parent.placement()->LocalToMaster(point_global, point_global_m);
-#if 0
-    std::cout<<"Exp-Local: "<<point_local[0]<<" "<<point_local[1]<<" "<<point_local[2]<<std::endl;
-    std::cout<<"Exp-Top: "<<point_global_m[0]<<" "<<point_global_m[1]<<" "<<point_global_m[2]<<std::endl;
+    parent.parent().placement()->LocalToMaster(point_global_m, point_global_t);
+    //#if 0
+    cout << "Exp-Loc: Name:     "        << placement()->GetName()
+	 << endl << "         Placement:" << placementPath()
+	 << endl << "         Path:     " << path()
+	 << endl << "         " << point_local[0] << " " << point_local[1] << " " << point_local[2] << endl;
+    printPlace(placement());
+    cout << "Exp-Par:" << parent.placement()->GetName()  << " " << parent.placementPath() << endl << "        "
+	 << " " << point_global[0] << " " << point_global[1] << " "<<point_global[2] << endl;
+    printPlace(parent.placement());
+    cout << "Exp-Env:" << parent.parent().placement()->GetName() << " " 
+	 << parent.parent().placementPath() << endl << "        "
+	 << " " << point_global_m[0] << " " << point_global_m[1] << " "<<point_global_m[2] << endl;
+    printPlace(parent.parent().placement());
+    cout << "Exp-Top:" << parent.parent().parent().placement()->GetName() << " "
+	 << parent.parent().parent().placementPath() << endl << "        "
+	 << " " << point_global_t[0] << " " << point_global_t[1] << " "<<point_global_t[2] << endl;
+    printPlace(parent.parent().parent().placement());
 
     // Now this should be equivalent
     Position global, global_w, global_r, global_p, local(point_local[0],point_local[1],point_local[2]);
     this->localToParent(local,global);
-    std::cout<<"Det-Par: " << global.x   << " " << global.y   << " " << global.z   << std::endl;
+    check_points("Det-Par: ", global, point_global);
+
     this->localToWorld(local,global_w);
-    std::cout<<"Det-Top: " << global_w.x << " " << global_w.y << " " << global_w.z << std::endl;
+    check_points("Det-Top: ", global_w, point_global_t);
+
     const_cast<TPCModule*>(this)->setReference(parent);
     this->localToReference(local,global_r);
-    std::cout<<"Det-Ref: " << global_r.x << " " << global_r.y << " " << global_r.z << std::endl;
+    check_points("Det-Ref: ", global_r, point_global);
+
     parent.setReference(*this);
-    parent.localToReference(local,global_p);
-    std::cout<<"Loc-Ref: " << global_p.x << " " << global_p.y << " " << global_p.z << std::endl;
-#endif
+    parent.localToReference(global,global_p);
+    cout << "Loc-Ref: " << global_p.x << " " << global_p.y << " " << global_p.z << endl;
+    //#endif
     // Need to check if the results are the same....
 
-    std::vector<double> center;
+    vector<double> center;
     center.push_back(point_global_m[0]);
     center.push_back(point_global_m[1]);
     return center;
@@ -183,9 +222,9 @@ namespace DD4hep {
     DetElement parent   = _data().parent;
     parent.placement()->MasterToLocal(point_global, point_global_m);
     placement()->MasterToLocal(point_global_m, point_local);
-  //   std::cout<<"Global: "<<point_global[0]<<" "<<point_global[1]<<" "<<point_global[2]<<std::endl;
-//     std::cout<<"Mother: "<<point_global_m[0]<<" "<<point_global_m[1]<<" "<<point_global_m[2]<<std::endl;
-//     std::cout<<"Local: "<<point_local[0]<<" "<<point_local[1]<<" "<<point_local[2]<<std::endl;
+    //   cout<<"Global: "<<point_global[0]<<" "<<point_global[1]<<" "<<point_global[2]<<endl;
+    //     cout<<"Mother: "<<point_global_m[0]<<" "<<point_global_m[1]<<" "<<point_global_m[2]<<endl;
+    //     cout<<"Local: "<<point_local[0]<<" "<<point_local[1]<<" "<<point_local[2]<<endl;
     
     //check if it is on that module
     bool onMod=volume().solid()->Contains(point_local);
-- 
GitLab