PFlow fix/improvement in the HF region for phase2

CommonTools/BaseParticlePropagator/src/BaseParticlePropagator.cc

@@ -87,7 +87,6 @@ bool BaseParticlePropagator::propagate() {
       solution = tminus < 0 ? tplus : tminus;
       if (solution < 0.)
         return false;
-      zProp = particle_.Z() + particle_.Pz() * solution;
       success = 2;
     } else {
       success = 1;
@@ -390,7 +389,6 @@ bool BaseParticlePropagator::propagateToClosestApproach(double x0, double y0, bo
   if (dist2 > dist1) {
     particle_.setVertex(vertex1);
     particle_.setMomentum(momentum1.X(), momentum1.Y(), momentum1.Z(), momentum1.E());
-    dist2 = dist1;
   }
 
   // Done
@@ -506,11 +504,10 @@ bool BaseParticlePropagator::propagateToHcalEntrance(bool first) {
 
 bool BaseParticlePropagator::propagateToVFcalEntrance(bool first) {
   //
-  // Propagation to VFCAL entrance
-  // TODO: include proper geometry
-  // Geometry taken from DAQ TDR Chapter 13
+  // Propagation to VFCAL (HF) entrance
+  // Geometry taken from xml: Geometry/CMSCommonData/data/cms.xml
 
-  setPropagationConditions(400.0, 1110.0, first);
+  setPropagationConditions(400.0, 1114.95, first);
   propDir = 0;
   bool done = propagate();
   propDir = 1;

DataFormats/ParticleFlowReco/interface/PFTrajectoryPoint.h

@@ -52,10 +52,15 @@ namespace reco {
       HOLayer = 8,
       /// VFcal(HF) front face
       VFcalEntrance = 9,
-
-      NLayers = 10
+      // Number of valid layers
+      NLayers = 10,
+      // Unknown
+      Unknown = -1
     };
 
+    static const std::array<std::string, NLayers> layerTypeNames;
+    static LayerType layerTypeByName(const std::string& name);
+
     /// default constructor. Set variables at default dummy values
     PFTrajectoryPoint();
 

DataFormats/ParticleFlowReco/src/PFTrajectoryPoint.cc

@@ -1,8 +1,22 @@
 #include "DataFormats/ParticleFlowReco/interface/PFTrajectoryPoint.h"
 #include <ostream>
+#include <algorithm>
 
 using namespace reco;
 
+// To be kept in synch with the enumerator definitions in PFTrajectoryPoint.h file
+// Don't consider "Unknown" and "NLayers"
+std::array<std::string, PFTrajectoryPoint::NLayers> const PFTrajectoryPoint::layerTypeNames{{"ClosestApproach",
+                                                                                             "BeamPipeOrEndVertex",
+                                                                                             "PS1",
+                                                                                             "PS2",
+                                                                                             "ECALEntrance",
+                                                                                             "ECALShowerMax",
+                                                                                             "HCALEntrance",
+                                                                                             "HCALExit",
+                                                                                             "HOLayer",
+                                                                                             "VFcalEntrance"}};
+
 PFTrajectoryPoint::PFTrajectoryPoint() : isTrackerLayer_(false), detId_(-1), layer_(-1) {}
 
 PFTrajectoryPoint::PFTrajectoryPoint(int detId,
@@ -25,6 +39,15 @@ PFTrajectoryPoint::PFTrajectoryPoint(const PFTrajectoryPoint& other)
 
 PFTrajectoryPoint::~PFTrajectoryPoint() {}
 
+PFTrajectoryPoint::LayerType PFTrajectoryPoint::layerTypeByName(const std::string& name) {
+  auto it = std::find(layerTypeNames.begin(), layerTypeNames.end(), name);
+  if (it == layerTypeNames.end()) {
+    return Unknown;  // better this or throw()?
+  }
+  int index = std::distance(layerTypeNames.begin(), it);
+  return LayerType(index);
+}
+
 bool PFTrajectoryPoint::operator==(const reco::PFTrajectoryPoint& other) const {
   if (posxyz_ == other.posxyz_ && momentum_ == other.momentum_)
     return true;

HLTrigger/Configuration/python/customizeHLTforCMSSW.py

@@ -46,7 +46,7 @@ def synchronizeHCALHLTofflineRun3on2018data(process):
     # this function bring back the Run3 menu to a Run2-2018 like meny, for testing in data 2018
 
     #----------------------------------------------------------------------------------------------------------
-    # adapt threshold for HB  - in 2018 only one depth 
+    # adapt threshold for HB  - in 2018 only one depth
 
     for producer in producers_by_type(process, "PFClusterProducer"):
         if producer.seedFinder.thresholdsByDetector[0].detector.value() == 'HCAL_BARREL1':
@@ -164,10 +164,51 @@ def customiseFor2017DtUnpacking(process):
 
     return process
 
+# for PFBlockProducer/Algo change to enable both track-HCAL and track-HF links
+hltPFBlockLinkDefPrePhase2 = cms.VPSet(
+    cms.PSet(  linkType = cms.string( "PS1:ECAL" ),
+               useKDTree = cms.bool( True ),
+               linkerName = cms.string( "PreshowerAndECALLinker" )
+    ),
+    cms.PSet(  linkType = cms.string( "PS2:ECAL" ),
+               useKDTree = cms.bool( True ),
+               linkerName = cms.string( "PreshowerAndECALLinker" )
+    ),
+    cms.PSet(  linkType = cms.string( "TRACK:ECAL" ),
+               useKDTree = cms.bool( True ),
+               linkerName = cms.string( "TrackAndECALLinker" )
+    ),
+    cms.PSet(  linkType = cms.string( "TRACK:HCAL" ),
+               useKDTree = cms.bool( True ),
+               linkerName = cms.string( "TrackAndHCALLinker" ),
+               trajectoryLayerEntrance = cms.string("HCALEntrance"), # added
+               trajectoryLayerExit = cms.string("HCALExit") # added
+    ),
+    cms.PSet(  linkType = cms.string( "ECAL:HCAL" ),
+               useKDTree = cms.bool( False ),
+               linkerName = cms.string( "ECALAndHCALLinker" )
+    ),
+    cms.PSet(  linkType = cms.string( "HFEM:HFHAD" ),
+               useKDTree = cms.bool( False ),
+               linkerName = cms.string( "HFEMAndHFHADLinker" )
+    )
+)
+def customiseFor28442(process):
+    if hasattr(process,'hltParticleFlowBlock'):
+        process.hltParticleFlowBlock.linkDefinitions = hltPFBlockLinkDefPrePhase2
+    if hasattr(process,'hltParticleFlowBlockForTaus'):
+        process.hltParticleFlowBlockForTaus.linkDefinitions = hltPFBlockLinkDefPrePhase2
+    if hasattr(process,'hltParticleFlowBlockReg'):
+        process.hltParticleFlowBlockReg.linkDefinitions = hltPFBlockLinkDefPrePhase2
+    if hasattr(process,'hltParticleFlowBlockPPOnAA'):
+        process.hltParticleFlowBlockPPOnAA.linkDefinitions = hltPFBlockLinkDefPrePhase2
+    return process
+
 # CMSSW version specific customizations
 def customizeHLTforCMSSW(process, menuType="GRun"):
 
     # add call to action function in proper order: newest last!
     # process = customiseFor12718(process)
+    process = customiseFor28442(process)
 
     return process

RecoParticleFlow/PFProducer/interface/KDTreeLinkerBase.h

@@ -1,6 +1,7 @@
 #ifndef KDTreeLinkerBase_h
 #define KDTreeLinkerBase_h
 
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
@@ -17,6 +18,7 @@ using BlockElt2BlockEltMap = std::map<reco::PFBlockElement *, BlockEltSet>;
 
 class KDTreeLinkerBase {
 public:
+  KDTreeLinkerBase(const edm::ParameterSet &conf) {}
   virtual ~KDTreeLinkerBase() {}
 
   void setTargetType(const reco::PFBlockElement::Type &tgt) { _targetType = tgt; }
@@ -81,11 +83,14 @@ class KDTreeLinkerBase {
   // substracting (adding) 2Pi. This field define the threshold of this operation.
   float phiOffset_ = 0.25;
 
+  // rechit with fraction this value will be ignored in KDTreeLinker
+  const float cutOffFrac = 1E-4;
+
   // Debug boolean. Not used until now.
   bool debug_ = false;
 };
 
 #include "FWCore/PluginManager/interface/PluginFactory.h"
-typedef edmplugin::PluginFactory<KDTreeLinkerBase *()> KDTreeLinkerFactory;
+typedef edmplugin::PluginFactory<KDTreeLinkerBase *(const edm::ParameterSet &)> KDTreeLinkerFactory;
 
 #endif /* !KDTreeLinkerBase_h */

RecoParticleFlow/PFProducer/interface/PFAlgo.h

@@ -55,6 +55,9 @@ class PFAlgo {
   /// constructor
   PFAlgo(double nSigmaECAL,
          double nSigmaHCAL,
+         double nSigmaHFEM,
+         double nSigmaHFHAD,
+         std::vector<double> resolHF_square,
          PFEnergyCalibration& calibration,
          PFEnergyCalibrationHF& thepfEnergyCalibrationHF,
          const edm::ParameterSet& pset);
@@ -153,10 +156,12 @@ class PFAlgo {
                          reco::TrackRef& trackRef);
 
   //Looks for a HF-associated element in the block and produces a PFCandidate from it with HF_EM and/or HF_HAD calibrations
-  void createCandidateHF(const reco::PFBlock& block,
-                         const reco::PFBlockRef& blockref,
-                         const edm::OwnVector<reco::PFBlockElement>& elements,
-                         ElementIndices& inds);
+  void createCandidatesHF(const reco::PFBlock& block,
+                          reco::PFBlock::LinkData& linkData,
+                          const edm::OwnVector<reco::PFBlockElement>& elements,
+                          std::vector<bool>& active,
+                          const reco::PFBlockRef& blockref,
+                          ElementIndices& inds);
 
   void createCandidatesHCAL(const reco::PFBlock& block,
                             reco::PFBlock::LinkData& linkData,
@@ -213,6 +218,11 @@ class PFAlgo {
 
   double nSigmaHCAL(double clusterEnergy, double clusterEta) const;
 
+  double hfEnergyResolution(double clusterEnergy) const;
+
+  double nSigmaHFEM(double clusterEnergy) const;
+  double nSigmaHFHAD(double clusterEnergy) const;
+
   std::unique_ptr<reco::PFCandidateCollection> pfCandidates_;
   // the post-HF-cleaned candidates
   reco::PFCandidateCollection pfCleanedCandidates_;
@@ -237,6 +247,13 @@ class PFAlgo {
   /// number of sigma to judge energy excess in HCAL
   const double nSigmaHCAL_;
 
+  /// number of sigma to judge energy excess in HF
+  const double nSigmaHFEM_;
+  const double nSigmaHFHAD_;
+
+  // HF resolution
+  const std::vector<double> resolHF_square_;
+
   PFEnergyCalibration& calibration_;
   PFEnergyCalibrationHF& thepfEnergyCalibrationHF_;
 
@@ -312,6 +329,11 @@ class PFAlgo {
   bool useVertices_ = false;
 
   edm::Handle<reco::MuonCollection> muonHandle_;
+
+  // Named constants
+  const double nSigmaEConstHCAL = 100.;
+  const double nSigmaEConstHFEM = 100.;
+  const double nSigmaEConstHFHAD = 100.;
 };
 
 #endif

RecoParticleFlow/PFProducer/plugins/PFProducer.cc

@@ -33,7 +33,7 @@
 
 #include "TFile.h"
 
-/**\class PFProducer 
+/**\class PFProducer
 \brief Producer for particle flow reconstructed particles (PFCandidates)
 
 This producer makes use of PFAlgo, the particle flow algorithm.
@@ -120,6 +120,9 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig)
                              iConfig.getParameter<std::vector<double>>("calibHF_b_EMHAD")),
       pfAlgo_(iConfig.getParameter<double>("pf_nsigma_ECAL"),
               iConfig.getParameter<double>("pf_nsigma_HCAL"),
+              iConfig.getParameter<double>("pf_nsigma_HFEM"),
+              iConfig.getParameter<double>("pf_nsigma_HFHAD"),
+              iConfig.getParameter<std::vector<double>>("resolHF_square"),
               pfEnergyCalibration_,
               pfEnergyCalibrationHF_,
               iConfig) {
@@ -407,6 +410,8 @@ void PFProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
   // number of sigmas for neutral energy detection
   desc.add<double>("pf_nsigma_ECAL", 0.0);
   desc.add<double>("pf_nsigma_HCAL", 1.0);
+  desc.add<double>("pf_nsigma_HFEM", 1.0);
+  desc.add<double>("pf_nsigma_HFHAD", 1.0);
 
   // ECAL/HCAL PF cluster calibration : take it from global tag ?
   desc.add<bool>("useCalibrationsFromDB", true);
@@ -448,5 +453,9 @@ void PFProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
   desc.add<std::vector<double>>("calibHF_b_HADonly", {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.});
   desc.add<std::vector<double>>("calibHF_b_EMHAD", {1., 1., 1., 1., 1., 1., 1., 1., 1., 1.});
 
+  // resolution parameters for HF: EPJC 53(2008)139, doi:10.1140/epjc/s10052-007-0459-4
+  desc.add<std::vector<double>>("resolHF_square", {2.799 * 2.799, 0.114 * 0.114, 0.0 * 0.0})
+      ->setComment("HF resolution - stochastic, constant, noise term squares");
+
   descriptions.add("particleFlow", desc);
 }

RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc

@@ -9,7 +9,7 @@
 // It is used in PFBlockAlgo.cc in the function links().
 class KDTreeLinkerPSEcal : public KDTreeLinkerBase {
 public:
-  KDTreeLinkerPSEcal();
+  KDTreeLinkerPSEcal(const edm::ParameterSet &conf);
   ~KDTreeLinkerPSEcal() override;
 
   // With this method, we create the list of psCluster that we want to link.
@@ -70,8 +70,8 @@ class KDTreeLinkerPSEcal : public KDTreeLinkerBase {
 // construct it when calling the factory
 DEFINE_EDM_PLUGIN(KDTreeLinkerFactory, KDTreeLinkerPSEcal, "KDTreePreshowerAndECALLinker");
 
-KDTreeLinkerPSEcal::KDTreeLinkerPSEcal()
-    : KDTreeLinkerBase(), resPSpitch_(0.19), resPSlength_(6.1), ps1ToEcal_(1.072), ps2ToEcal_(1.057) {}
+KDTreeLinkerPSEcal::KDTreeLinkerPSEcal(const edm::ParameterSet &conf)
+    : KDTreeLinkerBase(conf), resPSpitch_(0.19), resPSlength_(6.1), ps1ToEcal_(1.072), ps2ToEcal_(1.057) {}
 
 KDTreeLinkerPSEcal::~KDTreeLinkerPSEcal() { clear(); }
 
@@ -100,7 +100,7 @@ void KDTreeLinkerPSEcal::insertFieldClusterElt(reco::PFBlockElement *ecalCluster
     const reco::PFRecHitRef &rh = fraction[rhit].recHitRef();
     double fract = fraction[rhit].fraction();
 
-    if ((rh.isNull()) || (fract < 1E-4))
+    if ((rh.isNull()) || (fract < cutOffFrac))
       continue;
 
     const reco::PFRecHit &rechit = *rh;

RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc

@@ -9,7 +9,7 @@
 // It is used in PFBlockAlgo.cc in the function links().
 class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
 public:
-  KDTreeLinkerTrackEcal();
+  KDTreeLinkerTrackEcal(const edm::ParameterSet &conf);
   ~KDTreeLinkerTrackEcal() override;
 
   // With this method, we create the list of psCluster that we want to link.
@@ -58,7 +58,7 @@ class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
 // construct it when calling the factory
 DEFINE_EDM_PLUGIN(KDTreeLinkerFactory, KDTreeLinkerTrackEcal, "KDTreeTrackAndECALLinker");
 
-KDTreeLinkerTrackEcal::KDTreeLinkerTrackEcal() : KDTreeLinkerBase() {}
+KDTreeLinkerTrackEcal::KDTreeLinkerTrackEcal(const edm::ParameterSet &conf) : KDTreeLinkerBase(conf) {}
 
 KDTreeLinkerTrackEcal::~KDTreeLinkerTrackEcal() { clear(); }
 
@@ -84,7 +84,7 @@ void KDTreeLinkerTrackEcal::insertFieldClusterElt(reco::PFBlockElement *ecalClus
     const reco::PFRecHitRef &rh = fraction[rhit].recHitRef();
     double fract = fraction[rhit].fraction();
 
-    if ((rh.isNull()) || (fract < 1E-4))
+    if ((rh.isNull()) || (fract < cutOffFrac))
       continue;
 
     const reco::PFRecHit &rechit = *rh;