Skip to content
Snippets Groups Projects
Commit 423a16dd authored by Walter Lampl's avatar Walter Lampl
Browse files

Merge branch 'newClusterMoment' into 'master'

New topo-cluster moment & data size reduction for forward topo-towers

See merge request atlas/athena!42086
parents ace94e63 7ba6e9c9
No related branches found
No related tags found
No related merge requests found
Showing
with 262 additions and 189 deletions
......@@ -110,7 +110,8 @@ AODMoments=[#"LATERAL"
,"EM_PROBABILITY"
#,"PTD"
,"BadChannelList"
,#"LATERAL"
#,"LATERAL"
,"SECOND_TIME"
]
if jobproperties.CaloRecFlags.doExtendedClusterMoments.get_Value():
......@@ -215,8 +216,6 @@ if jobproperties.CaloRecFlags.doCaloFwdTopoTower.get_Value():
AuxListItem+="."+moment
CaloClusterItemList+=[AuxListItem]
CaloAODList+=CaloClusterItemList
# E4' cells
......
......@@ -234,7 +234,7 @@ class CaloClusterTopoGetter ( Configured ) :
,"AVG_TILE_Q"
,"PTD"
,"MASS"
,"EM_PROBABILITY"
,"SECOND_TIME"
]
doDigiTruthFlag = False
......
# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
......@@ -135,6 +135,7 @@ def getTopoMoments(configFlags):
,"AVG_TILE_Q"
,"PTD"
,"MASS"
, "SECOND_TIME"
]
# Disable for now, as broken on MC
......
......@@ -24,6 +24,7 @@ def TowersFromClustersDict(clusterBuilderName = 'TowerFromClusterTool',
cellEnergyThreshold = 0.,
applyLCW = False,
buildCombinedSignal = False,
removeSamplingData = True,
clusterRange = 5.):
''' Configuration dictionary for tower-to-cluster converter
'''
......@@ -38,7 +39,8 @@ def TowersFromClustersDict(clusterBuilderName = 'TowerFromClusterTool',
'PrepareLCW' : applyLCW, ### (control) prepare (and apply) LCW
'DoCellIndexCheck' : doCellIndexCheck, ### (control) check cell hash indices
'BuildCombinedTopoSignal' : buildCombinedSignal, ### (control) build combined topo-cluster/topo-tower container
'TopoClusterRange' : clusterRange, ### (control) range for topo-cluster in combined mode
'TopoClusterRange' : clusterRange, ### (control) range for topo-cluster in combined mode
'RemoveSamplingData' : removeSamplingData, ### (control) remove all sampling data from tower
}
return configDict
......@@ -146,15 +148,15 @@ def MakeTowersFromClusters(towerMakerName = 'CaloTowerBuilderAlg', #
clusterMoments.TwoGaussianNoise = jobproperties.CaloTopoClusterFlags.doTwoGaussianNoise()
clusterMoments.MinBadLArQuality = 4000
clusterMoments.MomentsNames = [
"CENTER_LAMBDA",
#"CENTER_MAG",
# "CENTER_LAMBDA",
# "CENTER_MAG",
"LONGITUDINAL",
#"FIRST_ENG_DENS",
#"ENG_FRAC_MAX",
# "FIRST_ENG_DENS",
# "ENG_FRAC_MAX",
"ENG_FRAC_EM",
#"PTD",
# "PTD",
"SIGNIFICANCE",
"ENG_POS"
# "ENG_POS",
]
from IOVDbSvc.CondDB import conddb
......
......@@ -58,7 +58,7 @@ caloTowerMerger.TopoClusterRange = caloTowerDict['TopoClusterRange']
caloTowerMerger.TopoClusterContainerKey = caloTowerDict['CaloTopoClusterContainerKey'] #### caloTowerAlgo.CaloFwdTopoTowerBuilder.CaloTopoClusterContainerKey
caloTowerMerger.TopoTowerContainerKey = caloTowerAlgo.TowersOutputName
caloTowerMerger.TopoSignalContainerKey = 'CaloCalTopoSignals'
caloTowerMerger.OutputLevel = Lvl.DEBUG
## caloTowerMerger.OutputLevel = Lvl.DEBUG
topSequence+=caloTowerAlgo
topSequence+=caloTowerMerger
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
//-----------------------------------------------------------------------
......@@ -12,6 +12,7 @@
//
// Author List:
// Sven Menke
// Peter Loch
//
//-----------------------------------------------------------------------
......@@ -39,78 +40,98 @@
#include <limits>
#include <sstream>
#include <map>
#include <string>
#include <cstdio>
using CLHEP::deg;
using CLHEP::cm;
// Known moments
namespace {
//FIXME, somehow make sure these names are in sync with the xAOD variable names
struct MomentName
{
const char* name;
xAOD::CaloCluster::MomentType mom;
};
// Must be sorted by name.
MomentName moment_names[] = {
{ "AVG_LAR_Q", xAOD::CaloCluster::AVG_LAR_Q },
{ "AVG_TILE_Q", xAOD::CaloCluster::AVG_TILE_Q },
{ "BADLARQ_FRAC", xAOD::CaloCluster::BADLARQ_FRAC },
{ "BAD_CELLS_CORR_E", xAOD::CaloCluster::BAD_CELLS_CORR_E },
{ "CELL_SIGNIFICANCE", xAOD::CaloCluster::CELL_SIGNIFICANCE },
{ "CELL_SIG_SAMPLING", xAOD::CaloCluster::CELL_SIG_SAMPLING },
{ "CENTER_LAMBDA", xAOD::CaloCluster::CENTER_LAMBDA },
{ "CENTER_MAG", xAOD::CaloCluster::CENTER_MAG },
{ "CENTER_X", xAOD::CaloCluster::CENTER_X },
{ "CENTER_Y", xAOD::CaloCluster::CENTER_Y },
{ "CENTER_Z", xAOD::CaloCluster::CENTER_Z },
{ "DELTA_ALPHA", xAOD::CaloCluster::DELTA_ALPHA },
{ "DELTA_PHI", xAOD::CaloCluster::DELTA_PHI },
{ "DELTA_THETA", xAOD::CaloCluster::DELTA_THETA },
{ "ENG_BAD_CELLS", xAOD::CaloCluster::ENG_BAD_CELLS },
{ "ENG_BAD_HV_CELLS", xAOD::CaloCluster::ENG_BAD_HV_CELLS },
{ "ENG_FRAC_CORE", xAOD::CaloCluster::ENG_FRAC_CORE },
{ "ENG_FRAC_EM", xAOD::CaloCluster::ENG_FRAC_EM },
{ "ENG_FRAC_MAX", xAOD::CaloCluster::ENG_FRAC_MAX },
{ "ENG_POS", xAOD::CaloCluster::ENG_POS },
{ "FIRST_ENG_DENS", xAOD::CaloCluster::FIRST_ENG_DENS },
{ "FIRST_ETA", xAOD::CaloCluster::FIRST_ETA },
{ "FIRST_PHI", xAOD::CaloCluster::FIRST_PHI },
{ "ISOLATION", xAOD::CaloCluster::ISOLATION },
{ "LATERAL", xAOD::CaloCluster::LATERAL },
{ "LONGITUDINAL", xAOD::CaloCluster::LONGITUDINAL },
{ "MASS", xAOD::CaloCluster::MASS },
{ "N_BAD_CELLS", xAOD::CaloCluster::N_BAD_CELLS },
{ "N_BAD_HV_CELLS", xAOD::CaloCluster::N_BAD_HV_CELLS },
{ "N_BAD_CELLS_CORR", xAOD::CaloCluster::N_BAD_CELLS_CORR },
{ "PTD", xAOD::CaloCluster::PTD },
{ "SECOND_ENG_DENS", xAOD::CaloCluster::SECOND_ENG_DENS },
{ "SECOND_LAMBDA", xAOD::CaloCluster::SECOND_LAMBDA },
{ "SECOND_R", xAOD::CaloCluster::SECOND_R },
{ "SIGNIFICANCE", xAOD::CaloCluster::SIGNIFICANCE },
};
MomentName* moment_names_end =
moment_names + sizeof(moment_names)/sizeof(moment_names[0]);
#if 0
bool operator< (const std::string& v, const MomentName& m)
{
return strcmp (v.c_str(), m.name) < 0;
}
#endif
bool operator< (const MomentName& m, const std::string& v)
{
return strcmp (m.name, v.c_str()) < 0;
}
// name -> enum translator
std::map<std::string,xAOD::CaloCluster::MomentType> momentNameToEnumMap = {
{ "AVG_LAR_Q", xAOD::CaloCluster::AVG_LAR_Q },
{ "AVG_TILE_Q", xAOD::CaloCluster::AVG_TILE_Q },
{ "BADLARQ_FRAC", xAOD::CaloCluster::BADLARQ_FRAC },
{ "BAD_CELLS_CORR_E", xAOD::CaloCluster::BAD_CELLS_CORR_E },
{ "CELL_SIGNIFICANCE", xAOD::CaloCluster::CELL_SIGNIFICANCE },
{ "CELL_SIG_SAMPLING", xAOD::CaloCluster::CELL_SIG_SAMPLING },
{ "CENTER_LAMBDA", xAOD::CaloCluster::CENTER_LAMBDA },
{ "CENTER_MAG", xAOD::CaloCluster::CENTER_MAG },
{ "CENTER_X", xAOD::CaloCluster::CENTER_X },
{ "CENTER_Y", xAOD::CaloCluster::CENTER_Y },
{ "CENTER_Z", xAOD::CaloCluster::CENTER_Z },
{ "DELTA_ALPHA", xAOD::CaloCluster::DELTA_ALPHA },
{ "DELTA_PHI", xAOD::CaloCluster::DELTA_PHI },
{ "DELTA_THETA", xAOD::CaloCluster::DELTA_THETA },
{ "ENG_BAD_CELLS", xAOD::CaloCluster::ENG_BAD_CELLS },
{ "ENG_BAD_HV_CELLS", xAOD::CaloCluster::ENG_BAD_HV_CELLS },
{ "ENG_FRAC_CORE", xAOD::CaloCluster::ENG_FRAC_CORE },
{ "ENG_FRAC_EM", xAOD::CaloCluster::ENG_FRAC_EM },
{ "ENG_FRAC_MAX", xAOD::CaloCluster::ENG_FRAC_MAX },
{ "ENG_POS", xAOD::CaloCluster::ENG_POS },
{ "FIRST_ENG_DENS", xAOD::CaloCluster::FIRST_ENG_DENS },
{ "FIRST_ETA", xAOD::CaloCluster::FIRST_ETA },
{ "FIRST_PHI", xAOD::CaloCluster::FIRST_PHI },
{ "ISOLATION", xAOD::CaloCluster::ISOLATION },
{ "LATERAL", xAOD::CaloCluster::LATERAL },
{ "LONGITUDINAL", xAOD::CaloCluster::LONGITUDINAL },
{ "MASS", xAOD::CaloCluster::MASS },
{ "N_BAD_CELLS", xAOD::CaloCluster::N_BAD_CELLS },
{ "N_BAD_HV_CELLS", xAOD::CaloCluster::N_BAD_HV_CELLS },
{ "N_BAD_CELLS_CORR", xAOD::CaloCluster::N_BAD_CELLS_CORR },
{ "PTD", xAOD::CaloCluster::PTD },
{ "SECOND_ENG_DENS", xAOD::CaloCluster::SECOND_ENG_DENS },
{ "SECOND_LAMBDA", xAOD::CaloCluster::SECOND_LAMBDA },
{ "SECOND_R", xAOD::CaloCluster::SECOND_R },
{ "SECOND_TIME", xAOD::CaloCluster::SECOND_TIME },
{ "SIGNIFICANCE", xAOD::CaloCluster::SIGNIFICANCE },
{ "EM_PROBABILITY", xAOD::CaloCluster::EM_PROBABILITY }
};
// enum -> name translator
std::map<xAOD::CaloCluster::MomentType,std::string> momentEnumToNameMap = {
{ xAOD::CaloCluster::AVG_LAR_Q, "AVG_LAR_Q" },
{ xAOD::CaloCluster::AVG_TILE_Q, "AVG_TILE_Q" },
{ xAOD::CaloCluster::BADLARQ_FRAC, "BADLARQ_FRAC" },
{ xAOD::CaloCluster::BAD_CELLS_CORR_E, "BAD_CELLS_CORR_E" },
{ xAOD::CaloCluster::CELL_SIGNIFICANCE, "CELL_SIGNIFICANCE"},
{ xAOD::CaloCluster::CELL_SIG_SAMPLING, "CELL_SIG_SAMPLING"},
{ xAOD::CaloCluster::CENTER_LAMBDA, "CENTER_LAMBDA" },
{ xAOD::CaloCluster::CENTER_MAG, "CENTER_MAG" },
{ xAOD::CaloCluster::CENTER_X, "CENTER_X" },
{ xAOD::CaloCluster::CENTER_Y, "CENTER_Y" },
{ xAOD::CaloCluster::CENTER_Z, "CENTER_Z" },
{ xAOD::CaloCluster::DELTA_ALPHA, "DELTA_ALPHA" },
{ xAOD::CaloCluster::DELTA_PHI, "DELTA_PHI" },
{ xAOD::CaloCluster::DELTA_THETA, "DELTA_THETA" },
{ xAOD::CaloCluster::ENG_BAD_CELLS, "ENG_BAD_CELLS" },
{ xAOD::CaloCluster::ENG_BAD_HV_CELLS, "ENG_BAD_HV_CELLS" },
{ xAOD::CaloCluster::ENG_FRAC_CORE, "ENG_FRAC_CORE" },
{ xAOD::CaloCluster::ENG_FRAC_EM, "ENG_FRAC_EM" },
{ xAOD::CaloCluster::ENG_FRAC_MAX, "ENG_FRAC_MAX" },
{ xAOD::CaloCluster::ENG_POS, "ENG_POS" },
{ xAOD::CaloCluster::FIRST_ENG_DENS, "FIRST_ENG_DENS" },
{ xAOD::CaloCluster::FIRST_ETA, "FIRST_ETA" },
{ xAOD::CaloCluster::FIRST_PHI, "FIRST_PHI" },
{ xAOD::CaloCluster::ISOLATION, "ISOLATION" },
{ xAOD::CaloCluster::LATERAL, "LATERAL" },
{ xAOD::CaloCluster::LONGITUDINAL, "LONGITUDINAL" },
{ xAOD::CaloCluster::MASS, "MASS" },
{ xAOD::CaloCluster::N_BAD_CELLS, "N_BAD_CELLS" },
{ xAOD::CaloCluster::N_BAD_HV_CELLS, "N_BAD_HV_CELLS" },
{ xAOD::CaloCluster::N_BAD_CELLS_CORR, "N_BAD_CELLS_CORR" },
{ xAOD::CaloCluster::PTD, "PTD" },
{ xAOD::CaloCluster::SECOND_ENG_DENS, "SECOND_ENG_DENS" },
{ xAOD::CaloCluster::SECOND_LAMBDA, "SECOND_LAMBDA" },
{ xAOD::CaloCluster::SECOND_R, "SECOND_R" },
{ xAOD::CaloCluster::SECOND_TIME, "SECOND_TIME" },
{ xAOD::CaloCluster::SIGNIFICANCE, "SIGNIFICANCE" },
{ xAOD::CaloCluster::EM_PROBABILITY, "EM_PROBABILITY" }
};
}
//###############################################################################
CaloClusterMomentsMaker::CaloClusterMomentsMaker(const std::string& type,
......@@ -134,117 +155,102 @@ CaloClusterMomentsMaker::CaloClusterMomentsMaker(const std::string& type,
// Name(s) of Moments to calculate
declareProperty("MomentsNames",m_momentsNames);
// Name(s) of Moments which can be stored on the AOD - all others go to ESD
// m_momentsNamesAOD.push_back(std::string("FIRST_PHI"));
// m_momentsNamesAOD.push_back(std::string("FIRST_ETA"));
// m_momentsNamesAOD.push_back(std::string("SECOND_R"));
// m_momentsNamesAOD.push_back(std::string("SECOND_LAMBDA"));
// m_momentsNamesAOD.push_back(std::string("CENTER_LAMBDA"));
// m_momentsNamesAOD.push_back(std::string("FIRST_ENG_DENS"));
// m_momentsNamesAOD.push_back(std::string("ENG_BAD_CELLS"));
// m_momentsNamesAOD.push_back(std::string("N_BAD_CELLS"));
//declareProperty("AODMomentsNames",m_momentsNamesAOD);
// maximum allowed angle between shower axis and the vector pointing
// to the shower center from the IP in degrees. This property is need
// Maximum allowed angle between shower axis and the vector pointing
// to the shower center from the IP in degrees. This property is needed
// to protect against cases where all significant cells are in one sampling
// and the shower axis can not be defined from them
// and the shower axis can thus not be defined.
declareProperty("MaxAxisAngle",m_maxAxisAngle);
declareProperty("MinRLateral",m_minRLateral);
declareProperty("MinLLongitudinal",m_minLLongitudinal);
declareProperty("MinBadLArQuality",m_minBadLArQuality);
// use 2-gaussian noise for Tile
// Use 2-gaussian noise for Tile
declareProperty("TwoGaussianNoise",m_twoGaussianNoise);
declareProperty("LArHVFraction",m_larHVFraction,"Tool Handle for LArHVFraction");
/// Not used anymore (with xAOD), but required when configured from COOL.
// Not used anymore (with xAOD), but required when configured from COOL.
declareProperty("AODMomentsNames",m_momentsNamesAOD);
// Use weighting of neg. clusters option?
declareProperty("WeightingOfNegClusters", m_absOpt);
}
//###############################################################################
StatusCode CaloClusterMomentsMaker::initialize()
{
m_calculateSignificance = false;
m_calculateIsolation = false;
// translate all moment names specified in MomentsNames property to moment enums,
// check that they are all valid and there are no repeating names
for(const auto& name: m_momentsNames) {
const MomentName* it =
std::lower_bound (moment_names, moment_names_end, name);
if (it != moment_names_end) {
m_validMoments.push_back (it->mom);
switch (it->mom) {
case xAOD::CaloCluster::SIGNIFICANCE:
case xAOD::CaloCluster::CELL_SIGNIFICANCE:
m_calculateSignificance = true;
break;
case xAOD::CaloCluster::ISOLATION:
m_calculateIsolation = true;
break;
case xAOD::CaloCluster::ENG_BAD_HV_CELLS:
m_calculateLArHVFraction = true;
default:
break;
}
}
else {
msg(MSG::ERROR) << "Moment " << name
<< " is not a valid Moment name and will be ignored! "
<< "Valid names are:";
int count = 0;
for (const MomentName& m : moment_names)
msg() << ((count++)==0?" ":", ") << m.name;
msg() << endmsg;
}
}
// loop list of requested moments
std::string::size_type nstr(0); int nmom(0);
for ( const auto& mom : m_momentsNames ) {
// check if moment is known (enumerator available)
auto fmap(momentNameToEnumMap.find(mom));
if ( fmap != momentNameToEnumMap.end() ) {
// valid moment found
nstr = std::max(nstr,mom.length()); ++nmom;
if ( fmap->second == xAOD::CaloCluster::SECOND_TIME ) {
// special flag for second moment of cell times - this moment is not
// calculated in this tool! Do not add to internal (!) valid moments list.
// Its value is available from xAOD::CaloCluster::secondTime()!
m_secondTime = true;
} else if ( fmap->second == xAOD::CaloCluster::EM_PROBABILITY ) {
ATH_MSG_WARNING( mom << " not calculated in this tool - misconfiguration?" );
} else {
// all other valid moments
m_validMoments.push_back(fmap->second);
// flag some special requests
switch (fmap->second) {
case xAOD::CaloCluster::SIGNIFICANCE:
case xAOD::CaloCluster::CELL_SIGNIFICANCE:
m_calculateSignificance = true;
break;
case xAOD::CaloCluster::ISOLATION:
m_calculateIsolation = true;
break;
case xAOD::CaloCluster::ENG_BAD_HV_CELLS:
m_calculateLArHVFraction = true;
default:
break;
} // set special processing flags
} // moment calculated with this tool
} else {
ATH_MSG_ERROR( "Moment name " << mom << " not known; known moments are:" );
char buffer[128]; std::string::size_type lstr(nstr);
// determine field size
for ( auto fmom : momentNameToEnumMap ) { lstr = std::max(lstr,fmom.first.length()); }
// print available moments
for ( auto fmom : momentNameToEnumMap ) {
sprintf(buffer,"moment name: %-*.*s - enumerator: %i",(int)lstr,(int)lstr,fmom.first.c_str(),(int)fmom.second);
ATH_MSG_INFO(buffer);
}
return StatusCode::FAILURE;
} // found unknown moment name
} // loop configured moment names
// sort and remove duplicates, order is not required for any of the code below
// but still may be useful property
// sort and remove duplicates
std::sort(m_validMoments.begin(), m_validMoments.end());
m_validMoments.erase(std::unique(m_validMoments.begin(),
m_validMoments.end()),
m_validMoments.end());
/*
// translate moment names in AODMomentsNames property into set of enums,
// only take valid names which are also in MomentsNames property
m_momentsAOD.reserve(m_momentsNamesAOD.size());
for(const auto& name: m_momentsNamesAOD) {
const MomentName* it =
std::lower_bound (moment_names, moment_names_end, name);
if (it != moment_names_end) {
if (std::find(m_validMoments.begin(), m_validMoments.end(), it->mom)
!= m_validMoments.end())
{
m_momentsAOD.push_back(it->mom);
}
}
m_validMoments.erase(std::unique(m_validMoments.begin(),m_validMoments.end()),m_validMoments.end());
// print configured moments
ATH_MSG_INFO( "Construct and save " << nmom << " cluster moments: " );
char buffer[128];
for ( auto menum : m_validMoments ) {
sprintf(buffer,"moment name: %-*.*s - enumerator: %i",(int)nstr,(int)nstr,momentEnumToNameMap.at(menum).c_str(),(int)menum);
ATH_MSG_INFO( buffer );
}
*/
if ( m_secondTime ) {
auto fmom(momentNameToEnumMap.find("SECOND_TIME"));
sprintf(buffer,"moment name: %-*.*s - enumerator: %i (save only)",(int)nstr,(int)nstr,fmom->first.c_str(),(int)fmom->second);
ATH_MSG_INFO( buffer );
}
// retrieve CaloCell ID server
CHECK(detStore()->retrieve(m_calo_id,"CaloCell_ID"));
// retrieve the calo depth tool
CHECK(m_caloDepthTool.retrieve());
if (m_calculateSignificance) {
ATH_CHECK(m_noiseCDOKey.initialize());
}
// retrieve specific servers and tools for selected processes
if (m_calculateSignificance) { ATH_CHECK(m_noiseCDOKey.initialize()); }
if (m_calculateLArHVFraction) { ATH_CHECK(m_larHVFraction.retrieve()); } else { m_larHVFraction.disable(); }
if (m_calculateLArHVFraction) {
ATH_CHECK(m_larHVFraction.retrieve());
}
else {
m_larHVFraction.disable();
}
return StatusCode::SUCCESS;
}
StatusCode CaloClusterMomentsMaker::finalize()
......@@ -256,6 +262,7 @@ StatusCode CaloClusterMomentsMaker::finalize()
namespace CaloClusterMomentsMaker_detail {
struct cellinfo {
double x;
double y;
......@@ -903,7 +910,7 @@ CaloClusterMomentsMaker::execute(const EventContext& ctx,
}
}
}
// normalize moments and copy to Cluster Moment Store
size_t size= m_validMoments.size();
for (size_t iMoment = 0; iMoment != size; ++iMoment) {
......@@ -913,9 +920,11 @@ CaloClusterMomentsMaker::execute(const EventContext& ctx,
if ( moment == xAOD::CaloCluster::FIRST_PHI )
myMoments[iMoment] = CaloPhiRange::fix(myMoments[iMoment]);
theCluster->insertMoment(moment,myMoments[iMoment]);
}
}
}
} // loop on moments for cluster
} // check on requested moments
// check on second moment of time if requested
if ( m_secondTime ) { theCluster->insertMoment(xAOD::CaloCluster::SECOND_TIME,theCluster->secondTime()); }
} // loop on clusters
return StatusCode::SUCCESS;
}
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
//Dear emacs, this is -*-c++-*-
......@@ -8,7 +8,8 @@
/**
* @class CaloClusterMomentsMaker
* @author Sven Menke <menke@mppmu.mpg.de>
* @date 28-February-2005
* @author Peter Loch <loch@physics.arizona.edu>
* @date 23-March-2021
* @brief Calculate moments for CaloCluster objects
*
* This is a CaloClusterCollectionProcessor which can be plugged into a
......@@ -25,7 +26,10 @@
* Note that only cells with positive energy are used in this definition.
* Common variables to calculate first and second moments of are
* \f$\phi\f$, \f$\eta\f$, and radial and longitudinal distances from
* the shower axis and the shower center, respectively. */
* the shower axis and the shower center, respectively.
*
* @since 23-March-2021: second moment of cell time distribution is calculated
*/
#include "GaudiKernel/ToolHandle.h"
......@@ -133,6 +137,10 @@ class CaloClusterMomentsMaker: public AthAlgTool, virtual public CaloClusterColl
* @brief if set to true use abs E value of cells to calculate
* cluster moments */
bool m_absOpt;
/**
* @brief Retreive second moment of cell times */
bool m_secondTime = { false };
};
#endif // CALOCLUSTERMOMENTSMAKER_H
......@@ -76,6 +76,7 @@ CaloTopoTowerFromClusterMaker::CaloTopoTowerFromClusterMaker(const std::string&
declareProperty("DoCellIndexCheck", m_doCellIndexCheck, "Check cell hash indices for consistency");
declareProperty("BuildCombinedTopoSignal", m_buildCombinedSignal, "Build topo-clusters and topo-towers");
declareProperty("TopoClusterRange", m_clusterRange, "Rapidity range for using topo-clusters in combined signal mode");
declareProperty("RemoveSamplingData", m_removeSamplingData, "Remove the associated sampling data");
}
StatusCode CaloTopoTowerFromClusterMaker::initialize()
......@@ -183,6 +184,7 @@ StatusCode CaloTopoTowerFromClusterMaker::initialize()
ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("BuildCombinedTopoSignal .... %s", blu[m_buildCombinedSignal].c_str()) );
ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("TopoClusterRange ........... %.2f", m_clusterRange) );
ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("ExcludedSamplings .......... %zu (number of)",m_excludedSamplingsName.size()) );
ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("RemoveSamplingData ......... %s", blu[m_removeSamplingData].c_str()) );
return StatusCode::SUCCESS;
}
......@@ -268,8 +270,11 @@ StatusCode CaloTopoTowerFromClusterMaker::execute(const EventContext& ctx,
clptr->addCellLink(lptr); // transfer cell links to CaloCluster
clptr->setClusterSize(csize); // set the cluster size spec
CaloRec::Helpers::calculateKine(clptr,false); // calculate kinematics and other signals from cells
clptr->setEta0(m_towerGeometrySvc->towerEta(ipc)); // save the tower center eta
clptr->setPhi0(m_towerGeometrySvc->towerPhi(ipc)); // save the tower center phi
if ( m_removeSamplingData ) { // remove sampling data and invalidate tower center
clptr->clearSamplingData(); clptr->setEta0(0.); clptr->setPhi0(0.);
} else { // keep sampling data and valid tower center
clptr->setEta0(m_towerGeometrySvc->towerEta(ipc)); clptr->setPhi0(m_towerGeometrySvc->towerPhi(ipc));
}
} else {
delete lptr;
}
......
......@@ -69,6 +69,7 @@ private:
bool m_buildCombinedSignal = { false }; ///< Build topo-clusters within given @f$ y @f$ range, else topo-towers
double m_energyThreshold; ///< Cell energy threshold, default is set in @c m_energyThresholdDef
double m_clusterRange; ///< Range where topo-clusters are used when <tt>m_buildCombinedSignal = true</tt>
bool m_removeSamplingData = { true }; ///< Remove sampling data for towers
/// @}
/// @name Constants and parameters
......
......@@ -56,6 +56,7 @@ namespace xAOD {
DEFINE_ACCESSOR( SIGNIFICANCE );
DEFINE_ACCESSOR( PTD );
DEFINE_ACCESSOR( MASS );
DEFINE_ACCESSOR( SECOND_TIME );
DEFINE_ACCESSOR( CELL_SIGNIFICANCE );
DEFINE_ACCESSOR( CELL_SIG_SAMPLING );
DEFINE_ACCESSOR( AVG_LAR_Q );
......
......@@ -33,7 +33,8 @@ namespace xAOD {
: IParticle(other),
m_samplingPattern(other.samplingPattern()),
m_cellLinks(nullptr),
m_recoStatus(other.m_recoStatus) {
m_recoStatus(other.m_recoStatus),
m_secondTime(other.m_secondTime) {
setSignalState(other.signalState());
this->makePrivateStore(other);
#if !(defined(SIMULATIONBASE) || defined(XAOD_ANALYSIS))
......@@ -58,6 +59,7 @@ namespace xAOD {
m_recoStatus=other.m_recoStatus;
setSignalState(other.signalState());
m_samplingPattern=other.m_samplingPattern;
m_secondTime = other.m_secondTime;
#if !(defined(SIMULATIONBASE) || defined(XAOD_ANALYSIS))
const CaloClusterCellLink* links=other.getCellLinks();
......@@ -922,6 +924,15 @@ namespace xAOD {
return true;
}
void CaloCluster_v1::setSecondTime(CaloCluster_v1::flt_t stime) { m_secondTime = stime; }
CaloCluster_v1::flt_t CaloCluster_v1::secondTime() const {
if ( m_secondTime < 0. ) {
double stime(0.); return this->retrieveMoment(SECOND_TIME,stime) ? stime : 0.;
} else {
return m_secondTime;
}
}
#if !(defined(SIMULATIONBASE) || defined(XAOD_ANALYSIS))
size_t CaloCluster_v1::size() const {
......
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#ifndef SIMULATIONBASE
......@@ -12,6 +12,8 @@
#include "CaloGeoHelpers/CaloPhiRange.h"
#include "CaloEvent/CaloPrefetch.h"
#include <cmath>
namespace {
......@@ -37,6 +39,7 @@ struct CellAccum
PresenceInSample(),
theNewEnergy(0),
theNewTime(0),
theNewSecondTime(0),
theNewEta(0),
theNewPhi(0),
phi0(-999),
......@@ -55,6 +58,7 @@ struct CellAccum
bool PresenceInSample[CaloSampling::Unknown];
double theNewEnergy;
double theNewTime;
double theNewSecondTime;
double theNewEta;
double theNewPhi;
double phi0;
......@@ -160,9 +164,10 @@ void accumCell (const CaloClusterCellLink::const_iterator& cellIt, CellAccum& ac
// Is time defined?
if ( cell.provenance() & pmask ) {
// keep the sign of weight for the time norm in case a cell is removed
double theTimeNorm = fabsweight * theEnergy * theEnergyNoWeight;
accum.theNewTime += theTimeNorm * cell.time();
accum.timeNorm += theTimeNorm;
double theTimeNorm = fabsweight * theEnergy * theEnergyNoWeight;
accum.theNewTime += theTimeNorm * cell.time();
accum.theNewSecondTime += theTimeNorm * cell.time() * cell.time();
accum.timeNorm += theTimeNorm;
}
}
}
......@@ -251,10 +256,17 @@ void CaloClusterKineHelper::calculateKine(xAOD::CaloCluster* clu, const bool use
clu->setE(accum.theNewEnergy);
clu->setM(0.0);
if ( timeNorm != 0. )
if ( timeNorm != 0. ) {
clu->setTime(accum.theNewTime/timeNorm);
else
clu->setTime(0);
// note that standard deviation of cell time distribution is >= 0.
// (no particular accommodation for <x^2> < <x><x>!)
// clu->setSecondTime(std::sqrt(std::max(accum.theNewSecondTime/timeNorm-clu->time()*clu->time(),0.)));
// consistency with other second moments: store variance!
clu->setSecondTime(accum.theNewSecondTime/timeNorm-clu->time()*clu->time());
} else {
clu->setTime(0.);
clu->setSecondTime(0.);
}
//Set Sampling pattern:
......
......@@ -6,6 +6,7 @@
<field name="m_signalState" transient="true" />
<field name="m_cellLinks" transient="true" />
<field name="m_recoStatus" transient="true" />
<field name="m_secondTime" transient="true" />
</class>
<read sourceClass="xAOD::CaloCluster_v1" version="[1-]"
targetClass="xAOD::CaloCluster_v1" source="" target="m_signalState" >
......@@ -19,6 +20,11 @@
m_cellLinks = 0;
]]>
</read>
<read sourceClass="xAOD::CaloCluster_v1" version="[1-]"
targetClass="xAOD::CaloCluster_v1" source="" target="m_secondTime" >
<![CDATA[
m_secondTime = -999;
]]>
<class name="xAOD::CaloClusterContainer_v1"
id="24997BCA-3F6A-4530-8826-822EE9FC3C08" />
<typedef name="xAOD::CaloCluster" />
......
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#ifndef XAODCALOEVENT_VERSIONS_CALOCLUSTER_V1_H
#define XAODCALOEVENT_VERSIONS_CALOCLUSTER_V1_H
......@@ -14,7 +14,7 @@ extern "C" {
// xAOD include(s):
#include "xAODBase/IParticle.h"
//Defintion of CaloSamples (enum)
//Definition of CaloSamples (enum)
#include "CaloGeoHelpers/CaloSampling.h"
#include "xAODCaloEvent/CaloClusterBadChannelData.h"
......@@ -45,8 +45,11 @@ namespace xAOD {
/// Description of a calorimeter cluster
/// @author Attila Krasznahorkay <Attila.Krasznahorkay@cern.ch>
/// @author Walter Lampl <Walter.Lampl@cern.ch>
/// @author Peter Loch <Peter.Loch@cern.ch>
///
///
/// @since 23-March-2021: added methods to set and retrieve second
/// moment of cell time distribution (persistified
/// as a cluster moment)
class CaloCluster_v1 : public IParticle {
friend class ::CaloClusterChangeSignalState;
......@@ -103,8 +106,8 @@ namespace xAOD {
DELTA_PHI = 301,
/// Angular shower axis deviation (\f$\theta\f$) from IP-to-Center
DELTA_THETA = 302,
/// Angular shower axis deviation from IP-to-Center
DELTA_ALPHA = 303,
/// Angular shower axis deviation (\f$\Delta\alpha\f$) from IP-to-Center
DELTA_ALPHA = 303,
CENTER_X = 401, ///< Cluster Centroid (\f$x\f$)
CENTER_Y = 402, ///< Cluster Centroid (\f$y\f$)
CENTER_Z = 403, ///< Cluster Centroid (\f$z\f$)
......@@ -153,6 +156,8 @@ namespace xAOD {
DM_WEIGHT = 903, ///< Dead-material weight (E_dm/E_ooc)
/// Confidence Level of a tile calorimeter cluster to be noise
TILE_CONFIDENCE_LEVEL = 904,
/// Second moment of cell time distribution in cluster
SECOND_TIME = 910,
VERTEX_FRACTION = 1000, /**< Vertex fraction of this cluster wrt. primary vertex of the event. Calculated in CaloRec/CaloClusterVertexFractionMaker.cxx */
NVERTEX_FRACTION = 1001, /**< slightly updated vertex fraction more pile up independent (similar to nJVF) */
......@@ -442,6 +447,16 @@ namespace xAOD {
void setTime(flt_t);
/// Access cluster time
flt_t time() const;
/// Set second moment of cell timing distribution
void setSecondTime(flt_t stime);
/// Access second moment of cell timing distribution
///
/// For clusters read from persistent storage, this method returns the value
/// stored for the @c SECOND_TIME moment.
///
/// @retval 0 if (1) moment is not available, (2) the cluster time could not be calculated,
/// or (3) the cluster has only one cell or all cells have exactly the same time.
flt_t secondTime() const;
/// Access to sampling pattern (one bit per sampling) (Method may be removed later)
unsigned samplingPattern() const;
......@@ -589,13 +604,16 @@ namespace xAOD {
/// Current signal state
State m_signalState;
///Unique ptr to cell links. For cluster building
/// Unique ptr to cell links. For cluster building
/// transient only , holds cells owned by the cluster if non-nullptr
std::unique_ptr<CaloClusterCellLink> m_cellLinks;
///Reco status (transient only)
/// Reco status (transient only)
CaloRecoStatus m_recoStatus;
/// Second cell time moment (transient only)
double m_secondTime = { -1. };
unsigned sampVarIdx(const CaloSample) const;
float getSamplVarFromAcc(const Accessor<std::vector<float > >& acc,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment