SimplifiedGeometryPropagator/src/LayerNavigator.cc

0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/LayerNavigator.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0003
0004 #include <vector>
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Geometry.h"
0009 #include "FastSimulation/SimplifiedGeometryPropagator/interface/BarrelSimplifiedGeometry.h"
0010 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ForwardSimplifiedGeometry.h"
0011 #include "FastSimulation/SimplifiedGeometryPropagator/interface/LayerNavigator.h"
0012 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Trajectory.h"
0013 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0014
0015 /**
0016 // find the next layer that the particle will cross
0017 //
0018 // motivation for a new algorithm
0019 //    - old algorithm is flawed
0020 //    - the new algorithm allows to put material and instruments on any plane perpendicular to z, or on any cylinder with the z-axis as axis
0021 //    - while the old algorith, with the requirement of nested layers, forbids the introduction of long narrow cylinders, required for a decent simulation of material in front of HF
0022 //
0023 // definitions
0024 //    the geometry is described by 2 sets of layers:
0025 //    - forward layers:
0026 //          flat layers, perpendicular to the z-axis, positioned at a given z
0027 //          these layers have material / instruments between a given materialMinR and materialMaxR
0028 //          no 2 forward layers should have the same z-position
0029 //    - barrel layers:
0030 //          cylindrically shaped layers, with the z-axis as axis, infinitely long
0031 //          these layers have material / instruments for |z| < materialMaxAbsZ
0032 //          no 2 barrel layers should have the same radius
0033 //    - forward(barrel) layers are ordered according to increasing z (r)
0034
0035 // principle
0036 //    - neutral particles follow a straight trajectory
0037 //    - charged particles follow a helix-shaped trajectory:
0038 //          constant speed along the z-axis
0039 //          circular path in the x-y plane
0040 //    => the next layer that the particle will cross is among the following 3 layers
0041 //    - closest forward layer with
0042 //         - z >(<) particle.z() for particles moving in the positive(negative) direction
0043 //    - closest barrel layer with r < particle.r
0044 //    - closest barrel layer with r > particle.r
0045
0046 // algorithm
0047 //    - find the 3 candidate layers
0048 //    - find the earliest positive intersection time for each of the 3 candidate layers
0049 //    - move the particle to the earliest intersection time
0050 //    - select and return the layer with the earliest positive intersection time
0051 //
0052 // notes
0053 //    - the implementation of the algorithm can probably be optimised, e.g.
0054 //       - one can probably gain time in moveToNextLayer if LayerNavigator is aware of the candidate layers of the previous call to moveToNextLayer
0055 //       - for straight tracks, the optimal strategy to find the next layer might be very different
0056 **/
0057
0058 const std::string fastsim::LayerNavigator::MESSAGECATEGORY = "FastSimulation";
0059
0060 fastsim::LayerNavigator::LayerNavigator(const fastsim::Geometry& geometry)
0061     : geometry_(&geometry),
0062       nextBarrelLayer_(nullptr),
0063       previousBarrelLayer_(nullptr),
0064       nextForwardLayer_(nullptr),
0065       previousForwardLayer_(nullptr) {
0066   ;
0067 }
0068
0069 bool fastsim::LayerNavigator::moveParticleToNextLayer(fastsim::Particle& particle,
0070                                                       const fastsim::SimplifiedGeometry*& layer) {
0071   LogDebug(MESSAGECATEGORY) << "   moveToNextLayer called";
0072
0073   // if the layer is provided, the particle must be on it
0074   if (layer) {
0075     if (!particle.isOnLayer(layer->isForward(), layer->index())) {
0076       throw cms::Exception("FastSimulation") << "If layer is provided, particle must be on layer."
0077                                              << "\n   Layer: " << *layer << "\n   Particle: " << particle;
0078     }
0079   }
0080
0081   // magnetic field at the current position of the particle
0082   // considered constant between the layers
0083   double magneticFieldZ =
0084       layer ? layer->getMagneticFieldZ(particle.position()) : geometry_->getMagneticFieldZ(particle.position());
0085   LogDebug(MESSAGECATEGORY) << "   magnetic field z component:" << magneticFieldZ;
0086
0087   // particle moves inwards?
0088   bool particleMovesInwards =
0089       particle.momentum().X() * particle.position().X() + particle.momentum().Y() * particle.position().Y() < 0;
0090
0091   //
0092   //  update nextBarrelLayer and nextForwardLayer
0093   //
0094
0095   ////////////
0096   // first time method is called
0097   ////////////
0098   if (!layer) {
0099     LogDebug(MESSAGECATEGORY) << "      called for first time";
0100
0101     // find the narrowest barrel layers with
0102     // layer.r > particle.r (the closest layer with layer.r < particle.r will then be considered, too)
0103     // assume barrel layers are ordered with increasing r
0104     for (const auto& layer : geometry_->barrelLayers()) {
0105       if (particle.isOnLayer(false, layer->index()) ||
0106           std::abs(layer->getRadius() - particle.position().Rho()) < 1e-2) {
0107         if (particleMovesInwards) {
0108           nextBarrelLayer_ = layer.get();
0109           break;
0110         } else {
0111           continue;
0112         }
0113       }
0114
0115       if (particle.position().Pt() < layer->getRadius()) {
0116         nextBarrelLayer_ = layer.get();
0117         break;
0118       }
0119
0120       previousBarrelLayer_ = layer.get();
0121     }
0122
0123     //  find the forward layer with smallest z with
0124     //  layer.z > particle z (the closest layer with layer.z < particle.z will then be considered, too)
0125     for (const auto& layer : geometry_->forwardLayers()) {
0126       if (particle.isOnLayer(true, layer->index()) || std::abs(layer->getZ() - particle.position().Z()) < 1e-3) {
0127         if (particle.momentum().Z() < 0) {
0128           nextForwardLayer_ = layer.get();
0129           break;
0130         } else {
0131           continue;
0132         }
0133       }
0134
0135       if (particle.position().Z() < layer->getZ()) {
0136         nextForwardLayer_ = layer.get();
0137         break;
0138       }
0139
0140       previousForwardLayer_ = layer.get();
0141     }
0142   }
0143   ////////////
0144   // last move worked, let's update
0145   ////////////
0146   else {
0147     LogDebug(MESSAGECATEGORY) << "      ordinary call";
0148
0149     // barrel layer was hit
0150     if (layer == nextBarrelLayer_) {
0151       if (!particleMovesInwards) {
0152         previousBarrelLayer_ = nextBarrelLayer_;
0153         nextBarrelLayer_ = geometry_->nextLayer(nextBarrelLayer_);
0154       }
0155     } else if (layer == previousBarrelLayer_) {
0156       if (particleMovesInwards) {
0157         nextBarrelLayer_ = previousBarrelLayer_;
0158         previousBarrelLayer_ = geometry_->previousLayer(previousBarrelLayer_);
0159       }
0160     }
0161     // forward layer was hit
0162     else if (layer == nextForwardLayer_) {
0163       if (particle.momentum().Z() > 0) {
0164         previousForwardLayer_ = nextForwardLayer_;
0165         nextForwardLayer_ = geometry_->nextLayer(nextForwardLayer_);
0166       }
0167     } else if (layer == previousForwardLayer_) {
0168       if (particle.momentum().Z() < 0) {
0169         nextForwardLayer_ = previousForwardLayer_;
0170         previousForwardLayer_ = geometry_->previousLayer(previousForwardLayer_);
0171       }
0172     }
0173
0174     // reset layer
0175     layer = nullptr;
0176   }
0177
0178   ////////////
0179   // move particle to first hit with one of the enclosing layers
0180   ////////////
0181
0182   LogDebug(MESSAGECATEGORY) << "   particle between BarrelLayers: "
0183                             << (previousBarrelLayer_ ? previousBarrelLayer_->index() : -1) << "/"
0184                             << (nextBarrelLayer_ ? nextBarrelLayer_->index() : -1)
0185                             << " (total: " << geometry_->barrelLayers().size() << ")"
0186                             << "\n   particle between ForwardLayers: "
0187                             << (previousForwardLayer_ ? previousForwardLayer_->index() : -1) << "/"
0188                             << (nextForwardLayer_ ? nextForwardLayer_->index() : -1)
0189                             << " (total: " << geometry_->forwardLayers().size() << ")";
0190
0191   // calculate and store some variables related to the particle's trajectory
0192   std::unique_ptr<fastsim::Trajectory> trajectory = Trajectory::createTrajectory(particle, magneticFieldZ);
0193
0194   // push back all possible candidates
0195   std::vector<const fastsim::SimplifiedGeometry*> layers;
0196   if (nextBarrelLayer_) {
0197     layers.push_back(nextBarrelLayer_);
0198   }
0199   if (previousBarrelLayer_) {
0200     layers.push_back(previousBarrelLayer_);
0201   }
0202
0203   if (particle.momentum().Z() > 0) {
0204     if (nextForwardLayer_) {
0205       layers.push_back(nextForwardLayer_);
0206     }
0207   } else {
0208     if (previousForwardLayer_) {
0209       layers.push_back(previousForwardLayer_);
0210     }
0211   }
0212
0213   // calculate time until each possible intersection
0214   // -> pick layer that is hit first
0215   double deltaTimeC = -1;
0216   for (auto _layer : layers) {
0217     double tempDeltaTime =
0218         trajectory->nextCrossingTimeC(*_layer, particle.isOnLayer(_layer->isForward(), _layer->index()));
0219     LogDebug(MESSAGECATEGORY) << "   particle crosses layer " << *_layer << " in time " << tempDeltaTime;
0220     if (tempDeltaTime > 0 && (layer == nullptr || tempDeltaTime < deltaTimeC || deltaTimeC < 0)) {
0221       layer = _layer;
0222       deltaTimeC = tempDeltaTime;
0223     }
0224   }
0225
0226   // if particle decays on the way to the next layer, stop propagation there and return
0227   double properDeltaTimeC = deltaTimeC / particle.gamma();
0228   if (!particle.isStable() && properDeltaTimeC > particle.remainingProperLifeTimeC()) {
0229     // move particle in space, time and momentum until it decays
0230     deltaTimeC = particle.remainingProperLifeTimeC() * particle.gamma();
0231
0232     trajectory->move(deltaTimeC);
0233     particle.position() = trajectory->getPosition();
0234     particle.momentum() = trajectory->getMomentum();
0235
0236     particle.setRemainingProperLifeTimeC(0.);
0237
0238     // particle no longer is on a layer
0239     particle.resetOnLayer();
0240     LogDebug(MESSAGECATEGORY) << "    particle about to decay. Will not be moved all the way to the next layer.";
0241     return false;
0242   }
0243
0244   if (layer) {
0245     // move particle in space, time and momentum so it is on the next layer
0246     trajectory->move(deltaTimeC);
0247     particle.position() = trajectory->getPosition();
0248     particle.momentum() = trajectory->getMomentum();
0249
0250     if (!particle.isStable())
0251       particle.setRemainingProperLifeTimeC(particle.remainingProperLifeTimeC() - properDeltaTimeC);
0252
0253     // link the particle to the layer
0254     particle.setOnLayer(layer->isForward(), layer->index());
0255     LogDebug(MESSAGECATEGORY) << "    moved particle to layer: " << *layer;
0256   } else {
0257     // particle no longer is on a layer
0258     particle.resetOnLayer();
0259   }
0260
0261   // return true / false if propagations succeeded /failed
0262   LogDebug(MESSAGECATEGORY) << "    success: " << bool(layer);
0263   return layer;
0264 }