Merge pull request cms-sw#120 from cerati/canonize-prop

Canonize propagation and update functions

Config.cc

@@ -48,6 +48,13 @@ namespace Config
   seedOpts  seedInput    = simSeeds;
   cleanOpts seedCleaning = noCleaning; 
 
+  bool             finding_requires_propagation_to_hit_pos;
+  PropagationFlags finding_inter_layer_pflags;
+  PropagationFlags finding_intra_layer_pflags;
+  PropagationFlags backward_fit_pflags;
+  PropagationFlags forward_fit_pflags;
+  PropagationFlags seed_fit_pflags;
+
   bool  useCMSGeom = false;
   bool  readCmsswTracks = false;
 

Config.h

@@ -15,6 +15,43 @@ class TrackerInfo;
   #define CUDA_CALLABLE 
 #endif
 
+
+//------------------------------------------------------------------------------
+
+enum PropagationFlagsEnum
+{
+  PF_none              = 0,
+
+  PF_use_param_b_field = 0x1,
+  PF_apply_material    = 0x2
+};
+
+struct PropagationFlags
+{
+  union
+  {
+    struct
+    {
+      bool use_param_b_field  : 1;
+      bool apply_material     : 1;
+      // Could add: bool use_trig_approx  -- now Config::useTrigApprox = true
+      // Could add: int  n_iter : 8       -- now Config::Niter = 5
+    };
+
+    unsigned int _raw_;
+
+  };
+
+  PropagationFlags() : _raw_(0) {}
+
+  PropagationFlags(int pfe) :
+    use_param_b_field       ( pfe & PF_use_param_b_field),
+    apply_material          ( pfe & PF_apply_material)
+  {}
+};
+
+//------------------------------------------------------------------------------
+
 // Enum for input seed options
 enum seedOpts {simSeeds, cmsswSeeds, findSeeds};
 typedef std::map<std::string,seedOpts> seedOptsMap;
@@ -27,6 +64,8 @@ typedef std::map<std::string,cleanOpts> cleanOptsMap;
 enum matchOpts {trkParamBased, hitBased, labelBased};
 typedef std::map<std::string, matchOpts> matchOptsMap;
 
+//------------------------------------------------------------------------------
+
 namespace Config
 {
   extern TrackerInfo TrkInfo;
@@ -182,18 +221,26 @@ namespace Config
   constexpr float seed_d0cut   = 0.5f; // 5mm
   extern bool cf_seeding;
 
-  // Config for propagation
+  // Config for propagation - could/should enter into PropagationFlags?!
   constexpr int Niter    =  5;
   constexpr int NiterSim = 10; // Can make more steps due to near volume misses.
   constexpr bool useTrigApprox = true;
 
-  // Config for Bfield
+  // PropagationFlags as used during finding and fitting. Defined for each Geom in its plugin.
+  extern bool             finding_requires_propagation_to_hit_pos;
+  extern PropagationFlags finding_inter_layer_pflags;
+  extern PropagationFlags finding_intra_layer_pflags;
+  extern PropagationFlags backward_fit_pflags;
+  extern PropagationFlags forward_fit_pflags;
+  extern PropagationFlags seed_fit_pflags;
+
+  // Config for Bfield. Note: for now the same for CMS-2017 and CylCowWLids.
   constexpr float Bfield = 3.8112;
   constexpr float mag_c1 = 3.8114;
   constexpr float mag_b0 = -3.94991e-06;
   constexpr float mag_b1 = 7.53701e-06;
   constexpr float mag_a  = 2.43878e-11;
-  
+
   // Config for seeding as well... needed bfield
   constexpr float maxCurvR = (100 * minSimPt) / (sol * Bfield); // in cm
 

Geoms/CMS-2017.cc

@@ -12,12 +12,20 @@ namespace
   void Create_CMS_2017(TrackerInfo& ti, bool verbose)
   {
     Config::nTotalLayers     = 18 + 2 * 27;
+
     Config::useCMSGeom       = true;
     Config::nlayers_per_seed = 4;
     Config::maxCandsPerSeed  = 6;  // GC said 3 is enough ???
     Config::maxHolesPerCand  = 12; // should be reduced
     Config::chi2Cut          = 30.0;
 
+    Config::finding_requires_propagation_to_hit_pos = true;
+    Config::finding_inter_layer_pflags = PropagationFlags(PF_apply_material);
+    Config::finding_intra_layer_pflags = PropagationFlags(PF_none);
+    Config::backward_fit_pflags        = PropagationFlags(PF_use_param_b_field | PF_apply_material);
+    Config::forward_fit_pflags         = PropagationFlags(PF_use_param_b_field | PF_apply_material);
+    Config::seed_fit_pflags            = PropagationFlags(PF_none);
+
     ti.set_eta_regions(0.9, 1.7, 2.45, false);
     ti.create_layers(18, 27, 27);
 

Geoms/CylCowWLids.cc

@@ -211,6 +211,13 @@ namespace
   {
     Config::nTotalLayers     = 10 + 2 * 9;
 
+    Config::finding_requires_propagation_to_hit_pos = false;
+    Config::finding_inter_layer_pflags = PropagationFlags(PF_none);
+    Config::finding_intra_layer_pflags = PropagationFlags(PF_none);
+    Config::backward_fit_pflags        = PropagationFlags(PF_use_param_b_field);
+    Config::forward_fit_pflags         = PropagationFlags(PF_use_param_b_field);
+    Config::seed_fit_pflags            = PropagationFlags(PF_none);
+
     CylCowWLidsCreator creator(ti);
 
     creator.FillTrackerInfo();

Matrix.h

@@ -87,6 +87,7 @@ inline void sincos4(const float x, float& sin, float& cos)
    cos  = 1.f - 0.5f*x2 + 0.04166667f*x2*x2;
    sin  = x - 0.16666667f*x*x2;
 }
+
 //==============================================================================
 
 // This ifdef needs to be changed to something like "use matriplex" and/or

Propagation.cc

@@ -234,13 +234,13 @@ void HelixState::propagateErrors(const HelixState& in, float totalDistance, bool
 // each step travels for a path length equal to the safe step between the current position and the nearest object.
 
 TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom,
-                                     const bool useParamBfield)
+                                     const PropagationFlags pflags)
 {
   bool debug = true;
 
-  const HelixState hsin(inputState,useParamBfield);
+  const HelixState hsin(inputState, pflags.use_param_b_field);
   TrackState result(inputState);
-  HelixState hsout(result,useParamBfield);
+  HelixState hsout(result, pflags.use_param_b_field);
 
 #ifdef CHECKSTATEVALID
   if (!hsout.state.valid) {
@@ -341,17 +341,19 @@ TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom
 
 // XXMT4K The following is only used in buildtest.cc, maybe obsolete?
 
-// Propagate to the next obj
-// each step travels for a path length equal to the safe step between the current position and the nearest object.
-TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const bool useParamBfield)
+// Propagate to the next obj. Eeach step travels for a path length equal to
+// the safe step between the current position and the nearest object.
+
+TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom,
+                                 const PropagationFlags pflags)
 {
   bool debug = true;
 
   const VUSolid* target = geom.Layer(layer);
 
-  const HelixState hsin(inputState,useParamBfield);
+  const HelixState hsin(inputState, pflags.use_param_b_field);
   TrackState result(inputState);
-  HelixState hsout(result,useParamBfield);
+  HelixState hsout(result, pflags.use_param_b_field);
 
 #ifdef CHECKSTATEVALID
   if (!hsout.state.valid) {
@@ -405,13 +407,13 @@ TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometr
 // each step travels for a path lenght equal to delta r between the current position and the target radius. 
 // for track with pT>=1 GeV this converges to the correct path lenght in <5 iterations
 // derivatives need to be updated at each iteration
-TrackState propagateHelixToR(TrackState inputState, float r, const bool useParamBfield)
+TrackState propagateHelixToR(TrackState inputState, float r, const PropagationFlags pflags)
 {
   bool debug = false;
 
-  const HelixState hsin(inputState,useParamBfield);
+  const HelixState hsin(inputState, pflags.use_param_b_field);
   TrackState result(inputState);
-  HelixState hsout(result,useParamBfield);
+  HelixState hsout(result, pflags.use_param_b_field);
 
 #ifdef CHECKSTATEVALID
   if (!hsout.state.valid) {
@@ -470,8 +472,8 @@ TrackState propagateHelixToR(TrackState inputState, float r, const bool useParam
   return hsout.state;
 }
 
-TrackState propagateHelixToZ(TrackState inputState, float zout, const bool useParamBfield) {
-
+TrackState propagateHelixToZ(TrackState inputState, float zout, const PropagationFlags pflags)
+{
   TrackState result = inputState;
 
   const float z = inputState.z();
@@ -485,7 +487,7 @@ TrackState propagateHelixToZ(TrackState inputState, float zout, const bool usePa
   const float cosP = std::cos(phi);
   const float sinP = std::sin(phi);
 
-  const float k = inputState.charge*100.f/(-Config::sol*(useParamBfield?Config::BfieldFromZR(z,inputState.posR()):Config::Bfield));
+  const float k = inputState.charge*100.f/(-Config::sol*(pflags.use_param_b_field ? Config::BfieldFromZR(z,inputState.posR()) : Config::Bfield));
 
   const float s = (zout - z)/cosT;
   const float angPath = s*sinT*ipT/k;

Propagation.h

@@ -8,16 +8,16 @@
 // assuming radial direction (i.e. origin at (0,0))
 TrackState propagateLineToR(const TrackState& inputState, float r);
 
-TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom, const bool useParamBfield = false);
-TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const bool useParamBfield = false);
+TrackState propagateHelixToNextSolid(TrackState inputState, const Geometry& geom, const PropagationFlags pflags);
+TrackState propagateHelixToLayer(TrackState inputState, int layer, const Geometry& geom, const PropagationFlags pflags);
 
 // helix propagation in steps along helix trajectory. 
 // each step travels for a path lenght equal to delta r between the current position and the target radius. 
 // for track with pT>=1 GeV this converges to the correct path lenght in <5 iterations
 // derivatives need to be updated at each iteration
-TrackState propagateHelixToR(TrackState inputState, float r, const bool useParamBfield = false);
+TrackState propagateHelixToR(TrackState inputState, float r, const PropagationFlags pflags);
 
 //
-TrackState propagateHelixToZ(TrackState inputState, float z, const bool useParamBfield = false);
+TrackState propagateHelixToZ(TrackState inputState, float z, const PropagationFlags pflags);
 
 #endif

Simulation.cc

@@ -89,14 +89,16 @@ void setupTrackByToyMC(SVector3& pos, SVector3& mom, SMatrixSym66& covtrk,
   // XXMT4M - This should become while not in outer layer (or at least propState.state == invalid)
   // Try the invalid thingy first ... but would be good to also know what layers are final.
 
+  const PropagationFlags pflags(PF_use_param_b_field);
+
   for (int ihit = 0; ihit < Config::nTotHit; ++ihit)
   {
     dprintf("\n================================================================================\n");
     dprintf("=== Going for next hit %d from layer %d, xyrzphi=(%f,%f,%f,%f,%f)\n",
             ihit, simLayer, tmpState.x(), tmpState.y(), tmpState.posR(), tmpState.z(), tmpState.posPhi());
     dprintf("--------------------------------------------------------------------------------\n");
 
-    auto propState = propagateHelixToNextSolid(tmpState,geom,true);
+    auto propState = propagateHelixToNextSolid(tmpState, geom, pflags);
 
     float initX   = propState.parameters.At(0);
     float initY   = propState.parameters.At(1);

TrackerInfo.cc

@@ -21,8 +21,8 @@ void LayerInfo::set_selection_limits(float p1, float p2, float q1, float q2)
 void LayerInfo::set_r_hole_range(float rh1, float rh2)
 {
   m_has_r_range_hole = true;
-  m_hole_r2_min = rh1 * rh1;
-  m_hole_r2_max = rh2 * rh2;
+  m_hole_r_min = rh1;
+  m_hole_r_max = rh2;
 }
 
 

TrackerInfo.h

@@ -33,7 +33,7 @@ struct WSR_Result
 class LayerInfo
 {
 private:
-  bool  is_in_r2_hole(float r2) const { return r2 > m_hole_r2_min && r2 < m_hole_r2_max; }
+  bool  is_in_r_hole_no_check(float r) const { return r > m_hole_r_min && r < m_hole_r_max; }
 
 public:
   enum  LayerType_e { Undef = -1, Barrel = 0, EndCapPos = 1, EndCapNeg = 2 };
@@ -49,7 +49,7 @@ class LayerInfo
 
   bool          m_is_outer         = false;
   bool          m_has_r_range_hole = false;
-  float         m_hole_r2_min, m_hole_r2_max; // This could be turned into std::function when needed.
+  float         m_hole_r_min, m_hole_r_max; // This could be turned into std::function when needed.
 
   // Selection limits
   float         m_q_bin; // > 0 - bin width, < 0 - number of bins
@@ -83,8 +83,7 @@ class LayerInfo
   bool  is_within_r_limits(float r) const { return r > m_rin  && r < m_rout; }
   bool  is_within_q_limits(float q) const { return is_barrel() ? is_within_z_limits(q) : is_within_r_limits(q); }
 
-  bool  is_in_rsqr_hole(float r2)       const { return m_has_r_range_hole ? is_in_r2_hole(r2) : false; }
-  bool  is_in_xy_hole(float x, float y) const { return m_has_r_range_hole ? is_in_r2_hole(x*x + y*y): false; }
+  bool  is_in_r_hole      (float r) const { return m_has_r_range_hole ? is_in_r_hole_no_check(r) : false; }
 
   WSR_Result is_within_z_sensitive_region(float z, float dz) const
   {
@@ -100,9 +99,8 @@ class LayerInfo
     {
       if (m_has_r_range_hole)
       {
-        const float r2 = r*r;
-        if (r2 < m_hole_r2_max - dr && r2 > m_hole_r2_min + dr)  return WSR_Result(WSR_Outside, true);
-        if (r2 < m_hole_r2_max + dr && r2 > m_hole_r2_min - dr ) return WSR_Result(WSR_Edge,    true);
+        if (r < m_hole_r_max - dr && r > m_hole_r_min + dr)  return WSR_Result(WSR_Outside, true);
+        if (r < m_hole_r_max + dr && r > m_hole_r_min - dr ) return WSR_Result(WSR_Edge,    true);
       }
       return WSR_Result(WSR_Inside, false);
     }

buildtest.cc

@@ -176,14 +176,16 @@ void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_
   const auto& evt_lay_hits(ev.layerHits_);
   const auto& segLayMap(segmentMap[ilayer]);
 
+  const PropagationFlags pflags(PF_none);
+
   dprint("processing candidate with nHits=" << tkcand.nFoundHits());
 #ifdef LINEARINTERP
-  TrackState propState = propagateHelixToR(updatedState,ev.geom_.Radius(ilayer));
+  TrackState propState = propagateHelixToR(updatedState, ev.geom_.Radius(ilayer), pflags);
 #else
 #ifdef TBB
 #error "Invalid combination of options (thread safety)"
 #endif
-  TrackState propState = propagateHelixToLayer(updatedState,ilayer,ev.geom_);
+  TrackState propState = propagateHelixToLayer(updatedState, ilayer,ev.geom_, pflags);
 #endif // LINEARINTERP
 #ifdef CHECKSTATEVALID
   if (!propState.valid) {
@@ -221,7 +223,7 @@ void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_
     const float deltaR = maxR - minR;
     dprint("min, max, delta: " << minR << ", " << maxR << ", " << deltaR);
     const TrackState propStateMin = propState;
-    const TrackState propStateMax = propagateHelixToR(updatedState,maxR);
+    const TrackState propStateMax = propagateHelixToR(updatedState, maxR, pflags);
 #ifdef CHECKSTATEVALID
     if (!propStateMax.valid) {
       return;