From 1b5d6b05a3278b15fc55f809d2eaac3d2a1310bf Mon Sep 17 00:00:00 2001
From: Astrid Munnich <astrid.muennich@cern.ch>
Date: Thu, 10 May 2012 12:35:38 +0000
Subject: [PATCH] First implementation of digi and reco toy chanin for TPC. Not
 correctly functioning yet, as it needs WorldToLocal trafo which is in the
 works.

---
 DDExamples/ILDExDet/compact/ILDExTPC.xml      |   6 +-
 DDExamples/ILDExDet/src/TPCModule.cpp         |  20 +--
 DDExamples/ILDExReco/CMakeLists.txt           |   6 +
 .../ILDExReco/CalculateTPCPadBinning.cpp      |  50 +++++++
 DDExamples/ILDExReco/GenerateTPCHits.cpp      |  82 +++++++++++
 DDExamples/ILDExReco/TPCDigitisation.cpp      | 126 ++++++++---------
 DDExamples/ILDExReco/TPCHitReco.cpp           | 129 ++++++++++++++++++
 DDExamples/ILDExReco/src/GearTPC.cpp          |  13 +-
 8 files changed, 353 insertions(+), 79 deletions(-)
 create mode 100644 DDExamples/ILDExReco/CalculateTPCPadBinning.cpp
 create mode 100644 DDExamples/ILDExReco/GenerateTPCHits.cpp
 create mode 100644 DDExamples/ILDExReco/TPCHitReco.cpp

diff --git a/DDExamples/ILDExDet/compact/ILDExTPC.xml b/DDExamples/ILDExDet/compact/ILDExTPC.xml
index f89aef800..63558e20c 100644
--- a/DDExamples/ILDExDet/compact/ILDExTPC.xml
+++ b/DDExamples/ILDExDet/compact/ILDExTPC.xml
@@ -113,15 +113,15 @@
 
     <readouts>
         <readout name="PadLayout0">
-          <segmentation type="ProjectiveCylinder"  phiBins="4.0" thetaBins="5.0"/> 
+          <segmentation type="ProjectiveCylinder"  phiBins="240" thetaBins="21"/> 
           <id>system:6</id>
         </readout>
        <readout name="PadLayout1">
-          <segmentation type="ProjectiveCylinder"  phiBins="9.0" thetaBins="9.0"/> 
+          <segmentation type="ProjectiveCylinder"  phiBins="260" thetaBins="21"/> 
           <id>system:7</id>
         </readout>
        <readout name="PadLayout2">
-          <segmentation type="ProjectiveCylinder"  phiBins="10.0" thetaBins="6.0"/> 
+          <segmentation type="ProjectiveCylinder"  phiBins="280" thetaBins="21"/> 
           <id>system:8</id>
         </readout>
     </readouts>
diff --git a/DDExamples/ILDExDet/src/TPCModule.cpp b/DDExamples/ILDExDet/src/TPCModule.cpp
index 6bff9df85..35efc0dff 100644
--- a/DDExamples/ILDExDet/src/TPCModule.cpp
+++ b/DDExamples/ILDExDet/src/TPCModule.cpp
@@ -58,7 +58,7 @@ namespace DD4hep {
   
   double TPCModule::getRowHeight(int row)const {
     if(row>getNRows())
-      throw OutsideGeometryException("Requested row not on module querried!");
+      throw OutsideGeometryException("getRowHeight: Requested row not on module querried!");
     //all rows are the same for FixedAnglePadLayout=ProjectiveCylinder 
     Tube       tube=volume().solid();
     double module_height= tube->GetRmax()-tube->GetRmin();
@@ -67,13 +67,13 @@ namespace DD4hep {
   
   int TPCModule::getRowNumber(int pad)const {
     if(pad>getNPads())
-      throw OutsideGeometryException("Requested pad not on module querried!");
+      throw OutsideGeometryException("getRowNumber: Requested pad not on module querried!");
     return pad / (getNPads()/getNRows());
   }
   
   double TPCModule::getPadPitch(int pad)const {
     if(pad>getNPads())
-      throw OutsideGeometryException("Requested pad not on module querried!");
+      throw OutsideGeometryException("getPadPitch: Requested pad not on module querried!");
     int row=getRowNumber(pad);
     Tube tube=volume().solid();
     double pad_radius=tube->GetRmin()+(row+0.5)*getRowHeight(0);
@@ -84,15 +84,15 @@ namespace DD4hep {
   
   int TPCModule::getPadNumber(int pad)const {
    if(pad>getNPads())
-      throw OutsideGeometryException("Requested pad not on module querried!");
+      throw OutsideGeometryException("getPadNumber: Requested pad not on module querried!");
     return pad % (getNPads()/getNRows());
   }
 
   int TPCModule::getPadIndex(int row,int padNr)const {
     if(padNr>(getNPads()/getNRows()))
-      throw OutsideGeometryException("Requested pad not on module querried!");
+      throw OutsideGeometryException("getPadIndex: Requested pad not on module querried!");
     if(row>getNRows())
-      throw OutsideGeometryException("Requested row not on module querried!");
+      throw OutsideGeometryException("getPadIndex: Requested row not on module querried!");
     return padNr + row*(getNPads()/getNRows());
   }
 
@@ -100,7 +100,7 @@ namespace DD4hep {
     //if on edge their is no neighbour, should throw an exception
     int row=getRowNumber(pad);
     if(getPadNumber(pad)==getNPadsInRow(row)-1)
-      throw OutsideGeometryException("Requested pad is on right edge and has no right neighbour!");
+      throw OutsideGeometryException("getRightNeighbour: Requested pad is on right edge and has no right neighbour!");
     // if not on edge
     return pad + 1;
   }
@@ -108,14 +108,14 @@ namespace DD4hep {
   int TPCModule::getLeftNeighbour(int pad)const {
     //if on edge their is no neighbour, should throw an exception
     if(getPadNumber(pad)==0)
-      throw OutsideGeometryException("Requested pad is on left edge and has no left neighbour!");
+      throw OutsideGeometryException("getLeftNeighbour: Requested pad is on left edge and has no left neighbour!");
     // if not on edge
     return pad - 1;
   }
 
   std::vector<double>  TPCModule::getPadCenter (int pad) const {
     if(pad>getNPads())
-      throw OutsideGeometryException("Requested pad not on module querried!");
+      throw OutsideGeometryException("getPadCenter: Requested pad not on module querried!");
     int row=getRowNumber(pad);
     Tube tube=volume().solid();
     double pad_radius=tube->GetRmin()+(row+0.5)*getRowHeight(0);
@@ -170,7 +170,7 @@ namespace DD4hep {
     //check if it is on that module
     bool onMod=volume().solid()->Contains(point_local);
     if(!onMod)
-      throw OutsideGeometryException("Requested point not on module querried!");
+      throw OutsideGeometryException("getNearestPad: Requested point not on module querried!");
     Tube tube=volume().solid();
     double module_width= tube->GetPhi2()-tube->GetPhi1();
     double radius=sqrt(point_local[0]*point_local[0]+point_local[1]*point_local[1]);
diff --git a/DDExamples/ILDExReco/CMakeLists.txt b/DDExamples/ILDExReco/CMakeLists.txt
index 2c79501f8..efd81d6ab 100644
--- a/DDExamples/ILDExReco/CMakeLists.txt
+++ b/DDExamples/ILDExReco/CMakeLists.txt
@@ -7,6 +7,12 @@ include_directories( ${CMAKE_SOURCE_DIR}/DDCore/include
 
 add_executable(ILDExReco ILDExReco.cpp src/GearTPC.cpp)
 add_executable(TPCDigitisation TPCDigitisation.cpp src/GearTPC.cpp)
+add_executable(TPCHitReco TPCHitReco.cpp src/GearTPC.cpp)
+add_executable(CalculateTPCPadBinning CalculateTPCPadBinning.cpp src/GearTPC.cpp)
+add_executable(GenerateTPCHits GenerateTPCHits.cpp src/GearTPC.cpp)
 target_link_libraries(ILDExReco DD4hepCore ILDEx ${ROOT_LIBRARIES})
 target_link_libraries(TPCDigitisation DD4hepCore ILDEx ${ROOT_LIBRARIES})
+target_link_libraries(TPCHitReco DD4hepCore ILDEx ${ROOT_LIBRARIES})
+target_link_libraries(CalculateTPCPadBinning DD4hepCore ILDEx ${ROOT_LIBRARIES})
+target_link_libraries(GenerateTPCHits DD4hepCore ILDEx ${ROOT_LIBRARIES})
 
diff --git a/DDExamples/ILDExReco/CalculateTPCPadBinning.cpp b/DDExamples/ILDExReco/CalculateTPCPadBinning.cpp
new file mode 100644
index 000000000..6f1220bc2
--- /dev/null
+++ b/DDExamples/ILDExReco/CalculateTPCPadBinning.cpp
@@ -0,0 +1,50 @@
+//  Copyright 2012 __MyCompanyName__. All rights reserved.
+//
+//====================================================================
+//  Demonstrator application for using TPC functionality 
+//  to calculate the binning of pads on modules
+//--------------------------------------------------------------------
+//
+//  Author     : A.Muennich
+//
+//====================================================================
+
+#include "DD4hep/LCDD.h"
+#include <iostream>
+#include <vector>
+#include <string>
+#include <TFile.h>
+#include <TTree.h>
+#include "GearTPC.h"
+#include "TGeoTube.h"
+
+using namespace std;
+using namespace DD4hep;
+using namespace Geometry;
+
+
+int main(int argc,char** argv)  {
+  
+  // instanciate geometry from copact file
+  LCDD& lcdd = LCDD::getInstance();
+  lcdd.fromCompact(argv[1]);
+  // get the TPC
+  GearTPC tpc(lcdd.detector("TPC"));
+ 
+  // for 5x10 pads
+  double pad_height=10;
+  double pad_width=5;
+  std::vector<TPCModule> mymods=tpc.getModules(0);
+  for(int m=0;m<mymods.size();m++)
+    {
+      TPCModule mymod=tpc.getModule(m,0);
+      Tube       tube=mymod.volume().solid();
+      double module_height= tube->GetRmax()-tube->GetRmin();
+      double module_angle= tube->GetPhi2()-tube->GetPhi1();
+      double pad_angle=pad_width/tube->GetRmin()*180/M_PI;
+      cout << "Mod "<<m<<" r= " << module_height <<" , phi= "<< module_angle
+	   <<"\t Binning theta= "<<module_height/pad_height
+	   <<"   phi= "<<module_angle/pad_angle<<endl;
+     }
+  return 0;
+}
diff --git a/DDExamples/ILDExReco/GenerateTPCHits.cpp b/DDExamples/ILDExReco/GenerateTPCHits.cpp
new file mode 100644
index 000000000..5ef89e893
--- /dev/null
+++ b/DDExamples/ILDExReco/GenerateTPCHits.cpp
@@ -0,0 +1,82 @@
+//  Copyright 2012 __MyCompanyName__. All rights reserved.
+//
+//====================================================================
+//  Fake hit generator to be replaced by G4 simulation
+//--------------------------------------------------------------------
+//
+//  Author     : A.Muennich
+//
+//====================================================================
+
+#include "DD4hep/LCDD.h"
+#include <iostream>
+#include <vector>
+#include <string>
+#include <TFile.h>
+#include <TTree.h>
+#include <TRandom.h>
+#include "GearTPC.h"
+
+using namespace std;
+using namespace DD4hep;
+using namespace Geometry;
+
+
+int main(int argc,char** argv)  {
+  
+  // instanciate geometry from copact file
+  LCDD& lcdd = LCDD::getInstance();
+  lcdd.fromCompact(argv[1]);
+  // get the TPC
+  GearTPC tpc(lcdd.detector("TPC"));
+ 
+
+  //Output file containing pads that are hit
+  TFile *out_file = new TFile("TPCSimulation_Output.root","RECREATE");
+  //vectors for output file
+  int EvID;
+  std::vector<double> xPos;
+  std::vector<double> yPos;
+  std::vector<double> zPos;
+  std::vector<double> charge;
+  //Tree structure
+  TTree *out_tree = new TTree("SimTree","SimTree");
+  out_tree->Branch("EvID",&EvID);
+  out_tree->Branch("xPos",&xPos) ;
+  out_tree->Branch("yPos",&yPos) ;
+  out_tree->Branch("zPos",&zPos) ;
+  out_tree->Branch("charge",&charge);
+
+  TRandom *rndm = new TRandom();
+  //event loop
+  int NEVENTS=10;
+  int NPOINTS=100;
+  for(int i=0;i<NEVENTS;i++)
+    {
+      //random direction
+      double phi=i*2*M_PI/NEVENTS;
+      double theta=40*M_PI/180+i*(50*M_PI/180)/NEVENTS;
+      //TPC half
+      int sign_z=1;
+   //    if(i%2)
+//  	sign_z=-1;
+      //point loop
+      for (int p=0;p<NPOINTS;p++)
+	{
+	  double radius=tpc.getInnerRadius()+(tpc.getOuterRadius()-tpc.getInnerRadius())*p/NPOINTS;
+	  xPos.push_back(radius*cos(phi));
+	  yPos.push_back(radius*sin(phi));
+	  zPos.push_back(sign_z*radius/tan(theta));
+	  charge.push_back(rndm->Uniform(5,10));
+	  //	  cout<<i<<"\t"<<p<<"\t"<<radius<<"\t"<<theta*180/M_PI<<"\t"<<sign_z*radius/tan(theta)<<"\t"<<tan(theta)<<endl;
+	}
+      EvID=i;
+      out_tree->Fill();
+      xPos.clear();
+      yPos.clear();
+      zPos.clear();
+      charge.clear();
+    }
+  out_file->Write(); 
+  return 0;
+}
diff --git a/DDExamples/ILDExReco/TPCDigitisation.cpp b/DDExamples/ILDExReco/TPCDigitisation.cpp
index 3a4830e80..e55e6f94e 100644
--- a/DDExamples/ILDExReco/TPCDigitisation.cpp
+++ b/DDExamples/ILDExReco/TPCDigitisation.cpp
@@ -49,80 +49,82 @@ int main(int argc,char** argv)  {
   out_tree->Branch("chargeDeposition",&chargeDeposition);
   
   //read in File with space points of energy deposits
-  //Event Loop
+  TFile *input=new TFile("TPCSimulation_Output.root","READ");
+  int EvID_in;
+  std::vector<double> *xPos=0;
+  std::vector<double> *yPos=0;
+  std::vector<double> *zPos=0;
+  std::vector<double> *charge=0;
+  TTree *sim_data = (TTree*)gDirectory->Get("SimTree");
+  // Set branch addresses.
+  sim_data->SetBranchAddress("EvID",&EvID_in);			 
+  sim_data->SetBranchAddress("xPos",&xPos) ;		 
+  sim_data->SetBranchAddress("yPos",&yPos) ;		 
+  sim_data->SetBranchAddress("zPos",&zPos) ;		 
+  sim_data->SetBranchAddress("charge",&charge) ;		 
   
-  double spacepoint[20][4];
-  //test fill
-  for(int p=0;p<10;p++)
-    {
-      spacepoint[p][0]=0;
-      spacepoint[p][1]=tpc.getInnerRadius()+(tpc.getOuterRadius()-tpc.getInnerRadius())*p/10;
-      spacepoint[p][2]=tpc.getMaxDriftLength()*p/10;
-      spacepoint[p][3]=1;
-    }
-  for(int p=10;p<20;p++)
-    {
-      spacepoint[p][0]=spacepoint[p-10][0];
-      spacepoint[p][1]=spacepoint[p-10][1];
-      spacepoint[p][2]=spacepoint[p-10][2]+10;
-      spacepoint[p][3]=4;
-    }
-
   //container to hold the count of pads hit
   //map<pair<moduleID,padID>,vector<z,Q>
   std::map< std::pair<int,int> , std::pair<double,double>  > padmap;
   map< std::pair<int,int> , std::pair<double,double> >::iterator it;
-  //loop the spacepoints x,y,z,E
-  for(int p=0;p<20;p++)
+
+  //Event Loop
+  for(int i=0;i<sim_data->GetEntries();i++)
     {
-      double x=spacepoint[p][0];
-      double y=spacepoint[p][1];
-      double z=spacepoint[p][2];
-      double E=spacepoint[p][3];
+      //pad loop
+      sim_data->GetEntry(i);
+      //loop the spacepoints x,y,z,E
+      for(int p=0;p<xPos->size();p++)
+	{
+	  double x=(*xPos)[p];
+	  double y=(*yPos)[p];
+	  double z=(*zPos)[p];
+	  double E=(*charge)[p];
 
-      int endplate=0;
-      if(z<0)
-	endplate=1;
-      //check if point is over module. If not move to next point
-      if(!tpc.isInsideModule(x,y,endplate))
-	continue;
-      TPCModule mymod=tpc.getNearestModule(x,y,endplate);
-      int modID=mymod.getID();
-      int padID= mymod.getNearestPad(x,y);
-      cout<<p<<"\t"<<spacepoint[p][0]<<"\t"<<spacepoint[p][1]<<"\t"<<spacepoint[p][2]<<"\t"<<modID<<"\t"<<padID<<endl;  
-      std::pair<int,int> ID_pair=make_pair(modID,padID);;
+	  int endplate=0;
+	  if(z<0)
+	    endplate=1;
+	  //check if point is over module. If not move to next point
+	  bool measured=tpc.isInsideModule(x,y,endplate);
+	  if(!measured)
+	    continue;
+	  TPCModule mymod=tpc.getNearestModule(x,y,endplate);
+	  int modID=mymod.getID();
+	  int padID= mymod.getNearestPad(x,y);
+	  std::pair<int,int> ID_pair=make_pair(modID,padID);
+	  
+	  //check if that pad has already been hit
+	  it = padmap.find(ID_pair);
+	  if(it!=padmap.end())
+	    {
+	      double new_z=(it->second.first*it->second.second+z*E)/(it->second.second+E);
+	      double new_E=it->second.second+E;
+	      padmap[ID_pair]=make_pair(new_z,new_E);
+	    }
+	  else //create a new entry
+	    {
+	      std::pair<double,double>  hit=make_pair(z,E);
+	      padmap[ID_pair]=hit;
+	    }
+	}//point loop
 
-      //check if that pad has already been hit
-      it = padmap.find(ID_pair);
-      if(it!=padmap.end())
-	{
-	  double new_z=(it->second.first*it->second.second+z*E)/(it->second.second+E);
-	  double new_E=it->second.second+E;
-	  padmap[ID_pair]=make_pair(new_z,new_E);
-	}
-      else //create a new entry
+      for ( it=padmap.begin() ; it != padmap.end(); it++ )
 	{
-	  std::pair<double,double>  hit=make_pair(z,E);
-	  padmap[ID_pair]=hit;
+	  moduleIDs.push_back(it->first.first);
+	  padIDs.push_back(it->first.second);
+	  zPositions.push_back(it->second.first);
+	  chargeDeposition.push_back(it->second.second);
 	}
+      EvID=EvID_in;
+      out_tree->Fill();
+      //clean up for next event
+      moduleIDs.clear();
+      padIDs.clear();
+      zPositions.clear();
+      chargeDeposition.clear();
+      padmap.clear();
     }
 
-  for ( it=padmap.begin() ; it != padmap.end(); it++ )
-    {
-      cout << it->first.first << "\t" <<(*it).first.second << "\t" << (*it).second.first <<"\t" << (*it).second.second << endl;
-      moduleIDs.push_back(it->first.first);
-      padIDs.push_back(it->first.second);
-      zPositions.push_back(it->second.first);
-      chargeDeposition.push_back(it->second.second);
-    }
-  out_tree->Fill();
-  //clean up for next event
-  moduleIDs.clear();
-  padIDs.clear();
-  zPositions.clear();
-  chargeDeposition.clear();
-
-
   //event loop ends, write file
   out_file->Write();
 
diff --git a/DDExamples/ILDExReco/TPCHitReco.cpp b/DDExamples/ILDExReco/TPCHitReco.cpp
new file mode 100644
index 000000000..a34717f1d
--- /dev/null
+++ b/DDExamples/ILDExReco/TPCHitReco.cpp
@@ -0,0 +1,129 @@
+//  Copyright 2012 __MyCompanyName__. All rights reserved.
+//
+//====================================================================
+//  Demonstrator application for using TPC functionality 
+//  for hit reconstruction
+//--------------------------------------------------------------------
+//
+//  Author     : A.Muennich
+//  Reads in pads with charge and reconstructs 3D space points.
+//
+//====================================================================
+
+#include "DD4hep/LCDD.h"
+#include <iostream>
+#include <vector>
+#include <string>
+#include <TFile.h>
+#include <TTree.h>
+#include "GearTPC.h"
+
+
+using namespace std;
+using namespace DD4hep;
+using namespace Geometry;
+
+
+int main(int argc,char** argv)  {
+  
+  // instanciate geometry from copact file
+  LCDD& lcdd = LCDD::getInstance();
+  lcdd.fromCompact(argv[1]);
+  // get the TPC
+  GearTPC tpc(lcdd.detector("TPC"));
+ 
+  //Output file containing the hits
+  TFile *out_file = new TFile("TPCHitReco_Output.root","RECREATE");
+  //vectors for output file
+  int EvID;
+  std::vector<double> xPos;
+  std::vector<double> yPos;
+  std::vector<double> zPos;
+  std::vector<double> charge;
+  //Tree structure
+  TTree *out_tree = new TTree("HitTree","HItTree");
+  out_tree->Branch("EvID",&EvID);
+  out_tree->Branch("xPos",&xPos) ;
+  out_tree->Branch("yPos",&yPos) ;
+  out_tree->Branch("zPos",&zPos) ;
+  out_tree->Branch("charge",&charge);
+  
+  //read in File with space points of energy deposits
+  TFile *input=new TFile("TPCDigitisation_Output.root","READ");
+  int EvID_in;
+  std::vector<int> *moduleIDs=0;
+  std::vector<int> *padIDs=0;
+  std::vector<double> *zPositions=0;
+  std::vector<double> *chargeDeposition=0;
+  
+  TTree *pad_data = (TTree*)gDirectory->Get("PadTree");
+  // Set branch addresses.
+  pad_data->SetBranchAddress("EvID",&EvID_in);			 
+  pad_data->SetBranchAddress("moduleIDs",&moduleIDs) ;		 
+  pad_data->SetBranchAddress("padIDs",&padIDs) ;			 
+  pad_data->SetBranchAddress("zPositions",&zPositions) ;		 
+  pad_data->SetBranchAddress("chargeDeposition",&chargeDeposition);
+
+  //Fill pads again in map for easier search of neighbours
+  // on map per endplate
+  //map<pair<moduleID,padID>,vector<z,Q>
+  std::map< std::pair<int,int> , std::pair<double,double>  > padmap_EP1;
+  map< std::pair<int,int> , std::pair<double,double> >::iterator it1;
+  std::map< std::pair<int,int> , std::pair<double,double>  > padmap_EP2;
+  map< std::pair<int,int> , std::pair<double,double> >::iterator it2;
+    
+  //Event Loop
+  for(int i=0;i<pad_data->GetEntries();i++)
+    {
+      //pad loop
+      pad_data->GetEntry(i);
+
+      for(int p=0;p<moduleIDs->size();p++)
+	{
+	  std::pair<int,int> ID_pair=make_pair((*moduleIDs)[p],(*padIDs)[p]);
+	  std::pair<double,double>  hit=make_pair((*zPositions)[p],(*chargeDeposition)[p]);
+	  if((*zPositions)[p] > 0)
+	    padmap_EP1[ID_pair]=hit;
+	  else
+	    padmap_EP2[ID_pair]=hit;
+	}
+   
+      //hit reco goes per module per row
+      //positive endplate
+      for ( it1=padmap_EP1.begin() ; it1 != padmap_EP1.end(); it1++ )
+	{
+	  TPCModule mymod=tpc.getModule(it1->first.first,0);
+	  cout << it1->first.first << "\t" <<(*it1).first.second <<" "<<mymod.getNPads()<< endl;
+	  std::vector<double> center=mymod.getPadCenter(it1->first.second);
+	   
+	  xPos.push_back(center[0]);
+	  yPos.push_back(center[1]);
+	  zPos.push_back(it1->second.first);
+	  charge.push_back(it1->second.second);
+	}
+     //negative endplate
+      for ( it2=padmap_EP2.begin() ; it2 != padmap_EP2.end(); it2++ )
+	{
+	  std::vector<double> center=tpc.getModule(it2->first.first,1).getPadCenter(it2->first.second);
+	  
+	  xPos.push_back(center[0]);
+	  yPos.push_back(center[1]);
+	  zPos.push_back(it2->second.first);
+	  charge.push_back(it2->second.second);
+	}
+      EvID=EvID_in;
+      out_tree->Fill();
+      //clean up for next event
+      xPos.clear();
+      yPos.clear();
+      zPos.clear();
+      charge.clear();
+      padmap_EP1.clear();
+      padmap_EP2.clear();
+    }
+  
+  //event loop ends, write file
+  out_file->Write();
+
+  return 0;
+}
diff --git a/DDExamples/ILDExReco/src/GearTPC.cpp b/DDExamples/ILDExReco/src/GearTPC.cpp
index fb63f0183..051fa753c 100644
--- a/DDExamples/ILDExReco/src/GearTPC.cpp
+++ b/DDExamples/ILDExReco/src/GearTPC.cpp
@@ -106,19 +106,22 @@ namespace DD4hep {
     TGeoManager *geoManager = ep.volume()->GetGeoManager();
     TGeoNode *mynode=geoManager->FindNode(c0,c1,zpos);
     Double_t point[3];
+    Double_t point_mother[3];
     Double_t point_local[3];
     point[0]=c0;
     point[1]=c1;
     point[2]=zpos;
     //FIXME: careful: master is mother not global=world, input is in world coordinates
-    ep.placements()[0]->MasterToLocal(point, point_local);
+    //ep.parent.placement()->MasterToLocal(point, point_mother);
+    ep.placement()->MasterToLocal(point, point_local);
+   
     
     bool onMod=false;
     std::map<std::string,DetElement>::const_iterator it;
     for ( it=ep.children().begin() ; it != ep.children().end(); it++ )
       {
 	Double_t point_local_node[3];
- 	it->second.placements()[0]->MasterToLocal(point_local, point_local_node);
+ 	it->second.placement()->MasterToLocal(point_local, point_local_node);
 	onMod=it->second.volume().solid()->Contains(point_local_node);
 	if(onMod)
 	  {
@@ -136,13 +139,15 @@ namespace DD4hep {
     TGeoManager *geoManager = ep.volume()->GetGeoManager();
     TGeoNode *mynode=geoManager->FindNode(c0,c1,zpos);
     Double_t point[3];
+    Double_t point_mother[3];
     Double_t point_local[3];
     point[0]=c0;
     point[1]=c1;
     point[2]=zpos;
    //FIXME: careful: master is mother not global=world, input is in world coordinates
-    ep.placements()[0]->MasterToLocal(point, point_local);
-    geoManager->SetCurrentPoint(point_local);
+    //  ep.parent.placement()->MasterToLocal(point, point_mother);
+    ep.placement()->MasterToLocal(point, point_local);
+   
     bool onMod=false;
     std::map<std::string,DetElement>::const_iterator it;
     //check if any of the modules contains that point
-- 
GitLab