Update of Calo Cluster Efficiency algorithm

Re-define the selection of "reconstructible" MC particles as particles that have deposited >90% of its energy in ECAL digits and those digits have >90% of contribution from the same MC particle. Then a MC particle is "reconstructed" if there exists a cluster with the same seedID with reconstructed energy >90% of the MC particle deposited energy.

This MR goes together with lhcb/Moore!945 as the function needs two new arguments when called:

• inputRelationsDigit: MC to Digits table
• caloDigits: ECAL Digits

And also with LHCb!3245 (merged).

This upgrade is ready to be tested.

Calo/CaloAssociators/src/CaloClusterEfficiency.cpp

@@ -45,14 +45,17 @@ using MC2CellTable = LHCb::CaloFuture2MC::MC2ClusterTable;
 namespace LHCb::Calo::Algorithms {
 
   class ClusterEfficiency
-      : public Gaudi::Functional::Consumer<void( const Cell2MCTable&, const LHCb::MCParticles&,
-                                                 const LHCb::Event::Calo::v2::Clusters&, const DeCalorimeter& ),
+      : public Gaudi::Functional::Consumer<void( const Cell2MCTable&, const LHCb::CaloFuture2MC::DigitTable&,
+                                                 const LHCb::MCParticles&, const LHCb::Event::Calo::v2::Clusters&,
+                                                 const LHCb::Event::Calo::Digits&, const DeCalorimeter& ),
                                            LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, DeCalorimeter>> {
   public:
-    using base_type = Gaudi::Functional::Consumer<void( const Cell2MCTable&, const LHCb::MCParticles&,
-                                                        const LHCb::Event::Calo::v2::Clusters&, const DeCalorimeter& ),
-                                                  LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, DeCalorimeter>>;
-    using KeyValue  = typename base_type::KeyValue;
+    using base_type =
+        Gaudi::Functional::Consumer<void( const Cell2MCTable&, const LHCb::CaloFuture2MC::DigitTable&,
+                                          const LHCb::MCParticles&, const LHCb::Event::Calo::v2::Clusters&,
+                                          const LHCb::Event::Calo::Digits&, const DeCalorimeter& ),
+                                    LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, DeCalorimeter>>;
+    using KeyValue = typename base_type::KeyValue;
 
     // ==========================================================================
     /// standard constructor
@@ -61,16 +64,19 @@ namespace LHCb::Calo::Algorithms {
                     pSvc,
                     {
                         KeyValue{"inputRelations", "Relations/" + LHCb::CaloClusterLocation::Default},
+                        KeyValue{"inputRelationsDigit", "Relations/" + LHCb::CaloDigitLocation::Default},
                         KeyValue{"inputMCParticles", LHCb::MCParticleLocation::Default},
                         KeyValue{"recoClusters", LHCb::Event::Calo::ClusterLocation::Default},
+                        KeyValue{"caloDigits", LHCb::CaloDigitLocation::Ecal},
                         KeyValue{"Detector", DeCalorimeterLocation::Ecal},
                     }} {}
 
     // ============================================================================
     // algorithm execution
     // ============================================================================
-    void operator()( const Cell2MCTable& table, const LHCb::MCParticles& mcparts,
-                     const LHCb::Event::Calo::v2::Clusters& clusters, const DeCalorimeter& calo ) const override {
+    void operator()( const Cell2MCTable& table, const LHCb::CaloFuture2MC::DigitTable& tableDigits,
+                     const LHCb::MCParticles& mcparts, const LHCb::Event::Calo::v2::Clusters& clusters,
+                     const LHCb::Event::Calo::Digits& digits, const DeCalorimeter& calo ) const override {
 
       const auto cluster_index = clusters.index();
       if ( !cluster_index ) {
@@ -87,9 +93,11 @@ namespace LHCb::Calo::Algorithms {
       unsigned int nreconstructed   = 0;
       bool         is_reconstructed;
 
-      // create MCParticle -> CaloCluster relations table from inverse one
+      // NEW: create MCParticle -> CaloDigit relations table from inverse one
       // LHCb::RelationWeighted1D(table, inverse_tag) creates inverse table
-      auto mc2cellTable = MC2CellTable( table, MC2CellTable::inverse_tag{} );
+      auto mc2cellTable = MC2CellTable( tableDigits, MC2CellTable::inverse_tag{} );
+      // create MCParticle -> CaloCluster relations table from inverse one
+      auto mc2clusterTable = MC2CellTable( table, MC2CellTable::inverse_tag{} );
 
       // loop over all mc particles
       for ( const auto& mcp : mcparts ) {
@@ -110,23 +118,43 @@ namespace LHCb::Calo::Algorithms {
         // check endVertex, also fills info into ntuple
         if ( not check_endVtx( mcp, tuple ) ) continue;
 
-        // find largest E deposited by mc particle in an ECAL cluster ("weight")
-        // use inverse MCParticle -> CellID table for this
-        // TODO: find E deposited in ECAL, not in a given cluster, to have a
-        // definition of reconstructible independent of the reconstruction itself!
-        const auto& [cell_id, weight] = largestWeight( mc2cellTable.relations( mcp ) );
-        if ( !cell_id ) continue;
+        // NEW find E deposited in ECAL from the mcp and find total E deposited in ECAL
+        // for the digits concerning the mcp (both are independent from the reconstructed clusters).
+        float e_deposited   = 0.0; // energy from the mcp deposited in ecal
+        float e_calo_digits = 0.0; // energy of the calo digits where the mcp has deposited energy
+
+        // Iterate through all the digits of the mcp
+        for ( const auto& [_, cell_id, weight] : mc2cellTable.relations( mcp ) ) {
+          e_deposited += weight; // accumulate weighted energy from the mcp
+          for ( const auto& digit : digits ) {
+            if ( digit.cellID() == cell_id ) { // find the same digit from the ecal
+              e_calo_digits += digit.energy(); // accumulate energy from the calo digits
+            }
+          }
+        }
+
         if ( msgLevel( MSG::DEBUG ) ) debug() << "MCParticle has desposited E in ECAL" << endmsg;
 
-        // reco'ble: fraction of MCParticle E deposited in cluster > m_minMCFraction
-        float f_mcp = weight / mcp->momentum().e();
+        // NEW: check that the mcp has deposited most of its energy in the calo
+        float f_mcp = e_deposited / mcp->momentum().e();
         if ( f_mcp < m_minMCfraction ) continue;
+
+        // NEW: check that the major contribution to the calo digits is from the mcp
+        float f_e = e_deposited / e_calo_digits;
+        if ( f_e < 0.9 ) continue;
+
+        // NEW: if the mcp has deposited >90% of its energy in the calo and
+        // there are no other mcps contributing to the same digits in more than 10%
+        // then mcp is reconstructible.
         nreconstructible++;
         if ( msgLevel( MSG::DEBUG ) ) debug() << "MCParticle is reco'ble" << endmsg;
 
         // fill tuple
         auto sc = tuple->column( "max_fMCP", f_mcp );
 
+        // find cell_id of the seed in the mcp
+        const auto& [cell_id, weight] = largestWeight( mc2clusterTable.relations( mcp ) );
+
         // find cluster matching MCParticle from cell_id
         auto it = cluster_index.find( cell_id ); // returns Iterator (see CaloClusters_v2)
         if ( it == cluster_index.end() ) {
@@ -144,7 +172,8 @@ namespace LHCb::Calo::Algorithms {
         auto cl_et   = eT_( &cl );
 
         // fraction of cluster energy from this MCParticle
-        float f_cl = ( cl_e > 0. ) ? weight / cl_e : 0.;
+        // NEW: e_deposited (energy from the mcp deposited in calo) instead of weight (energy of the cluster seed)
+        float f_cl = ( cl_e > 0. ) ? e_deposited / cl_e : 0.;
 
         // fill tuple
         sc &= tuple->column( "f_cl", f_cl );
@@ -243,4 +272,4 @@ namespace LHCb::Calo::Algorithms {
 
 // ============================================================================
 // The END
-// ============================================================================
+// ============================================================================
\ No newline at end of file