CaloFuture/CaloFutureMoniDst/src/CaloFutureEMuPIDMon.cpp

@@ -53,9 +53,6 @@ private:
   Gaudi::Property<float> m_maxEHcalE{this, "maxEHcalE", 1.e10,
                                      "maximal ProtoParticle::additionalInfo::CaloHcalE for electron id [MeV/c]"};
 
-  Gaudi::Property<float> m_minPrsE{this, "minPrsE", 15.,
-                                   "minimal ProtoParticle::additionalInfo::CaloPrsE for electron id [MeV/c]"};
-
   Gaudi::Property<unsigned long int> m_nEventMin{this, "nEventMin", 200, "minimal number of events to check"};
 
   Gaudi::Property<bool> m_muonLoose{this, "useIsMuonLoose", true,
@@ -95,7 +92,6 @@ StatusCode CaloFutureEMuPIDMon::initialize() {
   bool old_m_splitSides = m_splitSides;
   m_splitSides          = false;
 
-  hBook1( "prsem", "E Prs", 0.5, 25.5, 25 );
   hBook1( "ecalem", "E Ecal", 40., 2040., 50 );
   hBook1( "hcalem", "E Hcal", 100., 5100., 50 );
 
@@ -124,18 +120,15 @@ void CaloFutureEMuPIDMon::operator()( const Input& particles ) const {
     if ( track->type() < LHCb::Track::Types::Long ) continue;
     if ( ( track->pt() < m_minPt || track->pt() > m_maxPt ) && !m_uncut ) continue;
 
-    const bool inprs  = ( proto->info( LHCb::ProtoParticle::additionalInfo::InAccPrs, double( false ) ) != 0 );
     const bool inecal = ( proto->info( LHCb::ProtoParticle::additionalInfo::InAccEcal, double( false ) ) != 0 );
     const bool inhcal = ( proto->info( LHCb::ProtoParticle::additionalInfo::InAccHcal, double( false ) ) != 0 );
 
-    const float prse  = (float)proto->info( LHCb::ProtoParticle::additionalInfo::CaloPrsE, -1 * Gaudi::Units::GeV );
     const float ecale = (float)proto->info( LHCb::ProtoParticle::additionalInfo::CaloEcalE, -1 * Gaudi::Units::GeV );
     const float hcale = (float)proto->info( LHCb::ProtoParticle::additionalInfo::CaloHcalE, -1 * Gaudi::Units::GeV );
 
     // -----------------------------------------------------------------
     // electron histograms
     // e.g.:  "pt>0.5&&inecal==1&&(hcale<=0||hcale>0&&inhcal==1&&hcale<1000.)&&rdlle>4"
-    //   or:  "pt>250&&inecal==1&&prse>15"
     // -----------------------------------------------------------------
 
     do {
@@ -143,7 +136,6 @@ void CaloFutureEMuPIDMon::operator()( const Input& particles ) const {
       float rdlle = (float)proto->info( LHCb::ProtoParticle::additionalInfo::RichDLLe, -9999. );
       if ( rdlle < m_RichDLLe && !m_uncut ) break;            // ----
       if ( inhcal && hcale > m_maxEHcalE && !m_uncut ) break; // ----
-      if ( prse < m_minPrsE && !m_uncut ) break;              // ----
 
       const auto            hypos      = proto->calo();
       const LHCb::CaloHypo* m_electron = NULL;
@@ -174,7 +166,6 @@ void CaloFutureEMuPIDMon::operator()( const Input& particles ) const {
                                   : false ) );
 
     if ( ismuon || m_uncut ) {
-      if ( inprs ) hFill1( "prsem", prse );
       if ( inecal ) hFill1( "ecalem", ecale );
       if ( inhcal ) hFill1( "hcalem", hcale );
     }

CaloFuture/CaloFutureMoniDst/src/CaloFutureProtoElectronMonitor.cpp

@@ -51,7 +51,6 @@ private:
   Gaudi::Property<float>                           m_eOpMin{this, "HistoEoPMin", 0.};
   Gaudi::Property<float>                           m_eOpMax{this, "HistoEoPMax", 3.};
   Gaudi::Property<int>                             m_eOpBin{this, "HistoEoPBin", 100};
-  Gaudi::Property<float>                           m_prsCut{this, "PrsCut", 50. * Gaudi::Units::MeV};
   Gaudi::Property<bool>                            m_pairing{this, "ElectronPairing", false};
   Gaudi::Property<std::vector<LHCb::Track::Types>> m_tracks{
       this, "TrackTypes", {LHCb::Track::Types::Long, LHCb::Track::Types::Downstream}};
@@ -152,8 +151,6 @@ void CaloFutureProtoElectronMonitor::operator()( const Input& protos, const Dete
     const double         et = momentum.pt();
     if ( e < m_eFilter ) continue;
     if ( et < m_etFilter ) continue;
-    const double ePrs = proto->info( LHCb::ProtoParticle::additionalInfo::CaloPrsE, 0. );
-    if ( m_prsCut > 0 && ePrs < m_prsCut ) continue;
 
     // Retrieve the seed id
     auto id = LHCb::Detector::Calo::CellID();
@@ -173,7 +170,6 @@ void CaloFutureProtoElectronMonitor::operator()( const Input& protos, const Dete
       hFill2( id, "5", position->x(), position->y(), e );
     }
     hFill1( id, "6", eOp );
-    hFill1( id, "7", ePrs );
 
     // perform electron pairing
     if ( !m_pairing ) continue;
@@ -200,8 +196,6 @@ void CaloFutureProtoElectronMonitor::operator()( const Input& protos, const Dete
       const double         et2 = momentum2.pt();
       if ( e2 < m_eFilter ) continue;
       if ( et2 < m_etFilter ) continue;
-      const double ePrs2 = proto2->info( LHCb::ProtoParticle::additionalInfo::CaloPrsE, 0. );
-      if ( m_prsCut > 0 && ePrs2 < m_prsCut ) continue;
 
       // Combine hypo
       LHCb::Calo::Momentum momentumSum( hypo );

CaloFuture/CaloFutureMoniDst/src/ChargedPIDsChecker.cpp

@@ -175,6 +175,7 @@ namespace LHCb::Calo::Algorithms {
       sc &= tuple->column( base + "BremHypoMatch", brems ? brems->BremHypoMatch() : -1.f );
       sc &= tuple->column( base + "BremHypoEnergy", brems ? brems->BremHypoEnergy() : -1.f );
       sc &= tuple->column( base + "BremHypoDeltaX", brems ? brems->BremHypoDeltaX() : -1000.f );
+      sc &= tuple->column( base + "BremBendingCorrection", brems ? brems->BremBendingCorrection() : 1.f );
       sc &= tuple->column( base + "BremTrackBasedEnergy", brems ? brems->BremTrackBasedEnergy() : -1.f );
       sc &= tuple->column( base + "BremEnergy", brems ? brems->BremEnergy() : 0.f );
       sc &= tuple->column( base + "BremPIDe", brems ? brems->BremPIDe() : 0.f );

CaloFuture/CaloFuturePIDs/src/CaloChargedBremAlg.cpp

@@ -41,11 +41,11 @@ namespace LHCb::Calo {
 
   private:
     // bremsstrahlung recovery options
-    Gaudi::Property<BremMethod> m_bremmethod{this, "BremMethod", BremMethod::TrackBased, "Brem recovery method"};
-    Gaudi::Property<float>      m_bremchi2_max{this, "BremChi2Max", 4.,
+    Gaudi::Property<BremMethod> m_bremmethod{this, "BremMethod", BremMethod::Mixed, "Brem recovery method"};
+    Gaudi::Property<float>      m_bremchi2_max{this, "BremChi2Max", 6.,
                                           "Maximum chi2 of cluster for track to use brem recovery (any 'BremMethod')"};
     Gaudi::Property<float>      m_bremchi2_thr{
-        this, "BremChi2Threshold", 1.5,
+        this, "BremChi2Threshold", 1.0,
         "Threshold of chi2 of cluster when to switch to track-based brem recovery for 'BremMethod::Mixed'"};
 
     // DLL parametrization locations
@@ -138,13 +138,15 @@ namespace LHCb::Calo {
       auto pid = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>();
 
       // cluster-based brem recovery
-      auto brem_info = getBestBrem<BremMatch, BremDeltaX>( rel_brems.scalar()[idx], hypoclusters );
+      auto brem_info = getBestBrem<BremMatch, BremDeltaX, BremBendingCorr>( rel_brems.scalar()[idx], hypoclusters );
       auto brem_chi2 = brem_info ? brem_info->match.chi2 : -1.f;
       auto brem_dx   = brem_info ? brem_info->dx : 0.f;
+      auto brem_bc   = brem_info ? brem_info->bc : 1.f;
       pid.field<ChargedInfoTag::BremHypoIndex>().set( brem_info ? brem_info->match.index : -1 );
       pid.field<ChargedInfoTag::BremHypoMatch>().set( brem_chi2 );
       pid.field<ChargedInfoTag::BremHypoEnergy>().set( brem_info ? brem_info->energy : -1.f );
       pid.field<ChargedInfoTag::BremHypoDeltaX>().set( brem_dx );
+      pid.field<ChargedInfoTag::BremBendingCorrection>().set( brem_bc );
 
       // track-based brem recovery
       auto brem_e_trackbased = getRelationValue<BremEnergy, float>( rel_breme.scalar()[idx] );

CaloFuture/CaloFuturePIDs/src/CaloChargedInfoUtils.h

@@ -58,6 +58,7 @@ namespace LHCb::Calo {
         Match match;
         float energy;
         float dx;
+        float bc;
       };
 
       // helpers for DLL getter
@@ -100,7 +101,8 @@ namespace LHCb::Calo {
     }
 
     // above one, but for brem with additional info
-    template <typename Chi2Type, typename DeltaXType, typename ProxyType, typename ClusterView>
+    template <typename Chi2Type, typename DeltaXType, typename BendingCorrectionType, typename ProxyType,
+              typename ClusterView>
     inline std::optional<EnergyDeltaXMatch> getBestBrem( ProxyType const& proxy, ClusterView const& clusters ) {
       auto match = getMinChiSqMatch<ProxyType, Chi2Type>( proxy );
       if ( !( clusters && match.has_value() ) ) return std::nullopt;
@@ -109,20 +111,21 @@ namespace LHCb::Calo {
       auto const cluster  = clusters.find( idx );
       if ( cluster == clusters.end() ) return std::nullopt;
       match->index = idx;
-      return EnergyDeltaXMatch{*match, cluster->energy(), relation.template get<DeltaXType>().cast()};
+      return EnergyDeltaXMatch{*match, cluster->energy(), relation.template get<DeltaXType>().cast(),
+                               relation.template get<BendingCorrectionType>().cast()};
     }
 
     // to get brem energy, optimized between track-based and cluster method
     inline std::optional<float> getBremEnergy( BremMethod method, std::optional<EnergyDeltaXMatch> brem,
-                                               std::optional<float> brem_trackbased, float threshold_chi2 = 1.5,
-                                               float max_chi2 = 4 ) {
+                                               std::optional<float> brem_trackbased, float threshold_chi2 = 1.,
+                                               float max_chi2 = 6. ) {
       bool add_brem = brem.has_value() && ( brem->match.chi2 <= max_chi2 );
       switch ( method ) {
       case BremMethod::TrackBased:
-        return ( add_brem && brem_trackbased ) ? std::optional<float>( brem_trackbased ) : std::nullopt;
+        return ( add_brem && brem_trackbased ) ? std::optional<float>( *brem_trackbased ) : std::nullopt;
       case BremMethod::Mixed:
         return ( add_brem && brem_trackbased )
-                   ? std::optional<float>( ( brem->match.chi2 < threshold_chi2 ) ? brem->energy : brem_trackbased )
+                   ? std::optional<float>( brem->match.chi2 < threshold_chi2 ? brem->energy : *brem_trackbased )
                    : std::nullopt;
       case BremMethod::Cluster:
         return add_brem ? std::optional<float>( brem->energy ) : std::nullopt;

CaloFuture/CaloFuturePIDs/src/SelectiveBremMatchAlg.cpp

@@ -204,8 +204,9 @@ namespace LHCb::Calo {
     entable.reserve( tracksincalo.size() );
 
     // track state from where to extrapolate (linearly) to calo from
-    auto const state_locs = extrapolation_statelocs_brem( *tracksincalo.from() );
-    if ( !state_locs.has_value() ) {
+    auto const state_locs_brem    = extrapolation_statelocs_brem( *tracksincalo.from() );
+    auto const state_loc_electron = extrapolation_stateloc( *tracksincalo.from() );
+    if ( !state_locs_brem.has_value() ) {
       throw GaudiException( "Not a valid track type for brem recovery.", "LHCb::Event::Enum::Track::Type",
                             StatusCode::FAILURE );
     }
@@ -246,10 +247,18 @@ namespace LHCb::Calo {
       // propagate the two 'edge' states (first and last before magnet) to calo
       // NB: only using last-state covariance (also only calculated)
       auto track     = trackincalo.from();
-      auto ref_first = track.state( state_locs.value().first );
-      auto ref_last  = track.state( state_locs.value().second );
+      auto ref_first = track.state( state_locs_brem.value().first );
+      auto ref_last  = track.state( state_locs_brem.value().second );
       if ( !propagateToCaloForBrem<MyState>( state_first, state_last, ref_first, ref_last, plane_at_calo ) ) continue;
 
+      // information for bending correction
+      std::optional<float> full_dx = std::nullopt;
+      if ( state_loc_electron.has_value() ) {
+        auto ref_ele = track.state( state_loc_electron.value() );
+        full_dx      = {std::abs( state_first.x() -
+                             ( ref_ele.x().cast() + ref_ele.tx().cast() * ( state_first.z() - ref_ele.z().cast() ) ) )};
+      }
+
       // look for CellIDs from first to last state
       auto       dx         = state_last.x() - state_first.x();
       const auto cell_first = getClosestCellID( calo, state_first.position(), state_first.slopes(), other_planes );
@@ -327,10 +336,13 @@ namespace LHCb::Calo {
           continue;
         }
 
+        // bending correction
+        float bending_correction = full_dx.has_value() ? 1.f - std::abs( dx ) * rdx / full_dx.value() : 1.f;
+
         // only now push proper results
         auto cellid_v = simd_t::int_v{(int)LHCb::Detector::Calo::Index{caloobj->cellID()}};
         auto relation = trtable.emplace_back<SIMDWrapper::InstructionSet::Scalar>();
-        relation.set( track.indices(), cellid_v, chi2, teststat_dx );
+        relation.set( track.indices(), cellid_v, chi2, teststat_dx, bending_correction );
       }
     }
 

Phys/DaVinciInterfaces/dict/DaVinciInterfacesDict.h

@@ -23,7 +23,6 @@
 // DaVinciKernel
 // ============================================================================
 #include "Kernel/IBTaggingTool.h"
-#include "Kernel/IBremAdder.h"
 #include "Kernel/ICaloParticleMaker.h"
 #include "Kernel/IChangePIDTool.h"
 #include "Kernel/ICheckOverlap.h"

Phys/DaVinciInterfaces/dict/DaVinciInterfacesDict.xml

@@ -12,7 +12,6 @@
 
 <!-- DaVinciInterfaces -->
 
-  <class name = "IBremAdder"/>
   <class name = "IBTaggingTool"/>
   <class name = "ICaloParticleMaker"/>
   <class name = "IChangePIDTool"/>

Phys/DaVinciInterfaces/include/Kernel/IFunctionalBremAdder.h

@@ -20,7 +20,7 @@
 #include "DetDesc/IGeometryInfo.h"
 
 /**
- * Interface for BremStrahlung correction to electron particle
+ * Interface for Bremsstrahlung correction to electron particle
  *
  *  @author Olivier Deschamps
  *  @date   2006-10-25
@@ -37,14 +37,5 @@ struct GAUDI_API IFunctionalBremAdder : extend_interfaces<IAlgTool> {
                              IGeometryInfo const& geometry, bool force = false ) const = 0;
 
   // helper methods
-  virtual bool               hasBrem( const LHCb::Particle& particle ) const     = 0;
-  virtual LHCb::CaloMomentum bremMomentum( const LHCb::Particle& particle, const LHCb::Particle::Range& photons,
-                                           IGeometryInfo const& geometry ) const = 0;
-  virtual std::pair<unsigned int, LHCb::CaloMomentum>       bremNbAndMomentum( const LHCb::Particle&        particle,
-                                                                               const LHCb::Particle::Range& photons,
-                                                                               IGeometryInfo const& geometry ) const = 0;
-  virtual std::pair<LHCb::CaloMomentum, LHCb::CaloMomentum> bremMomenta( const LHCb::Particle&        p1,
-                                                                         const LHCb::Particle&        p2,
-                                                                         const LHCb::Particle::Range& photons,
-                                                                         IGeometryInfo const& geometry ) const = 0;
+  virtual bool hasBrem( const LHCb::Particle& particle ) const = 0;
 };

Phys/DaVinciInterfaces/include/Kernel/IBremAdder.h → Phys/DaVinciInterfaces/include/Kernel/ISelectiveBremAdder.h

@@ -10,36 +10,23 @@
 \*****************************************************************************/
 #pragma once
 
-#include "CaloUtils/CaloMomentum.h"
-#include "Event/CaloHypo.h"
+#include "DetDesc/IGeometryInfo.h"
 #include "Event/Particle.h"
-
 #include "GaudiKernel/IAlgTool.h"
-
 #include <string>
 
-#include "DetDesc/IGeometryInfo.h"
-
 /**
- * Interface for BremStrahlung correction to electron particle
- *
- *  @author Olivier Deschamps
- *  @date   2006-10-25
+ * Interface for bremsstrahlung correction to electron particle
  */
-struct GAUDI_API IBremAdder : extend_interfaces<IAlgTool> {
+struct GAUDI_API ISelectiveBremAdder : extend_interfaces<IAlgTool> {
 
-  DeclareInterfaceID( IBremAdder, 5, 0 );
+  DeclareInterfaceID( ISelectiveBremAdder, 4, 0 );
 
-  virtual bool addBrem( LHCb::Particle* particle, IGeometryInfo const& geometry, bool force = false )    = 0;
-  virtual bool removeBrem( LHCb::Particle* particle, IGeometryInfo const& geometry, bool force = false ) = 0;
+  // Add Brem
+  virtual bool addBrem( LHCb::Particle& particle, bool force = false ) const = 0;
   // Add Brem on particle pair (removing overlap)
-  virtual bool addBrem2Pair( LHCb::Particle* p1, LHCb::Particle* p2, IGeometryInfo const& geometry,
-                             bool force = false ) = 0;
+  virtual bool addBrem2Pair( LHCb::Particle& p1, LHCb::Particle& p2, bool force = false ) const = 0;
 
-  virtual bool                     hasBrem( LHCb::Particle const* particle )                                     = 0;
-  virtual const LHCb::CaloMomentum bremMomentum( LHCb::Particle const* particle, IGeometryInfo const& geometry ) = 0;
-  virtual const std::pair<unsigned int, LHCb::CaloMomentum> bremNbAndMomentum( const LHCb::Particle* particle,
-                                                                               IGeometryInfo const&  geometry )   = 0;
-  virtual const std::pair<LHCb::CaloMomentum, LHCb::CaloMomentum>
-  bremMomenta( LHCb::Particle const* p1, LHCb::Particle const* p2, IGeometryInfo const& geometry ) = 0;
+  // helper methods
+  virtual bool hasBrem( const LHCb::Particle& particle ) const = 0;
 };

Phys/DaVinciInterfaces/src/DaVinciInterfaces.cpp

@@ -17,7 +17,6 @@
 // ============================================================================
 // local
 // ============================================================================
-#include "Kernel/IBremAdder.h"
 #include "Kernel/ICheckOverlap.h"
 #include "Kernel/IDVAlgorithm.h"
 #include "Kernel/IDecayTreeFit.h"

Phys/DaVinciNeutralTools/CMakeLists.txt

@@ -15,8 +15,8 @@ Phys/DaVinciNeutralTools
 
 gaudi_add_module(DaVinciNeutralTools
     SOURCES
-        src/BremAdder.cpp
         src/FunctionalBremAdder.cpp
+        src/SelectiveBremAdder.cpp
     LINK
         Gaudi::GaudiAlgLib
         Gaudi::GaudiKernel

Phys/DaVinciNeutralTools/src/FunctionalBremAdder.cpp

@@ -9,8 +9,6 @@
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
 
-#include "CaloUtils/CaloMomentum.h"
-#include "CaloUtils/CaloParticle.h"
 #include "Event/StateParameters.h"
 
 #include "CaloDet/DeCalorimeter.h"
@@ -41,14 +39,6 @@ public:
   bool addBrem2Pair( LHCb::Particle& p1, LHCb::Particle& p2, const LHCb::Particle::Range& photons,
                      IGeometryInfo const& geometry, bool force = false ) const override;
   bool hasBrem( const LHCb::Particle& particle ) const override;
-  LHCb::CaloMomentum bremMomentum( const LHCb::Particle& particle, const LHCb::Particle::Range& photons,
-                                   IGeometryInfo const& geometry ) const override;
-  std::pair<unsigned int, LHCb::CaloMomentum>       bremNbAndMomentum( const LHCb::Particle&        particle,
-                                                                       const LHCb::Particle::Range& photons,
-                                                                       IGeometryInfo const&         geometry ) const override;
-  std::pair<LHCb::CaloMomentum, LHCb::CaloMomentum> bremMomenta( const LHCb::Particle& p1, const LHCb::Particle& p2,
-                                                                 const LHCb::Particle::Range& photons,
-                                                                 IGeometryInfo const&         geometry ) const override;
 
 protected:
   bool brem4particle( LHCb::Particle& particle, const std::vector<const LHCb::CaloHypo*>& brems,
@@ -156,35 +146,6 @@ FunctionalBremAdder::bremLists( const LHCb::Particle& p1, const LHCb::Particle&
   return {brems1, bb2};
 }
 
-//---- get the brem as CaloMomentum
-std::pair<unsigned int, LHCb::CaloMomentum>
-FunctionalBremAdder::bremNbAndMomentum( const LHCb::Particle& particle, const LHCb::Particle::Range& photons,
-                                        IGeometryInfo const& geometry ) const {
-  const std::vector<const LHCb::CaloHypo*>& brems = getBrem( particle, photons, geometry );
-  LHCb::CaloMomentum                        m     = LHCb::CaloMomentum( brems );
-  m.setReferencePoint( firstOrigin( brems, particle, geometry ) );
-  return {brems.size(), m};
-}
-
-LHCb::CaloMomentum FunctionalBremAdder::bremMomentum( const LHCb::Particle&        particle,
-                                                      const LHCb::Particle::Range& photons,
-                                                      IGeometryInfo const&         geometry ) const {
-  return bremNbAndMomentum( particle, photons, geometry ).second;
-}
-
-//--- get the brem as CaloMomentum pair
-std::pair<LHCb::CaloMomentum, LHCb::CaloMomentum>
-FunctionalBremAdder::bremMomenta( const LHCb::Particle& p1, const LHCb::Particle& p2,
-                                  const LHCb::Particle::Range& photons, IGeometryInfo const& geometry ) const {
-
-  const auto& [b1, b2]  = bremLists( p1, p2, photons, geometry );
-  LHCb::CaloMomentum m1 = LHCb::CaloMomentum( b1 );
-  LHCb::CaloMomentum m2 = LHCb::CaloMomentum( b2 );
-  m1.setReferencePoint( firstOrigin( b1, p1, geometry ) );
-  m2.setReferencePoint( firstOrigin( b2, p2, geometry ) );
-  return {m1, m2};
-}
-
 //============================================================================= PRIVATE METHODS
 bool FunctionalBremAdder::brem4particle( LHCb::Particle& particle, const std::vector<const LHCb::CaloHypo*>& brems,
                                          IGeometryInfo const& geometry, bool force ) const {

Phys/DaVinciNeutralTools/src/SelectiveBremAdder.cpp

+/*****************************************************************************\
+* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the GNU General Public      *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+*                                                                             *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization  *
+* or submit itself to any jurisdiction.                                       *
+\*****************************************************************************/
+
+#include "Gaudi/Accumulators.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/Point3DTypes.h"
+#include "GaudiKernel/SymmetricMatrixTypes.h"
+#include "Kernel/IParticle2State.h"
+#include "Kernel/ISelectiveBremAdder.h"
+#include <functional>
+
+/**
+ *  Adds brem corrections to particles using energy
+ *  (should be output of SelectiveBremMatchAlg)
+ *  stored in ProtoParticle
+ */
+
+// helper functions/classes
+namespace {
+
+  // brem info (needed for overlap check)
+  struct BremInfo {
+    int    inbrem;
+    int    index;
+    double chi2;
+    double energy;
+    double bc;
+    // for cluster overlap check
+    friend bool operator==( BremInfo const& lhs, BremInfo const& rhs ) {
+      return ( lhs.index == rhs.index ) && ( lhs.energy > 0. ) && ( rhs.energy > 0. );
+    }
+  };
+
+  // 4-momentum covariance
+  // -- this assumes one scales the 3-momentum, which itself has a error associated to it
+  // -- brem correction should not affect the direction (uncertainties)
+  // -- ignoring mass terms in covariance calculation
+  Gaudi::SymMatrix4x4 getMomCov( Gaudi::LorentzVector const& mom0, Gaudi::SymMatrix4x4 const& cov0, double frac,
+                                 double cov_f ) {
+    Gaudi::SymMatrix4x4 ppt;
+    ppt( 0, 0 ) = mom0.Px() * mom0.Px();
+    ppt( 0, 1 ) = mom0.Px() * mom0.Py();
+    ppt( 0, 2 ) = mom0.Px() * mom0.Pz();
+    ppt( 0, 3 ) = mom0.Px() * mom0.E();
+    ppt( 1, 1 ) = mom0.Py() * mom0.Py();
+    ppt( 1, 2 ) = mom0.Py() * mom0.Pz();
+    ppt( 1, 3 ) = mom0.Py() * mom0.E();
+    ppt( 2, 2 ) = mom0.Pz() * mom0.Pz();
+    ppt( 2, 3 ) = mom0.Pz() * mom0.E();
+    ppt( 3, 3 ) = mom0.E() * mom0.E();
+    return cov_f * ppt + std::pow( frac, 2 ) * cov0;
+  }
+
+  // main covariance updater
+  void updateCovariances( LHCb::Particle& particle, Gaudi::LorentzVector const& mom0, double frac,
+                          double relative_momentum_err ) {
+    // adapt with relative resolution estimate (if electron, ignore mass term)
+    auto momCov = getMomCov( mom0, particle.momCovMatrix(), frac, relative_momentum_err );
+    particle.setMomCovMatrix( momCov );
+    // scale position-momentum correlation accordingly
+    particle.setPosMomCovMatrix( frac * particle.posMomCovMatrix() );
+  }
+
+} // namespace
+
+class SelectiveBremAdder : public extends<GaudiTool, ISelectiveBremAdder> {
+public:
+  using extends::extends;
+
+  bool addBrem( LHCb::Particle& particle, bool force = false ) const override;
+  bool addBrem2Pair( LHCb::Particle& p1, LHCb::Particle& p2, bool force = false ) const override;
+  bool hasBrem( const LHCb::Particle& particle ) const override;
+
+protected:
+  bool brem4particle( LHCb::Particle& particle, std::optional<BremInfo> const& breminfo, bool force = false ) const;
+
+private:
+  // properties
+  Gaudi::Property<int>    m_overlap_method{this, "OverlapAssignmentMethod", 0,
+                                        "For di-electron (pair) in case of overlap, determine assignment: 0 is to "
+                                        "one with lowest matching chi2; 1 is with highest (pre-brem) momentum track"};
+  Gaudi::Property<bool>   m_use_bending_corr{this, "UseBendingCorrection", true};
+  Gaudi::Property<bool>   m_update_covs_brem{this, "UpdateCovariancesBrem", true};
+  Gaudi::Property<bool>   m_update_covs_nobrem{this, "UpdateCovariancesNoBrem", true};
+  double                  m_relative_mom_cov_inbrem;
+  Gaudi::Property<double> m_relative_mom_inbrem_err{
+      this,
+      "RelativeMomentumScalingErrorInBrem",
+      0.25,
+      [this]( auto& ) { m_relative_mom_cov_inbrem = std::pow( m_relative_mom_inbrem_err, 2 ); },
+      Gaudi::Details::Property::ImmediatelyInvokeHandler{true},
+      "Update convariances with this relative momentum scaling error (in brem acceptance)"};
+  double                  m_relative_mom_cov_outbrem;
+  Gaudi::Property<double> m_relative_mom_outbrem_err{
+      this,
+      "RelativeMomentumScalingErrorNotInBrem",
+      0.35,
+      [this]( auto& ) { m_relative_mom_cov_outbrem = std::pow( m_relative_mom_outbrem_err, 2 ); },
+      Gaudi::Details::Property::ImmediatelyInvokeHandler{true},
+      "Update convariances with this relative momentum scaling error (not in brem acceptance)"};
+
+  // other members
+  ToolHandle<IParticle2State> m_p2s{this, "Particle2State", "Particle2State"};
+
+  // functions
+  std::optional<BremInfo> getBremInfo( LHCb::Particle const& proto ) const;
+
+  // Counters
+  mutable Gaudi::Accumulators::Counter<>     m_revBeamCorrCounter{this, "Revert brem correction"};
+  mutable Gaudi::Accumulators::StatCounter<> m_nAddPhotons{this, "Nb photons added to single electrons"};
+  mutable Gaudi::Accumulators::StatCounter<> m_nAddPairPhotons{this, "Nb photons added to electron pair"};
+  mutable Gaudi::Accumulators::StatCounter<> m_nAddPairSharedPhotons{this, "Nb photons shared by electron pair"};
+  mutable Gaudi::Accumulators::StatCounter<> m_DeltaE{this, "Delta(E)"};
+};
+
+//============================================================================= INTERFACED METHODS
+
+// ---- Add brem to particle
+bool SelectiveBremAdder::addBrem( LHCb::Particle& particle, bool force ) const {
+  auto brem    = getBremInfo( particle );
+  auto success = brem4particle( particle, brem, force );
+  m_nAddPhotons += hasBrem( particle );
+  return success;
+}
+
+//---- Add brem(s) to a pair of particles removing overlaps
+bool SelectiveBremAdder::addBrem2Pair( LHCb::Particle& p1, LHCb::Particle& p2, bool force ) const {
+  auto b1 = getBremInfo( p1 );
+  auto b2 = getBremInfo( p2 );
+  if ( b1.has_value() && b1 == b2 ) {
+    m_nAddPairPhotons += 1;
+    m_nAddPairSharedPhotons += 1;
+    switch ( m_overlap_method ) {
+    case ( 0 ):
+      ( b1->chi2 < b2->chi2 ? b2 : b1 )->energy = 0.;
+      break;
+    case ( 1 ):
+      ( p1.p() > p2.p() ? b2 : b1 )->energy = 0.;
+      break;
+    }
+    return brem4particle( p1, b1, force ) && brem4particle( p2, b2, force );
+  } else {
+    auto success = brem4particle( p1, b1, force ) && brem4particle( p2, b2, force );
+    m_nAddPairPhotons += hasBrem( p1 );
+    m_nAddPairPhotons += hasBrem( p2 );
+    return success;
+  }
+}
+
+//---- indicate whether the particle has a brem added or not
+bool SelectiveBremAdder::hasBrem( const LHCb::Particle& particle ) const {
+  return ( particle.info( LHCb::Particle::additionalInfo::HasBremAdded, 0. ) > 0. );
+}
+
+//============================================================================= PRIVATE METHODS
+bool SelectiveBremAdder::brem4particle( LHCb::Particle& particle, std::optional<BremInfo> const& breminfo,
+                                        bool force ) const {
+  // boolean 'force' : force re-adding brem
+
+  //-- consistency checks
+  const auto* proto = particle.proto();
+  if ( !proto ) return false;
+  if ( !proto->track() ) return false;
+
+  //-- check if bremsstrahlung has already been added
+  if ( hasBrem( particle ) && !force ) {
+    if ( msgLevel( MSG::DEBUG ) ) debug() << " The bremsstrahlung has already been added - nothing to add" << endmsg;
+    return false;
+  }
+
+  // revert the particle momentum  (as in Charged CombinedParticle)
+  if ( force ) {
+    ++m_revBeamCorrCounter;
+    auto state = &proto->track()->firstState();
+    m_p2s->state2Particle( *state, particle ).ignore();
+  }
+
+  // ----- start bremsstrahlung recovery -----
+  particle.eraseInfo( LHCb::Particle::additionalInfo::HasBremAdded );
+  particle.addInfo( LHCb::Particle::additionalInfo::HasBremAdded,
+                    breminfo.has_value() ? ( ( breminfo->energy > 0. ) ? 1. : 0. ) : 0. );
+
+  // only continue if one is in brem acceptance
+  const Gaudi::LorentzVector mom = particle.momentum();
+  if ( !breminfo.has_value() ) {
+    if ( m_update_covs_nobrem ) updateCovariances( particle, mom, 1., m_relative_mom_cov_outbrem );
+    return true;
+  }
+
+  // only continue if one has brem
+  if ( !hasBrem( particle ) ) {
+    if ( m_update_covs_nobrem ) updateCovariances( particle, mom, 1., m_relative_mom_cov_inbrem );
+    return true;
+  }
+
+  //----- add brem momentum
+  // keep momentum direction and particle mass
+  auto newenergy = mom.E() + ( m_use_bending_corr ? breminfo->bc * breminfo->energy : breminfo->energy );
+  auto momfrac   = std::sqrt( newenergy * newenergy - mom.M2() ) / mom.P();
+  auto newmom    = Gaudi::LorentzVector( momfrac * mom.Px(), momfrac * mom.Py(), momfrac * mom.Pz(), newenergy );
+  particle.setMomentum( newmom );
+  m_DeltaE += newenergy - mom.E();
+
+  // update covariance matrices
+  if ( m_update_covs_brem ) updateCovariances( particle, mom, momfrac, m_relative_mom_cov_inbrem );
+
+  return true;
+}
+
+std::optional<BremInfo> SelectiveBremAdder::getBremInfo( LHCb::Particle const& particle ) const {
+  //-- consistency checks
+  const auto* proto = particle.proto();
+  if ( !proto ) return std::nullopt;
+
+  //-- get info from ProtoParticle
+  auto inbrem = (int)proto->info( LHCb::ProtoParticle::additionalInfo::InAccBrem, 0. );
+  if ( !inbrem ) return std::nullopt;
+
+  auto energy = proto->info( LHCb::ProtoParticle::additionalInfo::CaloBremEnergy, 0. );
+  auto index  = (int)proto->info( LHCb::ProtoParticle::additionalInfo::CaloBremHypoIndex, -1. );
+  auto chi2   = proto->info( LHCb::ProtoParticle::additionalInfo::CaloBremMatch, -1. );
+  auto bc     = proto->info( LHCb::ProtoParticle::additionalInfo::CaloBremBendingCorr, 1. );
+  return std::optional<BremInfo>{{inbrem, index, chi2, energy, bc}};
+}
+
+// Declaration of the Tool Factory
+DECLARE_COMPONENT( SelectiveBremAdder )

Phys/DecayTreeFitter/src/lib/RecoTrack.cpp

@@ -64,9 +64,6 @@ namespace DecayTreeFitter {
           particle().momCovMatrix()( 3, 3 ) * std::pow( particle().momentum().E() / momentumFromParticle, 2 );
       m_bremEnergyCov = momentumError2FromParticle - momentumError2FromTrack;
       m_bremEnergy    = momentumFromParticle - momentumFromTrack;
-      // if the correction is too small or unphysical, ignore it. the units are MeV.
-      if ( !( m_bremEnergyCov > 1 * Gaudi::Units::MeV * Gaudi::Units::MeV && m_bremEnergy > 1 * Gaudi::Units::MeV ) )
-        m_bremEnergy = m_bremEnergyCov = 0;
       if ( vtxverbose >= 1 ) {
         std::cout << "Fitting an electron. (" << momentumFromTrack << " +/- " << std::sqrt( momentumError2FromTrack )
                   << ") --> (" << momentumFromParticle << " +/- " << std::sqrt( momentumError2FromParticle ) << ")"

Phys/ExtraInfoTools/src/ConeVariablesForEW.cpp

@@ -313,21 +313,18 @@ StatusCode ConeVariablesForEW::ChargedCone( const LHCb::Particle* seed, const LH
           if ( rcut == 0. ) {
 
             // Extra Electron
-            double prsE    = 50.;
             double eCalEoP = .10;
             double hCalEoP = .05;
-            if ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloPrsE, -1. ) > prsE ) {
-              if ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloEcalE, -1. ) / track->p() > eCalEoP ) {
-                if ( ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloHcalE, -1. ) > 0 ) &&
-                     ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloHcalE, -1. ) / track->p() < hCalEoP ) ) {
-                  if ( track->pt() < minPtE ) {
-                    minPtE = track->pt();
-                    minQE  = track->charge();
-                  }
-                  if ( track->pt() > maxPtE ) {
-                    maxPtE = track->pt();
-                    maxQE  = track->charge();
-                  }
+            if ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloEcalE, -1. ) / track->p() > eCalEoP ) {
+              if ( ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloHcalE, -1. ) > 0 ) &&
+                   ( proto->info( LHCb::ProtoParticle::additionalInfo::CaloHcalE, -1. ) / track->p() < hCalEoP ) ) {
+                if ( track->pt() < minPtE ) {
+                  minPtE = track->pt();
+                  minQE  = track->charge();
+                }
+                if ( track->pt() > maxPtE ) {
+                  maxPtE = track->pt();
+                  maxQE  = track->charge();
                 }
               }
             }

Phys/LoKiProtoParticles/include/LoKi/ProtoParticles.h

@@ -1386,18 +1386,9 @@ namespace LoKi {
     /// @var PP_CaloNeutralSpd
     const LoKi::ProtoParticles::Info PP_CaloNeutralSpd( LHCb::ProtoParticle::CaloNeutralSpd, -1 );
     // ========================================================================
-    /// @var PP_CaloNeutralPrs
-    const LoKi::ProtoParticles::Info PP_CaloNeutralPrs( LHCb::ProtoParticle::CaloNeutralPrs, -1.e+10 );
-    // ========================================================================
     /// @var PP_CaloNeutralEcal
     const LoKi::ProtoParticles::Info PP_CaloNeutralEcal( LHCb::ProtoParticle::CaloNeutralEcal, -1.e+10 );
     // ========================================================================
-    /// @var PP_CaloSpdE
-    const LoKi::ProtoParticles::Info PP_CaloSpdE( LHCb::ProtoParticle::CaloSpdE, -1.e+10 );
-    // ========================================================================
-    /// @var PP_CaloPrsE
-    const LoKi::ProtoParticles::Info PP_CaloPrsE( LHCb::ProtoParticle::CaloPrsE, -1.e+10 );
-    // ========================================================================
     /// @var PP_CaloEcalE
     const LoKi::ProtoParticles::Info PP_CaloEcalE( LHCb::ProtoParticle::CaloEcalE, -1.e+10 );
     // ========================================================================

Phys/LoKiProtoParticles/python/LoKiProtoParticles/functions.py

@@ -285,15 +285,9 @@ PP_ClusterMass = PP_INFO(LHCb.ProtoParticle.ClusterMass, -1.e+10)
 
 ## @see LoKi::Cuts::PP_CaloNeutralSpd
 PP_CaloNeutralSpd = PP_INFO(LHCb.ProtoParticle.CaloNeutralSpd, -1)
-## @see LoKi::Cuts::PP_CaloNeutralPrs
-PP_CaloNeutralPrs = PP_INFO(LHCb.ProtoParticle.CaloNeutralPrs, -1.e+10)
 ## @see LoKi::Cuts::PP_CaloNeutralEcal
 PP_CaloNeutralEcal = PP_INFO(LHCb.ProtoParticle.CaloNeutralEcal, -1.e+10)
 
-## @see LoKi::Cuts::PP_CaloSpdE
-PP_CaloSpdE = PP_INFO(LHCb.ProtoParticle.CaloSpdE, -1.e+10)
-## @see LoKi::Cuts::PP_CaloPesE
-PP_CaloPrsE = PP_INFO(LHCb.ProtoParticle.CaloPrsE, -1.e+10)
 ## @see LoKi::Cuts::PP_CaloEcalE
 PP_CaloEcalE = PP_INFO(LHCb.ProtoParticle.CaloEcalE, -1.e+10)
 ## @see LoKi::Cuts::PP_CaloHcalE

Phys/LoKiProtoParticles/src/ProtoParticles.cpp

@@ -688,16 +688,15 @@ bool LoKi::ProtoParticles::HasDetector::hasCaloDLL( const LHCb::ProtoParticle* p
 }
 // ============================================================================
 bool LoKi::ProtoParticles::HasDetector::hasSpdInfo( const LHCb::ProtoParticle* p ) const {
-  return p && ( //
-                  p->hasInfo( LHCb::ProtoParticle::CaloChargedSpd ) ||
-                  p->hasInfo( LHCb::ProtoParticle::CaloNeutralSpd ) || p->hasInfo( LHCb::ProtoParticle::CaloSpdE ) ||
-                  p->hasInfo( LHCb::ProtoParticle::CaloDepositID ) || p->hasInfo( LHCb::ProtoParticle::PhotonID ) );
+  return p &&
+         ( //
+             p->hasInfo( LHCb::ProtoParticle::CaloChargedSpd ) || p->hasInfo( LHCb::ProtoParticle::CaloNeutralSpd ) ||
+             p->hasInfo( LHCb::ProtoParticle::CaloDepositID ) || p->hasInfo( LHCb::ProtoParticle::PhotonID ) );
 }
 // ============================================================================
 bool LoKi::ProtoParticles::HasDetector::hasPrsInfo( const LHCb::ProtoParticle* p ) const {
   return p && ( //
                   p->hasInfo( LHCb::ProtoParticle::PrsPIDe ) || p->hasInfo( LHCb::ProtoParticle::CaloChargedPrs ) ||
-                  p->hasInfo( LHCb::ProtoParticle::CaloNeutralPrs ) || p->hasInfo( LHCb::ProtoParticle::CaloPrsE ) ||
                   p->hasInfo( LHCb::ProtoParticle::CaloDepositID ) || p->hasInfo( LHCb::ProtoParticle::PhotonID ) );
 }
 // ============================================================================
@@ -755,10 +754,8 @@ bool LoKi::ProtoParticles::HasDetector::hasCaloInfo( const LHCb::ProtoParticle*
              p->hasInfo( LHCb::ProtoParticle::CaloDepositID ) || p->hasInfo( LHCb::ProtoParticle::ShowerShape ) ||
              p->hasInfo( LHCb::ProtoParticle::ClusterMass ) ||
              //
-             p->hasInfo( LHCb::ProtoParticle::CaloNeutralSpd ) || p->hasInfo( LHCb::ProtoParticle::CaloNeutralPrs ) ||
-             p->hasInfo( LHCb::ProtoParticle::CaloNeutralEcal ) ||
+             p->hasInfo( LHCb::ProtoParticle::CaloNeutralSpd ) || p->hasInfo( LHCb::ProtoParticle::CaloNeutralEcal ) ||
              //
-             p->hasInfo( LHCb::ProtoParticle::CaloSpdE ) || p->hasInfo( LHCb::ProtoParticle::CaloPrsE ) ||
              p->hasInfo( LHCb::ProtoParticle::CaloEcalE ) || p->hasInfo( LHCb::ProtoParticle::CaloHcalE ) ||
              //
              p->hasInfo( LHCb::ProtoParticle::CaloEcalChi2 ) || p->hasInfo( LHCb::ProtoParticle::CaloBremChi2 ) ||