Skip to content

GitLab

Commits (89)
Expand all Hide whitespace changes
Inline Side-by-side
Calorimeter/CaloRec/CaloRec/CaloTopoClusterMaker.h
View file @ 32cd521a
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
//Dear emacs, this is -*-c++-*-
......@@ -66,7 +66,7 @@ class CaloTopoClusterMaker: public AthAlgTool, virtual public CaloClusterCollect
private:
inline bool passCellTimeCut(const CaloCell*) const;
inline bool passCellTimeCut(const CaloCell*, float) const;
const CaloDetDescrManager* m_calo_dd_man;
......@@ -154,9 +154,8 @@ class CaloTopoClusterMaker: public AthAlgTool, virtual public CaloClusterCollect
seed level */
float m_seedThresholdOnEtorAbsEt;
/**
* threshold used for timing cut. Implemented as |seed_cell_time|<m_seedThresholdOnTAbs. No such cut on neighbouring cells.*/
* threshold used for timing cut on seed cells. Implemented as |seed_cell_time|<m_seedThresholdOnTAbs. No such cut on neighboring cells.*/
float m_seedThresholdOnTAbs;
......@@ -283,6 +282,11 @@ class CaloTopoClusterMaker: public AthAlgTool, virtual public CaloClusterCollect
*/
bool m_seedCutsInT;
/**
* if set to true, seed cells failing the time cut are also excluded from cluster at all
*/
bool m_cutOOTseed;
/**
* @brief if set to true use 2-gaussian noise description for
......
Calorimeter/CaloRec/python/CaloClusterTopoGetter.py
View file @ 32cd521a
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
#
# $Id: CaloClusterTopoGetter.py,v 1.10 2009-05-19 09:41:18 menke Exp $
......@@ -370,6 +370,8 @@ class CaloClusterTopoGetter ( Configured ) :
TopoMaker.SeedThresholdOnEorAbsEinSigma = 4.0
#timing
TopoMaker.SeedCutsInT = jobproperties.CaloTopoClusterFlags.doTimeCut()
TopoMaker.CutOOTseed = jobproperties.CaloTopoClusterFlags.extendTimeCut() and jobproperties.CaloTopoClusterFlags.doTimeCut()
# note E or AbsE
#
......
Calorimeter/CaloRec/python/CaloTopoClusterFlags.py
View file @ 32cd521a
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
#
# $Id: CaloTopoClusterFlags.py,v 1.5 2009-05-04 16:23:04 lochp Exp $
......@@ -116,6 +116,13 @@ class doTimeCut(JobProperty):
allowedTypes=['bool']
StoredValue=False
class extendTimeCut(JobProperty):
"""
"""
statusOn=True
allowedTypes=['bool']
StoredValue=False
# add the flags container to the top container
jobproperties.add_Container(CaloTopoClusterFlags)
......@@ -138,6 +145,7 @@ list_jobproperties = [
,doTreatEnergyCutAsAbsolute
,doMomentsfromAbs
,doTimeCut
,extendTimeCut
]
for i in list_jobproperties:
......
Calorimeter/CaloRec/src/CaloTopoClusterMaker.cxx
View file @ 32cd521a
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
//-----------------------------------------------------------------------
......@@ -99,6 +99,7 @@ CaloTopoClusterMaker::CaloTopoClusterMaker(const std::string& type,
m_cellCutsInAbsE (true),
m_clusterEtorAbsEtCut ( 0.*MeV),
m_seedCutsInT (false),
m_cutOOTseed (false),
m_twogaussiannoise (false),
m_treatL1PredictedCellsAsGood (true),
m_minSampling (0),
......@@ -147,6 +148,8 @@ CaloTopoClusterMaker::CaloTopoClusterMaker(const std::string& type,
//do Seed cuts on Time
declareProperty("SeedCutsInT",m_seedCutsInT);
//exclude out-of-time seeds from neighbouring and cell stage
declareProperty("CutOOTseed",m_cutOOTseed);
// Noise Sigma
declareProperty("NoiseSigma",m_noiseSigma);
......@@ -523,8 +526,17 @@ StatusCode CaloTopoClusterMaker::execute(xAOD::CaloClusterContainer* clusColl)
bool passedSeedAndTimeCut = (passedSeedCut && (!m_seedCutsInT || passCellTimeCut(pCell))); //time cut is applied here
if ( passedCellCut || passedNeighborCut || passedSeedAndTimeCut ) {
bool passTimeCut_seedCell = (!m_seedCutsInT || passCellTimeCut(pCell,m_seedThresholdOnTAbs));
bool passedSeedAndTimeCut = (passedSeedCut && passTimeCut_seedCell); //time cut is applied here
bool passedNeighborAndTimeCut = passedNeighborCut;
if(m_cutOOTseed && passedSeedCut && !passTimeCut_seedCell) passedNeighborAndTimeCut=false; //exclude Out-Of-Time seeds from neighbouring stage as well (if required)
bool passedCellAndTimeCut = passedCellCut;
if(m_cutOOTseed && passedSeedCut && !passTimeCut_seedCell) passedCellAndTimeCut=false; //exclude Out-Of-Time seeds from cluster (if required)
if ( passedCellAndTimeCut || passedNeighborAndTimeCut || passedSeedAndTimeCut ) {
const CaloDetDescrElement* dde = pCell->caloDDE();
IdentifierHash hashid = dde ? dde->calo_hash() : m_calo_id->calo_cell_hash(pCell->ID());
CaloTopoTmpClusterCell *tmpClusterCell =
......@@ -548,7 +560,7 @@ StatusCode CaloTopoClusterMaker::execute(xAOD::CaloClusterContainer* clusColl)
}
#endif
HashCell hashCell(tmpClusterCell);
if ( passedNeighborCut || passedSeedAndTimeCut ) {
if ( passedNeighborAndTimeCut || passedSeedAndTimeCut ) {
HashCluster *tmpCluster =
new (tmpclus_pool.allocate()) HashCluster (tmplist_pool);
tmpClusterCell->setCaloTopoTmpHashCluster(tmpCluster);
......@@ -813,7 +825,7 @@ void CaloTopoClusterMaker::getClusterSize(){
inline bool CaloTopoClusterMaker::passCellTimeCut(const CaloCell* pCell) const
inline bool CaloTopoClusterMaker::passCellTimeCut(const CaloCell* pCell,float threshold) const
{
// get the cell time to cut on (the same as in CaloEvent/CaloCluster.h)
......@@ -833,7 +845,7 @@ inline bool CaloTopoClusterMaker::passCellTimeCut(const CaloCell* pCell) const
// Is time defined?
if(pCell->provenance() & pmask)
{
return std::abs(pCell->time())<m_seedThresholdOnTAbs;
return std::abs(pCell->time())<threshold;
}
}
}
......
Database/IOVDbTPCnv/IOVDbTPCnv/IOVPayloadContainerCnv_p1.h
View file @ 32cd521a
......@@ -48,6 +48,9 @@ private:
// Attribute types
std::map<std::string, unsigned int> m_attributeTypeMap;
std::vector<std::string> m_attributeTypes;
unsigned int m_objIndexOffset[IOVPayloadContainer_p1::ATTR_TIME_STAMP+1];
};
#endif
Database/IOVDbTPCnv/src/IOVPayloadContainerCnv_p1.cxx
View file @ 32cd521a
......@@ -23,6 +23,10 @@ IOVPayloadContainerCnv_p1::persToTrans(const IOVPayloadContainer_p1* persObj,
<< endreq;
}
for( auto & offset : m_objIndexOffset ) {
offset = 1; // this is just the initial value before first use
}
IOVPayloadContainer::payloadVec& payloadVec = transObj->m_payloadVec;
// Make sure transient container is empty - may be reused
......@@ -135,7 +139,7 @@ IOVPayloadContainerCnv_p1::persToTrans(const IOVPayloadContainer_p1* persObj,
}
}
}
// std::ostringstream attrStr;
// attrList.toOutputStream(attrStr);
......@@ -560,53 +564,85 @@ IOVPayloadContainerCnv_p1::fillAttributeData(const IOVPayloadContainer_p1* persO
coral::AttributeList& attrList,
MsgStream & log)
{
/*
this offset calculation solves the problem reported in
ATR-22116 trying to get past the limitation introduced by
AttrListIndexes::m_objIndex being of type short, which doesn't
cover the full length of the type-wise data vectors
It is assumed that when reading the persistent object entries,
the objIndex starts at 0, and increases by 1 each read for a given type.
One must also not call fillAttributeData multiple times for the
same (AttrListIndexes index), nor call them out of order
* There is one offset per type, initialized to 1
* In the first read of a type the offset is set to 0
* In all later reads the offset is not changed unless the objIndex equals 0
- if objIndex is fixed to be unsigned int, then objIndex should never be 0 in later reads
- if objIndex is unsigned short, then objIndex==0 only when it runs into the 65536 boundary
and the offset then gets increased by this amount
So this will work also when objIndex is integer
*/
unsigned int objIndex = index.objIndex();
unsigned int & offset = m_objIndexOffset[index.typeIndex()];
if(offset == 1) {
offset = 0;
} else {
if(objIndex == 0) {
offset += 65536;
}
}
// Fill persistent object with attribute data - must use switch
// for all possible types
switch (index.typeIndex()) {
case IOVPayloadContainer_p1::ATTR_BOOL:
attrList[name].setValue(persObj->m_bool[index.objIndex()]);
attrList[name].setValue(persObj->m_bool[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_CHAR:
attrList[name].setValue(persObj->m_char[index.objIndex()]);
attrList[name].setValue(persObj->m_char[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_UNSIGNED_CHAR:
attrList[name].setValue(persObj->m_unsignedChar[index.objIndex()]);
attrList[name].setValue(persObj->m_unsignedChar[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_SHORT:
attrList[name].setValue(persObj->m_short[index.objIndex()]);
attrList[name].setValue(persObj->m_short[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_UNSIGNED_SHORT:
attrList[name].setValue(persObj->m_unsignedShort[index.objIndex()]);
attrList[name].setValue(persObj->m_unsignedShort[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_INT:
attrList[name].setValue(persObj->m_int[index.objIndex()]);
attrList[name].setValue(persObj->m_int[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_UNSIGNED_INT:
attrList[name].setValue(persObj->m_unsignedInt[index.objIndex()]);
attrList[name].setValue(persObj->m_unsignedInt[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_LONG:
attrList[name].setValue(persObj->m_long[index.objIndex()]);
attrList[name].setValue(persObj->m_long[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_UNSIGNED_LONG:
attrList[name].setValue(persObj->m_unsignedLong[index.objIndex()]);
attrList[name].setValue(persObj->m_unsignedLong[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_LONG_LONG:
attrList[name].setValue(persObj->m_longLong[index.objIndex()]);
attrList[name].setValue(persObj->m_longLong[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_UNSIGNED_LONG_LONG:
attrList[name].setValue(persObj->m_unsignedLongLong[index.objIndex()]);
attrList[name].setValue(persObj->m_unsignedLongLong[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_FLOAT:
attrList[name].setValue(persObj->m_float[index.objIndex()]);
attrList[name].setValue(persObj->m_float[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_DOUBLE:
attrList[name].setValue(persObj->m_double[index.objIndex()]);
attrList[name].setValue(persObj->m_double[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_LONG_DOUBLE:
// attrList[name].setValue(persObj->m_longDouble[index.objIndex()]);
// attrList[name].setValue(persObj->m_longDouble[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_STRING:
attrList[name].setValue(persObj->m_string[index.objIndex()]);
case IOVPayloadContainer_p1::ATTR_STRING:
attrList[name].setValue(persObj->m_string[objIndex + offset]);
break;
case IOVPayloadContainer_p1::ATTR_BLOB:
log << MSG::ERROR
......@@ -615,14 +651,14 @@ IOVPayloadContainerCnv_p1::fillAttributeData(const IOVPayloadContainer_p1* persO
return;
case IOVPayloadContainer_p1::ATTR_DATE:
{
coral::TimeStamp::ValueType ns( persObj->m_date[index.objIndex()] );
coral::TimeStamp::ValueType ns( persObj->m_date[objIndex + offset] );
attrList[name].setValue( coral::Date(coral::TimeStamp(ns).time()) );
break;
}
case IOVPayloadContainer_p1::ATTR_TIME_STAMP:
{
coral::TimeStamp::ValueType ns =
coral::TimeStamp::ValueType( persObj->m_timeStamp[index.objIndex()] );
coral::TimeStamp::ValueType( persObj->m_timeStamp[objIndex + offset] );
attrList[name].setValue( coral::TimeStamp(ns) );
break;
}
......
Event/EventOverlay/EventOverlayJobTransforms/share/CaloOverlay_jobOptions.py
View file @ 32cd521a
......@@ -79,6 +79,7 @@ if DetFlags.overlay.Tile_on():
ServiceMgr.ByteStreamAddressProviderSvc.TypeNames += [ "TileDigitsContainer/TileDigitsCnt"]
ServiceMgr.ByteStreamAddressProviderSvc.TypeNames += [ "TileL2Container/TileL2Cnt"]
ServiceMgr.ByteStreamAddressProviderSvc.TypeNames += [ "TileLaserObject/TileLaserObj"]
ServiceMgr.ByteStreamCnvSvc.ROD2ROBmap = [ "-1" ]
#--------------------
InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt
View file @ 32cd521a
......@@ -27,6 +27,7 @@ atlas_depends_on_subdirs( PUBLIC
InnerDetector/InDetDetDescr/InDetReadoutGeometry
InnerDetector/InDetDetDescr/SCT_ModuleDistortions
InnerDetector/InDetRawEvent/InDetSimData
MagneticField/MagFieldInterfaces
Simulation/Tools/AtlasCLHEP_RandomGenerators )
# External dependencies:
......@@ -39,7 +40,7 @@ atlas_add_component( SCT_Digitization
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesSvcLib InDetIdentifier InDetReadoutGeometry InDetSimData AtlasCLHEP_RandomGenerators )
LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesSvcLib InDetIdentifier InDetReadoutGeometry InDetSimData AtlasCLHEP_RandomGenerators MagFieldInterfaces PathResolver )
# Install files from the package:
atlas_install_headers( SCT_Digitization )
......
InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.h 0 → 100644
View file @ 32cd521a
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
#define SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
//-----------------------------------------------
// 2020.8.12 Implementation in Athena by Manabu Togawa
// 2017.7.24 update for the case of negative VD (type-P)
// 2017.7.17 updated
// 2016.4.3 for ICM (induced charge model) by Taka Kondo (KEK)
//-----------------------------------------------
// Athena
#include "AthenaKernel/IAtRndmGenSvc.h"
#include "AthenaKernel/MsgStreamMember.h"
#include "Identifier/IdentifierHash.h"
#include "InDetConditionsSummaryService/ISiliconConditionsSvc.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "MagFieldInterfaces/IMagFieldSvc.h"
#include "PathResolver/PathResolver.h"
#include "SiPropertiesSvc/ISiPropertiesSvc.h"
// Gaudi
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/ServiceHandle.h"
// Random number
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGaussZiggurat.h" // for RandGaussZiggurat
#include "CLHEP/Random/RandPoisson.h"
#include "CLHEP/Units/SystemOfUnits.h"
// C++ Standard Library
#include <array>
#include <memory>
#include <mutex>
#include <string>
#include <utility>
#include <vector>
class SCT_InducedChargeModel {
public:
enum EFieldModel {FlatDiodeModel=0, FEMsolutions=1, UniformE=2};
enum TransportStep {NTransportSteps=100};
enum FEMNums {NRamoPoints=81, NEFieldPoints=17, NDepthPoints=115};
enum InducedStrips {StartStrip=-2, EndStrip=+2, Offset=-StartStrip, NStrips=EndStrip+Offset+1};
struct SCT_InducedChargeModelData {
float m_VD; // full depletion voltage [Volt] negative for type-P
float m_VB; // applied bias voltage [Volt]
float m_T; // temperature
float m_depletion_depth;
Amg::Vector3D m_magneticField;
EFieldModel m_EFieldModel;
const InDetDD::SiDetectorElement* m_element;
std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_ExValue;
std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_EyValue;
SCT_InducedChargeModelData(const float vdepl,
const float vbias,
const InDetDD::SiDetectorElement* element,
const Amg::Vector3D& magneticField, // in kTesla
const float bulk_depth,
const EFieldModel model,
const ServiceHandle<ISiliconConditionsSvc> siConditionsSvc) {
m_VD = vdepl; // full depletion voltage [Volt] negative for type-P
m_VB = vbias; // applied bias voltage [Volt]
m_element = element;
m_magneticField = magneticField;
//------------ find delepletion deph for model=0 and 1 -------------
m_depletion_depth = bulk_depth;
// for type N (before type inversion)
if (m_VD >= 0.) {
if (m_VB < m_VD) m_depletion_depth = std::sqrt(m_VB/m_VD) * bulk_depth;
} else {
// for type P
if (m_VB <= std::abs(m_VD)) m_depletion_depth = 0.;
}
m_EFieldModel = model;
m_T = siConditionsSvc->temperature(m_element->identifyHash()) + Gaudi::Units::STP_Temperature;
}
};
SCT_InducedChargeModel(size_t maxHash, EFieldModel model=FEMsolutions);
SCT_InducedChargeModelData*
setWaferData(const float vdepl,
const float vbias,
const InDetDD::SiDetectorElement* element,
const ServiceHandle<MagField::IMagFieldSvc> magFieldSvc,
const ServiceHandle<ISiliconConditionsSvc> siConditionsSvc) const;
void setRandomEngine(CLHEP::HepRandomEngine* rndmEngine) { m_rndmEngine = rndmEngine; };
void setEField(SCT_InducedChargeModelData& data) const;
void transport(const SCT_InducedChargeModelData& data,
const bool isElectron,
const double x0, const double y0,
double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
const IdentifierHash hashId,
const ServiceHandle<ISiPropertiesSvc> siPropertiesSvc) const;
void holeTransport(const SCT_InducedChargeModelData& data,
const double x0, const double y0,
double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
const IdentifierHash hashId,
const ServiceHandle<ISiPropertiesSvc> siPropertiesSvc) const;
void electronTransport(const SCT_InducedChargeModelData& data,
const double x0, const double y0,
double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
const IdentifierHash hashId,
const ServiceHandle<ISiPropertiesSvc> siPropertiesSvc) const;
private:
void loadICMParameters();
bool getVxVyD(const SCT_InducedChargeModelData& data,
const bool isElectron,
const double x, const double y, double& vx, double& vy, double& D,
const IdentifierHash hashId,
const ServiceHandle<ISiPropertiesSvc> siPropertiesSvc) const;
double induced(const SCT_InducedChargeModelData& data,
const int istrip, const double x, const double y) const;
void getEField(const SCT_InducedChargeModelData& data,
const double x, const double y, double& Ex, double& Ey) const;
size_t getFEMIndex(SCT_InducedChargeModelData& data) const;
/// Log a message using the Athena controlled logging system
MsgStream& msg(MSG::Level lvl) const { return m_msg << lvl; }
/// Check whether the logging system is active at the provided verbosity level
bool msgLvl(MSG::Level lvl) const { return m_msg.get().level() <= lvl; }
//-------- parameters for e, h transport --------------------------------
static const double s_kB; // [m^2*kg/s^2/K]
static const double s_e; // [Coulomb]
// Private message stream member
mutable Athena::MsgStreamMember m_msg;
CLHEP::HepRandomEngine* m_rndmEngine; //!< Random number generation engine
//------parameters given externally by jobOptions ------------------
EFieldModel m_EFieldModel; // 0 (flat diode model), 1 (FEM solusions), 2 (uniform E)
double m_transportTimeStep = 0.50; // one step side in time [nsec]
double m_transportTimeMax = m_transportTimeStep*NTransportSteps; // maximun tracing time [nsec]
//------parameters mostly fixed but can be changed externally ------------
double m_bulk_depth = 0.0285; // in [cm]
double m_strip_pitch = 0.0080; // in [cm]
double m_y_origin_min = 0.;
//---------- arrays of FEM analysis -----------------------------------
// Storages for Ramo potential and Electric field.
// In strip pitch directions :
// Ramo potential : 80 divisions (81 points) with 5 um intervals from 40-440 um.
// Electric field : 16 divisions (17 points) with 2.5 um intervals from 0-40 um.
// In sensor depth directions (for both potential and Electric field):
// 114 divisions (115 points) with 2.5 nm intervals for 285 um.
std::vector<std::array<std::array<double, NDepthPoints>, NRamoPoints>> m_PotentialValue;
std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_ExValue;
std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_EyValue;
// Cache of SCT_InducedChargeModelData for each wafer.
// Assuming wafer parameters do not change during a job.
// This should be almost always satisfied.
mutable std::vector<std::unique_ptr<SCT_InducedChargeModelData>> m_data;
// To protect concurrent access to m_data
mutable std::mutex m_mutex;
static const std::vector<float> s_VFD0;
static const std::vector<std::string> s_VFD0str;
};
#endif // SCT_DIGITIZATION_SCTINDUCEDCHARGEDMODEL_H
InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
View file @ 32cd521a
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "SCT_SurfaceChargesGenerator.h"
#include <cmath>
#include "InDetIdentifier/SCT_ID.h"
// DD
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "InDetReadoutGeometry/SCT_ModuleSideDesign.h"
......@@ -91,9 +94,11 @@ parent)
m_siConditionsSvc("SCT_SiliconConditionsSvc", name),
m_siPropertiesSvc("SCT_SiPropertiesSvc", name),
m_radDamageSvc("SCT_RadDamageSummarySvc", name),
m_magFieldSvc("AtlasFieldSvc", name),
m_element(0),
m_rndmEngine(0),
m_rndmEngineName("SCT_Digitization") {
m_rndmEngineName("SCT_Digitization"),
m_doInducedChargeModel(false) {
declareInterface< ISCT_SurfaceChargesGenerator >(this);
declareProperty("FixedTime", m_tfix = -999); // !< fixed
......@@ -107,6 +112,7 @@ parent)
declareProperty("SmallStepLength", m_smallStepLength = 5);
declareProperty("SiConditionsSvc", m_siConditionsSvc);
declareProperty("SiPropertiesSvc", m_siPropertiesSvc);
declareProperty("MagFieldSvc", m_magFieldSvc);
// declareProperty("rndmEngineName",m_rndmEngineName="SCT_Digitization");
declareProperty("doDistortions", m_doDistortions,
"Simulation of module distortions");
......@@ -124,6 +130,8 @@ parent)
"Tool to retrieve SCT distortions");
declareProperty("SCT_RadDamageSummarySvc", m_radDamageSvc);
declareProperty("isOverlay", m_isOverlay=false);
declareProperty("m_doInducedChargeModel", m_doInducedChargeModel,
"Flag for Induced Charge Model"); //
}
// Destructor:
......@@ -152,6 +160,9 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
// Get ISCT_ModuleDistortionsTool
ATH_CHECK(m_distortionsTool.retrieve());
// Get IMagFieldSvc
ATH_CHECK(m_magFieldSvc.retrieve());
if (m_doTrapping) {
///////////////////////////////////////////////////
// -- Get Radiation Damage Service
......@@ -236,6 +247,13 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
}
///////////////////////////////////////////////////
// Induced Charge Module. M.Togawa
if ( m_doInducedChargeModel ){
const SCT_ID* sct_id;
ATH_CHECK(detStore()->retrieve(sct_id, "SCT_ID"));
m_InducedChargeModel = std::make_unique<SCT_InducedChargeModel>(sct_id->wafer_hash_max());
}
m_smallStepLength *= CLHEP::micrometer;
m_tSurfaceDrift *= CLHEP::ns;
......@@ -295,7 +313,7 @@ float SCT_SurfaceChargesGenerator::DriftTime(float zhit) const {
float denominator = m_depletionVoltage + m_biasVoltage - (2.0 * zhit * m_depletionVoltage / m_thickness);
if (denominator <= 0.0) {
if (m_biasVoltage >= m_depletionVoltage) { // Should not happen
if(!m_isOverlay) {
if (!m_isOverlay) {
ATH_MSG_ERROR("DriftTime: negative argument X for log(X) " << zhit);
}
return -1.0;
......@@ -496,6 +514,15 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
float yhit = xPhi;
float zhit = xDep;
SCT_InducedChargeModel::SCT_InducedChargeModelData* data = nullptr;
if ( m_doInducedChargeModel ){ // Setting magnetic field for the ICM.
data = m_InducedChargeModel->setWaferData(m_depletionVoltage/CLHEP::volt,
m_biasVoltage/CLHEP::volt,
m_element,
m_magFieldSvc,
m_siConditionsSvc);
}
if (m_doDistortions) {
if (m_isBarrel) {// Only apply disortions to barrel
// modules
......@@ -583,9 +610,12 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
if (tdrift > 0.0) {
float x1 = xhit + StepX * dstep;// (static_cast<float>(istep)+0.5) ;
float y1 = yhit + StepY * dstep;// (static_cast<float>(istep)+0.5) ;
y1 += m_tanLorentz * zReadout; // !< Taking into account the magnetic
// field
float diffusionSigma = DiffusionSigma(zReadout);
float diffusionSigma = 0;
if ( !m_doInducedChargeModel ){
diffusionSigma = DiffusionSigma(zReadout);
y1 += m_tanLorentz * zReadout; // !< Taking into account the magnetic field
} // Should be treated in Induced Charge Model.
for (int i = 0; i < m_numberOfCharges; ++i) {
float rx = CLHEP::RandGaussZiggurat::shoot(m_rndmEngine);
......@@ -733,6 +763,57 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
} // m_doTrapping==true
if (!m_doRamo) {
if ( m_doInducedChargeModel ){ // Induced Charge Model
// Charges storages for 50 ns. 0.5 ns steps.
double Q_m2[SCT_InducedChargeModel::NTransportSteps]={0};
double Q_m1[SCT_InducedChargeModel::NTransportSteps]={0};
double Q_00[SCT_InducedChargeModel::NTransportSteps]={0};
double Q_p1[SCT_InducedChargeModel::NTransportSteps]={0};
double Q_p2[SCT_InducedChargeModel::NTransportSteps]={0};
const double mm2cm = 0.1; // For mm -> cm conversion
// Unit for y and z : mm -> cm in SCT_InducedChargeModel
m_InducedChargeModel->holeTransport(*data,
y0*mm2cm, z0*mm2cm,
Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
m_hashId, m_siPropertiesSvc);
m_InducedChargeModel->electronTransport(*data,
y0*mm2cm, z0*mm2cm,
Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
m_hashId, m_siPropertiesSvc);
for (int it=0; it<SCT_InducedChargeModel::NTransportSteps; it++){
if ( Q_00[it] == 0.0) continue;
double ICM_time = (it+0.5)*0.5 + timeOfFlight;
double Q_new[SCT_InducedChargeModel::NStrips]={0};
Q_new[0] = Q_m2[it];
Q_new[1] = Q_m1[it];
Q_new[2] = Q_00[it];
Q_new[3] = Q_p1[it];
Q_new[4] = Q_p2[it];
for (int strip = SCT_InducedChargeModel::StartStrip; strip <= SCT_InducedChargeModel::EndStrip; strip++) {
double ystrip = y1 + strip * stripPitch;
SiLocalPosition position(m_element->hitLocalToLocal(x1, ystrip));
if (m_design->inActiveArea(position)) {
inserter(SiSurfaceCharge(position,
SiCharge(q1 * Q_new[strip+SCT_InducedChargeModel::Offset],
ICM_time, hitproc,
trklink)));
}
}
}
} // End Induced Charge Model
if ( !m_doInducedChargeModel ){ // SCT_Digitization
SiLocalPosition position(m_element->hitLocalToLocal(xd,
yd));
if (m_design->inActiveArea(position)) {
......@@ -772,6 +853,9 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiHit &phit, const
// "<<position<<" of the element is out of active area,
// charge = "<<q1) ;
}
} // End of SCT_Digitization
} // end of loop on charges
}
}
......@@ -852,7 +936,7 @@ void SCT_SurfaceChargesGenerator::setVariables() {
m_electronHolePairsPerEnergy = m_siPropertiesSvc->getSiProperties(m_hashId).electronHolePairsPerEnergy();
// get sensor thickness and tg lorentz from SiDetectorDesign
if(m_design) {
if (m_design) {
m_thickness = m_design->thickness();
} else {
m_thickness = 0.;
......
InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h
View file @ 32cd521a
......@@ -35,9 +35,13 @@
#include "SCT_Digitization/ISCT_SurfaceChargesGenerator.h"
#include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
#include "SCT_InducedChargeModel.h"
#include "GaudiKernel/ToolHandle.h"
#include <iostream>
#include "MagFieldInterfaces/IMagFieldSvc.h"
// Charges and hits
class SiHit;
#include "Identifier/IdentifierHash.h"
......@@ -87,7 +91,12 @@ private:
void setFixedTime(float fixedTime) {m_tfix = fixedTime;}
void setCosmicsRun(bool cosmicsRun) {m_cosmicsRun = cosmicsRun;}
void setComTimeFlag(bool useComTime) {m_useComTime = useComTime;}
void setRandomEngine(CLHEP::HepRandomEngine *rndmEngine) {m_rndmEngine = rndmEngine;}
void setRandomEngine(CLHEP::HepRandomEngine *rndmEngine) {
m_rndmEngine = rndmEngine;
if ( m_doInducedChargeModel ){
m_InducedChargeModel->setRandomEngine( m_rndmEngine );
}
}
void setDetectorElement(const InDetDD::SiDetectorElement *ele) {m_element = ele; setVariables();}
/** create a list of surface charges from a hit */
......@@ -162,6 +171,7 @@ private:
ServiceHandle<ISiliconConditionsSvc> m_siConditionsSvc;
ServiceHandle<ISiPropertiesSvc> m_siPropertiesSvc;
ServiceHandle<ISCT_RadDamageSummarySvc> m_radDamageSvc;
ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc;
const InDetDD::SiDetectorElement * m_element;
CLHEP::HepRandomEngine * m_rndmEngine; //!< Random Engine
......@@ -180,6 +190,10 @@ private:
double m_center;
double m_tanLorentz;
bool m_isBarrel;
// For Induced Charge Module, M.Togawa
std::unique_ptr<SCT_InducedChargeModel> m_InducedChargeModel;
bool m_doInducedChargeModel; //!< Flag to use Induced Charge Model
};
#endif // SCT_SURFACECHARGESGENERATOR_H
......
InnerDetector/InDetEventCnv/InDetEventAthenaPool/InDetEventAthenaPool/InDetEventAthenaPoolCnvDict.h
View file @ 32cd521a
/*
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef INDETEVENTATHENAPOOLDICT_H
......@@ -18,6 +18,7 @@
#include "InDetEventAthenaPool/InDetRawDataCollection_p1.h"
#include "InDetEventAthenaPool/InDetRawDataContainer_p1.h"
#include "InDetEventAthenaPool/InDetRawDataContainer_p2.h"
#include "InDetEventAthenaPool/InDetRawDataContainer_p3.h"
#include "InDetEventAthenaPool/SCT_RawDataContainer_p1.h"
#include "InDetEventAthenaPool/SCT_RawDataContainer_p2.h"
#include "InDetEventAthenaPool/SCT_RawDataContainer_p3.h"
......
InnerDetector/InDetEventCnv/InDetEventAthenaPool/InDetEventAthenaPool/InDetRawDataContainer_p3.h 0 → 100644
View file @ 32cd521a
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef INDETRAWDATACONTAINER_P3_H
#define INDETRAWDATACONTAINER_P3_H
/*
For the TRT, the channel Ids are dropped from the persistent container and so
the p2 version "std::vector<InDetRawData_p2> m_rawdata;" is replaced with
the p3 version "std::vector<unsigned int> m_rawdata;".
*/
#include <vector>
#include <string>
#include "InDetEventAthenaPool/InDetRawDataCollection_p1.h"
class InDetRawDataContainer_p3
{
public:
/// Default constructor
InDetRawDataContainer_p3 ();
friend class TRT_LoLumRawDataContainerCnv_p3;
private:
std::vector<InDetRawDataCollection_p1> m_collections;
std::vector<unsigned int> m_rawdata;
};
// inlines
inline
InDetRawDataContainer_p3::InDetRawDataContainer_p3 () {}
#endif
InnerDetector/InDetEventCnv/InDetEventAthenaPool/InDetEventAthenaPool/selection.xml
View file @ 32cd521a
......@@ -29,6 +29,7 @@
<class name="std::vector<SCT3_RawData_p4>" />
<class name="InDetRawDataContainer_p1" id="DA76970C-E019-43D2-B2F9-25660DCECD9D" />
<class name="InDetRawDataContainer_p2" id="7138342E-0A80-4A32-A387-2842A01C2539" />
<class name="InDetRawDataContainer_p3" id="D0313948-9BC8-415F-BE58-7BA8178F93CD" />
<class name="SCT_RawDataContainer_p1" id="8E13963E-13E5-4D10-AA8B-73F00AFF8FA8" />
<class name="SCT_RawDataContainer_p2" id="D1258125-2CBA-476E-8578-E09D54F477E1" />
<class name="SCT_RawDataContainer_p3" id="5FBC8D4D-7B4D-433A-8487-0EA0C870CBDB" />
......
InnerDetector/InDetEventCnv/InDetEventAthenaPool/src/TRT_LoLumRawDataContainerCnv_p3.cxx 0 → 100644
View file @ 32cd521a
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "InDetEventAthenaPool/InDetRawDataCollection_p1.h"
#include "InDetRawData/TRT_RDO_Container.h"
#include "InDetIdentifier/TRT_ID.h"
#include "TRT_LoLumRawDataContainerCnv_p3.h"
#include "TRT_RDO_Elements.h"
#include "MsgUtil.h"
void TRT_LoLumRawDataContainerCnv_p3::transToPers(const TRT_RDO_Container* transCont, InDetRawDataContainer_p3* persCont, MsgStream &log)
{
// The transient model has a container holding collections (up to 14912 TRT detector elements/modules)
// and the collections hold channels (e.g. up to 24 straws in the TRT endcaps). There are a total of 350848 TRT straws.
// The persistent model flattens this so that the persistent container has two vectors: 1) all collections, 2) all RDO
// The persistent collections then only maintain indexes into the container's vector of all channels.
// So here we loop over all collection and add their channels to the container's vector, saving the indexes in the collection.
// Hard-code the number of collections their ids and their sizes to generate collections and channels for the whole detector
const unsigned int dummy_digit=0;
const unsigned int trt_number_of_collections=14912;
const unsigned int trt_number_of_channels=350848;
const unsigned int trt_channel_id_increment=0x20u;
std::vector<unsigned int> trt_collection_id; // There are 14912 elements/modules in the TRT
std::vector<unsigned int> trt_collection_size; // Each element has a specific number of straws
set_vectors_of_collection_ids_and_size(trt_collection_id,
trt_collection_size,
trt_number_of_collections);
unsigned int trt_collection_index=0; // module index (0,14911), initialized first collection
unsigned int trt_channel_index=0; // straw index (0,350847)
unsigned int trt_channel_id; // straw id (32-bit), initialized to the first channel in the first collection
unsigned int tcoll_id; // transient collection id from transCont
unsigned int tchan_id; // transient channel id from transCont
unsigned int tchan_word; // transient channel word from transCont
typedef TRT_RDO_Container TRANS;
TRANS::const_iterator it_transColl = transCont->begin(); // The transient container has an incomplete list of collections
TRANS::const_iterator it_transCollEnd = transCont->end(); // and channels, we can access them with this iterator
persCont->m_collections.resize(trt_number_of_collections);
persCont->m_rawdata.resize(trt_number_of_channels);
// Loop over all existing transient collections, add missing collections and missing channels
for (; it_transColl != it_transCollEnd; it_transColl++) {
const TRT_RDO_Collection& collection = (**it_transColl);
tcoll_id = collection.identify().get_identifier32().get_compact();
// create collections and channels until we find a match to a transient container
while ( trt_collection_id[trt_collection_index] != tcoll_id ) {
// The transient collection is missing; create the persistent collection with empty digits
trt_channel_id = trt_collection_id[trt_collection_index];
InDetRawDataCollection_p1& pcollection = persCont->m_collections[trt_collection_index];
pcollection.m_id = 0;
pcollection.m_hashId = 0; // no hashId for this empty collection
pcollection.m_begin = 0;
pcollection.m_end = 0;
// Fill all the channels with with empty digits.
for (unsigned int i=0; i<trt_collection_size[trt_collection_index]; ++i) {
persCont->m_rawdata[trt_channel_index] = dummy_digit;
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
}
trt_collection_index++;
}
// Here we have a match to a transient collection
// create the persistent collection just as we did above
trt_channel_id = trt_collection_id[trt_collection_index];
InDetRawDataCollection_p1& pcollection = persCont->m_collections[trt_collection_index];
pcollection.m_id = 0;
pcollection.m_hashId = (unsigned int) collection.identifyHash(); // Required by the overlay alogrithm (not the reco algorithms)
pcollection.m_begin = 0;
pcollection.m_end = 0;
const unsigned int collection_end = trt_channel_index+trt_collection_size[trt_collection_index];
// Fill all the channels; some with existing digits, the others with empty digits.
for (unsigned int i = 0; i < collection.size(); ++i) {
// Get the channel id and word from the digitization transient object
// We will use the id to determine when to insert missing straws to create a container that represents the complete detector.
// We will only store the digit word.
const TRT_LoLumRawData* tchan0 = dynamic_cast<const TRT_LoLumRawData*>(collection[i]);
tchan_id = tchan0->identify().get_identifier32().get_compact(); // the i'th channel id in the transient collection
tchan_word = tchan0->getWord(); // the i'th channel word in the transient collection
while ( trt_channel_id != tchan_id) {
persCont->m_rawdata[trt_channel_index] = dummy_digit;
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
}
// Here we have a matching transient channel; write it.
persCont->m_rawdata[trt_channel_index]=tchan_word;
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
} // end transient channel loop
// write all other missing channels
while (trt_channel_index != collection_end) {
persCont->m_rawdata[trt_channel_index] = dummy_digit;
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
}
// increment to the next collection
trt_collection_index++;
} // end transient collection loop
// Finally, create any remaining missing collections
while ( trt_collection_index < trt_number_of_collections ) {
// The transient collection is missing; create the persistent collection with empty digits
trt_channel_id = trt_collection_id[trt_collection_index];
InDetRawDataCollection_p1& pcollection = persCont->m_collections[trt_collection_index];
pcollection.m_id = 0;
pcollection.m_hashId = 0; // no hashId for this empty collection
pcollection.m_begin = 0;
pcollection.m_end = 0;
// Fill all the channels with with empty digits.
for (unsigned int i=0; i<trt_collection_size[trt_collection_index]; ++i) {
persCont->m_rawdata[trt_channel_index] = dummy_digit;
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
}
trt_collection_index++;
}
MSG_DEBUG(log," Prepared " << persCont->m_collections.size() << "Collections");
MSG_DEBUG(log," *** Writing TRT_RDO_Container (TRT_LoLumRawData concrete type)");
trt_collection_id.clear();
trt_collection_size.clear();
} // end of transToPers
void TRT_LoLumRawDataContainerCnv_p3::persToTrans(const InDetRawDataContainer_p3* persCont, TRT_RDO_Container* transCont, MsgStream &log)
{
// The transient model has a container holding collections (up to 14912 TRT detector elements/modules)
// and the collections hold channels (e.g. up to 24 straws in the TRT endcaps). There are a total of 350848 TRT straws.
// The persistent model flattens this so that the persistent container has two vectors: 1) all collections, 2) all RDO
// The persistent collections then only maintain indexes into the container's vector of all channels.
// So here we loop over all collection and extract their channels from the vector.
// Hard-code the number of collections, their ids and their sizes to generate collections and channels for the whole detector
const unsigned int dummy_digit=0;
const unsigned int trt_number_of_collections=14912;
const unsigned int trt_channel_id_increment=0x20u;
std::vector<unsigned int> trt_collection_id; // There are 14912 elements/modules in the TRT
std::vector<unsigned int> trt_collection_size; // Each element has a specific number of straws
set_vectors_of_collection_ids_and_size(trt_collection_id,
trt_collection_size,
trt_number_of_collections);
unsigned int trt_channel_id;
unsigned int trt_channel_index=0;
unsigned int trt_channel_index_old;
unsigned int total_number_of_channels=0;
TRT_RDO_Collection* tcoll=0; // transient collection to be constructed
MSG_DEBUG(log," Reading " << persCont->m_collections.size() << "Collections");
if (persCont->m_collections.size() != trt_number_of_collections)
log << MSG::ERROR << "TRT_LoLumRawDataContainerCnv_p3::persToTrans expected 14912 collections but got " << persCont->m_collections.size() << ". We should be reading the whole detector!" << endmsg;
for (unsigned int trt_collection_index=0; trt_collection_index<trt_number_of_collections; ++trt_collection_index) {
// Create trans collection - in NOT owner of TRT_RDO_RawData (SG::VIEW_ELEMENTS); IDet collection don't have the Ownership policy c'tor
// count the number of non-dummy digits and skip to the next persistent collection if there are none.
unsigned int nchans = trt_collection_size[trt_collection_index];
trt_channel_index_old = trt_channel_index; // the beginning of this collection
unsigned int mchans = 0;
for (unsigned int ichan = 0; ichan < nchans; ++ichan) {
const unsigned int pword = persCont->m_rawdata[trt_channel_index];
if ( pword != dummy_digit ) mchans++;
trt_channel_index++;
}
if (!mchans) continue;
// This collection has "mchans" digits
total_number_of_channels += mchans;
// Create the transeint collection and fill with channels
tcoll = new TRT_RDO_Collection(IdentifierHash(trt_collection_index));
tcoll->setIdentifier(Identifier(trt_collection_id[trt_collection_index]));
tcoll->resize(mchans);
trt_channel_id = trt_collection_id[trt_collection_index]; // the id of the first channel in this collection
trt_channel_index = trt_channel_index_old; // go back to the beginning of this collection and process it again
unsigned int jchan=0; // counter for the non-dummy digits
for (unsigned int ichan = 0; ichan < nchans; ++ichan) {
const unsigned int pword = persCont->m_rawdata[trt_channel_index];
if ( pword == dummy_digit ) {
// advance to the next straw
trt_channel_id += trt_channel_id_increment;
trt_channel_index++;
continue; // don't write a dummy digit
}