diff --git a/CaloFuture/CaloFutureReco/src/FutureClusterSpreadTool.cpp b/CaloFuture/CaloFutureReco/src/FutureClusterSpreadTool.cpp index 12a6fbc3ce20aca0734d8a071074d0c03f429f66..87b673570b2939f6cb3d12a1ee47cd6c4841db50 100644 --- a/CaloFuture/CaloFutureReco/src/FutureClusterSpreadTool.cpp +++ b/CaloFuture/CaloFutureReco/src/FutureClusterSpreadTool.cpp @@ -65,7 +65,7 @@ public: // ========================================================================== private: // ========================================================================== - SpreadEstimator m_estimator; + LHCb::Calo::SpreadEstimator m_estimator; Gaudi::Property<std::string> m_detData{this, "Detector", LHCb::CaloFutureAlgUtils::DeCaloFutureLocation( name() )}; const DeCalorimeter* m_det = nullptr; // ========================================================================== diff --git a/CaloFuture/CaloFutureReco/src/SpreadEstimator.cpp b/CaloFuture/CaloFutureReco/src/SpreadEstimator.cpp index 9e7aaf652adfd37ada59d39f93a8154cecb67537..ca6469a1be3b4f8a14f89700bd470025cbb8ed26 100644 --- a/CaloFuture/CaloFutureReco/src/SpreadEstimator.cpp +++ b/CaloFuture/CaloFutureReco/src/SpreadEstimator.cpp @@ -35,209 +35,212 @@ namespace { } // ========================================================================== } // namespace -// ============================================================================ -/* standard/default constructor - * @param Det pointer to calorimeter detector - */ -// ============================================================================ -SpreadEstimator::SpreadEstimator( const DeCalorimeter* Det ) : m_detector( Det ) {} -// ============================================================================ -/* set new value for calorimeter - * @param Det pointer to calorimeter detector - */ -// ============================================================================ -void SpreadEstimator::setDetector( const DeCalorimeter* Det ) { m_detector = Det; } -// ============================================================================ -/* calculation of cluster spread - * - * Error codes: - * - * - 221 - invalid source of detector information - * - 222 - the seed cell was not found - * - 223 - strange combination - * - 224 - energy is not positive - * - * @param pointer to cluster object - * @return status code - */ -// ============================================================================ -StatusCode SpreadEstimator::operator()( LHCb::CaloCluster* cluster ) const { - // ignore trivial cases - if ( !cluster ) { return StatusCode::SUCCESS; } - if ( cluster->entries().empty() ) { return StatusCode::SUCCESS; } - /// detector? - if ( !m_detector ) { return StatusCode( 221 ); } - - auto et = LHCb::CaloDataFunctor::EnergyTransverse{m_detector}; - - const double x_prec = 0.2 * Gaudi::Units::mm; - const double x_prec2 = x_prec * x_prec; - const double e_prec = 0.2 * Gaudi::Units::MeV; - - /// - double covxx = 0; - double covyx = 0; - double covyy = 0; - double xmean = 0; - double ymean = 0; - double etot = 0; - - double x0 = 0; - double y0 = 0; - - unsigned int ncells = 0; - double eall = 0; - - double cellsize = -10; - bool first = true; - - for ( const auto& entry : cluster->entries() ) { - const LHCb::CaloDigit* digit = entry.digit(); - if ( !digit ) { continue; } - /// - const double fraction = entry.fraction(); - /// - const double energy = digit->e() * fraction; - /// - const Gaudi::XYZPoint& pos = m_detector->cellCenter( digit->cellID() ); - /// - if ( entry.status() & LHCb::CaloDigitStatus::SeedCell ) cellsize = m_detector->cellSize( digit->cellID() ); - /// - eall += energy; - /// - if ( energy <= 0 ) { continue; } // CONTINUE! - // - const double weight = energy; // ATTENTION! - if ( weight <= 0 ) { continue; } // CONTINUE! - /// - ++ncells; + +namespace LHCb::Calo { + // ============================================================================ + /* standard/default constructor + * @param Det pointer to calorimeter detector + */ + // ============================================================================ + SpreadEstimator::SpreadEstimator( const DeCalorimeter* Det ) : m_detector( Det ) {} + // ============================================================================ + /* set new value for calorimeter + * @param Det pointer to calorimeter detector + */ + // ============================================================================ + void SpreadEstimator::setDetector( const DeCalorimeter* Det ) { m_detector = Det; } + // ============================================================================ + /* calculation of cluster spread + * + * Error codes: + * + * - 221 - invalid source of detector information + * - 222 - the seed cell was not found + * - 223 - strange combination + * - 224 - energy is not positive + * + * @param pointer to cluster object + * @return status code + */ + // ============================================================================ + StatusCode SpreadEstimator::operator()( LHCb::CaloCluster* cluster ) const { + // ignore trivial cases + if ( !cluster ) { return StatusCode::SUCCESS; } + if ( cluster->entries().empty() ) { return StatusCode::SUCCESS; } + /// detector? + if ( !m_detector ) { return StatusCode( 221 ); } + + auto et = LHCb::CaloDataFunctor::EnergyTransverse{m_detector}; + + const double x_prec = 0.2 * Gaudi::Units::mm; + const double x_prec2 = x_prec * x_prec; + const double e_prec = 0.2 * Gaudi::Units::MeV; + /// - const double x = pos.x(); - const double y = pos.y(); - - /// accumulate the energy - etot += weight; - - // the first good cell? - if ( first ) { - x0 = x; // adjust the bias - y0 = y; // adjust the bias - first = false; // can skip since everything propto dx and/or dy - continue; // continue + double covxx = 0; + double covyx = 0; + double covyy = 0; + double xmean = 0; + double ymean = 0; + double etot = 0; + + double x0 = 0; + double y0 = 0; + + unsigned int ncells = 0; + double eall = 0; + + double cellsize = -10; + bool first = true; + + for ( const auto& entry : cluster->entries() ) { + const LHCb::CaloDigit* digit = entry.digit(); + if ( !digit ) { continue; } + /// + const double fraction = entry.fraction(); + /// + const double energy = digit->e() * fraction; + /// + const Gaudi::XYZPoint& pos = m_detector->cellCenter( digit->cellID() ); + /// + if ( entry.status() & LHCb::CaloDigitStatus::SeedCell ) cellsize = m_detector->cellSize( digit->cellID() ); + /// + eall += energy; + /// + if ( energy <= 0 ) { continue; } // CONTINUE! + // + const double weight = energy; // ATTENTION! + if ( weight <= 0 ) { continue; } // CONTINUE! + /// + ++ncells; + /// + const double x = pos.x(); + const double y = pos.y(); + + /// accumulate the energy + etot += weight; + + // the first good cell? + if ( first ) { + x0 = x; // adjust the bias + y0 = y; // adjust the bias + first = false; // can skip since everything propto dx and/or dy + continue; // continue + } + + const double dx = x - x0; // temporary shift + const double dy = y - y0; // temporary shift + + xmean += dx * weight; + ymean += dy * weight; + // + covxx += dx * dx * weight; + covyx += dy * dx * weight; + covyy += dy * dy * weight; + /// } - const double dx = x - x0; // temporary shift - const double dy = y - y0; // temporary shift - - xmean += dx * weight; - ymean += dy * weight; - // - covxx += dx * dx * weight; - covyx += dy * dx * weight; - covyy += dy * dy * weight; + // strange combinations + if ( ncells <= 0 ) { return StatusCode( 223 ); } + // energy is not positive + if ( eall <= 1 * Gaudi::Units::MeV ) { return StatusCode( 224 ); } + // energy is not positive + if ( etot <= 1 * Gaudi::Units::MeV ) { return StatusCode( 224 ); } + // seed cell was not found + if ( cellsize <= 0 ) { return StatusCode( 222 ); } /// - } - // strange combinations - if ( ncells <= 0 ) { return StatusCode( 223 ); } - // energy is not positive - if ( eall <= 1 * Gaudi::Units::MeV ) { return StatusCode( 224 ); } - // energy is not positive - if ( etot <= 1 * Gaudi::Units::MeV ) { return StatusCode( 224 ); } - // seed cell was not found - if ( cellsize <= 0 ) { return StatusCode( 222 ); } - /// - - const double uniform = cellsize * cellsize / 12.; - - xmean /= etot; - ymean /= etot; - - covxx /= etot; - covyx /= etot; - covyy /= etot; - - covxx -= xmean * xmean; - covyx -= ymean * xmean; - covyy -= ymean * ymean; - - xmean += x0; // shift back - ymean += y0; // shift back - - LHCb::CaloPosition::Center& center = cluster->position().center(); - center( LHCb::CaloPosition::Index::X ) = xmean; - center( LHCb::CaloPosition::Index::Y ) = ymean; - - /// could do nothing else for 1 cell "clusters" - if ( 1 == ncells ) { - covxx = uniform; - covyy = uniform; - covyx = 0.0; - } else { - const double trace = covxx + covyy; - const double diff = covxx - covyy; - const double disc = std::sqrt( diff * diff + 4 * ( covyx * covyx ) ); - - // eigen values: - const double lambda1 = 0.5 * ( trace - disc ); - const double lambda2 = 0.5 * ( trace + disc ); - - // minimal eigenvalue - const double lambdaMin = std::min( lambda1, lambda2 ) / uniform; - const double covMin = std::min( covxx, covyy ) / uniform; - - const double s_Cut = 0.5; - const unsigned int s_Cells = 4; - - if ( s_Cells >= ncells || // small amout of cells wth e>0 - s_Cut >= lambdaMin || // small eigenvalue - s_Cut >= covMin ) // small eigenvalue - { - // ======================================================================= - // the matrix must be modified a bit: - const double eT = round( et( cluster ), e_prec ); - m_ratio += (float)lambdaMin; - m_energy += (float)( eT / Gaudi::Units::GeV ); - m_cells += ncells; - - // construct the matrix of eigen vectors - - const double diff1 = covxx - lambda1; - const double diff2 = covxx - lambda2; - - const LHCb::Math::abs_less<double> absLess = LHCb::Math::abs_less<double>(); - - const double phi = - // the matrix is numerically diagonal? => trivial eigenvectors - absLess( covyx, x_prec2 ) ? 0.0 : - // find "the best non-zero" eigenvector - absLess( diff1, diff2 ) ? std::atan2( covyx, diff2 ) : std::atan2( diff1, -covyx ); - - const double cphi = std::cos( phi ); - const double sphi = std::sin( phi ); - - const double s2 = sphi * sphi; - const double c2 = cphi * cphi; - const double cs = cphi * sphi; - - const double newLambda1 = std::max( lambda1, uniform ); - const double newLambda2 = std::max( lambda2, uniform ); - - // recalculate the matrix using new eigenvalues - covxx = c2 * newLambda1 + s2 * newLambda2; - covyx = cs * ( newLambda1 - newLambda2 ); - covyy = s2 * newLambda1 + c2 * newLambda2; + const double uniform = cellsize * cellsize / 12.; + + xmean /= etot; + ymean /= etot; + + covxx /= etot; + covyx /= etot; + covyy /= etot; + + covxx -= xmean * xmean; + covyx -= ymean * xmean; + covyy -= ymean * ymean; + + xmean += x0; // shift back + ymean += y0; // shift back + + LHCb::CaloPosition::Center& center = cluster->position().center(); + center( LHCb::CaloPosition::Index::X ) = xmean; + center( LHCb::CaloPosition::Index::Y ) = ymean; + + /// could do nothing else for 1 cell "clusters" + if ( 1 == ncells ) { + covxx = uniform; + covyy = uniform; + covyx = 0.0; + } else { + const double trace = covxx + covyy; + const double diff = covxx - covyy; + const double disc = std::sqrt( diff * diff + 4 * ( covyx * covyx ) ); + + // eigen values: + const double lambda1 = 0.5 * ( trace - disc ); + const double lambda2 = 0.5 * ( trace + disc ); + + // minimal eigenvalue + const double lambdaMin = std::min( lambda1, lambda2 ) / uniform; + const double covMin = std::min( covxx, covyy ) / uniform; + + const double s_Cut = 0.5; + const unsigned int s_Cells = 4; + + if ( s_Cells >= ncells || // small amout of cells wth e>0 + s_Cut >= lambdaMin || // small eigenvalue + s_Cut >= covMin ) // small eigenvalue + { + // ======================================================================= + // the matrix must be modified a bit: + const double eT = round( et( cluster ), e_prec ); + m_ratio += (float)lambdaMin; + m_energy += (float)( eT / Gaudi::Units::GeV ); + m_cells += ncells; + + // construct the matrix of eigen vectors + + const double diff1 = covxx - lambda1; + const double diff2 = covxx - lambda2; + + const LHCb::Math::abs_less<double> absLess = LHCb::Math::abs_less<double>(); + + const double phi = + // the matrix is numerically diagonal? => trivial eigenvectors + absLess( covyx, x_prec2 ) ? 0.0 : + // find "the best non-zero" eigenvector + absLess( diff1, diff2 ) ? std::atan2( covyx, diff2 ) : std::atan2( diff1, -covyx ); + + const double cphi = std::cos( phi ); + const double sphi = std::sin( phi ); + + const double s2 = sphi * sphi; + const double c2 = cphi * cphi; + const double cs = cphi * sphi; + + const double newLambda1 = std::max( lambda1, uniform ); + const double newLambda2 = std::max( lambda2, uniform ); + + // recalculate the matrix using new eigenvalues + covxx = c2 * newLambda1 + s2 * newLambda2; + covyx = cs * ( newLambda1 - newLambda2 ); + covyy = s2 * newLambda1 + c2 * newLambda2; + } } - } - LHCb::CaloPosition::Spread& spread = cluster->position().spread(); - spread( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = covxx; - spread( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::X ) = covyx; - spread( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = covyy; + LHCb::CaloPosition::Spread& spread = cluster->position().spread(); + spread( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = covxx; + spread( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::X ) = covyx; + spread( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = covyy; - /// - return StatusCode::SUCCESS; -} + /// + return StatusCode::SUCCESS; + } +} // namespace LHCb::Calo // ============================================================================ // The End // ============================================================================ diff --git a/CaloFuture/CaloFutureReco/src/SpreadEstimator.h b/CaloFuture/CaloFutureReco/src/SpreadEstimator.h index 55dbf05e36b5ced386a91bcb94c7a98fd9255dbd..2a3018e638380eecaff2b76dbaff2906d9a99106 100644 --- a/CaloFuture/CaloFutureReco/src/SpreadEstimator.h +++ b/CaloFuture/CaloFutureReco/src/SpreadEstimator.h @@ -26,71 +26,73 @@ // ============================================================================ class DeCalorimeter; // from CaloDet package // ============================================================================ -/** @class SpreadEstimator SpreadEstimator.h CaloFutureUtils/SpreadEstimator.h - * - * simple helper class for estimation of spread of - * CaloCluster object. - * - * - * @author Vanya Belyaev Ivan.Belyaev@itep.ru - * @date 22/11/2001 - */ -class SpreadEstimator { - // ========================================================================== -public: - // ========================================================================== - /** standard/default constructor - * @param Det pointer to calorimeter detector - */ - SpreadEstimator( const DeCalorimeter* Det = nullptr ); - // ========================================================================== - /** calculation of cluster spread +namespace LHCb::Calo { + /** @class SpreadEstimator SpreadEstimator.h CaloFutureUtils/SpreadEstimator.h * - * Error codes: - * - 221 - invalid source of detector information + * simple helper class for estimation of spread of + * CaloCluster object. * - * @param cluster pointer to cluster object - * @return status code - */ - StatusCode operator()( LHCb::CaloCluster* cluster ) const; - // ========================================================================== - /** calculate spread for cluster - * @param cluster pointer to cluster - * @return status code - */ - StatusCode calculateSpread( LHCb::CaloCluster* cluster ) const { return ( *this )( cluster ); } - // ========================================================================== -public: - // ========================================================================== - /** set new value for calorimeter - * @param Det pointer to calorimeter detector - */ - void setDetector( const DeCalorimeter* Det ); - // ========================================================================== - /** simple accessor to DeCalorimeter object - * @return pointer to detector + * + * @author Vanya Belyaev Ivan.Belyaev@itep.ru + * @date 22/11/2001 */ - [[deprecated]] const DeCalorimeter* detector() const { return m_detector; } - // ========================================================================== - /// get the counter for problematic cases - const StatEntity& invalidRatio() const { return m_ratio; } - /// get the counter for problematic cases - const StatEntity& invalidCells() const { return m_cells; } - /// get the counter for problematic cases - const StatEntity& invalidEnergy() const { return m_energy; } - // ========================================================================== -private: - // ========================================================================== - /// the detector elemenet - const DeCalorimeter* m_detector; // the detector elemenet - // ========================================================================== - /// counter of invalid cells - mutable StatEntity m_cells; // counter for invalid cells - /// counter of invalid size ratio - mutable StatEntity m_ratio; // counter for invalid size ratio - /// counter of invalid energy - mutable StatEntity m_energy; // counter for invalid energy -}; + class SpreadEstimator { + // ========================================================================== + public: + // ========================================================================== + /** standard/default constructor + * @param Det pointer to calorimeter detector + */ + SpreadEstimator( const DeCalorimeter* Det = nullptr ); + // ========================================================================== + /** calculation of cluster spread + * + * Error codes: + * - 221 - invalid source of detector information + * + * @param cluster pointer to cluster object + * @return status code + */ + StatusCode operator()( LHCb::CaloCluster* cluster ) const; + // ========================================================================== + /** calculate spread for cluster + * @param cluster pointer to cluster + * @return status code + */ + StatusCode calculateSpread( LHCb::CaloCluster* cluster ) const { return ( *this )( cluster ); } + // ========================================================================== + public: + // ========================================================================== + /** set new value for calorimeter + * @param Det pointer to calorimeter detector + */ + void setDetector( const DeCalorimeter* Det ); + // ========================================================================== + /** simple accessor to DeCalorimeter object + * @return pointer to detector + */ + [[deprecated]] const DeCalorimeter* detector() const { return m_detector; } + // ========================================================================== + /// get the counter for problematic cases + const StatEntity& invalidRatio() const { return m_ratio; } + /// get the counter for problematic cases + const StatEntity& invalidCells() const { return m_cells; } + /// get the counter for problematic cases + const StatEntity& invalidEnergy() const { return m_energy; } + // ========================================================================== + private: + // ========================================================================== + /// the detector elemenet + const DeCalorimeter* m_detector; // the detector elemenet + // ========================================================================== + /// counter of invalid cells + mutable StatEntity m_cells; // counter for invalid cells + /// counter of invalid size ratio + mutable StatEntity m_ratio; // counter for invalid size ratio + /// counter of invalid energy + mutable StatEntity m_energy; // counter for invalid energy + }; +} // namespace LHCb::Calo // ============================================================================ // The END // ============================================================================