Allow the user to change the material of the surrounding volume from

air to whatever is required. Is an optional parameter to the telescope
constructor.

Author: simonspa

README.md

@@ -11,6 +11,7 @@ This small collection of scripts provides a simple interface for the simulation
 * Error propagation of track extrapolation uncertainty using GBL
 * Includes scattering in material, estimated via the PDG Highland formula
 * Automatically accounts for air between the telescope planes
+* Allows to exchange the volume material considered for scattering (Air, vacuum, Helium...)
 * Automatic placement of the thin scatterers at correct positions
 * Planes ordered automatically in `z` for correct GBL trajectory building
 * Radiation length for some common materials are defined in `utils/materials.h`
@@ -69,4 +70,10 @@ This small collection of scripts provides a simple interface for the simulation
   double MIM26 = 50e-3 / X0_Si + 50e-3 / X0_Kapton;
   ```
 
-* The resolution is always given as intrinsic resolution of the respective sensor in units of micrometer.
+* The resolution should always be given as intrinsic resolution of the respective sensor in units of micrometer.
+
+* The constructor of the telescope class takes the radiation length of the surrounding volume as optional parameter:
+
+`telescope(std::vector<gblsim::plane> planes, double beam_energy, double material = X0_Air);`
+
+It defaults to the radiation length of dry air but can be replaced with other materials or with vacuum (`X0 = 0`) for comparison.

telescope/assembly.cc

@@ -25,7 +25,8 @@ plane::plane(double position, double material, bool measurement, double resoluti
   m_position(position) {}
 
 
-telescope::telescope(std::vector<gblsim::plane> planes, double beam_energy) :
+telescope::telescope(std::vector<gblsim::plane> planes, double beam_energy, double material) :
+  m_volumeMaterial(material),
   m_listOfPoints(),
   m_listOfLabels()
 {
@@ -63,24 +64,32 @@ telescope::telescope(std::vector<gblsim::plane> planes, double beam_energy) :
     // Let's first add the air:
     double plane_distance = pl->m_position - oldpos;
     LOG(logDEBUG2) << "Distance to next plane: " << plane_distance;
+    double distance = 0;
 
-    // Propagate [mm] from 0 to 0.21 = 0.5 - 1/sqrt(12)
-    double distance = 0.21 * plane_distance; arclength += distance;
+    // Check if a volume scatterer with radiation length != 0 has been defined:
+    if(m_volumeMaterial > 0.0) {
+      // Propagate [mm] from 0 to 0.21 = 0.5 - 1/sqrt(12)
+      distance = 0.21 * plane_distance; arclength += distance;
 
-    // Add air:
-    m_listOfPoints.push_back(getPoint(distance,getScatterer(beam_energy,0.5*plane_distance/X0_Air,total_materialbudget)));
-    LOG(logDEBUG3) << "Added air scat at " << arclength;
+      // Add volume scatterer:
+      m_listOfPoints.push_back(getPoint(distance,getScatterer(beam_energy,0.5*plane_distance/m_volumeMaterial,total_materialbudget)));
+      LOG(logDEBUG3) << "Added volume scat at " << arclength;
 
-    // Propagate [mm] 0.58 = from 0.21 to 0.79 = 0.5 + 1/sqrt(12)
-    distance = 0.58 * plane_distance; arclength += distance;
+      // Propagate [mm] 0.58 = from 0.21 to 0.79 = 0.5 + 1/sqrt(12)
+      distance = 0.58 * plane_distance; arclength += distance;
 
-    // Factor 0.5 for the air as it is split into two scatterers:
-    m_listOfPoints.push_back(getPoint(distance,getScatterer(beam_energy,0.5*plane_distance/X0_Air,total_materialbudget)));
-    LOG(logDEBUG3) << "Added air scat at " << arclength;
+      // Factor 0.5 for the volume as it is split into two scatterers:
+      m_listOfPoints.push_back(getPoint(distance,getScatterer(beam_energy,0.5*plane_distance/m_volumeMaterial,total_materialbudget)));
+      LOG(logDEBUG3) << "Added volume scat at " << arclength;
 
-    // Propagate [mm] from 0 to 0.21 = 0.5 - 1/sqrt(12)
-    distance = 0.21 * plane_distance; arclength += distance;
-    LOG(logDEBUG) << "Added air scatterers.";
+      // Propagate [mm] from 0 to 0.21 = 0.5 - 1/sqrt(12)
+      distance = 0.21 * plane_distance; arclength += distance;
+      LOG(logDEBUG) << "Added volume scatterers.";
+    }
+    else {
+      // No volume scatterer defined (vacuum), simply propagate to the next plane:
+      distance = plane_distance;
+    }
 
     if(pl->m_measurement) {
       m_listOfPoints.push_back(getPoint(distance,pl->m_resolution,getScatterer(beam_energy,pl->m_materialbudget,total_materialbudget)));
@@ -111,11 +120,13 @@ double telescope::getTotalMaterialBudget(std::vector<gblsim::plane> planes) {
     total_materialbudget += p->m_materialbudget;
   }
 
-  // Add the air as scattering material:
-  double total_distance = (planes.back().m_position - planes.front().m_position);
-  LOG(logDEBUG2) << "Adding x/X0=" << (total_distance/X0_Air) << " (air)";
-  total_materialbudget += total_distance/X0_Air;
-
+  if(m_volumeMaterial > 0.0) {
+    // Add the air as scattering material:
+    double total_distance = (planes.back().m_position - planes.front().m_position);
+    LOG(logDEBUG2) << "Adding x/X0=" << (total_distance/m_volumeMaterial) << " (air)";
+    total_materialbudget += total_distance/m_volumeMaterial;
+  }
+  
   LOG(logDEBUG) << "Total track material budget x/X0=" << total_materialbudget;
   return total_materialbudget;
 }

telescope/assembly.h

@@ -1,4 +1,5 @@
 #include "GblTrajectory.h"
+#include "materials.h"
 
 namespace gblsim {
   class plane {
@@ -25,7 +26,7 @@ namespace gblsim {
 
   class telescope {
   public:
-    telescope(std::vector<gblsim::plane> planes, double beam_energy);
+    telescope(std::vector<gblsim::plane> planes, double beam_energy, double material = X0_Air);
 
     // Return the trajectory
     gbl::GblTrajectory getTrajectory();
@@ -35,6 +36,9 @@ namespace gblsim {
 
     void printLabels();
   private:
+    // Radiationlength of the material of the surrounding volume, defaults to dry air:
+    double m_volumeMaterial;
+    
     double getTotalMaterialBudget(std::vector<gblsim::plane> planes);
     std::vector<gbl::GblPoint> m_listOfPoints;
     std::vector<int> m_listOfLabels;