Muon Tracking Geometry

Working with control version git (but in parallel with svn using the git svn tool).

Repository from svn
Just using the asetup X.Y.Z,here the SVNROOT variable is already set. Anyway, it should be
SVNROOT=svn+ssh://svn.cern.ch/reps/atlasoff
And obtaining the code with
git svn clone svn+ssh://duarte@svn.cern.ch/reps/atlasoff/MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry -s
Note that initially I did just: git svn clone MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry -s assuming the previous export variable, but appears the error: E: 'trunk' is not a complete URL and a separate URL is not specified.

Note also that the whole process lasts some time due that git is localizing and keeping track of all branches and tags. The svn and git link explains how to work with a git repository pushing to a remote svn. The repository has been cloned in the lxplus path: /afs/cern.ch/user/d/duarte/work/private/MTG/CODE_19_2_0/MuonTrackingGeometry.

Apart from the atlasoff repository, a remote repository has been put in place in github, just to push the commits and have a double copy of the work. The forked repository can be found at MuonTrackingGeometry github repo.

Muon Tracking Geometry package validation

Problem: Geometry defined in Geant4 presents important discrepancies with the tracking geometry (a optimized and non-detailed geometry to perform tracking).
Tentative Strategy:
  • ....

Setting-up the tools

  • Athena release: 19.2.0. Needs some fixes in the following packages:
Package cmt version
Generators/ParticleGenerator 00-00-61
Simulation/G4Atlas/G4AtlasApps 00-07-63
Simulation/ISF/ISF_Config 00-00-54
Simulation/ISF/ISF_Core/ISF_Services 00-03-07
Simulation/ISF/ISF_Geant4/ISF_Geant4Tools 00-02-00
Tracking/TrkDetDescr/TrkVolumes 01-01-03
Tracking/TrkExtrapolation//TrkExSTEP_Propagator 01-00-29
Tracking/TrkExtrapolation/TrkExTools 03-00-34-branch

Generator/ParticleGenerator
===================================================================
--- src/ParticleKinematics.cxx	(revisión: 618359)
+++ src/ParticleKinematics.cxx	(copia de trabajo)
@@ -217,7 +217,8 @@
 
     const HepPDT::ParticleData* particle 
       = p_particleTable->particle(HepPDT::ParticleID(abs( m_pdg )));
-    m_mass = particle->mass().value();
+ 
+    m_mass = particle ? particle->mass().value() : 0.;
 
     m_time = getValue(kTime);

Simulation/G4Atlas/G4AtlasApps
Index: Simulation/G4Atlas/G4AtlasApps/python/SimAtlasKernel.py
===================================================================
--- Simulation/G4Atlas/G4AtlasApps/python/SimAtlasKernel.py	(revisión: 592906)
+++ Simulation/G4Atlas/G4AtlasApps/python/SimAtlasKernel.py	(copia de trabajo)
@@ -266,6 +266,11 @@
           from G4AtlasApps import AtlasG4Eng,PyG4Atlas
           actions = AtlasG4Eng.G4Eng.menu_UserActions()
 
+          # PhysicsValidationUserAction
+          AtlasG4Eng.G4Eng.log.verbose('ISF_AtlasSimSkeleton::do_UserActions add ISF_G4PhysicsValidation')
+          ISF_G4PhysicsValidation = PyG4Atlas.UserAction('ISF_Geant4Tools','ToolSvc.ISFG4PhysicsValidationUserAction',['BeginOfEvent','EndOfEvent','BeginOfRun','EndOfRun','Step'])
+          actions.add_UserAction(ISF_G4PhysicsValidation)
+
           # TrackProcessorUserAction
           AtlasG4Eng.G4Eng.log.verbose('ISF_AtlasSimSkeleton::do_UserActions add ISF_G4TrackProcessor')
           ISF_G4TrackProcessor = PyG4Atlas.UserAction('ISF_Geant4Tools','ToolSvc.ISFG4TrackProcessorUserAction',['BeginOfEvent','EndOfEvent','BeginOfRun','EndOfRun','Step'])
@@ -636,7 +641,7 @@
 
         from G4AtlasApps.SimFlags import simFlags
         if not simFlags.ISFRun:
-            from atlas_utilities import G4SimTimer
+            from atlas_utilities import G4SimTimer           
             actions.add_UserAction(G4SimTimer)
             
         actions.add_UserAction(G4TrackCounter)

Simulation/ISF/ISF_Config
Index: Simulation/ISF/ISF_Config/python/ISF_MainConfig.py
===================================================================
--- Simulation/ISF/ISF_Config/python/ISF_MainConfig.py	(revision 615726)
+++ Simulation/ISF/ISF_Config/python/ISF_MainConfig.py	(working copy)
@@ -107,33 +107,33 @@
                                         'vertY: constant 0.0',
                                         'vertZ: constant 0.0',
                                         't: constant 0.0',
-                                        'eta: flat -4.0 4.0',
-                                        'phi: flat  0 6.28318',
-                                        'pt: constant 50000']))
+                                        'eta: flat 0.0 1.0',
+                                        'phi: flat  0.3 0.5',
+                                        'p: constant 50000']))
     return getInput_GenericGenerator(name, **kwargs)
 
 ############## Input: pions ###############
 def getInput_pions(name="ISF_Input_pions", **kwargs):
     kwargs.setdefault('orders', sorted(['pdgcode: sequence 211 -211',
-                                        'vertX: constant 0.0',
+                                        'vertX: constant 1140.0',
                                         'vertY: constant 0.0',
                                         'vertZ: constant 0.0',
                                         't: constant 0.0',
-                                        'eta: flat -4.0 4.0',
-                                        'phi: flat  0 6.28318',
-                                        'pt: constant 50000']))
+                                        'eta: flat -0.02 0.02',
+                                        'phi: flat  1.55 1.57',
+                                        'p: constant 50000']))
     return getInput_GenericGenerator(name, **kwargs)
 
 ############## Input: photons ###############
 def getInput_photons(name="ISF_Input_photons", **kwargs):
     kwargs.setdefault('orders', sorted(['pdgcode: fixed 22',
-                                        'vertX: constant 0.0',
+                                        'vertX: constant 1140.0',
                                         'vertY: constant 0.0',
                                         'vertZ: constant 0.0',
                                         't: constant 0.0',
-                                        'eta: flat -4.0 4.0',
-                                        'phi: flat  0 6.28318',
-                                        'pt: constant 50000']))
+                                        'eta: flat  0.0 0.0',
+                                        'phi: flat  0 0.',
+                                        'p: constant 10000']))
     return getInput_GenericGenerator(name, **kwargs)
 
 
@@ -161,9 +161,9 @@
                                         'vertY: constant 0.0',
                                         'vertZ: constant 0.0',
                                         't: constant 0.0',
-                                        'eta: flat -5.0 5.0',
-                                        'phi: flat  0 6.28318',
-                                        'pt: constant 10000']))
+                                        'eta: flat -4.0 4.0',
+                                        'phi: flat  0. 6.28316',
+                                        'p: constant 10000']))
     return getInput_GenericGenerator(name, **kwargs)
 
 
@@ -285,9 +285,10 @@
 
     SimKernel.BeamPipeSimulationSelectors = [ getPublicTool('ISF_DefaultParticleKillerSelector') ] 
     SimKernel.IDSimulationSelectors       = [ getPublicTool('ISF_DefaultFatrasSelector') ] 
-    SimKernel.CaloSimulationSelectors     = [ getPublicTool('ISF_MuonFatrasSelector'),
-                                                getPublicTool('ISF_EtaGreater5ParticleKillerSimSelector'),
-                                                getPublicTool('ISF_DefaultFastCaloSimSelector')  ]
+    SimKernel.CaloSimulationSelectors       = [ getPublicTool('ISF_DefaultFatrasSelector') ] 
+    #SimKernel.CaloSimulationSelectors     = [ getPublicTool('ISF_MuonFatrasSelector'),
+    #                                            getPublicTool('ISF_EtaGreater5ParticleKillerSimSelector'),
+    #                                            getPublicTool('ISF_DefaultFastCaloSimSelector')  ]
     SimKernel.MSSimulationSelectors       = [ getPublicTool('ISF_DefaultFatrasSelector') ]
     SimKernel.CavernSimulationSelectors   = [ getPublicTool('ISF_DefaultParticleKillerSelector') ]
     # set the simFlags accordingly (TODO: is this even needed?)
Index: Simulation/ISF/ISF_Config/share/jobOptions_ConfGetter.py
===================================================================
--- Simulation/ISF/ISF_Config/share/jobOptions_ConfGetter.py	(revision 615726)
+++ Simulation/ISF/ISF_Config/share/jobOptions_ConfGetter.py	(working copy)
@@ -52,8 +52,8 @@
 #ISF_Flags.OverrideInputFiles = [ '/afs/cern.ch/atlas/groups/fatras/fatras_input/singlepart_mu_pt20_evgen_merged.root' ]
 #ISF_Flags.OverrideInputFiles = [ '/afs/cern.ch/atlas/groups/fatras/fatras_input/singlepart_pi_pt5_evgen_merged.root' ]
 
-#from ISF_FatrasServices.FatrasTuning import FatrasTuningFlags
-#FatrasTuningFlags.MomCutOffSec = 50.
+from ISF_FatrasServices.FatrasTuning import FatrasTuningFlags
+FatrasTuningFlags.MomCutOffSec = 1000.
 
 from G4AtlasApps.SimFlags import simFlags
 simFlags.ReleaseGeoModel = False    # Comment this out when not making making DCube plots
@@ -80,7 +80,38 @@
 
 include('ISF_Config/ISF_ConfigJobInclude.py')
 
+from AGDD2Geo.AGDD2GeoConf import AGDD2GeoSvc
+AGDD2GeoSvc=AGDD2GeoSvc()
+AGDD2GeoSvc.PrintSections=False
+AGDD2GeoSvc.OverrideConfiguration = True
+AGDD2GeoSvc.Locked = False
+AGDD2GeoSvc.DisableSections = False
+AGDD2GeoSvc.Volumes += ["ECT_Toroids"]
+AGDD2GeoSvc.Volumes += ["BAR_Toroid"]
+AGDD2GeoSvc.Volumes += ["Feet"]
+AGDD2GeoSvc.Volumes += ["RailAssembly"]
+AGDD2GeoSvc.Volumes += ["JFSH_Shield"]
+AGDD2GeoSvc.Volumes += ["JTSH_Shield"]
+AGDD2GeoSvc.Volumes += ["JDSH_Shield"]
+AGDD2GeoSvc.Volumes += ["pp2"]
+AGDD2GeoSvc.Volumes += ["servicesAtZ0"]
+AGDD2GeoSvc.Volumes += ["MBAP_AccessPlatform"]
+AGDD2GeoSvc.Volumes += ["MBWH_BigWheels"]
+AGDD2GeoSvc.Volumes += ["SADL_CalorimeterSaddle"]
+AGDD2GeoSvc.Volumes += ["TBWH_BigWheels"]
+AGDD2GeoSvc.Volumes += ["TGC3_BigWheels"]
+AGDD2GeoSvc.Volumes += ["TGC1_BigWheels"]
+AGDD2GeoSvc.Volumes += ["MDTRail"]
+AGDD2GeoSvc.Volumes += ["HFTruckRail"]
+theApp.CreateSvc += ["AGDD2GeoSvc"]
+svcMgr += AGDD2GeoSvc
 
+from TrkExTools.AtlasExtrapolator import AtlasExtrapolator
+AtlasExtrapolator = AtlasExtrapolator()
+AtlasExtrapolator.TrackingGeometrySvc = 'AtlasTrackingGeometrySvc'
+#AtlasExtrapolator.OutputLevel = VERBOSE
+ToolSvc += AtlasExtrapolator
+
 # the particle broker THist stream
 if ISF_Flags.ValidationMode:
   if not hasattr(ServiceMgr, 'THistSvc'):
@@ -91,6 +122,22 @@
   ServiceMgr.THistSvc.Output += ["ISFFatras DATAFILE='ISFFatras.root' OPT='RECREATE'"]
   ServiceMgr.THistSvc.Output += ["ISFG4SimKernel DATAFILE='ISFG4SimKernel.root' OPT='RECREATE'"]
 
-from AthenaCommon.ConfigurationShelve import saveToAscii
-saveToAscii("config.txt")
+#from AthenaCommon.ConfigurationShelve import saveToAscii
+#saveToAscii("config.txt")
 
+
+#from LArTrackingGeometry.LArTrackingGeometryConf import LAr__LArVolumeBuilder
+#LArVolumeBuilder = LAr__LArVolumeBuilder()
+#ToolSvc += LArVolumeBuilder
+#ToolSvc.LArVolumeBuilder.OutputLevel = VERBOSE
+
+#from TileTrackingGeometry.TileTrackingGeometryConf import Tile__TileVolumeBuilder
+#TileVolumeBuilder = Tile__TileVolumeBuilder()
+#ToolSvc += TileVolumeBuilder
+#ToolSvc.TileVolumeBuilder.OutputLevel = VERBOSE
+
+#from CaloTrackingGeometry.ConfiguredCaloTrackingGeometryBuilder import ConfiguredCaloTrackingGeometryBuilder
+#CaloTrackingGeometryBuilder = ConfiguredCaloTrackingGeometryBuilder(name='CaloTrackingGeometryBuilder')
+#ToolSvc += CaloTrackingGeometryBuilder
+#ToolSvc.CaloTrackingGeometryBuilder.OutputLevel = VERBOSE

Simulation/ISF/ISF_Core/ISF_Services
Index: Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx
===================================================================
--- Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx	(revision 615511)
+++ Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx	(working copy)
@@ -162,10 +162,12 @@
     }
 
     // setup for validation mode
+
     if ( m_validationOutput) {
 
       // retrieve the histogram service
       if ( m_thistSvc.retrieve().isSuccess() ) {
+        /*  this is crashing the job
         ATH_CHECK( registerPosValTree( "push_position",
                                        "push() particle positions",
                                        m_t_pushPosition) );
@@ -178,6 +180,7 @@
         ATH_CHECK( registerPosValTree( "muonExit_pos",
                                        "MuonExitLayer positions",
                                        m_t_entryLayerPos[ISF::fAtlasMuonExit]  ) );
+	*/
       }
 
       // error when trying to retrieve the THistSvc
@@ -417,7 +420,7 @@
       m_simSelectorSet.insert( m_simSelector[curGeoID].begin(), m_simSelector[curGeoID].end() );
     }
 
-    ATH_MSG_DEBUG("Number of unique SimulationSelctors registered: "
+    ATH_MSG_DEBUG("Number of unique SimulationSelectors registered: "
                   << m_simSelectorSet.size() );
   }
 
@@ -441,6 +444,8 @@
   ISFParticleContainer::iterator particleIterEnd = initContainer.end();
   for ( ; particleIter != particleIterEnd; ++particleIter) {
 
+    std::cout <<"looping over input particles:"<< std::endl;
+
     // identify the geoID of the particle
     m_geoIDSvcQuick->identifyAndRegNextGeoID(**particleIter);
     // the geoID at this point better makes sense :)
@@ -609,6 +614,7 @@
   ISF::EntryLayer layer = m_entryLayerToolQuick->registerParticle( *particle, layerInput);
 
   // ---> if validation mode: fill the corresponding entry layer ROOT tree
+  /*
   if ( m_validationOutput) {
     // fill the push() position TTree
     fillPosValTree( m_t_pushPosition, *particle);
@@ -616,6 +622,7 @@
     if ( validEntryLayer(layer) )
       fillPosValTree( m_t_entryLayerPos[layer], *particle);
   }
+  */
   // <--- end validation output
 
   // validation mode: check whether the particle position corresponds to the GeoID given

Simulation/ISF/ISF_Geant4/ISF_Geant4Tools
isfgeant4tools.diff: Differences between rev. copy in table above and the local copy (Too many changes to be put explicitly.)

Tracking/TrkDetDescr/TrkVolumes
Index: Tracking/TrkDetDescr/TrkVolumes/TrkVolumes/SubtractedVolumeBounds.h
===================================================================
--- Tracking/TrkDetDescr/TrkVolumes/TrkVolumes/SubtractedVolumeBounds.h	(revision 618906)
+++ Tracking/TrkDetDescr/TrkVolumes/TrkVolumes/SubtractedVolumeBounds.h	(working copy)
@@ -101,7 +101,7 @@
 
  inline bool SubtractedVolumeBounds::inside(const Amg::Vector3D &pos, double tol) const
  { 
-   return (m_outer->inside(pos,tol) && !m_inner->inside(pos,0.) );
+   return (m_outer->inside(pos,tol) && !m_inner->inside(pos,-tol) );
  }
 
  inline Volume* SubtractedVolumeBounds::outer() const { return m_outer; }

Tracking/TrkExtrapolation/TrkExSTEP_Propagator
There is no explicit changes between local copy and rev.

Tracking/TrkExtrapolation/TrkExTools
Index: Tracking/TrkExtrapolation/TrkExTools/src/Extrapolator.cxx
===================================================================
--- Tracking/TrkExtrapolation/TrkExTools/src/Extrapolator.cxx	(revision 617387)
+++ Tracking/TrkExtrapolation/TrkExTools/src/Extrapolator.cxx	(working copy)
@@ -1153,7 +1153,9 @@
   Amg::Vector3D gp = parm.position();
   if ( vol && vol->inside(gp,m_tolerance) ) {
     staticVol = vol; 
+    std::cout << vol->volumeName() <<" assigned as input static at:"<< gp.perp()<<","<<gp.z()<<","<< staticVol->geometrySignature()<< std::endl;
   } else {
+    if (vol) std::cout << vol->volumeName() <<"not correctly assigned:"<< gp.perp()<<","<<gp.z()<< std::endl;
     staticVol =  m_navigator->trackingGeometry()->lowestStaticTrackingVolume(gp);    
     const Trk::TrackingVolume* nextStatVol = 0;
     if ( m_navigator->atVolumeBoundary(currPar,staticVol,dir,nextStatVol,m_tolerance) && nextStatVol != staticVol ) 
@@ -1525,7 +1527,9 @@
      // check missing volume boundary
      if (nextPar && !(m_currentDense->inside(nextPar->position(),m_tolerance)
 		      || m_navigator->atVolumeBoundary(nextPar,m_currentDense,dir,assocVol,m_tolerance) ) ) {
-       ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() );
+       ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() <<":inside volume:"<<
+		      m_currentDense->inside(nextPar->position(),m_tolerance)<<" at volume boundary:"<<
+		      m_navigator->atVolumeBoundary(nextPar,m_currentDense,dir,assocVol,m_tolerance) );
        if ( m_currentDense->zOverAtimesRho() != 0.) {     
          ATH_MSG_DEBUG( "  [!] ERROR: trying to recover: repeat the propagation step in"<< m_highestVolume->volumeName() );
          m_currentDense = m_highestVolume;
Index: Tracking/TrkExtrapolation/TrkExTools/src/TimedExtrapolator.cxx
===================================================================
--- Tracking/TrkExtrapolation/TrkExTools/src/TimedExtrapolator.cxx	(revision 617388)
+++ Tracking/TrkExtrapolation/TrkExTools/src/TimedExtrapolator.cxx	(working copy)
@@ -784,7 +784,10 @@
      // check missing volume boundary
      if (nextPar && !(m_currentDense->inside(nextPar->position(),m_tolerance) 
 		      || m_navigator->atVolumeBoundary(nextPar,m_currentDense,dir,assocVol,m_tolerance) ) ) {
-       ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() );
+       //ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() );
+       ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() <<":inside volume:"<<
+		      m_currentDense->inside(nextPar->position(),m_tolerance)<<" at volume boundary:"<<
+		      m_navigator->atVolumeBoundary(nextPar,m_currentDense,dir,assocVol,m_tolerance) );
      }
      //  if ( m_currentDense->zOverAtimesRho() != 0.) {     
      //    ATH_MSG_DEBUG( "  [!] ERROR: trying to recover: repeat the propagation step in"<< m_highestVolume->volumeName() );
@@ -906,15 +909,16 @@
 	 if ( dIter != m_denseVols.end() ) {
 	   currVol = (*dIter).first;
 	   nextVol = ( (*dIter).first->boundarySurfaces())[index].getPtr()->attachedVolume(*nextPar,dir);
-           // the boundary orientation is not reliable
-	   Amg::Vector3D tp = nextPar->position()+2*m_tolerance*dir*nextPar->momentum().normalized();
-           if (!nextVol || !nextVol->inside(tp,m_tolerance) ) {   // search for dense volumes
+           // the boundary orientation is not reliable : step few microns in to see if good choice
+           // to be removed when boundaries smart enough ...
+	   Amg::Vector3D tp = nextPar->position()+10*m_tolerance*dir*nextPar->momentum().normalized();
+           if (!nextVol || !nextVol->inside(tp,0.) ) {   // search for dense volumes
 	     m_currentDense =  m_highestVolume;
 	     if (m_dense && !m_denseVols.size()) m_currentDense = m_currentStatic;
 	     else {
 	       for (unsigned int i=0;i<m_denseVols.size(); i++) {
 		 const Trk::TrackingVolume* dVol = m_denseVols[i].first;
-		 if ( dVol->inside(tp,m_tolerance)  && dVol->zOverAtimesRho()!=0. ){
+		 if ( dVol->inside(tp,0.)  && dVol->zOverAtimesRho()!=0. ){
 		   m_currentDense = dVol;
 		   ATH_MSG_DEBUG( "  [+] Next dense volume found: '" << m_currentDense->volumeName() << "'."); 
 		   break;
@@ -1331,9 +1335,12 @@
       nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
       return currPar;
     }
-    ATH_MSG_DEBUG( "  [+] current static volume resolved by navigator at boundary: "<< positionOutput(currPar->position())
-		   <<m_currentStatic->volumeName()<<"->"<< nextVol->volumeName()<<" with tolerance " << m_tolerance );
-    m_currentStatic = nextVol;
+    // boundary orientation not reliable
+    if ( nextVol->inside(currPar->position()+0.01*dir*currPar->momentum().normalized(),0.) ) {
+      ATH_MSG_DEBUG( "  [+] current static volume resolved by navigator at boundary: "<< positionOutput(currPar->position())
+		     <<m_currentStatic->volumeName()<<"->"<< nextVol->volumeName()<<" with tolerance " << m_tolerance );
+      m_currentStatic = nextVol;
+    }
   }
 
   // current frame volume known-retrieve geoID
@@ -1592,6 +1599,8 @@
   // resolve the use of dense volumes
   m_dense = (m_currentStatic->geometrySignature()==Trk::MS && m_useMuonMatApprox ) || (m_currentStatic->geometrySignature()!=Trk::MS && m_useDenseVolumeDescription );
 
+  // std::cout <<"current geometry signature:"<< m_currentStatic->geometrySignature() <<","<< m_useMuonMatApprox<<","<<m_useDenseVolumeDescription<< std::endl;
+
   // reset remaining counters
   m_currentDense = m_dense ?  m_currentStatic : m_highestVolume;
   m_navigBoundaries.clear(); 
@@ -1689,18 +1698,26 @@
      Amg::Vector3D nextPos = currPar->position()+dir*currPar->momentum().normalized()*m_trSurfs[sols[is]].second;
      //Amg::Vector3D halfStep = nextPos - 0.5*step*dir*currPar->momentum().normalized();
 
-     // check missing volume boundary 
-     if ( !(m_currentDense->inside(nextPos,m_tolerance) ) ) {
-       ATH_MSG_DEBUG( "  [!] ERROR: missing volume boundary for volume"<< m_currentDense->volumeName() );
-       // new search
-       m_currentDense =  m_highestVolume;
-       for (unsigned int i=0;i<m_denseVols.size(); i++) {
-	 const Trk::TrackingVolume* dVol = m_denseVols[i].first;
-	 if ( dVol->inside(nextPos,m_tolerance)  && dVol->zOverAtimesRho()!=0. ) m_currentDense = dVol; 
+     // check missing volume boundary :
+     if ( !(m_currentDense->inside(nextPos,m_tolerance)) && step>m_tolerance ) {
+       ATH_MSG_DEBUG( "  [!] WARNING: missing volume boundary for volume"<< m_currentDense->volumeName()<<" after step:"<< step );
+       unsigned il=0;
+       for ( ; il<int(fabs(step)/0.001); il++) {
+         if ( m_currentDense->inside(nextPos-il*0.001*dir*currPar->momentum().normalized(),0.) ) {
+	   ATH_MSG_DEBUG("last seen inside volume before:"<< il<<" microns");
+           break;
+	 }
        }
-       if (m_dense && m_currentDense==m_highestVolume) m_currentDense = m_currentStatic;
-
-       ATH_MSG_DEBUG( "  [!] new search for dense volume : "<< m_currentDense->volumeName() );
+       if (il>10) {     // reassign current dense - use 10 microns tolerance ( limited by distance calculation ?)
+         // new search
+         m_currentDense =  m_highestVolume;
+         for (unsigned int i=0;i<m_denseVols.size(); i++) {
+	   const Trk::TrackingVolume* dVol = m_denseVols[i].first;
+	   if ( dVol->inside(nextPos,m_tolerance)  && dVol->zOverAtimesRho()!=0. ) m_currentDense = dVol; 
+         }
+         if (m_dense && m_currentDense==m_highestVolume) m_currentDense = m_currentStatic;
+	 ATH_MSG_DEBUG( "  [!] new search for dense volume : "<< m_currentDense->volumeName() );
+       }
      }
      
      double tDelta = step/beta;
@@ -1726,6 +1743,9 @@
 
      //std::cout << "looping over intersections:"<<is<<","<< m_trSurfs[sols[is]].second<<","<<step << ","<< tDelta<<","<<mDelta << std::endl;
 
+     if (fr*mDelta>0 && m_currentDense->averageZ()>0) ATH_MSG_VERBOSE( "collecting material from dense volume: "<< m_currentDense->volumeName()<<
+	     ":thickness[X0]"<<mDelta<<" at position[R,z]="<<nextPos.perp()<<","<<nextPos.z() );  
+
      if (fr<1.) { // decay or material interaction during the step
 
        int process = frT < frM ? timeLim.process : m_path.process;
@@ -1772,7 +1792,7 @@
        // use global coordinates to retrieve attached volume (just for static!)
        nextVol = (m_currentStatic->boundarySurfaces())[index].getPtr()->attachedVolume(nextPar->position(),nextPar->momentum(),dir);
        // double check the next volume
-       if ( nextVol && !(nextVol->inside(nextPar->position()+0.01*nextPar->momentum().normalized(),m_tolerance) ) ) {
+       if ( nextVol && !(nextVol->inside(nextPar->position()+0.01*nextPar->momentum().normalized(),0.) ) ) {
 	 ATH_MSG_DEBUG( "  [!] WARNING: wrongly assigned static volume ?"<< m_currentStatic->volumeName()<<"->" << nextVol->volumeName() );
 	 nextVol = m_navigator->trackingGeometry()->lowestStaticTrackingVolume(nextPar->position()+0.01*nextPar->momentum().normalized());
 	 if (nextVol) ATH_MSG_DEBUG( "  new search yields: "<< nextVol->volumeName() );          
@@ -1794,6 +1814,7 @@
 	   ATH_MSG_DEBUG( "  [+] Crossing position is         - at " <<  positionOutput(nextPar->position()) );
 	   if (!destVol && m_currentStatic->geometrySignature()!=nextVol->geometrySignature())
 	     { nextGeoID=nextVol->geometrySignature(); return nextPar; }   
+	   m_currentStatic = nextVol;
 	 }     
 	 return transportToVolumeWithPathLimit(*nextPar, timeLim, dir, particle, nextGeoID, destVol);
        }
@@ -1844,23 +1865,29 @@
 	 currVol = (*dIter).first;
 
          if ( m_navigator->trackingGeometry()->atVolumeBoundary(nextPos,nextPar->momentum(), currVol,assocVol, dir,m_tolerance) ) {
-            if ( assocVol && assocVol->zOverAtimesRho()!= 0. ) m_currentDense = assocVol;
-            else if ( currVol->inside(nextPos+0.002*dir*nextPar->momentum().normalized()) ) m_currentDense = currVol;
-            else {
+	   if ( assocVol && assocVol->zOverAtimesRho()!= 0. && assocVol->inside(nextPos+0.002*dir*nextPar->momentum().normalized(),0.) ) {
+	     m_currentDense = assocVol;
+	     ATH_MSG_VERBOSE( " dense volume assigned by boundary association :"<< m_currentDense->volumeName() );
+	   } else if ( currVol->inside(nextPos+0.002*dir*nextPar->momentum().normalized(),0.) ) {
+	     m_currentDense = currVol;
+	     ATH_MSG_VERBOSE( " dense volume assigned by continuation (current) :"<< m_currentDense->volumeName() );
+	   } else {
 	      // new search
 	      m_currentDense =  m_highestVolume;
-	      if (m_useMuonMatApprox && !m_denseVols.size()) m_currentDense = m_currentStatic;
+	      if (m_currentStatic->geometrySignature()==Trk::MS && m_useMuonMatApprox && !m_denseVols.size()) m_currentDense = m_currentStatic;
 	      else {
 		for (unsigned int i=0;i<m_denseVols.size(); i++) {
 		  const Trk::TrackingVolume* dVol = m_denseVols[i].first;
-		  if ( dVol->inside(nextPos+0.002*dir*nextPar->momentum().normalized(),m_tolerance)  && dVol->zOverAtimesRho()!=0. ) m_currentDense = dVol; 
+		  if ( dVol->inside(nextPos+0.002*dir*nextPar->momentum().normalized(),0.) && dVol->zOverAtimesRho()!=0. ) m_currentDense = dVol; 
 		}
 	      }
+	      ATH_MSG_VERBOSE( " dense volume assigned by new search:"<< m_currentDense->volumeName() );           
 	    }
 	 }
        }
      } else {   // detached volume bounds - not relevant ?
-       
+       ATH_MSG_VERBOSE( " detached volume bouns: not relevant ? " );                  
+
      }
 
      throwIntoGarbageBin(nextPar);

Configuration

The

the default builds just endcap toroid for dead material and BME stations for active volumes, you can easily change the configuration ). The same script allows you to change the output level for different tools ( for debugging ) and to switch on the VP1 to visualize both geometries ( GeoModel under 'Geo', TrackingGeometry under 'TrkGeo', after loading plugin "All studies", when asked for the name of TrackingGeometrySvc and TrackingGeometryName, fill in : ISF_FatrasTrackingGeometrySvc and ISF_FatrasTrackingGeometry ).

The input tracks are shooted in ranges defined in mtg/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py.

The ntuples contain tree "particles" which gives you, among others, the amount of the material traversed by the track (X0) in a given subdetector (geoID=1 (ID), 3(Calo), 4(MS) ). This is the info used to plot the integrated material maps ( as function of eta and phi ).

Please try to produce these ntuples and make some rough comparisons ( for some of discrepancies you will see I have a fix already, some are not yet understood, we can discussed this next week over skype ).

Let me know if something does not work , or if I forgot to explain, ok ?

-- JordiDuarte - 02 Oct 2014

Topic attachments
I Attachment History Action Size Date Who Comment
Unknown file formatdiff isfgeant4tools.diff r1 manage 54.3 K 2014-10-05 - 14:06 JordiDuarte Diff between local copy and rev. copy of the Gean4Tools packages
Texttxt jobOptions_ConfGetter.py.txt r1 manage 8.3 K 2014-10-05 - 14:21 JordiDuarte Main configuration file for validation
Edit | Attach | Watch | Print version | History: r8 | r6 < r5 < r4 < r3 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r4 - 2014-10-05 - JordiDuarte
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Sandbox All webs login

  • Edit
  • Attach
This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2021 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback