Commit 982c63c0 authored by Christos Anastopoulos's avatar Christos Anastopoulos Committed by Edward Moyse
Browse files

Egamma Move helper methods , no needing any class member to anonymous namespace

parent a081765e
......@@ -6,7 +6,7 @@
#define EGAMMACALOTOOLS_EGAMMASTRIPSSHAPE_H
/// @class egammaStripsShape
/// @brief EM cluster shower shape calculations in 1st ECAL sampling
/// @brief EM cluster shower shape calculations in 1st ECAL sampling
/// Calculate the width in the strip layer around the eta,phi of
/// the hottest cell in the **middle layer**.
///
......@@ -14,85 +14,46 @@
/// @author Christos Anastopoulos
///
class CaloCellContainer;
class CaloDetDescrManager;
class LArEM_ID;
#include "xAODCaloEvent/CaloCluster.h"
#include "CaloGeoHelpers/CaloSampling.h"
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
#include "CaloGeoHelpers/CaloSampling.h"
#include "CaloIdentifier/CaloCell_ID.h"
#include "GaudiKernel/ToolHandle.h"
#include "egammaInterfaces/IegammaStripsShape.h"
#include "xAODCaloEvent/CaloCluster.h"
class egammaStripsShape : public AthAlgTool, virtual public IegammaStripsShape {
class egammaStripsShape
: public AthAlgTool
, virtual public IegammaStripsShape
{
public:
/** @brief Default constructor*/
egammaStripsShape(const std::string& type,
const std::string& name,
const IInterface* parent);
/** @brief Destructor*/
~egammaStripsShape();
/** @brief AlgTool initialize method.*/
StatusCode initialize() override;
/** @brief AlgTool finalize method */
StatusCode finalize() override;
/** @brief AlgTool main method */
virtual StatusCode execute(const xAOD::CaloCluster& cluster,
const CaloDetDescrManager& cmgr,
Info& info) const override final;
private:
/** @brief From the original (eta,phi) position, find the location
(sampling, barrel/end-cap, granularity) */
/** @brief set an array of energies,eta,phi in ~40 strips around max*/
void setArray(const xAOD::CaloCluster& cluster,
const CaloDetDescrManager& cmgr,
CaloSampling::CaloSample sam, double eta, double phi,
double deta, double dphi, double* enecell, double* etacell,
double* gracell, int* ncell) const;
/** @brief check index of seed in the array*/
void setIndexSeed(Info& info, const double* etacell, const double* gracell) const;
/** @brief set total width in strips*/
void setWstot(Info& info, double deta, const double* enecell, const double* etacell, const int* ncell) const;
/** @brief set fraction of energy in 2nd sampling*/
void setF2(Info& info, const double* enecell,const double eallsamples) const ;
/** @brief set energy in 3x1 and in 15x3 strips*/
void setEnergy(Info& info, const double* enecell) const;
/** @brief set asymmetry*/
void setAsymmetry(Info& info, const double* enecell) const;
/** @brief set width in three strips*/
void setWs3(Info& info, const xAOD::CaloCluster::CaloSample sam, const xAOD::CaloCluster& cluster,
const double* enecell, const double* etacell,const int* ncell) const;
/** @brief set difference between eta of max and eta of cells*/
double setDeltaEtaTrackShower(int nstrips,int ieta, const double* enecell) const;
/** @brief set width in 5 strips*/
void setWidths5(Info& info, const double* enecell) const;
/** @brief set energy of strip with maximum energy*/
void setEmax(Info& info, const double* enecell) const;
/** @brief set energy of the second local maximum*/
int setEmax2(Info& info, const double* enecell, const double* gracell, const int* ncell) const;
/** @brief set energy of strip with minimum energy*/
void setEmin(int ncsec1,Info& info, const double* enecell, const double* gracell, const int* ncell ) const;
/** @brief set M.S's valley*/
void setValley(Info& info, double* enecell) const;
/** @brief set fraction of energy outside shower core
(E(+/-3strips)-E(+/-1strips))/ E(+/-1strips) */
void setFside(Info& info, const double* enecell, const double* gracell, const int* ncell) const;
/** @brief set F1core*/
void setF1core(Info& info, const xAOD::CaloCluster& cluster) const;
// calculate quantities based on information in the strips in a region
// around the cluster.
/** @brief boolean to calculate all variables*/
Gaudi::Property<bool> m_ExecAllVariables {this,
"ExecAllVariables", true, "flag used by trigger"};
/** @brief Default constructor*/
egammaStripsShape(const std::string& type,
const std::string& name,
const IInterface* parent);
/** @brief Destructor*/
~egammaStripsShape() = default;
/** @brief AlgTool initialize method.*/
StatusCode initialize() override;
/** @brief AlgTool finalize method */
StatusCode finalize() override;
/** @brief AlgTool main method */
virtual StatusCode execute(const xAOD::CaloCluster& cluster,
const CaloDetDescrManager& cmgr,
Info& info) const override final;
private:
/** @brief boolean to calculate all variables*/
Gaudi::Property<bool> m_ExecAllVariables{ this,
"ExecAllVariables",
true,
"flag used by trigger" };
};
#endif
......@@ -11,11 +11,11 @@
#include <cmath>
#include <vector>
#include "GaudiKernel/SystemOfUnits.h"
#include "CaloConditions/CaloAffectedRegionInfoVec.h"
#include "CaloIdentifier/CaloCell_ID.h"
#include "CaloIdentifier/LArEM_ID.h"
#include "FourMomUtils/P4Helpers.h"
#include "GaudiKernel/SystemOfUnits.h"
#include "Identifier/HWIdentifier.h"
#include "StoreGate/ReadHandle.h"
#include "StoreGate/StoreGateSvc.h"
......@@ -44,6 +44,40 @@ isbadtilecell(CaloCellList& ccl,
}
return isbadtilecell;
}
bool
findCentralCell(const xAOD::CaloCluster* cluster, Identifier& cellCentrId)
{
bool thereIsACentrCell = false;
// LOOP OVER CLUSTER TO FIND THE CENTRAL CELL IN S2
xAOD::CaloCluster::const_cell_iterator cellIter = cluster->cell_begin();
xAOD::CaloCluster::const_cell_iterator cellIterEnd = cluster->cell_end();
float clusEta = cluster->eta();
float clusPhi = cluster->phi();
float energymax = -999999.;
for (; cellIter != cellIterEnd; cellIter++) {
const CaloCell* cell = (*cellIter);
if (!cell) {
continue;
}
float eta = cell->eta();
float phi = cell->phi();
float energy = cell->energy();
CaloSampling::CaloSample layer = cell->caloDDE()->getSampling();
if (fabs(eta - clusEta) < 0.025 &&
fabs(P4Helpers::deltaPhi(phi, clusPhi)) < 0.025 &&
(layer == CaloSampling::EMB2 || layer == CaloSampling::EME2) &&
(energy > energymax)) {
energymax = energy;
cellCentrId = cellIter->ID();
thereIsACentrCell = true;
}
}
return thereIsACentrCell;
}
}
egammaOQFlagsBuilder::egammaOQFlagsBuilder(const std::string& type,
......@@ -93,41 +127,6 @@ egammaOQFlagsBuilder::finalize()
return StatusCode::SUCCESS;
}
bool
egammaOQFlagsBuilder::findCentralCell(const xAOD::CaloCluster* cluster,
Identifier& cellCentrId) const
{
bool thereIsACentrCell = false;
// LOOP OVER CLUSTER TO FIND THE CENTRAL CELL IN S2
xAOD::CaloCluster::const_cell_iterator cellIter = cluster->cell_begin();
xAOD::CaloCluster::const_cell_iterator cellIterEnd = cluster->cell_end();
float clusEta = cluster->eta();
float clusPhi = cluster->phi();
float energymax = -999999.;
for (; cellIter != cellIterEnd; cellIter++) {
const CaloCell* cell = (*cellIter);
if (!cell){
continue;
}
float eta = cell->eta();
float phi = cell->phi();
float energy = cell->energy();
CaloSampling::CaloSample layer = cell->caloDDE()->getSampling();
if (fabs(eta - clusEta) < 0.025 &&
fabs(P4Helpers::deltaPhi(phi, clusPhi)) < 0.025 &&
(layer == CaloSampling::EMB2 || layer == CaloSampling::EME2) &&
(energy > energymax)) {
energymax = energy;
cellCentrId = cellIter->ID();
thereIsACentrCell = true;
}
}
return thereIsACentrCell;
}
bool
egammaOQFlagsBuilder::isCore(
const Identifier Id,
......@@ -188,8 +187,7 @@ egammaOQFlagsBuilder::execute(const EventContext& ctx,
// Find the central cell in the middle layer
Identifier cellCentrId;
bool foundCentralCell =
egammaOQFlagsBuilder::findCentralCell(cluster, cellCentrId);
bool foundCentralCell = findCentralCell(cluster, cellCentrId);
// Set timing bit
const double absEnergyGeV = fabs(cluster->e() * (1. / Gaudi::Units::GeV));
......
......@@ -102,9 +102,6 @@ private:
std::vector<IdentifierHash> findNeighbours(
const Identifier cellCentrId) const;
bool findCentralCell(const xAOD::CaloCluster* cluster,
Identifier& cellCentrId) const;
Gaudi::Property<double> m_QCellCut{ this, "QCellCut", 4000. };
Gaudi::Property<double> m_QCellHECCut{ this, "QCellHECCut", 60000. };
Gaudi::Property<double> m_QCellSporCut{ this, "QCellSporCut", 4000. };
......
......@@ -18,6 +18,76 @@
#include "StoreGate/ReadHandle.h"
namespace {
double
electronPhiResoA(double eta)
{
if (eta < 0.30)
return 0.000191492;
if (eta < 0.60)
return 9.35047e-05 + 0.000392766 * eta;
if (eta < 0.80)
return 0.000327201;
if (eta < 1.05)
return 0.000141755;
if (eta < 1.35)
return (-1.07475 + 1.15372 * eta) * 1e-3;
if (eta < 1.55)
return (-15.2133 + 11.2163 * eta) * 1e-3;
if (eta < 1.85)
return 0.00128452 - 0.00053016 * eta;
if (eta < 2.30)
return -0.000665622 + 0.00052136 * eta;
return 0.000327754;
}
double
electronPhiResoB(double eta)
{
if (eta < 0.65)
return 0.0285262 + 0.00985529 * eta;
if (eta < 1.04)
return -0.0690774 + 0.166424 * eta;
if (eta < 1.25)
return 0.0769113 + 0.0149434 * eta;
if (eta < 1.55)
return -0.407594 + 0.393218 * eta;
if (eta < 1.95)
return 0.415602 - 0.172824 * eta;
if (eta < 2.05)
return 0.0840844;
if (eta < 2.40)
return 0.187563 - 0.0472463 * eta;
return 0.0693652;
}
double
electronPhiResolution(double eta, double energy)
{
eta = fabs(eta);
return electronPhiResoA(eta) + electronPhiResoB(eta) / energy;
}
} // end of anonymous namespace
CaloCluster_OnTrackBuilder::CaloCluster_OnTrackBuilder(const std::string& t,
const std::string& n,
......@@ -32,7 +102,7 @@ CaloCluster_OnTrackBuilder::~CaloCluster_OnTrackBuilder() {}
StatusCode CaloCluster_OnTrackBuilder::initialize() {
ATH_MSG_INFO("Initializing CaloCluster_OnTrackBuilder");
......@@ -41,8 +111,8 @@ StatusCode CaloCluster_OnTrackBuilder::initialize() {
ATH_MSG_FATAL ( "Unable to retrieve the instance " << m_calosurf.name() << "... Exiting!" );
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
......@@ -54,31 +124,31 @@ Trk::CaloCluster_OnTrack* CaloCluster_OnTrackBuilder::buildClusterOnTrack(
const xAOD::CaloCluster* cluster, int charge) const {
ATH_MSG_DEBUG("Building Trk::CaloCluster_OnTrack");
if(!m_useClusterPhi && !m_useClusterEta && !m_useClusterEnergy){
ATH_MSG_WARNING("CaloCluster_OnTrackBuilder is configured incorrectly");
return nullptr;
ATH_MSG_WARNING("CaloCluster_OnTrackBuilder is configured incorrectly");
return nullptr;
}
if(!cluster) return nullptr;
const Trk::Surface* surface = getCaloSurface( cluster );
if(!surface) return nullptr;
const Trk::LocalParameters* lp =getClusterLocalParameters( cluster, surface, charge );
if (!lp){
delete surface;
return nullptr;
}
const Amg::MatrixX *em =getClusterErrorMatrix( cluster, surface, charge );
const Amg::MatrixX *em =getClusterErrorMatrix( cluster, surface, charge );
if (!em){
delete surface;
delete lp;
return nullptr;
}
Trk::CaloCluster_OnTrack* ccot = new Trk::CaloCluster_OnTrack( *lp, *em, *surface );
delete surface;
delete lp;
......@@ -87,50 +157,49 @@ Trk::CaloCluster_OnTrack* CaloCluster_OnTrackBuilder::buildClusterOnTrack(
if(ccot) {
ATH_MSG_DEBUG("Successful build of Trk::CaloCluster_OnTrack");
}
return ccot;
}
const Trk::Surface* CaloCluster_OnTrackBuilder::getCaloSurface( const xAOD::CaloCluster* cluster ) const
{
const Trk::Surface* destinationSurface = nullptr;
// Determine if we want to extrapolate to the barrel or endcap. If in the crack choose the
// detector with largest amount of energy in the second sampling layer
// Determine if we want to extrapolate to the barrel or endcap. If in the crack choose the
// detector with largest amount of energy in the second sampling layer
if ( xAOD::EgammaHelpers::isBarrel( cluster ) )
{
destinationSurface = m_calosurf->CreateUserSurface (CaloCell_ID::EMB2, 0. , cluster->eta() );
} else{
} else{
destinationSurface = m_calosurf->CreateUserSurface (CaloCell_ID::EME2, 0. , cluster->eta() );
}
return destinationSurface;
}
const Trk::LocalParameters* CaloCluster_OnTrackBuilder::getClusterLocalParameters( const xAOD::CaloCluster* cluster,
const Trk::LocalParameters* CaloCluster_OnTrackBuilder::getClusterLocalParameters( const xAOD::CaloCluster* cluster,
const Trk::Surface* surf,
int charge) const
{
Amg::Vector3D surfRefPoint = surf->globalReferencePoint();
//std::cout << "REFPOINT " << "[r,phi,z] = [ " << surfRefPoint.perp() << ", " << surfRefPoint.phi() << ", " << surfRefPoint.z() << " ]" <<std::endl;
double eta = cluster->eta();
double theta = 2*atan(exp(-eta)); // -log(tan(theta/2));
double tantheta = tan(theta);
double phi = cluster->phi();
double clusterQoverE = cluster->calE() !=0 ? (double)charge/cluster->calE() : 0;
Trk::LocalParameters* newLocalParameters(nullptr);
if ( xAOD::EgammaHelpers::isBarrel( cluster ) ){
//Two corindate in a cyclinder are
//Two corindate in a cyclinder are
//Trk::locRPhi = 0 (ie phi)
//Trk::locZ = 1(ie z)
double r = surfRefPoint.perp() ;
double z = tantheta == 0 ? 0. : r/tantheta;
double z = tantheta == 0 ? 0. : r/tantheta;
Trk::DefinedParameter locRPhi( r * phi , Trk::locRPhi ) ;
Trk::DefinedParameter locZ ( z , Trk::locZ ) ;
Trk::DefinedParameter qOverP ( clusterQoverE , Trk::qOverP ) ;
......@@ -139,12 +208,12 @@ const Trk::LocalParameters* CaloCluster_OnTrackBuilder::getClusterLocalParamet
if(m_useClusterEta)defPar.push_back( locZ ) ;
if(m_useClusterEnergy)defPar.push_back( qOverP );
newLocalParameters = new Trk::LocalParameters( defPar ) ;
} else{
} else{
//Local paramters of a disk are
//Trk::locR = 0
//Trk::locPhi = 1
double z = surfRefPoint.z();
double r = z*tantheta;
double r = z*tantheta;
Trk::DefinedParameter locR ( r , Trk::locR ) ;
Trk::DefinedParameter locPhi( phi , Trk::locPhi ) ;
Trk::DefinedParameter qOverP ( clusterQoverE , Trk::qOverP ) ;
......@@ -155,7 +224,7 @@ const Trk::LocalParameters* CaloCluster_OnTrackBuilder::getClusterLocalParamet
newLocalParameters = new Trk::LocalParameters( defPar ) ;
}
return newLocalParameters;
}
......@@ -173,21 +242,21 @@ const Amg::MatrixX* CaloCluster_OnTrackBuilder::getClusterErrorMatrix( const
double etaerr = 10; //10mm large error as currently we dont want to use this measurement
etaerr = etaerr < 1.e-10 ? 10. : pow(etaerr,2);
double energyerr = pow( 0.10 * sqrt(cluster->calE()*1e-3)*1000 ,-4 ) ;
double energyerr = pow( 0.10 * sqrt(cluster->calE()*1e-3)*1000 ,-4 ) ;
Amg::MatrixX covMatrix;
if ( xAOD::EgammaHelpers::isBarrel( cluster ) ){
//Two corindate in a cyclinder are
//Two corindate in a cyclinder are
//Trk::locRPhi = 0 (ie phi)
//Trk::locZ = 1(ie z)
//Trk::locZ = 1(ie z)
Amg::Vector3D surfRefPoint = surf->globalReferencePoint();
double r2 = pow(surfRefPoint[0],2 );
int indexCount(0);
if(m_useClusterPhi){
if(m_useClusterPhi){
++indexCount;
covMatrix( indexCount, indexCount ) = phierr * r2 ;
}
......@@ -199,13 +268,13 @@ const Amg::MatrixX* CaloCluster_OnTrackBuilder::getClusterErrorMatrix( const
++indexCount;
covMatrix( indexCount, indexCount ) = energyerr ;
}
} else{
} else{
//Local paramters of a disk are
//Trk::locR = 0
//Trk::locPhi = 1
int indexCount(0);
if(m_useClusterEta){
++indexCount;
covMatrix( indexCount, indexCount ) = etaerr ;
......@@ -223,91 +292,26 @@ const Amg::MatrixX* CaloCluster_OnTrackBuilder::getClusterErrorMatrix( const
const Amg::MatrixX *result= new Amg::MatrixX(covMatrix);
return result;
}
double CaloCluster_OnTrackBuilder::getClusterPhiError( const xAOD::CaloCluster* cluster ) const
double CaloCluster_OnTrackBuilder::getClusterPhiError( const xAOD::CaloCluster* cluster )
{
/** Error on theta = C(eta) mrad/sqrt(Energy) */
/** NOTE: Eta dependence of C will be implemented later. C=10mrad based on EG? note*/
/** NOTE: Eta dependence of C will be implemented later. C=10mrad based on EG? note*/
/** Note these should be take from EMError details*/
double clusterEnergy = cluster->calE()*1e-3;
/** Error on phi = C(eta) mrad/sqrt(Energy) */
double error = electronPhiResolution( cluster->eta() , clusterEnergy );
double error = electronPhiResolution( cluster->eta() , clusterEnergy );
return error * 1.1;
}
double CaloCluster_OnTrackBuilder::electronPhiResolution(double eta, double energy) const
{
eta = fabs( eta );
return electronPhiResoA( eta ) + electronPhiResoB( eta ) / energy;
}
double CaloCluster_OnTrackBuilder::electronPhiResoA(double eta) const
{
if (eta < 0.30)
return 0.000191492;
if (eta < 0.60)
return 9.35047e-05 + 0.000392766 * eta;
if (eta < 0.80)
return 0.000327201;
if (eta < 1.05)
return 0.000141755;
if (eta < 1.35)
return (-1.07475 + 1.15372 * eta) * 1e-3;
if (eta < 1.55)
return (-15.2133 + 11.2163 * eta) * 1e-3;
if (eta < 1.85)
return 0.00128452 - 0.00053016 * eta;
if (eta < 2.30)
return -0.000665622 + 0.00052136 * eta;
return 0.000327754;
}
double CaloCluster_OnTrackBuilder::electronPhiResoB(double eta) const
{
if (eta < 0.65)
return 0.0285262 + 0.00985529 * eta;
if (eta < 1.04)