Use BaseParticlePropagator also for propagation to HGCal for phase-2

cms-sw · Oct 16, 2020 · 5363318 · 5363318
1 parent 2ad8a7b
commit 5363318
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 98 deletions.
diff --git a/RecoTauTag/RecoTau/BuildFile.xml b/RecoTauTag/RecoTau/BuildFile.xml
@@ -23,7 +23,6 @@
 <use name="PhysicsTools/JetMCUtils"/>
 <use name="CommonTools/Utils"/>
 <use name="CommonTools/BaseParticlePropagator"/>
-<use name="TrackPropagation/RungeKutta"/>
 <use name="RecoLocalCalo/HGCalRecAlgos"/>
 <use name="roottmva"/>
 <use name="PhysicsTools/TensorFlow"/>

diff --git a/RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h b/RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h
@@ -38,10 +38,10 @@ class PositionAtECalEntranceComputer {
   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeo_esToken_;
   double bField_z_;
   bool isPhase2_;
-  MagneticField const* bField_;
   hgcal::RecHitTools recHitTools_;
   float hgcalFace_z_;
   static constexpr float ecalBarrelEndcapEtaBorder_ = 1.479;
+  static constexpr float hgcalHfEtaBorder_ = 3.0;
 };
 
 #endif  // RecoTauTag_RecoTau_PositionAtECalEntranceComputer_h
diff --git a/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc b/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc
@@ -5,78 +5,9 @@
 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
-#include "DataFormats/GeometrySurface/interface/Plane.h"
-#include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
 
 #include <cassert>
 
-//HGCal helper classes
-//MB: looks be copy of HGCal utils: L1Trigger/L1THGCalUtilities/plugins/ntuples/HGCalTriggerNtupleGen.cc
-namespace hgcal_helpers {
-  class SimpleTrackPropagator {
-  public:
-    SimpleTrackPropagator(MagneticField const* f) : field_(f), prod_(field_, alongMomentum), absz_target_(0) {
-      ROOT::Math::SMatrixIdentity id;
-      AlgebraicSymMatrix55 C(id);
-      //MB: To define uncertainty of starting point of trajectory propagation scale identity matrix created above by 0.001
-      C *= 0.001;
-      err_ = CurvilinearTrajectoryError(C);
-    }
-    void setPropagationTargetZ(const float& z);
-
-    bool propagate(const double& px,
-                   const double& py,
-                   const double& pz,
-                   const double& x,
-                   const double& y,
-                   const double& z,
-                   const float& charge,
-                   reco::Candidate::Point& output) const;
-
-  private:
-    SimpleTrackPropagator() : field_(nullptr), prod_(field_, alongMomentum), absz_target_(0) {}
-    const RKPropagatorInS& RKProp() const { return prod_.propagator; }
-    Plane::PlanePointer targetPlaneForward_, targetPlaneBackward_;
-    MagneticField const* field_;
-    CurvilinearTrajectoryError err_;
-    defaultRKPropagator::Product prod_;
-    float absz_target_;
-  };
-  void SimpleTrackPropagator::setPropagationTargetZ(const float& z) {
-    targetPlaneForward_ = Plane::build(Plane::PositionType(0, 0, std::abs(z)), Plane::RotationType());
-    targetPlaneBackward_ = Plane::build(Plane::PositionType(0, 0, -std::abs(z)), Plane::RotationType());
-    absz_target_ = std::abs(z);
-  }
-  bool SimpleTrackPropagator::propagate(const double& px,
-                                        const double& py,
-                                        const double& pz,
-                                        const double& x,
-                                        const double& y,
-                                        const double& z,
-                                        const float& charge,
-                                        reco::Candidate::Point& output) const {
-    typedef TrajectoryStateOnSurface TSOS;
-    GlobalPoint startingPosition(x, y, z);
-    GlobalVector startingMomentum(px, py, pz);
-    Plane::PlanePointer startingPlane = Plane::build(Plane::PositionType(x, y, z), Plane::RotationType());
-    TSOS startingStateP(
-        GlobalTrajectoryParameters(startingPosition, startingMomentum, charge, field_), err_, *startingPlane);
-    TSOS trackStateP;
-    if (pz > 0) {
-      trackStateP = RKProp().propagate(startingStateP, *targetPlaneForward_);
-    } else {
-      trackStateP = RKProp().propagate(startingStateP, *targetPlaneBackward_);
-    }
-    if (trackStateP.isValid()) {
-      output.SetXYZ(
-          trackStateP.globalPosition().x(), trackStateP.globalPosition().y(), trackStateP.globalPosition().z());
-      return true;
-    }
-    output.SetXYZ(0, 0, 0);
-    return false;
-  }
-}  // namespace hgcal_helpers
-
 PositionAtECalEntranceComputer::PositionAtECalEntranceComputer(edm::ConsumesCollector&& cc, bool isPhase2)
     : bField_esToken_(cc.esConsumes<MagneticField, IdealMagneticFieldRecord>()),
       caloGeo_esToken_(cc.esConsumes<CaloGeometry, CaloGeometryRecord>()),
@@ -92,8 +23,7 @@ PositionAtECalEntranceComputer::PositionAtECalEntranceComputer(edm::ConsumesColl
 PositionAtECalEntranceComputer::~PositionAtECalEntranceComputer() {}
 
 void PositionAtECalEntranceComputer::beginEvent(const edm::EventSetup& es) {
-  bField_ = &es.getData(bField_esToken_);
-  bField_z_ = bField_->inTesla(GlobalPoint(0., 0., 0.)).z();
+  bField_z_ = es.getData(bField_esToken_).inTesla(GlobalPoint(0., 0., 0.)).z();
   if (isPhase2_) {
     recHitTools_.setGeometry(es.getData(caloGeo_esToken_));
     hgcalFace_z_ = recHitTools_.getPositionLayer(1).z();  // HGCal 1st layer
@@ -104,36 +34,32 @@ reco::Candidate::Point PositionAtECalEntranceComputer::operator()(const reco::Ca
                                                                   bool& success) const {
   assert(bField_z_ != -1.);
   reco::Candidate::Point position;
+  BaseParticlePropagator propagator = BaseParticlePropagator(
+      RawParticle(particle->p4(),
+                  math::XYZTLorentzVector(particle->vertex().x(), particle->vertex().y(), particle->vertex().z(), 0.),
+                  particle->charge()),
+      0.,
+      0.,
+      bField_z_);
   if (!isPhase2_ || std::abs(particle->eta()) < ecalBarrelEndcapEtaBorder_) {  // ECal
-    BaseParticlePropagator propagator = BaseParticlePropagator(
-        RawParticle(particle->p4(),
-                    math::XYZTLorentzVector(particle->vertex().x(), particle->vertex().y(), particle->vertex().z(), 0.),
-                    particle->charge()),
-        0.,
-        0.,
-        bField_z_);
     propagator.propagateToEcalEntrance(false);
-    if (propagator.getSuccess() != 0) {
-      position = propagator.particle().vertex().Vect();
-      success = true;
-    } else {
+  } else {  // HGCal
+    if (std::abs(particle->vertex().z()) >= hgcalFace_z_) {
       success = false;
+      return position;
     }
-  } else {  // HGCal
+    propagator.setPropagationConditions(152.6, hgcalFace_z_, false);
+    propagator.propagate();
+  }
+  if (propagator.getSuccess() != 0) {
+    position = propagator.particle().vertex().Vect();
+    if (std::abs(position.eta()) > hgcalHfEtaBorder_) {
+      success = false;
+    } else {
+      success = true;
+    }
+  } else {
     success = false;
-    if (std::abs(particle->vertex().z()) >= hgcalFace_z_)
-      return position;
-
-    hgcal_helpers::SimpleTrackPropagator propagator(bField_);
-    propagator.setPropagationTargetZ(hgcalFace_z_);
-    success = propagator.propagate(particle->px(),
-                                   particle->py(),
-                                   particle->pz(),
-                                   particle->vertex().x(),
-                                   particle->vertex().y(),
-                                   particle->vertex().z(),
-                                   particle->charge(),
-                                   position);
   }
   return position;
 }