From 5bc1399bc69d6b4b7ed3d7c68c4ba52c88c471fb Mon Sep 17 00:00:00 2001
From: Simone Pagan Griso <simone.pagan.griso@cern.ch>
Date: Thu, 30 Apr 2015 10:47:17 +0200
Subject: [PATCH] 'Fix python to determine if is simulation'
(InDetPrepRawDataToxAOD-00-01-16)
---
.../InDetPrepRawDataToxAOD/cmt/requirements | 44 +
.../InDetPrepRawDataToxAOD/doc/TRT.txt | 37 +
.../share/InDetDxAOD.py | 219 ++++
.../src/PixelPrepDataToxAOD.cxx | 1066 +++++++++++++++++
.../src/PixelPrepDataToxAOD.h | 104 ++
.../src/SCT_PrepDataToxAOD.cxx | 590 +++++++++
.../src/SCT_PrepDataToxAOD.h | 87 ++
.../src/TRT_PrepDataToxAOD.cxx | 342 ++++++
.../src/TRT_PrepDataToxAOD.h | 66 +
.../InDetPrepRawDataToxAOD_entries.cxx | 18 +
.../InDetPrepRawDataToxAOD_load.cxx | 4 +
11 files changed, 2577 insertions(+)
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/cmt/requirements
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/doc/TRT.txt
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/share/InDetDxAOD.py
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.cxx
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.h
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.cxx
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.h
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.cxx
create mode 100644 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.h
create mode 100755 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_entries.cxx
create mode 100755 InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_load.cxx
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/cmt/requirements b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/cmt/requirements
new file mode 100644
index 00000000000..0fbff6f8bd2
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/cmt/requirements
@@ -0,0 +1,44 @@
+package InDetPrepRawDataToxAOD
+
+author Alex Alonso
+
+public
+use AtlasPolicy AtlasPolicy-*
+use GaudiInterface GaudiInterface-* External
+
+private
+
+use InDetPrepRawData InDetPrepRawData-* InnerDetector/InDetRecEvent
+
+use InDetIdentifier InDetIdentifier-* InnerDetector/InDetDetDescr
+
+use TRT_ConditionsServices TRT_ConditionsServices-* InnerDetector/InDetConditions
+use TRT_DriftFunctionTool TRT_DriftFunctionTool-* InnerDetector/InDetRecTools
+
+use TrkTruthData TrkTruthData-* Tracking/TrkEvent
+
+use AthenaBaseComps AthenaBaseComps-* Control
+use Identifier Identifier-* DetectorDescription
+
+use xAODTracking xAODTracking-* Event/xAOD
+
+use InDetRawData InDetRawData-* InnerDetector/InDetRawEvent
+
+use InDetSimData InDetSimData-* InnerDetector/InDetRawEvent
+use InDetSimEvent InDetSimEvent-* InnerDetector
+#use GeoPrimitives GeoPrimitives-* DetectorDescription
+#use EventPrimitives EventPrimitives-* Event
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+
+
+use InDetReadoutGeometry InDetReadoutGeometry-* InnerDetector/InDetDetDescr
+
+use AtlasHepMC AtlasHepMC-* External
+use AtlasCLHEP AtlasCLHEP-* External
+use AtlasROOT AtlasROOT-* External
+
+public
+apply_pattern component_library
+apply_pattern declare_joboptions files="*.py"
+library InDetPrepRawDataToxAOD *.cxx components/*.cxx
+
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/doc/TRT.txt b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/doc/TRT.txt
new file mode 100644
index 00000000000..4795afa2ef1
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/doc/TRT.txt
@@ -0,0 +1,37 @@
+This tools is inspired on:
+MuonSpectrometer/MuonCnv/MuonPrepRawDataToxAOD
+
+To run, you have to have flag:
+
+rec.doAOD=True
+InDetFlags.doxAOD.set_Value_and_Lock()
+
+Then, in: InnerDetector/InDetExample/InDetRecExample/share/WriteInDetAOD.py
+Add:
+if InDetFlags.doxAOD():
+ InDetAODList += ['xAOD::PrepRawDataContainer_v1#TRT_DriftCircles']
+ InDetAODList += ['xAOD::PrepRawDataAuxContainer_v1#TRT_DriftCirclesAux.']
+
+In: InnerDetector/InDetExample/InDetRecExample/share/InDetxAODCreator.py
+
+from TRT_ConditionsServices.TRT_ConditionsServicesConf import
+TRT_StrawNeighbourSvc
+TRTStrawNeighbourSvc=TRT_StrawNeighbourSvc()
+ServiceMgr += TRTStrawNeighbourSvc
+
+from TRT_ConditionsServices.TRT_ConditionsServicesConf import TRT_CalDbSvc
+TRTCalibDBSvc=TRT_CalDbSvc()
+ServiceMgr += TRTCalibDBSvc
+
+from InDetPrepRawDataToxAOD.InDetPrepRawDataToxAODConf import
+TRT_PrepDataToxAOD
+xAOD_TRT_PrepDataToxAOD = TRT_PrepDataToxAOD( name =
+"xAOD_TRT_PrepDataToxAOD")
+xAOD_TRT_PrepDataToxAOD.OutputLevel=DEBUG
+print "Add TRT xAOD PrepRawData:"
+print xAOD_TRT_PrepDataToxAOD
+topSequence += xAOD_TRT_PrepDataToxAOD
+
+
+
+
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/share/InDetDxAOD.py b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/share/InDetDxAOD.py
new file mode 100644
index 00000000000..4887e638ab8
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/share/InDetDxAOD.py
@@ -0,0 +1,219 @@
+#################
+### Steering options
+#################
+# Select active sub-systems
+dumpPixInfo=True
+dumpSctInfo=True
+dumpTrtInfo=False
+
+# Bytestream errors (for sub-systems who have implemented it)
+dumpBytestreamErrors=True
+
+# Force to do not dump truth info if set to False
+# (otherwise determined by autoconf below)
+dumpTruthInfo=True
+
+# Print settings for main tools
+printIdTrkDxAODConf = True
+
+## Autoconfiguration adjustements
+isIdTrkDxAODSimulation = False
+if (globalflags.DataSource == 'geant4'):
+ isIdTrkDxAODSimulation = True
+
+if ( 'dumpTruthInfo' in dir() ):
+ dumpTruthInfo = dumpTruthInfo and isIdTrkDxAODSimulation
+
+print "isIdTrkDxAODSimulation=", isIdTrkDxAODSimulation
+print "dumpTruthInfo=", dumpTruthInfo
+
+if InDetFlags.doSLHC():
+ dumpTrtInfo=False
+
+## Other settings
+# Prefix for decoration, if any
+prefixName = ""
+
+## More fine-tuning available for each tool/alg below (default value shown)
+
+#################
+### Setup tools
+#################
+if dumpTrtInfo:
+ from TRT_ConditionsServices.TRT_ConditionsServicesConf import TRT_StrawNeighbourSvc
+ TRTStrawNeighbourSvc=TRT_StrawNeighbourSvc()
+ ServiceMgr += TRTStrawNeighbourSvc
+
+ from TRT_ConditionsServices.TRT_ConditionsServicesConf import TRT_CalDbSvc
+ TRTCalibDBSvc=TRT_CalDbSvc()
+ ServiceMgr += TRTCalibDBSvc
+
+ from TRT_ToT_Tools.TRT_ToT_ToolsConf import TRT_ToT_dEdx
+ TRT_dEdx_Tool = TRT_ToT_dEdx(name="TRT_ToT_dEdx")
+ ToolSvc += TRT_dEdx_Tool
+
+if dumpPixInfo:
+ from AthenaCommon.AlgSequence import AlgSequence
+ topSequence = AlgSequence()
+ from PixelCalibAlgs.PixelCalibAlgsConf import PixelChargeToTConversion
+ PixelChargeToTConversionSetter = PixelChargeToTConversion(name = "PixelChargeToTConversionSetter")
+ topSequence += PixelChargeToTConversionSetter
+ if (printIdTrkDxAODConf):
+ print PixelChargeToTConversionSetter
+ print PixelChargeToTConversionSetter.properties()
+
+
+#################
+### Setup decorators tools
+#################
+
+if dumpTrtInfo:
+ from InDetPrepRawDataToxAOD.InDetPrepRawDataToxAODConf import TRT_PrepDataToxAOD
+ xAOD_TRT_PrepDataToxAOD = TRT_PrepDataToxAOD( name = "xAOD_TRT_PrepDataToxAOD")
+ ## Content steering Properties (default value shown as comment)
+ xAOD_TRT_PrepDataToxAOD.OutputLevel=INFO
+ xAOD_TRT_PrepDataToxAOD.UseTruthInfo = dumpTruthInfo
+ #xAOD_TRT_PrepDataToxAOD.WriteSDOs = True
+
+ topSequence += xAOD_TRT_PrepDataToxAOD
+ if (printIdTrkDxAODConf):
+ print xAOD_TRT_PrepDataToxAOD
+ print xAOD_TRT_PrepDataToxAOD.properties()
+
+if dumpSctInfo:
+ from InDetPrepRawDataToxAOD.InDetPrepRawDataToxAODConf import SCT_PrepDataToxAOD
+ xAOD_SCT_PrepDataToxAOD = SCT_PrepDataToxAOD( name = "xAOD_SCT_PrepDataToxAOD")
+ ## Content steering Properties (default value shown as comment)
+ xAOD_SCT_PrepDataToxAOD.OutputLevel=INFO
+ xAOD_SCT_PrepDataToxAOD.UseTruthInfo = dumpTruthInfo
+ xAOD_SCT_PrepDataToxAOD.WriteRDOinformation = False
+ #xAOD_SCT_PrepDataToxAOD.WriteSDOs = True
+ #xAOD_SCT_PrepDataToxAOD.WriteSiHits = True # if available
+
+ topSequence += xAOD_SCT_PrepDataToxAOD
+ if (printIdTrkDxAODConf):
+ print xAOD_SCT_PrepDataToxAOD
+ print xAOD_SCT_PrepDataToxAOD.properties()
+
+if dumpPixInfo:
+ from InDetPrepRawDataToxAOD.InDetPrepRawDataToxAODConf import PixelPrepDataToxAOD
+ xAOD_PixelPrepDataToxAOD = PixelPrepDataToxAOD( name = "xAOD_PixelPrepDataToxAOD")
+ ## Content steering Properties (default value shown as comment)
+ xAOD_PixelPrepDataToxAOD.OutputLevel = INFO
+ xAOD_PixelPrepDataToxAOD.UseTruthInfo = dumpTruthInfo
+ xAOD_PixelPrepDataToxAOD.WriteRDOinformation = False
+ xAOD_PixelPrepDataToxAOD.WriteNNinformation = False
+ #xAOD_PixelPrepDataToxAOD.WriteSDOs = True
+ #xAOD_PixelPrepDataToxAOD.WriteSiHits = True # if available
+ if InDetFlags.doSLHC():
+ xAOD_PixelPrepDataToxAOD.WriteNNinformation=False
+
+ topSequence += xAOD_PixelPrepDataToxAOD
+ if (printIdTrkDxAODConf):
+ print xAOD_PixelPrepDataToxAOD
+ print xAOD_PixelPrepDataToxAOD.properties()
+
+
+#################
+### Setup derivation framework
+#################
+from AthenaCommon import CfgMgr
+
+# DerivationJob is COMMON TO ALL DERIVATIONS
+DerivationFrameworkJob = CfgMgr.AthSequencer("MySeq2")
+
+# Set up stream auditor
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+if not hasattr(svcMgr, 'DecisionSvc'):
+ svcMgr += CfgMgr.DecisionSvc()
+svcMgr.DecisionSvc.CalcStats = True
+
+
+# Add the TSOS augmentation tool to the derivation framework
+augmentationTools=[]
+
+from DerivationFrameworkInDet.DerivationFrameworkInDetConf import DerivationFramework__TrackStateOnSurfaceDecorator
+DFTSOS = DerivationFramework__TrackStateOnSurfaceDecorator(name = "DFTrackStateOnSurfaceDecorator",
+ ContainerName = "InDetTrackParticles",
+ DecorationPrefix = prefixName,
+ StoreTRT = dumpTrtInfo,
+ StoreSCT = dumpSctInfo,
+ StorePixel = dumpPixInfo,
+ IsSimulation = isIdTrkDxAODSimulation,
+ OutputLevel = INFO)
+if dumpTrtInfo:
+ #Add tool to calculate TRT-based dEdx
+ DFTSOS.TRT_ToT_dEdx = TRT_dEdx_Tool
+
+ToolSvc += DFTSOS
+augmentationTools+=[DFTSOS]
+if (printIdTrkDxAODConf):
+ print DFTSOS
+ print DFTSOS.properties()
+
+# Add BS error augmentation tool
+if dumpBytestreamErrors:
+ from DerivationFrameworkInDet.DerivationFrameworkInDetConf import DerivationFramework__EventInfoBSErrDecorator
+ DFEI = DerivationFramework__EventInfoBSErrDecorator(name = "DFEventInfoBSErrDecorator",
+ ContainerName = "EventInfo",
+ DecorationPrefix = prefixName,
+ OutputLevel =INFO)
+ ToolSvc += DFEI
+ augmentationTools+=[DFEI]
+ if (printIdTrkDxAODConf):
+ print DFEI
+ print DFEI.properties()
+
+# Add the derivation job to the top AthAlgSeqeuence
+from DerivationFrameworkCore.DerivationFrameworkCoreConf import DerivationFramework__CommonAugmentation
+DerivationFrameworkJob += CfgMgr.DerivationFramework__CommonAugmentation("DFTSOS_KERN",
+ AugmentationTools = augmentationTools,
+ OutputLevel =INFO)
+
+topSequence += DerivationFrameworkJob
+if (printIdTrkDxAODConf):
+ print DerivationFrameworkJob
+ print DerivationFrameworkJob.properties()
+
+#################
+### Steer output file content
+#################
+from OutputStreamAthenaPool.MultipleStreamManager import MSMgr
+from D2PDMaker.D2PDHelpers import buildFileName
+from PrimaryDPDMaker.PrimaryDPDFlags import primDPD
+streamName = primDPD.WriteDAOD_IDTRKVALIDStream.StreamName
+fileName = buildFileName( primDPD.WriteDAOD_IDTRKVALIDStream )
+IDTRKVALIDStream = MSMgr.NewPoolRootStream( streamName, fileName )
+excludedAuxData = "-caloExtension.-cellAssociation.-clusterAssociation.-trackParameterCovarianceMatrices.-parameterX.-parameterY.-parameterZ.-parameterPX.-parameterPY.-parameterPZ.-parameterPosition"
+# Add generic event information
+IDTRKVALIDStream.AddItem("xAOD::EventInfo#*")
+IDTRKVALIDStream.AddItem("xAOD::EventAuxInfo#*")
+# Add track particles collection
+IDTRKVALIDStream.AddItem("xAOD::TrackParticleContainer#InDetTrackParticles")
+IDTRKVALIDStream.AddItem("xAOD::TrackParticleAuxContainer#InDetTrackParticlesAux."+excludedAuxData)
+# Add vertices
+IDTRKVALIDStream.AddItem("xAOD::VertexContainer#PrimaryVertices")
+IDTRKVALIDStream.AddItem("xAOD::VertexAuxContainer#PrimaryVerticesAux.-vxTrackAtVertex")
+# Add links and measurements
+IDTRKVALIDStream.AddItem("xAOD::TrackStateValidationContainer#*")
+IDTRKVALIDStream.AddItem("xAOD::TrackStateValidationAuxContainer#*")
+IDTRKVALIDStream.AddItem("xAOD::TrackMeasurementValidationContainer#*")
+IDTRKVALIDStream.AddItem("xAOD::TrackMeasurementValidationAuxContainer#*")
+# Add info about electrons and muons (are small containers)
+IDTRKVALIDStream.AddItem("xAOD::MuonContainer#Muons")
+IDTRKVALIDStream.AddItem("xAOD::MuonAuxContainer#MuonsAux.")
+IDTRKVALIDStream.AddItem("xAOD::ElectronContainer#Electrons")
+IDTRKVALIDStream.AddItem("xAOD::ElectronAuxContainer#ElectronsAux.")
+IDTRKVALIDStream.AddItem("xAOD::TrackParticleContainer#GSFTrackParticles")
+IDTRKVALIDStream.AddItem("xAOD::TrackParticleAuxContainer#GSFTrackParticlesAux."+excludedAuxData)
+# Add truth-related information
+if dumpTruthInfo:
+ IDTRKVALIDStream.AddItem("xAOD::TruthParticleContainer#*")
+ IDTRKVALIDStream.AddItem("xAOD::TruthParticleAuxContainer#TruthParticlesAux.-caloExtension")
+ IDTRKVALIDStream.AddItem("xAOD::TruthVertexContainer#*")
+ IDTRKVALIDStream.AddItem("xAOD::TruthVertexAuxContainer#*")
+ IDTRKVALIDStream.AddItem("xAOD::TruthEventContainer#*")
+ IDTRKVALIDStream.AddItem("xAOD::TruthEventAuxContainer#*")
+
+if (printIdTrkDxAODConf):
+ print IDTRKVALIDStream
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.cxx b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.cxx
new file mode 100644
index 00000000000..2824e53949a
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.cxx
@@ -0,0 +1,1066 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PixelPrepDataToxAOD.cxx
+// Implementation file for class PixelPrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#include "PixelPrepDataToxAOD.h"
+
+#include "InDetPrepRawData/PixelClusterContainer.h"
+
+#include "xAODTracking/TrackMeasurementValidation.h"
+#include "xAODTracking/TrackMeasurementValidationContainer.h"
+#include "xAODTracking/TrackMeasurementValidationAuxContainer.h"
+
+#include "Identifier/Identifier.h"
+#include "InDetIdentifier/PixelID.h"
+#include "InDetReadoutGeometry/PixelModuleDesign.h"
+
+
+#include "TrkTruthData/PRD_MultiTruthCollection.h"
+#include "HepMC/GenParticle.h"
+#include "HepMC/GenVertex.h"
+#include "InDetSimEvent/SiHit.h"
+#include "InDetSimData/InDetSimDataCollection.h"
+
+
+#include "TMath.h"
+#include "CLHEP/Geometry/Point3D.h"
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// Constructor with parameters:
+//
+/////////////////////////////////////////////////////////////////////
+PixelPrepDataToxAOD::PixelPrepDataToxAOD(const std::string &name, ISvcLocator *pSvcLocator) :
+ AthAlgorithm(name,pSvcLocator),
+ m_PixelHelper(0),
+ m_firstEventWarnings(true)
+{
+ // --- Steering and configuration flags
+
+ declareProperty("UseTruthInfo", m_useTruthInfo=false);
+ declareProperty("WriteSDOs", m_writeSDOs = true);
+ declareProperty("WriteSiHits", m_writeSiHits = true);
+ declareProperty("WriteNNinformation", m_writeNNinformation = true);
+ declareProperty("WriteRDOinformation", m_writeRDOinformation = true);
+
+ // --- Configuration keys
+ declareProperty("SiClusterContainer", m_clustercontainer = "PixelClusters");
+ declareProperty("MC_SDOs", m_SDOcontainer = "PixelSDO_Map");
+ declareProperty("MC_Hits", m_sihitContainer = "PixelHits");
+ declareProperty("PRD_MultiTruth", m_multiTruth = "PRD_MultiTruthPixel");
+
+ // --- Services and Tools
+
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Initialize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode PixelPrepDataToxAOD::initialize()
+{
+ ATH_CHECK( detStore()->retrieve(m_PixelHelper, "PixelID") );
+
+ //make sure we don't write what we don't have
+ if (not m_useTruthInfo) {
+ m_writeSDOs = false;
+ m_writeSiHits = false;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Execute method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode PixelPrepDataToxAOD::execute()
+{
+ //Mandatory. Require if the algorithm is scheduled.
+ const InDet::PixelClusterContainer* PixelClusterContainer = 0;
+ if( evtStore()->retrieve(PixelClusterContainer,m_clustercontainer).isFailure() ) {
+ ATH_MSG_ERROR("Cannot retrieve Pixel PrepDataContainer " << m_clustercontainer);
+ return StatusCode::FAILURE;
+ }
+
+ // Optional. Normally only available in Hits files -- samples need to digitised and Hits need to be copied for this to work
+ const SiHitCollection* sihitCollection = 0;
+ if((m_writeNNinformation||m_writeSiHits)&&m_useTruthInfo) {
+ if ( evtStore()->contains<SiHitCollection>(m_sihitContainer) ) {
+ ATH_CHECK(evtStore()->retrieve(sihitCollection, m_sihitContainer));
+ } else {
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("Si Hit cotainer no found (" << m_sihitContainer << "). Skipping it although requested.");
+ sihitCollection = 0;
+ }
+ }
+ }
+
+ // Optional. On RDO
+ const InDetSimDataCollection* sdoCollection = 0;
+ if(m_writeSDOs) {
+ if ( evtStore()->contains<InDetSimDataCollection>(m_SDOcontainer) ) {
+ ATH_CHECK(evtStore()->retrieve(sdoCollection, m_SDOcontainer));
+ } else {
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("SDO Collection not found (" << m_SDOcontainer << "). Skipping it although requested.");
+ sdoCollection = 0;
+ }
+ }
+ }
+
+ // Optional. On ESD and AOD
+ const PRD_MultiTruthCollection* prdmtColl = 0;
+ if (m_useTruthInfo) {
+ if ( evtStore()->contains<PRD_MultiTruthCollection>(m_multiTruth) ) {
+ ATH_CHECK(evtStore()->retrieve(prdmtColl, m_multiTruth));
+ } else {
+ ATH_MSG_WARNING("PRD Truth collection missing (" << m_multiTruth << "). Skipping it although requested.");
+ prdmtColl = 0;
+ }
+ }
+
+
+ // Create the xAOD container and its auxiliary store:
+ xAOD::TrackMeasurementValidationContainer* xaod = new xAOD::TrackMeasurementValidationContainer();
+ CHECK( evtStore()->record( xaod, m_clustercontainer ) );
+ xAOD::TrackMeasurementValidationAuxContainer* aux = new xAOD::TrackMeasurementValidationAuxContainer();
+ CHECK( evtStore()->record( aux, m_clustercontainer + "Aux." ) );
+ xaod->setStore( aux );
+
+ std::vector<unsigned int>* offsets = new std::vector<unsigned int>( m_PixelHelper->wafer_hash_max(), 0 );
+ CHECK( evtStore()->record( offsets, m_clustercontainer + "Offsets" ) );
+
+
+ // Loop over the container
+ unsigned int counter(0);
+
+ for( const auto& clusterCollection : *PixelClusterContainer){
+
+ //Fill Offset container
+ (*offsets)[clusterCollection->identifyHash()] = counter;
+
+ // skip empty collections
+ if( clusterCollection->empty() ) continue;
+
+ // loop over collection and convert to xAOD
+ for( auto& prd : *clusterCollection ){
+ ++counter;
+
+ Identifier clusterId = prd->identify();
+ if ( !clusterId.is_valid() ) {
+ ATH_MSG_WARNING("Pixel cluster identifier is not valid");
+ }
+
+ // create and add xAOD object
+ xAOD::TrackMeasurementValidation* xprd = new xAOD::TrackMeasurementValidation();
+ xaod->push_back(xprd);
+
+ //Set Identifier
+ xprd->setIdentifier( clusterId.get_compact() );
+
+ //Set Global Position
+ Amg::Vector3D gpos = prd->globalPosition();
+ xprd->setGlobalPosition(gpos.x(),gpos.y(),gpos.z());
+
+ //Set Local Position
+ const Amg::Vector2D& locpos = prd->localPosition();
+
+ // Set local error matrix
+ xprd->setLocalPosition( locpos.x(), locpos.y() );
+
+ const Amg::MatrixX& localCov = prd->localCovariance();
+ //std::cout << localCov << std::endl;
+ if(localCov.size() == 1){
+ //std::cout << "Size == 1" << std::endl;
+ xprd->setLocalPositionError( localCov(0,0), 0., 0. );
+ } else if(localCov.size() == 4){
+ //std::cout << "Size == 2" << std::endl;
+ xprd->setLocalPositionError( localCov(0,0), localCov(1,1), localCov(0,1) );
+ } else {
+ //std::cout << "Size == "<< localCov.size() << std::endl;
+ xprd->setLocalPositionError(0.,0.,0.);
+ }
+
+ // Set vector of hit identifiers
+ std::vector< uint64_t > rdoIdentifierList;
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ rdoIdentifierList.push_back( hitIdentifier.get_compact() );
+ //May want to addinformation about the individual hits here
+ }
+ xprd->setRdoIdentifierList(rdoIdentifierList);
+
+ //Add pixel cluster properties
+ xprd->auxdata<int>("bec") = m_PixelHelper->barrel_ec(clusterId) ;
+ xprd->auxdata<int>("layer") = m_PixelHelper->layer_disk(clusterId) ;
+ xprd->auxdata<int>("phi_module") = m_PixelHelper->phi_module(clusterId) ;
+ xprd->auxdata<int>("eta_module") = m_PixelHelper->eta_module(clusterId) ;
+
+ //xprd->auxdata<int>("col") = m_PixelHelper->eta_index(clusterId);
+ //xprd->auxdata<int>("row") = m_PixelHelper->phi_index(clusterId);
+ xprd->auxdata<int>("eta_pixel_index") = m_PixelHelper->eta_index(clusterId);
+ xprd->auxdata<int>("phi_pixel_index") = m_PixelHelper->phi_index(clusterId);
+
+
+ const InDet::SiWidth cw = prd->width();
+ xprd->auxdata<int>("sizePhi") = (int)cw.colRow()[0];
+ xprd->auxdata<int>("sizeZ") = (int)cw.colRow()[1];
+ xprd->auxdata<int>("size") = (int)prd->rdoList().size();
+
+ xprd->auxdata<float>("charge") = prd->totalCharge();
+ xprd->auxdata<int>("ToT") = prd->totalToT();
+ xprd->auxdata<int>("LVL1A") = prd->LVL1A();
+
+ xprd->auxdata<char>("isFake") = (char)prd->isFake();
+ xprd->auxdata<char>("gangedPixel") = (char)prd->gangedPixel();
+ xprd->auxdata<int>("isSplit") = (int)prd->isSplit();
+ xprd->auxdata<float>("splitProbability1") = prd->splitProbability1();
+ xprd->auxdata<float>("splitProbability2") = prd->splitProbability2();
+
+ // Need to add something to Add the NN splitting information
+ if(m_writeNNinformation)addNNInformation( xprd, prd, 7, 7);
+
+ // Add information for each contributing hit
+ if(m_writeRDOinformation) addRdoInformation(xprd, prd);
+
+
+
+ // Add the Detector element ID -- not sure if needed as we have the informations above
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ uint64_t detElementId(0);
+ if(de){
+ Identifier detId = de->identify();
+ if ( detId.is_valid() ) {
+ detElementId = detId.get_compact();
+ }
+ }
+ xprd->auxdata<uint64_t>("detectorElementID") = detElementId;
+
+ // Use the MultiTruth Collection to get a list of all true particle contributing to the cluster
+ if(prdmtColl){
+ std::vector<int> barcodes;
+ //std::pair<PRD_MultiTruthCollection::const_iterator,PRD_MultiTruthCollection::const_iterator>;
+ auto range = prdmtColl->equal_range(clusterId);
+ for (auto i = range.first; i != range.second; ++i) {
+ barcodes.push_back( i->second.barcode() );
+ }
+ xprd->auxdata< std::vector<int> >("truth_barcode") = barcodes;
+ }
+
+ // Use the SDO Collection to get a list of all true particle contributing to the cluster per readout element
+ // Also get the energy deposited by each true particle per readout element
+ if( sdoCollection ){
+ addSDOInformation( xprd, prd,sdoCollection);
+ }
+
+ // Now Get the most detailed truth from the SiHits
+ // Note that this could get really slow if there are a lot of hits and clusters
+ if( sihitCollection ){
+ if(m_writeSiHits)
+ addSiHitInformation( xprd, prd, sihitCollection);
+ if(m_writeNNinformation)
+ addNNTruthInfo( xprd, prd, sihitCollection );
+ }
+ }
+ }
+ ATH_MSG_DEBUG( " recorded PixelPrepData objects: size " << xaod->size() );
+
+ m_firstEventWarnings = false;
+
+ return StatusCode::SUCCESS;
+}
+
+
+void PixelPrepDataToxAOD::addSDOInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* prd,
+ const InDetSimDataCollection* sdoCollection ) const
+{
+ std::vector<int> sdo_word;
+ std::vector< std::vector< int > > sdo_depositsBarcode;
+ std::vector< std::vector< float > > sdo_depositsEnergy;
+ // find hit
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ auto pos = sdoCollection->find(hitIdentifier);
+ if( pos != sdoCollection->end() ) {
+ sdo_word.push_back( pos->second.word() ) ;
+
+ std::vector< int > sdoDepBC;
+ std::vector< float > sdoDepEnergy;
+ for( auto deposit : pos->second.getdeposits() ){
+ if(deposit.first){
+ sdoDepBC.push_back( deposit.first->barcode());
+ } else {
+ sdoDepBC.push_back( -1 );
+ }
+ sdoDepEnergy.push_back( deposit.second );
+ ATH_MSG_DEBUG(" SDO Energy Deposit " << deposit.second ) ;
+ }
+ sdo_depositsBarcode.push_back( sdoDepBC );
+ sdo_depositsEnergy.push_back( sdoDepEnergy );
+ }
+ }
+ xprd->auxdata< std::vector<int> >("sdo_words") = sdo_word;
+ xprd->auxdata< std::vector< std::vector<int> > >("sdo_depositsBarcode") = sdo_depositsBarcode;
+ xprd->auxdata< std::vector< std::vector<float> > >("sdo_depositsEnergy") = sdo_depositsEnergy;
+}
+
+
+
+void PixelPrepDataToxAOD::addSiHitInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* prd,
+ const SiHitCollection* sihitCollection) const
+
+
+{
+
+ std::vector<SiHit> matchingHits = findAllHitsCompatibleWithCluster( prd, sihitCollection );
+
+ int numHits = matchingHits.size();
+
+ std::vector<float> sihit_energyDeposit(numHits,0);
+ std::vector<float> sihit_meanTime(numHits,0);
+ std::vector<int> sihit_barcode(numHits,0);
+ std::vector<int> sihit_pdgid(numHits,0);
+
+ std::vector<float> sihit_startPosX(numHits,0);
+ std::vector<float> sihit_startPosY(numHits,0);
+ std::vector<float> sihit_startPosZ(numHits,0);
+
+ std::vector<float> sihit_endPosX(numHits,0);
+ std::vector<float> sihit_endPosY(numHits,0);
+ std::vector<float> sihit_endPosZ(numHits,0);
+
+ int hitNumber(0);
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ if(de){
+ for ( auto sihit : matchingHits ) {
+ sihit_energyDeposit[hitNumber] = sihit.energyLoss() ;
+ sihit_meanTime[hitNumber] = sihit.meanTime() ;
+ sihit_barcode[hitNumber] = sihit.particleLink().barcode() ;
+ if(sihit.particleLink().isValid()){
+ sihit_pdgid[hitNumber] = sihit.particleLink()->pdg_id();
+ }
+
+ // Convert Simulation frame into reco frame
+ const HepGeom::Point3D<double>& startPos=sihit.localStartPosition();
+
+ Amg::Vector2D pos= de->hitLocalToLocal( startPos.z(), startPos.y() );
+ sihit_startPosX[hitNumber] = pos[0];
+ sihit_startPosY[hitNumber] = pos[1];
+ sihit_startPosZ[hitNumber] = startPos.x();
+
+
+ const HepGeom::Point3D<double>& endPos=sihit.localEndPosition();
+ pos= de->hitLocalToLocal( endPos.z(), endPos.y() );
+ sihit_endPosX[hitNumber] = pos[0];
+ sihit_endPosY[hitNumber] = pos[1];
+ sihit_endPosZ[hitNumber] = endPos.x();
+ ++hitNumber;
+ }
+ }
+
+ xprd->auxdata<std::vector<float> >("sihit_energyDeposit") = sihit_energyDeposit;
+ xprd->auxdata<std::vector<float> >("sihit_meanTime") = sihit_meanTime;
+ xprd->auxdata<std::vector<int> >("sihit_barcode") = sihit_barcode;
+ xprd->auxdata<std::vector<int> >("sihit_pdgid") = sihit_pdgid;
+
+ xprd->auxdata<std::vector<float> >("sihit_startPosX") = sihit_startPosX;
+ xprd->auxdata<std::vector<float> >("sihit_startPosY") = sihit_startPosY;
+ xprd->auxdata<std::vector<float> >("sihit_startPosZ") = sihit_startPosZ;
+
+ xprd->auxdata<std::vector<float> >("sihit_endPosX") = sihit_endPosX;
+ xprd->auxdata<std::vector<float> >("sihit_endPosY") = sihit_endPosY;
+ xprd->auxdata<std::vector<float> >("sihit_endPosZ") = sihit_endPosZ;
+
+
+}
+
+
+
+
+
+
+std::vector<SiHit> PixelPrepDataToxAOD::findAllHitsCompatibleWithCluster( const InDet::PixelCluster* prd,
+ const SiHitCollection* collection) const
+{
+ ATH_MSG_VERBOSE( "Got " << collection->size() << " SiHits to look through" );
+ std::vector<SiHit> matchingHits;
+
+ // Check if we have detector element -- needed to find the local position of the SiHits
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ if(!de)
+ return matchingHits;
+
+ std::vector<const SiHit* > multiMatchingHits;
+
+ for ( const auto& siHit : *collection) {
+ // Check if it is a Pixel hit
+ if( !siHit.isPixel() )
+ continue;
+
+ //Check if it is on the correct module
+ Identifier clusterId = prd->identify();
+
+ if( m_PixelHelper->barrel_ec(clusterId) != siHit.getBarrelEndcap() ||
+ m_PixelHelper->layer_disk(clusterId)!= siHit.getLayerDisk() ||
+ m_PixelHelper->phi_module(clusterId)!= siHit.getPhiModule() ||
+ m_PixelHelper->eta_module(clusterId)!= siHit.getEtaModule() )
+ continue;
+
+ // Now we have all hits in the module that match lets check to see if they match the cluster
+ // Must be within +/- 1 hits of any hit in the cluster to be included
+ ATH_MSG_DEBUG("Hit is on the same module");
+
+ HepGeom::Point3D<double> averagePosition = siHit.localStartPosition() + siHit.localEndPosition();
+ averagePosition *= 0.5;
+ Amg::Vector2D pos = de->hitLocalToLocal( averagePosition.z(), averagePosition.y() );
+ InDetDD::SiCellId diode = de->cellIdOfPosition(pos);
+
+ for( const auto &hitIdentifier : prd->rdoList() ){
+
+ ATH_MSG_DEBUG("Truth Phi " << diode.phiIndex() << " Cluster Phi " << m_PixelHelper->phi_index( hitIdentifier ) );
+ ATH_MSG_DEBUG("Truth Eta " << diode.etaIndex() << " Cluster Eta " << m_PixelHelper->eta_index( hitIdentifier ) );
+
+ if( abs( int(diode.etaIndex()) - m_PixelHelper->eta_index( hitIdentifier ) ) <=1
+ && abs( int(diode.phiIndex()) - m_PixelHelper->phi_index( hitIdentifier ) ) <=1 )
+ {
+ multiMatchingHits.push_back(&siHit);
+ break;
+ }
+ }
+ }
+
+
+ //Now we will now make 1 SiHit for each true particle if the SiHits "touch" other
+ std::vector<const SiHit* >::iterator siHitIter = multiMatchingHits.begin();
+ std::vector<const SiHit* >::iterator siHitIter2 = multiMatchingHits.begin();
+ ATH_MSG_DEBUG( "Found " << multiMatchingHits.size() << " SiHit " );
+ for ( ; siHitIter != multiMatchingHits.end(); siHitIter++) {
+ const SiHit* lowestXPos = *siHitIter;
+ const SiHit* highestXPos = *siHitIter;
+
+
+ // We will merge these hits
+ std::vector<const SiHit* > ajoiningHits;
+ ajoiningHits.push_back( *siHitIter );
+
+ siHitIter2 = siHitIter+1;
+ while ( siHitIter2 != multiMatchingHits.end() ) {
+ // Need to come from the same truth particle
+
+ if( (*siHitIter)->particleLink().barcode() != (*siHitIter2)->particleLink().barcode() ){
+ ++siHitIter2;
+ continue;
+ }
+
+ // Check to see if the SiHits are compatible with each other.
+ if (fabs((highestXPos->localEndPosition().x()-(*siHitIter2)->localStartPosition().x()))<0.00005 &&
+ fabs((highestXPos->localEndPosition().y()-(*siHitIter2)->localStartPosition().y()))<0.00005 &&
+ fabs((highestXPos->localEndPosition().z()-(*siHitIter2)->localStartPosition().z()))<0.00005 )
+ {
+ highestXPos = *siHitIter2;
+ ajoiningHits.push_back( *siHitIter2 );
+ // Dont use hit more than once
+ siHitIter2 = multiMatchingHits.erase( siHitIter2 );
+ }else if (fabs((lowestXPos->localStartPosition().x()-(*siHitIter2)->localEndPosition().x()))<0.00005 &&
+ fabs((lowestXPos->localStartPosition().y()-(*siHitIter2)->localEndPosition().y()))<0.00005 &&
+ fabs((lowestXPos->localStartPosition().z()-(*siHitIter2)->localEndPosition().z()))<0.00005)
+ {
+ lowestXPos = *siHitIter2;
+ ajoiningHits.push_back( *siHitIter2 );
+ // Dont use hit more than once
+ siHitIter2 = multiMatchingHits.erase( siHitIter2 );
+ } else {
+ ++siHitIter2;
+ }
+ }
+
+ if( ajoiningHits.size() == 0){
+ ATH_MSG_WARNING("This should really never happen");
+ continue;
+ }else if(ajoiningHits.size() == 1){
+ // Copy Si Hit ready to return
+ matchingHits.push_back( *ajoiningHits[0] );
+ continue;
+ } else {
+ // Build new SiHit and merge information together.
+ ATH_MSG_DEBUG("Merging " << ajoiningHits.size() << " SiHits together." );
+
+
+ float energyDep(0);
+ float time(0);
+ for( auto& siHit : ajoiningHits){
+ energyDep += siHit->energyLoss();
+ time += siHit->meanTime();
+ }
+ time /= (float)ajoiningHits.size();
+
+ matchingHits.push_back( SiHit(lowestXPos->localStartPosition(),
+ highestXPos->localEndPosition(),
+ energyDep,
+ time,
+ (*siHitIter)->particleLink().barcode(),
+ 0, // 0 for pixel 1 for Pixel
+ (*siHitIter)->getBarrelEndcap(),
+ (*siHitIter)->getLayerDisk(),
+ (*siHitIter)->getEtaModule(),
+ (*siHitIter)->getPhiModule(),
+ (*siHitIter)->getSide() ) );
+ ATH_MSG_DEBUG("Finished Merging " << ajoiningHits.size() << " SiHits together." );
+
+ }
+ }
+
+
+ return matchingHits;
+
+}
+
+void PixelPrepDataToxAOD::addRdoInformation(xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* pixelCluster ) const
+{
+ ATH_MSG_VERBOSE( " Starting creating input from cluster " );
+
+
+ const std::vector<Identifier>& rdos = pixelCluster->rdoList();
+
+ const std::vector<float> chList = pixelCluster->chargeList();
+ const std::vector<int> totList = pixelCluster->totList();
+
+ // std::vector<int> rowList;
+ // std::vector<int> colList;
+ std::vector<int> etaIndexList;
+ std::vector<int> phiIndexList;
+
+ ATH_MSG_VERBOSE( "Number of RDOs: " << rdos.size() );
+
+ //Itererate over all elements hits in the cluster and fill the charge and tot matricies
+ std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
+ std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
+
+ ATH_MSG_VERBOSE(" Putting together the n. " << rdos.size() << " rdos into a matrix.");
+
+ for (; rdosBegin!= rdosEnd; ++rdosBegin)
+ {
+ Identifier rId = *rdosBegin;
+ // rowList.push_back( m_PixelHelper->phi_index(rId) );
+ // colList.push_back( m_PixelHelper->eta_index(rId) );
+ phiIndexList.push_back( m_PixelHelper->phi_index(rId) );
+ etaIndexList.push_back( m_PixelHelper->eta_index(rId) );
+ }//end iteration on rdos
+
+
+ // xprd->auxdata< std::vector<int> >("rdo_row") = rowList;
+ // xprd->auxdata< std::vector<int> >("rdo_col") = colList;
+ xprd->auxdata< std::vector<int> >("rdo_phi_pixel_index") = phiIndexList;
+ xprd->auxdata< std::vector<int> >("rdo_eta_pixel_index") = etaIndexList;
+ xprd->auxdata< std::vector<float> >("rdo_charge") = chList;
+ xprd->auxdata< std::vector<int> >("rdo_tot") = totList;
+
+}
+
+
+void PixelPrepDataToxAOD::addNNInformation(xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* pixelCluster,
+ const unsigned int sizeX, const unsigned int sizeY ) const
+{
+ ATH_MSG_VERBOSE( " Starting creating input from cluster " );
+
+ const InDetDD::SiDetectorElement* de = pixelCluster->detectorElement();
+ if (de==0) {
+ ATH_MSG_ERROR("Could not get detector element");
+ return;
+ }
+
+
+ const InDetDD::PixelModuleDesign* design(dynamic_cast<const InDetDD::PixelModuleDesign*>(&de->design()));
+ if (not design) {
+ ATH_MSG_WARNING("PixelModuleDesign was not retrieved in function 'addNNInformation'");
+ return;
+ }
+ const std::vector<Identifier>& rdos = pixelCluster->rdoList();
+
+ const std::vector<float>& chList = pixelCluster->chargeList();
+ const std::vector<int>& totList = pixelCluster->totList();
+
+ ATH_MSG_VERBOSE( "Number of RDOs: " << rdos.size() );
+ ATH_MSG_VERBOSE( "Number of charges: " << chList.size() );
+ ATH_MSG_VERBOSE( "Number of TOT: " << totList.size() );
+
+
+ //Calculate the centre of the cluster
+ int phiPixelIndexMin, phiPixelIndexMax, etaPixelIndexMin, etaPixelIndexMax;
+ InDetDD::SiCellId cellIdWeightedPosition= getCellIdWeightedPosition( pixelCluster, &phiPixelIndexMin, &phiPixelIndexMax, &etaPixelIndexMin, &etaPixelIndexMax);
+
+ if (!cellIdWeightedPosition.isValid())
+ {
+ ATH_MSG_WARNING( "Weighted position is on invalid CellID." );
+ }
+
+ int etaPixelIndexWeightedPosition=cellIdWeightedPosition.etaIndex();
+ int phiPixelIndexWeightedPosition=cellIdWeightedPosition.phiIndex();
+
+
+ ATH_MSG_DEBUG(" weighted pos phiPixelIndex: " << phiPixelIndexWeightedPosition << " etaPixelIndex: " << etaPixelIndexWeightedPosition );
+
+ // SiLocalPosition PixelModuleDesign::positionFromColumnRow(const int column, const int row) const;
+ //
+ // Given row and column index of diode, returns position of diode center
+ // ALTERNATIVE/PREFERED way is to use localPositionOfCell(const SiCellId & cellId) or
+ // rawLocalPositionOfCell method in SiDetectorElement.
+ // DEPRECATED (but used in numerous places)
+ //
+ // Comment by Hide (referring the original comment in the code) : 2015-02-04
+ // I automatically replaced column to etaPixelIndex and row to phiPixelIndex here. It was bofore:
+ // InDetDD::SiLocalPosition siLocalPosition( design->positionFromColumnRow(columnWeightedPosition,rowWeightedPosition) );
+ //
+ // Then I assume the argument of column/row in this function is in offline manner, not the real hardware column/row.
+ //
+ InDetDD::SiLocalPosition w = design->positionFromColumnRow(etaPixelIndexWeightedPosition,phiPixelIndexWeightedPosition);
+
+
+ double localEtaPixelIndexWeightedPosition = w.xEta();
+ double localPhiPixelIndexWeightedPosition = w.xPhi();
+
+ int centralIndexX=(sizeX-1)/2;
+ int centralIndexY=(sizeY-1)/2;
+
+
+
+ // Check to see if the cluster is too big for the NN
+
+ if (abs(phiPixelIndexWeightedPosition-phiPixelIndexMin)>centralIndexX ||
+ abs(phiPixelIndexWeightedPosition-phiPixelIndexMax)>centralIndexX)
+ {
+ ATH_MSG_DEBUG(" Cluster too large phiPixelIndexMin " << phiPixelIndexMin << " phiPixelIndexMax " << phiPixelIndexMax << " centralX " << centralIndexX);
+ //return;
+ }
+
+ if (abs(etaPixelIndexWeightedPosition-etaPixelIndexMin)>centralIndexY ||
+ abs(etaPixelIndexWeightedPosition-etaPixelIndexMax)>centralIndexY)
+ {
+ ATH_MSG_DEBUG(" Cluster too large etaPixelIndexMin" << etaPixelIndexMin << " etaPixelIndexMax " << etaPixelIndexMax << " centralY " << centralIndexY);
+ //return;
+ }
+
+ std::vector< std::vector<float> > matrixOfToT (sizeX, std::vector<float>(sizeY,0) );
+ std::vector< std::vector<float> > matrixOfCharge(sizeX, std::vector<float>(sizeY,0));
+ std::vector<float> vectorOfPitchesY(sizeY,0.4);
+
+
+ //Itererate over all elements hits in the cluster and fill the charge and tot matrices
+ std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
+ std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
+ auto charge = chList.begin();
+ auto tot = totList.begin();
+
+ ATH_MSG_VERBOSE(" Putting together the n. " << rdos.size() << " rdos into a matrix.");
+
+ for (; rdosBegin!= rdosEnd; ++rdosBegin)
+ {
+
+ Identifier rId = *rdosBegin;
+ int absphiPixelIndex = m_PixelHelper->phi_index(rId)-phiPixelIndexWeightedPosition + centralIndexX;
+ int absetaPixelIndex = m_PixelHelper->eta_index(rId)-etaPixelIndexWeightedPosition + centralIndexY;
+
+ ATH_MSG_VERBOSE( " Phi Index: " << m_PixelHelper->phi_index(rId) << " absphiPixelIndex: " << absphiPixelIndex << " eta Idx: " << m_PixelHelper->eta_index(rId) << " absetaPixelIndex: " << absetaPixelIndex << " charge " << *charge );
+
+ if (absphiPixelIndex <0 || absphiPixelIndex >= (int)sizeX)
+ {
+ ATH_MSG_DEBUG(" problem with index: " << absphiPixelIndex << " min: " << 0 << " max: " << sizeX);
+ continue;
+ }
+
+ if (absetaPixelIndex <0 || absetaPixelIndex >= (int)sizeY)
+ {
+ ATH_MSG_DEBUG(" problem with index: " << absetaPixelIndex << " min: " << 0 << " max: " << sizeY);
+ continue;
+ }
+
+ InDetDD::SiCellId cellId = de->cellIdFromIdentifier(*rdosBegin);
+ InDetDD::SiDiodesParameters diodeParameters = design->parameters(cellId);
+ float pitchY = diodeParameters.width().xEta();
+
+ if ( (not totList.empty()) && tot != totList.end()) {
+ matrixOfToT[absphiPixelIndex][absetaPixelIndex] =*tot;
+ ++tot;
+ } else matrixOfToT[absphiPixelIndex][absetaPixelIndex] = -1;
+
+ if ( (not chList.empty()) && charge != chList.end()){
+ matrixOfCharge[absphiPixelIndex][absetaPixelIndex]=*charge;
+ ++charge;
+ } else matrixOfCharge[absphiPixelIndex][absetaPixelIndex] = -1;
+
+ if (pitchY > 0.1)
+ {
+ vectorOfPitchesY[absetaPixelIndex]=pitchY;
+ }
+ }//end iteration on rdos
+
+
+ ATH_MSG_VERBOSE( " End RDO LOOP " );
+
+ // Using the centre of the module and beam spot calculate
+ // the incidence angles of the tracks
+ const Amg::Vector2D& prdLocPos = pixelCluster->localPosition();
+ InDetDD::SiLocalPosition centroid(prdLocPos);
+
+ Amg::Vector3D globalPos = de->globalPosition(centroid);
+ Amg::Vector3D trackDir = globalPos; // - beamSpotPosition;
+ trackDir.normalize();
+
+ Amg::Vector3D module_normal = de->normal();
+ Amg::Vector3D module_phiax = de->phiAxis();
+ Amg::Vector3D module_etaax = de->etaAxis();
+
+ // Calculate the phi incidence angle
+ float trkphicomp = trackDir.dot(module_phiax);
+ float trketacomp = trackDir.dot(module_etaax);
+ float trknormcomp = trackDir.dot(module_normal);
+ double bowphi = atan2(trkphicomp,trknormcomp);
+ double boweta = atan2(trketacomp,trknormcomp);
+ double tanl = de->getTanLorentzAnglePhi();
+ if(bowphi > TMath::Pi()/2) bowphi -= TMath::Pi();
+ if(bowphi < -TMath::Pi()/2) bowphi += TMath::Pi();
+ int readoutside = design->readoutSide();
+ double angle = atan(tan(bowphi)-readoutside*tanl);
+
+
+ // Calculate the theta incidence angle
+ ATH_MSG_VERBOSE( " Angle theta bef corr: " << boweta );
+ if (boweta>TMath::Pi()/2.) boweta-=TMath::Pi();
+ if (boweta<-TMath::Pi()/2.) boweta+=TMath::Pi();
+
+
+ ATH_MSG_VERBOSE(" Angle phi: " << angle << " theta: " << boweta );
+ ATH_MSG_VERBOSE(" PhiPixelIndexWeightedPosition: " << phiPixelIndexWeightedPosition << " EtaPixelIndexWeightedPosition: " << etaPixelIndexWeightedPosition );
+
+ // store the matrixOfToT in a vector
+ std::vector<float> vectorOfCharge(sizeX*sizeY,0);
+ std::vector<float> vectorOfToT(sizeX*sizeY,0);
+ int counter(0);
+ for (unsigned int u=0;u<sizeX;u++)
+ {
+ for (unsigned int s=0;s<sizeY;s++)
+ {
+ vectorOfToT[counter] = matrixOfToT[u][s];
+ vectorOfCharge[counter] = matrixOfCharge[u][s];
+ ++counter;
+ }
+ }
+
+ ATH_MSG_VERBOSE( "matrixOfToT converted in a std::vector<float> " );
+
+ ATH_MSG_VERBOSE( "... and saved " );
+ // Add information to xAOD
+ xprd->auxdata< int >("NN_sizeX") = sizeX;
+ xprd->auxdata< int >("NN_sizeY") = sizeY;
+
+ xprd->auxdata< float >("NN_phiBS") = angle;
+ xprd->auxdata< float >("NN_thetaBS") = boweta;
+
+ xprd->auxdata< std::vector<float> >("NN_matrixOfToT") = vectorOfToT;
+ xprd->auxdata< std::vector<float> >("NN_matrixOfCharge") = vectorOfCharge;
+ xprd->auxdata< std::vector<float> >("NN_vectorOfPitchesY") = vectorOfPitchesY;
+
+
+ xprd->auxdata< int >("NN_etaPixelIndexWeightedPosition") = etaPixelIndexWeightedPosition;
+ xprd->auxdata< int >("NN_phiPixelIndexWeightedPosition") = phiPixelIndexWeightedPosition;
+
+ xprd->auxdata< float >("NN_localEtaPixelIndexWeightedPosition") = localEtaPixelIndexWeightedPosition;
+ xprd->auxdata< float >("NN_localPhiPixelIndexWeightedPosition") = localPhiPixelIndexWeightedPosition;
+
+ ATH_MSG_VERBOSE( "NN training Written" );
+}
+
+void PixelPrepDataToxAOD::addNNTruthInfo( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* pixelCluster,
+ const SiHitCollection* sihitCollection ) const
+{
+
+ std::vector<SiHit> matchingHits = findAllHitsCompatibleWithCluster( pixelCluster, sihitCollection );
+
+ unsigned int numberOfSiHits = matchingHits.size();
+
+ std::vector<float> positionsX(numberOfSiHits,0);
+ std::vector<float> positionsY(numberOfSiHits,0);
+
+ std::vector<float> positions_indexX(numberOfSiHits,0);
+ std::vector<float> positions_indexY(numberOfSiHits,0);
+
+ std::vector<float> theta(numberOfSiHits,0);
+ std::vector<float> phi(numberOfSiHits,0);
+
+ std::vector<int> barcode(numberOfSiHits,0);
+ std::vector<int> pdgid(numberOfSiHits,0);
+ std::vector<float> chargeDep(numberOfSiHits,0);
+ std::vector<float> truep(numberOfSiHits,0);
+
+ std::vector<float> pathlengthX(numberOfSiHits,0);
+ std::vector<float> pathlengthY(numberOfSiHits,0);
+ std::vector<float> pathlengthZ(numberOfSiHits,0);
+
+ std::vector<int> motherBarcode(numberOfSiHits,0);
+ std::vector<int> motherPdgid(numberOfSiHits,0);
+
+
+
+ // Check if we have detector element -- needed to find the local position of the SiHits
+ const InDetDD::SiDetectorElement* de = pixelCluster->detectorElement();
+ if(!de)
+ return;
+
+ InDetDD::SiCellId cellIdWeightedPosition = getCellIdWeightedPosition( pixelCluster );
+
+ const InDetDD::PixelModuleDesign* design(dynamic_cast<const InDetDD::PixelModuleDesign*>(&de->design()));
+ if (not design) {
+ ATH_MSG_WARNING("PixelModuleDesign was not retrieved in function 'addNNTruthInfo'");
+ return;
+ }
+ // lorentz shift correction
+ double shift = de->getLorentzCorrection();
+
+ unsigned hitNumber(0);
+ for( auto& siHit : matchingHits ){
+
+ HepGeom::Point3D<double> averagePosition = (siHit.localStartPosition() + siHit.localEndPosition()) * 0.5;
+
+ ATH_MSG_VERBOSE("Truth Part X: " << averagePosition.y() << " shift " << shift << " Y: " << averagePosition.z() );
+
+ // position lorentz shift corrected
+ float YposC = averagePosition.y()-shift;
+
+ if (fabs(YposC)>design->width()/2 &&
+ fabs(averagePosition.y())<design->width()/2)
+ {
+ if (YposC>design->width()/2)
+ {
+ YposC=design->width()/2-1e-6;
+ } else if (YposC<-design->width()/2)
+ {
+ YposC=-design->width()/2+1e-6;
+ }
+ }
+
+ positionsX[hitNumber] = YposC;
+ positionsY[hitNumber] = averagePosition.z();
+
+ HepGeom::Point3D<double> deltaPosition = siHit.localEndPosition() - siHit.localStartPosition();
+
+ pathlengthX[hitNumber] = deltaPosition.y();
+ pathlengthY[hitNumber] = deltaPosition.z();
+ pathlengthZ[hitNumber] = deltaPosition.x();
+
+
+
+ InDetDD::SiLocalPosition siLocalTruthPosition(averagePosition.z(),YposC ) ;
+ InDetDD::SiCellId cellIdOfTruthPosition =design->cellIdOfPosition(siLocalTruthPosition);
+
+ int truthEtaIndex = cellIdOfTruthPosition.etaIndex();
+ int truthPhiIndex = cellIdOfTruthPosition.phiIndex();
+
+ InDetDD::SiDiodesParameters diodeParameters = design->parameters(cellIdOfTruthPosition);
+ double pitchY = diodeParameters.width().xEta();
+ double pitchX = diodeParameters.width().xPhi();
+
+ // pixel center
+ // SiLocalPosition PixelModuleDesign::positionFromColumnRow(const int column, const int row) const;
+ //
+ // Given row and column index of diode, returns position of diode center
+ // ALTERNATIVE/PREFERED way is to use localPositionOfCell(const SiCellId & cellId) or
+ // rawLocalPositionOfCell method in SiDetectorElement.
+ // DEPRECATED (but used in numerous places)
+ //
+ // Comment by Hide (referring the original comment in the code) : 2015-02-04
+ // I automatically replaced column to etaPixelIndex and row to phiPixelIndex here. It was bofore:
+ // InDetDD::SiLocalPosition siLocalPosition( design->positionFromColumnRow(truthColumn,truthRow) );
+ //
+ // Then I assume the argument of column/row in this function is in offline manner, not the real hardware column/row.
+ //
+ InDetDD::SiLocalPosition siLocalPositionCenter(design->positionFromColumnRow(truthEtaIndex,truthPhiIndex));
+ double pixelCenterY = siLocalPositionCenter.xEta();
+ double pixelCenterX = siLocalPositionCenter.xPhi();
+
+
+ // truth index
+ double truthIndexY = truthEtaIndex + (averagePosition.z() - pixelCenterY)/pitchY;
+ double truthIndexX = truthPhiIndex + (YposC - pixelCenterX)/pitchX;
+
+ positions_indexX[hitNumber] = truthIndexX - cellIdWeightedPosition.phiIndex();
+ positions_indexY[hitNumber] = truthIndexY - cellIdWeightedPosition.etaIndex();
+
+ HepGeom::Point3D<double> diffPositions = (siHit.localEndPosition() - siHit.localStartPosition());
+ double bowphi = atan2( diffPositions.y(), diffPositions.x() );
+
+
+ //Truth Track incident angle theta
+ theta[hitNumber] = atan2(diffPositions.z() ,diffPositions.x());
+ //Truth track incident angle phi -- correct for lorentz angle
+ float tanlorentz = de->getTanLorentzAnglePhi();
+ int readoutside = design->readoutSide();
+ phi[hitNumber] = atan(tan(bowphi)-readoutside*tanlorentz);
+
+ if(siHit.particleLink().isValid()){
+ barcode[hitNumber] = siHit.particleLink().barcode();
+
+ auto particle = siHit.particleLink();
+ pdgid[hitNumber] = particle->pdg_id();
+ truep[hitNumber] = particle->momentum().rho();
+ if ( particle->production_vertex() ){
+ auto vertex = particle->production_vertex();
+ if( vertex->particles_in_const_begin() != vertex->particles_in_const_end() ){
+ motherBarcode[hitNumber] = (*vertex->particles_in_const_begin())->barcode();
+ motherPdgid[hitNumber] = (*vertex->particles_in_const_begin())->pdg_id();
+ }
+ }
+ }
+ chargeDep[hitNumber] = siHit.energyLoss() ;
+
+ ++hitNumber;
+ }
+
+
+ xprd->auxdata< std::vector<float> >("NN_positionsX") = positionsX;
+ xprd->auxdata< std::vector<float> >("NN_positionsY") = positionsY;
+
+ xprd->auxdata< std::vector<float> >("NN_positions_indexX") = positions_indexX;
+ xprd->auxdata< std::vector<float> >("NN_positions_indexY") = positions_indexY;
+
+ xprd->auxdata< std::vector<float> >("NN_theta") = theta;
+ xprd->auxdata< std::vector<float> >("NN_phi") = phi;
+
+ xprd->auxdata< std::vector<int> >("NN_barcode") = barcode;
+ xprd->auxdata< std::vector<int> >("NN_pdgid") = pdgid;
+ xprd->auxdata< std::vector<float> >("NN_energyDep") = chargeDep;
+ xprd->auxdata< std::vector<float> >("NN_trueP") = truep;
+
+ xprd->auxdata< std::vector<int> >("NN_motherBarcode") = motherBarcode;
+ xprd->auxdata< std::vector<int> >("NN_motherPdgid") = motherPdgid;
+
+
+ xprd->auxdata< std::vector<float> >("NN_pathlengthX") = pathlengthX;
+ xprd->auxdata< std::vector<float> >("NN_pathlengthY") = pathlengthY;
+ xprd->auxdata< std::vector<float> >("NN_pathlengthZ") = pathlengthZ;
+
+
+}
+
+
+
+
+InDetDD::SiCellId PixelPrepDataToxAOD::getCellIdWeightedPosition( const InDet::PixelCluster* pixelCluster,
+ int *rphiPixelIndexMin,
+ int *rphiPixelIndexMax,
+ int *retaPixelIndexMin,
+ int *retaPixelIndexMax ) const
+{
+
+ const InDetDD::SiDetectorElement* de = pixelCluster->detectorElement();
+ if (de==0) {
+ ATH_MSG_ERROR("Could not get detector element");
+ return InDetDD::SiCellId();
+ }
+
+ const InDetDD::PixelModuleDesign* design(dynamic_cast<const InDetDD::PixelModuleDesign*>(&de->design()));
+ if (not design) {
+ ATH_MSG_WARNING("PixelModuleDesign was not retrieved in function 'getCellIdWeightedPosition'");
+ return InDetDD::SiCellId();
+ }
+ const std::vector<Identifier>& rdos = pixelCluster->rdoList();
+
+ ATH_MSG_VERBOSE( "Number of RDOs: " << rdos.size() );
+ const std::vector<float>& chList = pixelCluster->chargeList();
+
+ ATH_MSG_VERBOSE( "Number of charges: " << chList.size() );
+ std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
+ std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
+
+ auto charge = chList.begin();
+
+ InDetDD::SiLocalPosition sumOfWeightedPositions(0,0,0);
+ double sumOfCharge=0;
+
+ int phiPixelIndexMin = 99999;
+ int phiPixelIndexMax = -99999;
+ int etaPixelIndexMin = 99999;
+ int etaPixelIndexMax = -99999;
+
+ for (; rdosBegin!= rdosEnd; ++rdosBegin, ++charge)
+ {
+
+ Identifier rId = *rdosBegin;
+ int phiPixelIndex = m_PixelHelper->phi_index(rId);
+ int etaPixelIndex = m_PixelHelper->eta_index(rId);
+
+ ATH_MSG_VERBOSE(" Adding pixel phiPixelIndex: " << phiPixelIndex << " etaPixelIndex: " << etaPixelIndex << " charge: " << *charge );
+
+ // SiLocalPosition PixelModuleDesign::positionFromColumnRow(const int column, const int row) const;
+ //
+ // Given row and column index of diode, returns position of diode center
+ // ALTERNATIVE/PREFERED way is to use localPositionOfCell(const SiCellId & cellId) or
+ // rawLocalPositionOfCell method in SiDetectorElement.
+ // DEPRECATED (but used in numerous places)
+ //
+ // Comment by Hide (referring the original comment in the code): 2015-02-04
+ // I automatically replaced column to etaPixelIndex and row to phiPixelIndex here. It was bofore:
+ // InDetDD::SiLocalPosition siLocalPosition( design->positionFromColumnRow(column,row) );
+ //
+ // Then I assume the argument of column/row in this function is in offline manner, not the real hardware column/row.
+ //
+ InDetDD::SiLocalPosition siLocalPosition( design->positionFromColumnRow(etaPixelIndex,phiPixelIndex) );
+ ATH_MSG_VERBOSE ( "Local Position: Row = " << siLocalPosition.xRow() << ", Col = " << siLocalPosition.xColumn() );
+
+ sumOfWeightedPositions += (*charge)*siLocalPosition;
+ sumOfCharge += (*charge);
+
+ if (phiPixelIndex < phiPixelIndexMin)
+ phiPixelIndexMin = phiPixelIndex;
+
+ if (phiPixelIndex > phiPixelIndexMax)
+ phiPixelIndexMax = phiPixelIndex;
+
+ if (etaPixelIndex < etaPixelIndexMin)
+ etaPixelIndexMin = etaPixelIndex;
+
+ if (etaPixelIndex > etaPixelIndexMax)
+ etaPixelIndexMax = etaPixelIndex;
+
+ }
+ sumOfWeightedPositions /= sumOfCharge;
+
+ ATH_MSG_VERBOSE ( "Wighted position: Row = " << sumOfWeightedPositions.xRow() << ", Col = " << sumOfWeightedPositions.xColumn() );
+
+ if(rphiPixelIndexMin) *rphiPixelIndexMin = phiPixelIndexMin;
+ if(rphiPixelIndexMax) *rphiPixelIndexMax = phiPixelIndexMax;
+ if(retaPixelIndexMin) *retaPixelIndexMin = etaPixelIndexMin;
+ if(retaPixelIndexMax) *retaPixelIndexMax = etaPixelIndexMax;
+
+ //what you want to know is simple:
+ //just the phiPixelIndex and etaPixelIndex of this average position!
+
+ InDetDD::SiCellId cellIdWeightedPosition=design->cellIdOfPosition(sumOfWeightedPositions);
+
+
+ return cellIdWeightedPosition;
+
+}
+
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// Finalize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode PixelPrepDataToxAOD::finalize()
+{
+ return StatusCode::SUCCESS;
+}
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.h b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.h
new file mode 100644
index 00000000000..98d86b9aba7
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/PixelPrepDataToxAOD.h
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PixelPrepDataToxAOD.h
+// Header file for class PixelPrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#ifndef PIXELPREPDATATOXAOD_H
+#define PIXELPREPDATATOXAOD_H
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+#include "InDetSimEvent/SiHitCollection.h"
+#include "xAODTracking/TrackMeasurementValidation.h"
+
+
+#include <string>
+
+class PixelID;
+class SiHit;
+class InDetSimDataCollection;
+
+
+namespace InDet
+{
+ class PixelCluster;
+}
+
+namespace InDetDD
+{
+ class SiCellId;
+}
+
+
+
+class PixelPrepDataToxAOD : public AthAlgorithm {
+
+public:
+ // Constructor with parameters:
+ PixelPrepDataToxAOD(const std::string &name,ISvcLocator *pSvcLocator);
+
+ // Basic algorithm methods:
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+
+private:
+
+ void addSDOInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* prd,
+ const InDetSimDataCollection* sdoCollection ) const;
+
+
+ void addSiHitInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* prd,
+ const SiHitCollection* collection) const;
+
+ std::vector<SiHit> findAllHitsCompatibleWithCluster(const InDet::PixelCluster* prd,
+ const SiHitCollection* collection) const;
+
+
+ void addNNTruthInfo( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* prd,
+ const SiHitCollection* collection ) const;
+
+ void addNNInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* pixelCluster,
+ const unsigned int SizeX,
+ const unsigned int SizeY ) const;
+
+ void addRdoInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::PixelCluster* pixelCluster) const;
+
+
+
+ InDetDD::SiCellId getCellIdWeightedPosition( const InDet::PixelCluster* pixelCluster,
+ int *rrowMin = 0,
+ int *rrowMax = 0,
+ int *rcolMin = 0,
+ int *rcolMax = 0 ) const;
+
+
+ const PixelID *m_PixelHelper;
+ std::string m_clustercontainer;
+ std::string m_SDOcontainer;
+ std::string m_sihitContainer;
+ std::string m_multiTruth;
+
+ bool m_useTruthInfo;
+ bool m_writeSDOs;
+ bool m_writeSiHits;
+ bool m_writeNNinformation;
+ bool m_writeRDOinformation;
+
+ // -- Private members
+ bool m_firstEventWarnings;
+
+};
+
+
+#endif
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.cxx b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.cxx
new file mode 100644
index 00000000000..44b6d41c0ad
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.cxx
@@ -0,0 +1,590 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// SCT_PrepDataToxAOD.cxx
+// Implementation file for class SCT_PrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#include "SCT_PrepDataToxAOD.h"
+
+#include "InDetPrepRawData/SCT_ClusterContainer.h"
+
+#include "xAODTracking/TrackMeasurementValidationContainer.h"
+#include "xAODTracking/TrackMeasurementValidationAuxContainer.h"
+
+#include "Identifier/Identifier.h"
+#include "InDetIdentifier/SCT_ID.h"
+
+#include "InDetRawData/SCT_RDO_Container.h"
+#include "InDetRawData/SCT_RDO_Collection.h"
+
+#include "TrkTruthData/PRD_MultiTruthCollection.h"
+#include "HepMC/GenParticle.h"
+#include "InDetSimEvent/SiHit.h"
+#include "InDetSimData/InDetSimDataCollection.h"
+
+#include "CLHEP/Geometry/Point3D.h"
+
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// Constructor with parameters:
+//
+/////////////////////////////////////////////////////////////////////
+SCT_PrepDataToxAOD::SCT_PrepDataToxAOD(const std::string &name, ISvcLocator *pSvcLocator) :
+ AthAlgorithm(name,pSvcLocator),
+ m_incidentSvc("IncidentSvc", name),
+ m_SCTHelper(0),
+ m_firstEventWarnings(true)
+{
+ // --- Steering and configuration flags
+ declareProperty("UseTruthInfo", m_useTruthInfo=false);
+ declareProperty("WriteRDOinformation", m_writeRDOinformation =true);
+ declareProperty("WriteSDOs", m_writeSDOs = true);
+ declareProperty("WriteSiHits", m_writeSiHits = true);
+
+ // --- Configuration keys
+ declareProperty("SiClusterContainer", m_clustercontainer = "SCT_Clusters");
+ declareProperty("MC_SDOs", m_SDOcontainer = "SCT_SDO_Map");
+ declareProperty("MC_Hits", m_sihitContainer = "SCT_Hits");
+ declareProperty("PRD_MultiTruth", m_multiTruth = "PRD_MultiTruthSCT");
+
+ // --- Services and Tools
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Initialize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode SCT_PrepDataToxAOD::initialize()
+{
+
+ CHECK ( detStore()->retrieve(m_SCTHelper, "SCT_ID") );
+
+ CHECK ( m_incidentSvc.retrieve() );
+ // register to the incident service:
+ m_incidentSvc->addListener( this, IncidentType::EndEvent);
+
+ //make sure we don't write what we don't have
+ if (not m_useTruthInfo) {
+ m_writeSDOs = false;
+ m_writeSiHits = false;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Execute method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode SCT_PrepDataToxAOD::execute()
+{
+ const SCT_RDO_Container * rdoContainer(0);
+
+ // the cluster ambiguity map
+ if ( m_writeRDOinformation ) {
+ if(!evtStore()->contains<SCT_RDO_Container>("SCT_RDOs")){
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("RDO ASSOC: No SCT RDO container in StoreGate");
+ }
+ }
+ else {
+ if(evtStore()->retrieve(rdoContainer,"SCT_RDOs").isFailure()) {
+ ATH_MSG_WARNING( "Failed to retrieve SCT RDO container" );
+ }
+ }
+
+ if ( rdoContainer != 0){
+
+ // get all the RIO_Collections in the container
+
+ for( auto& collection: *rdoContainer ){
+
+ //get all the RDOs in the collection
+ for (auto& rdo : *collection) {
+
+ if ( !rdo) {
+ ATH_MSG_WARNING( "Null SCT RDO. Skipping it");
+ continue;
+ }
+
+ Identifier rdoId = rdo->identify();
+
+ m_IDtoRAWDataMap.insert( std::pair< Identifier, const SCT_RDORawData*>( rdoId, rdo ) );
+ } // collection
+ } // container
+ } // Have container;
+ }
+ ATH_MSG_DEBUG("Size of RDO map is "<<m_IDtoRAWDataMap.size());
+
+ // Mandatory. This is needed and required if this algorithm is scheduled.
+ const InDet::SCT_ClusterContainer* sctClusterContainer = 0;
+ if( evtStore()->retrieve(sctClusterContainer,m_clustercontainer).isFailure() ) {
+ ATH_MSG_ERROR("Cannot retrieve SCT PrepDataContainer " << m_clustercontainer);
+ return StatusCode::FAILURE;
+ }
+
+ // Optional. Normally only available in Hits files -- samples need to digitised and Hits need to be copied for this to work
+ const SiHitCollection* sihitCollection = 0;
+ if (m_writeSiHits) {
+ if (evtStore()->contains<SiHitCollection>(m_sihitContainer)) {
+ ATH_CHECK(evtStore()->retrieve(sihitCollection, m_sihitContainer));
+ } else {
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("SiHit collection not available (" << m_sihitContainer << "). Skipping although requested.");
+ sihitCollection = 0;
+ }
+ }
+ }
+
+ // Optional. On RDO
+ const InDetSimDataCollection* sdoCollection = 0;
+ if (m_writeSDOs) {
+ if ( evtStore()->contains<InDetSimDataCollection>(m_SDOcontainer) ) {
+ ATH_CHECK(evtStore()->retrieve(sdoCollection, m_SDOcontainer));
+ } else {
+ ATH_MSG_WARNING("SDO container not available (" << m_SDOcontainer << "). Skipping although requested.");
+ sdoCollection = 0;
+ }
+ }
+
+ // Optional. On ESD and AOD
+ const PRD_MultiTruthCollection* prdmtColl = 0;
+ if (m_useTruthInfo) {
+ if ( evtStore()->contains<PRD_MultiTruthCollection>(m_multiTruth) ) {
+ ATH_CHECK(evtStore()->retrieve(prdmtColl, m_multiTruth));
+ } else {
+ ATH_MSG_WARNING("MultiTruth container not available (" << m_multiTruth << "). Skipping although requested.");
+ prdmtColl = 0;
+ }
+ }
+
+
+ // Create the xAOD container and its auxiliary store:
+ xAOD::TrackMeasurementValidationContainer* xaod = new xAOD::TrackMeasurementValidationContainer();
+ CHECK( evtStore()->record( xaod, m_clustercontainer ) );
+ xAOD::TrackMeasurementValidationAuxContainer* aux = new xAOD::TrackMeasurementValidationAuxContainer();
+ CHECK( evtStore()->record( aux, m_clustercontainer + "Aux." ) );
+ xaod->setStore( aux );
+
+ std::vector<unsigned int>* offsets = new std::vector<unsigned int>( m_SCTHelper->wafer_hash_max(), 0 );
+ CHECK( evtStore()->record( offsets, m_clustercontainer + "Offsets" ) );
+
+ // Loop over the container
+ unsigned int counter(0);
+ for( const auto& clusterCollection : *sctClusterContainer){
+
+ //Fill Offset container
+ (*offsets)[clusterCollection->identifyHash()] = counter;
+
+ // skip empty collections
+ if( clusterCollection->empty() ) continue;
+
+
+ // loop over collection and convert to xAOD
+ for( auto& prd : *clusterCollection ){
+ ++counter;
+
+ Identifier clusterId = prd->identify();
+ if ( !clusterId.is_valid() ) {
+ ATH_MSG_WARNING("SCT cluster identifier is not valid!");
+ }
+
+ // create and add xAOD object
+ xAOD::TrackMeasurementValidation* xprd = new xAOD::TrackMeasurementValidation();
+ xaod->push_back(xprd);
+
+ //Set Identifier
+ xprd->setIdentifier( clusterId.get_compact() );
+
+ //Set Global Position
+ Amg::Vector3D gpos = prd->globalPosition();
+ xprd->setGlobalPosition(gpos.x(),gpos.y(),gpos.z());
+
+ //Set Local Position
+ const Amg::Vector2D& locpos = prd->localPosition();
+ float locY(0.);
+ float locX = locpos.x();
+ if ( !(std::isinf(locpos.y()) || std::isnan(locpos.y())) ){
+ if (locpos.y()>=1e-07)
+ locY=locpos.y();
+ } else {
+ locY = -9999.;
+ }
+
+ // Set local error matrix
+ xprd->setLocalPosition(locX,locY);
+
+ const Amg::MatrixX& localCov = prd->localCovariance();
+ if(localCov.size() == 1){
+ xprd->setLocalPositionError( localCov(0,0), 0., 0. );
+ } else if(localCov.size() == 4){
+ xprd->setLocalPositionError( localCov(0,0), localCov(1,1), localCov(0,1) );
+ } else {
+ xprd->setLocalPositionError(0.,0.,0.);
+ }
+
+ // Set vector of hit identifiers
+ std::vector< uint64_t > rdoIdentifierList;
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ rdoIdentifierList.push_back( hitIdentifier.get_compact() );
+ }
+ xprd->setRdoIdentifierList(rdoIdentifierList);
+
+ //Add SCT specific information
+ const InDet::SiWidth cw = prd->width();
+ xprd->auxdata<int>("size") = (int)cw.colRow()[0];
+
+ xprd->auxdata<int>("bec") = m_SCTHelper->barrel_ec(clusterId) ;
+ xprd->auxdata<int>("layer") = m_SCTHelper->layer_disk(clusterId) ;
+ xprd->auxdata<int>("phi_module") = m_SCTHelper->phi_module(clusterId) ;
+ xprd->auxdata<int>("eta_module") = m_SCTHelper->eta_module(clusterId) ;
+ xprd->auxdata<int>("side") = m_SCTHelper->side(clusterId) ;
+
+ // Add the Detector element ID -- not sure if needed as we have the informations above
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ uint64_t detElementId(0);
+ if(de){
+ Identifier detId = de->identify();
+ if ( detId.is_valid() ) {
+ detElementId = detId.get_compact();
+ }
+ }
+ xprd->auxdata<uint64_t>("detectorElementID") = detElementId;
+
+ //Add details about the individual hits
+ if(m_writeRDOinformation)
+ addRDOInformation(xprd, prd);
+
+
+ // Use the MultiTruth Collection to get a list of all true particle contributing to the cluster
+ if(prdmtColl){
+ std::vector<int> barcodes;
+ //std::pair<PRD_MultiTruthCollection::const_iterator,PRD_MultiTruthCollection::const_iterator>;
+ auto range = prdmtColl->equal_range(clusterId);
+ for (auto i = range.first; i != range.second; ++i) {
+ barcodes.push_back( i->second.barcode() );
+ }
+ xprd->auxdata< std::vector<int> >("truth_barcode") = barcodes;
+ }
+
+ // Use the SDO Collection to get a list of all true particle contributing to the cluster per readout element
+ // Also get the energy deposited by each true particle per readout element
+ if( sdoCollection ){
+ addSDOInformation( xprd, prd,sdoCollection);
+ }
+
+ // Now Get the most detailed truth from the SiHits
+ // Note that this could get really slow if there are a lot of hits and clusters
+ if( sihitCollection ){
+ addSiHitInformation( xprd, prd, sihitCollection);
+ }
+ }
+ }
+ ATH_MSG_DEBUG( " recorded SCT_PrepData objects: size " << xaod->size() );
+
+ m_firstEventWarnings = false; //disable one-time warnings
+
+ return StatusCode::SUCCESS;
+}
+
+
+
+void SCT_PrepDataToxAOD::addSDOInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::SCT_Cluster* prd,
+ const InDetSimDataCollection* sdoCollection ) const
+{
+ std::vector<int> sdo_word;
+ std::vector< std::vector< int > > sdo_depositsBarcode;
+ std::vector< std::vector< float > > sdo_depositsEnergy;
+ // find hit
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ auto pos = sdoCollection->find(hitIdentifier);
+ if( pos != sdoCollection->end() ) {
+ sdo_word.push_back( pos->second.word() ) ;
+
+ std::vector< int > sdoDepBC;
+ std::vector< float > sdoDepEnergy;
+ for( auto deposit : pos->second.getdeposits() ){
+ if(deposit.first){
+ sdoDepBC.push_back( deposit.first->barcode());
+ } else {
+ sdoDepBC.push_back( -1 );
+ }
+ sdoDepEnergy.push_back( deposit.second );
+ ATH_MSG_DEBUG(" SDO Energy Deposit " << deposit.second ) ;
+ }
+ sdo_depositsBarcode.push_back( sdoDepBC );
+ sdo_depositsEnergy.push_back( sdoDepEnergy );
+ }
+ }
+ xprd->auxdata< std::vector<int> >("sdo_words") = sdo_word;
+ xprd->auxdata< std::vector< std::vector<int> > >("sdo_depositsBarcode") = sdo_depositsBarcode;
+ xprd->auxdata< std::vector< std::vector<float> > >("sdo_depositsEnergy") = sdo_depositsEnergy;
+}
+
+
+
+void SCT_PrepDataToxAOD::addSiHitInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::SCT_Cluster* prd,
+ const SiHitCollection* sihitCollection) const
+
+
+{
+
+ std::vector<SiHit> matchingHits = findAllHitsCompatibleWithCluster( prd, sihitCollection );
+
+ int numHits = matchingHits.size();
+
+ std::vector<float> sihit_energyDeposit(numHits,0);
+ std::vector<float> sihit_meanTime(numHits,0);
+ std::vector<int> sihit_barcode(numHits,0);
+
+ std::vector<float> sihit_startPosX(numHits,0);
+ std::vector<float> sihit_startPosY(numHits,0);
+ std::vector<float> sihit_startPosZ(numHits,0);
+
+ std::vector<float> sihit_endPosX(numHits,0);
+ std::vector<float> sihit_endPosY(numHits,0);
+ std::vector<float> sihit_endPosZ(numHits,0);
+
+ int hitNumber(0);
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ if(de){
+ for ( auto sihit : matchingHits ) {
+ sihit_energyDeposit[hitNumber] = sihit.energyLoss() ;
+ sihit_meanTime[hitNumber] = sihit.meanTime() ;
+ sihit_barcode[hitNumber] = sihit.particleLink().barcode() ;
+
+ // Convert Simulation frame into reco frame
+ const HepGeom::Point3D<double>& startPos=sihit.localStartPosition();
+
+ Amg::Vector2D pos= de->hitLocalToLocal( startPos.z(), startPos.y() );
+ sihit_startPosX[hitNumber] = pos[0];
+ sihit_startPosY[hitNumber] = pos[1];
+ sihit_startPosZ[hitNumber] = startPos.x();
+
+
+ const HepGeom::Point3D<double>& endPos=sihit.localEndPosition();
+ pos= de->hitLocalToLocal( endPos.z(), endPos.y() );
+ sihit_endPosX[hitNumber] = pos[0];
+ sihit_endPosY[hitNumber] = pos[1];
+ sihit_endPosZ[hitNumber] = endPos.x();
+ ++hitNumber;
+ }
+ }
+
+ xprd->auxdata<std::vector<float> >("sihit_energyDeposit") = sihit_energyDeposit;
+ xprd->auxdata<std::vector<float> >("sihit_meanTime") = sihit_meanTime;
+ xprd->auxdata<std::vector<int> >("sihit_barcode") = sihit_barcode;
+
+ xprd->auxdata<std::vector<float> >("sihit_startPosX") = sihit_startPosX;
+ xprd->auxdata<std::vector<float> >("sihit_startPosY") = sihit_startPosY;
+ xprd->auxdata<std::vector<float> >("sihit_startPosZ") = sihit_startPosZ;
+
+ xprd->auxdata<std::vector<float> >("sihit_endPosX") = sihit_endPosX;
+ xprd->auxdata<std::vector<float> >("sihit_endPosY") = sihit_endPosY;
+ xprd->auxdata<std::vector<float> >("sihit_endPosZ") = sihit_endPosZ;
+
+
+}
+
+
+
+
+
+std::vector<SiHit> SCT_PrepDataToxAOD::findAllHitsCompatibleWithCluster( const InDet::SCT_Cluster* prd,
+ const SiHitCollection* collection) const
+{
+ ATH_MSG_VERBOSE( "Got " << collection->size() << " SiHits to look through" );
+ std::vector<SiHit> matchingHits;
+
+ // Check if we have detector element -- needed to find the local position of the SiHits
+ const InDetDD::SiDetectorElement* de = prd->detectorElement();
+ if(!de)
+ return matchingHits;
+
+ std::vector<const SiHit* > multiMatchingHits;
+
+ for ( const auto& siHit : *collection) {
+ // Check if it is a SCT hit
+ if( !siHit.isSCT() )
+ continue;
+
+ //Check if it is on the correct module
+ Identifier clusterId = prd->identify();
+
+ if( m_SCTHelper->barrel_ec(clusterId) != siHit.getBarrelEndcap() ||
+ m_SCTHelper->layer_disk(clusterId)!= siHit.getLayerDisk() ||
+ m_SCTHelper->phi_module(clusterId)!= siHit.getPhiModule() ||
+ m_SCTHelper->eta_module(clusterId)!= siHit.getEtaModule() ||
+ m_SCTHelper->side(clusterId) != siHit.getSide() )
+ continue;
+
+ // Now we have all hits in the module that match lets check to see if they match the cluster
+ // Must be within +/- 1 hits of any hit in the cluster to be included
+ ATH_MSG_DEBUG("Hit is on the same module");
+
+ HepGeom::Point3D<double> averagePosition = siHit.localStartPosition() + siHit.localEndPosition();
+ averagePosition *= 0.5;
+ Amg::Vector2D pos = de->hitLocalToLocal( averagePosition.z(), averagePosition.y() );
+ InDetDD::SiCellId diode = de->cellIdOfPosition(pos);
+
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ //if( abs( int(diode.etaIndex()) - m_SCTHelper->eta_index( hitIdentifier ) ) <=1 )
+ //if( abs( int(diode.phiIndex() - m_SCTHelper->phi_index( hitIdentifier ) ) <=1 )
+ ATH_MSG_DEBUG("Truth Strip " << diode.phiIndex() << " Cluster Strip " << m_SCTHelper->strip( hitIdentifier ) );
+
+ if( abs( int(diode.phiIndex()) - m_SCTHelper->strip( hitIdentifier ) ) <=1)
+ {
+ multiMatchingHits.push_back(&siHit);
+ break;
+ }
+ }
+ }
+
+
+ //Now we will now make 1 SiHit for each true particle if the SiHits "touch" other
+ std::vector<const SiHit* >::iterator siHitIter = multiMatchingHits.begin();
+ std::vector<const SiHit* >::iterator siHitIter2 = multiMatchingHits.begin();
+ ATH_MSG_DEBUG( "Found " << multiMatchingHits.size() << " SiHit " );
+ for ( ; siHitIter != multiMatchingHits.end(); siHitIter++) {
+ const SiHit* lowestXPos = *siHitIter;
+ const SiHit* highestXPos = *siHitIter;
+
+
+ // We will merge these hits
+ std::vector<const SiHit* > ajoiningHits;
+ ajoiningHits.push_back( *siHitIter );
+
+ siHitIter2 = siHitIter+1;
+ while ( siHitIter2 != multiMatchingHits.end() ) {
+ // Need to come from the same truth particle
+ if( (*siHitIter)->particleLink().barcode() != (*siHitIter2)->particleLink().barcode() ){
+ ++siHitIter2;
+ continue;
+ }
+
+ // Check to see if the SiHits are compatible with each other.
+ if (fabs((highestXPos->localEndPosition().x()-(*siHitIter2)->localStartPosition().x()))<0.00005 &&
+ fabs((highestXPos->localEndPosition().y()-(*siHitIter2)->localStartPosition().y()))<0.00005 &&
+ fabs((highestXPos->localEndPosition().z()-(*siHitIter2)->localStartPosition().z()))<0.00005 )
+ {
+ highestXPos = *siHitIter2;
+ ajoiningHits.push_back( *siHitIter2 );
+ // Dont use hit more than once
+ siHitIter2 = multiMatchingHits.erase( siHitIter2 );
+ }else if (fabs((lowestXPos->localStartPosition().x()-(*siHitIter2)->localEndPosition().x()))<0.00005 &&
+ fabs((lowestXPos->localStartPosition().y()-(*siHitIter2)->localEndPosition().y()))<0.00005 &&
+ fabs((lowestXPos->localStartPosition().z()-(*siHitIter2)->localEndPosition().z()))<0.00005)
+ {
+ lowestXPos = *siHitIter2;
+ ajoiningHits.push_back( *siHitIter2 );
+ // Dont use hit more than once
+ siHitIter2 = multiMatchingHits.erase( siHitIter2 );
+ } else {
+ ++siHitIter2;
+ }
+ }
+
+ if( ajoiningHits.size() == 0){
+ ATH_MSG_WARNING("This should really never happen");
+ continue;
+ }else if(ajoiningHits.size() == 1){
+ // Copy Si Hit ready to return
+ matchingHits.push_back( *ajoiningHits[0] );
+ continue;
+ } else {
+ // Build new SiHit and merge information together.
+ ATH_MSG_DEBUG("Merging " << ajoiningHits.size() << " SiHits together." );
+
+
+ float energyDep(0);
+ float time(0);
+ for( auto& siHit : ajoiningHits){
+ energyDep += siHit->energyLoss();
+ time += siHit->meanTime();
+ }
+ time /= (float)ajoiningHits.size();
+
+ matchingHits.push_back( SiHit(lowestXPos->localStartPosition(),
+ highestXPos->localEndPosition(),
+ energyDep,
+ time,
+ (*siHitIter)->particleLink().barcode(),
+ 1, // 0 for pixel 1 for SCT
+ (*siHitIter)->getBarrelEndcap(),
+ (*siHitIter)->getLayerDisk(),
+ (*siHitIter)->getEtaModule(),
+ (*siHitIter)->getPhiModule(),
+ (*siHitIter)->getSide() ) );
+ }
+ }
+
+
+ return matchingHits;
+
+}
+
+void SCT_PrepDataToxAOD::addRDOInformation(xAOD::TrackMeasurementValidation* xprd,
+ const InDet::SCT_Cluster* prd) const{
+
+
+ std::vector<int> strip;
+ std::vector<int> timebin;
+ std::vector<int> groupsize;
+
+ for( const auto &hitIdentifier : prd->rdoList() ){
+ auto result = m_IDtoRAWDataMap.find(hitIdentifier);
+ if( result != m_IDtoRAWDataMap.end() ){
+ const SCT_RDORawData *sctRdo = result->second;
+ const SCT3_RawData* rdo3 = dynamic_cast<const SCT3_RawData*>(sctRdo);
+ int tbin(-1);
+ int gs(-1);
+ if (rdo3!=0){
+ tbin = rdo3->getTimeBin();
+ gs = rdo3->getGroupSize();
+ }
+ timebin.push_back( tbin);
+ groupsize.push_back( gs);
+ strip.push_back(sctRdo->getStrip());
+ } else {
+ timebin.push_back( -1 );
+ strip.push_back( -1 );
+ groupsize.push_back( -1 );
+ }
+ }
+
+ xprd->auxdata< std::vector<int> >("rdo_strip") = strip;
+ xprd->auxdata< std::vector<int> >("rdo_timebin") = timebin;
+ xprd->auxdata< std::vector<int> >("rdo_groupsize") = groupsize;
+
+}
+
+void SCT_PrepDataToxAOD::handle(const Incident& inc) {
+
+ /// clear map of RDOs<->identifiers
+ if ( m_writeRDOinformation && inc.type() == IncidentType::EndEvent ){
+ ATH_MSG_VERBOSE("'EndEvent' incident caught. Refreshing Cache.");
+ m_IDtoRAWDataMap.clear();
+ }
+}
+
+
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// Finalize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode SCT_PrepDataToxAOD::finalize()
+{
+ return StatusCode::SUCCESS;
+}
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.h b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.h
new file mode 100644
index 00000000000..d9821e973d4
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/SCT_PrepDataToxAOD.h
@@ -0,0 +1,87 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// SCT_PrepDataToxAOD.h
+// Header file for class SCT_PrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#ifndef SCT_PREPDATATOXAOD_H
+#define SCT_PREPDATATOXAOD_H
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IIncidentListener.h"
+
+#include "InDetSimEvent/SiHitCollection.h"
+#include "xAODTracking/TrackMeasurementValidation.h"
+
+
+#include <string>
+
+class SCT_ID;
+class SiHit;
+class InDetSimDataCollection;
+class SCT_RDORawData;
+class Identifier;
+
+namespace InDet
+{
+ class SCT_Cluster;
+}
+
+
+class SCT_PrepDataToxAOD : public AthAlgorithm, virtual public IIncidentListener {
+
+public:
+ // Constructor with parameters:
+ SCT_PrepDataToxAOD(const std::string &name,ISvcLocator *pSvcLocator);
+
+ // Basic algorithm methods:
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+
+ virtual void handle(const Incident& inc);
+
+private:
+
+
+ void addSDOInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::SCT_Cluster* prd,
+ const InDetSimDataCollection* sdoCollection ) const;
+
+
+ void addSiHitInformation( xAOD::TrackMeasurementValidation* xprd,
+ const InDet::SCT_Cluster* prd,
+ const SiHitCollection* collection) const;
+
+ std::vector<SiHit> findAllHitsCompatibleWithCluster(const InDet::SCT_Cluster* prd,
+ const SiHitCollection* collection) const;
+
+ void addRDOInformation(xAOD::TrackMeasurementValidation*, const InDet::SCT_Cluster*) const;
+
+
+ ServiceHandle<IIncidentSvc> m_incidentSvc; //!< IncidentSvc to catch begin of event and end of envent
+
+ const SCT_ID *m_SCTHelper;
+ std::string m_clustercontainer;
+ std::string m_SDOcontainer;
+ std::string m_sihitContainer;
+ std::string m_multiTruth;
+
+ bool m_useTruthInfo;
+ bool m_writeRDOinformation;
+ bool m_writeSDOs;
+ bool m_writeSiHits;
+
+ std::map< Identifier, const SCT_RDORawData* > m_IDtoRAWDataMap;
+
+ // --- private members
+ bool m_firstEventWarnings;
+
+};
+
+
+#endif
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.cxx b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.cxx
new file mode 100644
index 00000000000..72e65dc383d
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.cxx
@@ -0,0 +1,342 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// TRT_PrepDataToxAOD.cxx
+// Implementation file for class TRT_PrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#include "TRT_PrepDataToxAOD.h"
+
+#include "InDetPrepRawData/TRT_DriftCircleContainer.h"
+//#include "EventPrimitives/EventPrimitivesHelpers.h"
+
+#include "xAODTracking/TrackMeasurementValidation.h"
+#include "xAODTracking/TrackMeasurementValidationContainer.h"
+#include "xAODTracking/TrackMeasurementValidationAuxContainer.h"
+
+
+#include "Identifier/Identifier.h"
+#include "InDetIdentifier/TRT_ID.h"
+
+
+#include "TrkSurfaces/Surface.h"
+#include "TRT_ConditionsServices/ITRT_StrawNeighbourSvc.h"
+#include "TRT_DriftFunctionTool/ITRT_DriftFunctionTool.h"
+
+#include "TrkTruthData/PRD_MultiTruthCollection.h"
+#include "InDetSimEvent/TRTUncompressedHitCollection.h"
+#include "InDetSimData/InDetSimDataCollection.h"
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// Constructor with parameters:
+//
+/////////////////////////////////////////////////////////////////////
+TRT_PrepDataToxAOD::TRT_PrepDataToxAOD(const std::string &name, ISvcLocator *pSvcLocator) :
+ AthAlgorithm(name,pSvcLocator),
+ m_driftFunctionTool("TRT_DriftFunctionTool"),
+ m_trtcaldbSvc("TRT_CalDbSvc", name),
+ m_neighbourSvc("TRT_StrawNeighbourSvc", name),
+ m_TRTStrawSummarySvc("InDetTRTStrawStatusSummarySvc",name),
+ m_TRTHelper(0),
+ m_firstEventWarnings(true)
+{
+ // --- Steering and configuration flags
+ declareProperty("UseTruthInfo", m_useTruthInfo=false);
+ declareProperty("WriteSDOs", m_writeSDOs = true);
+
+ // --- Configuration keys
+ declareProperty("DriftCircleContainer", m_driftcirclecontainer="TRT_DriftCircles");
+ declareProperty("MC_TRTUncompressedHit", m_SDOcontainer="TRT_SDO_Map");
+ declareProperty("PRD_MultiTruth", m_multiTruth="PRD_MultiTruthTRT");
+
+ // --- Services and Tools
+ declareProperty("TRTDriftFunctionTool", m_driftFunctionTool);
+ declareProperty("TRTCalDbSvc", m_trtcaldbSvc);
+ declareProperty("NeighbourSvc", m_neighbourSvc);
+ declareProperty("TRTStrawSummarySvc", m_TRTStrawSummarySvc);
+ declareProperty("NeighbourSvc", m_neighbourSvc);
+
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Initialize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode TRT_PrepDataToxAOD::initialize()
+{
+ ATH_MSG_DEBUG("Initialize");
+
+ // --- Retrieve services and tools
+ CHECK ( detStore()->retrieve(m_TRTHelper, "TRT_ID") );
+
+ CHECK ( m_neighbourSvc.retrieve() );
+
+ CHECK ( m_trtcaldbSvc.retrieve() );
+
+ CHECK ( m_TRTStrawSummarySvc.retrieve() );
+
+ CHECK( m_driftFunctionTool.retrieve() );
+
+ return StatusCode::SUCCESS;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Execute method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode TRT_PrepDataToxAOD::execute()
+{
+ //This is needed for the algorithm. If not there, it fails
+ const InDet::TRT_DriftCircleContainer* m_trtPrds = 0;
+ if( evtStore()->retrieve(m_trtPrds,m_driftcirclecontainer).isFailure() ) {
+ ATH_MSG_ERROR("Cannot retrieve TRT PrepDataContainer " << m_driftcirclecontainer);
+ return StatusCode::FAILURE;
+ }
+
+ //This is optional for the algorithm. If not there, just print a one-time warning
+ // On ESD
+ const PRD_MultiTruthCollection* m_prdmtColl = 0;
+ if (m_useTruthInfo) {
+ if ( evtStore()->contains<PRD_MultiTruthCollection>(m_multiTruth) ) {
+ if ( evtStore()->retrieve(m_prdmtColl, m_multiTruth).isFailure() ) {
+ ATH_MSG_ERROR("ERROR in retrieving PRD MultiTruth collection although available (" << m_multiTruth << ").");
+ return StatusCode::FAILURE;
+ }
+ } else {
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("PRD MultiTruth collection not available (" << m_multiTruth << "). Skipping this info although requested.");
+ m_prdmtColl = 0;
+ }
+ }
+ }
+
+ //This is optional for the algorithm. If not there, just print a one-time warning
+ // On RDO
+ const InDetSimDataCollection* m_sdoCollection = 0;
+ if (m_useTruthInfo && m_writeSDOs) {
+ if ( evtStore()->contains<InDetSimDataCollection>(m_SDOcontainer) ) {
+ if ( evtStore()->retrieve(m_sdoCollection, m_SDOcontainer).isFailure() ) {
+ ATH_MSG_ERROR("ERROR in retrieving SDO container despite being available. Collection = " << m_SDOcontainer);
+ return StatusCode::FAILURE;
+ }
+ } else {
+ if (m_firstEventWarnings) {
+ ATH_MSG_WARNING("SDO Collection not available (" << m_SDOcontainer << "). Skipping this info although requested.");
+ m_sdoCollection = 0;
+ }
+ }
+ }
+
+
+ // Create the xAOD container and its auxiliary store:
+ xAOD::TrackMeasurementValidationContainer* xaod = new xAOD::TrackMeasurementValidationContainer();
+ CHECK( evtStore()->record( xaod, m_driftcirclecontainer ) );
+ xAOD::TrackMeasurementValidationAuxContainer* aux = new xAOD::TrackMeasurementValidationAuxContainer();
+ CHECK( evtStore()->record( aux, m_driftcirclecontainer + "Aux." ) );
+ xaod->setStore( aux );
+
+ std::vector<unsigned int>* offsets = new std::vector<unsigned int>( m_TRTHelper->straw_layer_hash_max() , 0 );
+ CHECK( evtStore()->record( offsets, m_driftcirclecontainer + "Offsets" ) );
+
+ InDet::TRT_DriftCircleContainer::const_iterator it = m_trtPrds->begin();
+ InDet::TRT_DriftCircleContainer::const_iterator it_end = m_trtPrds->end();
+ unsigned int counter(0);
+ for( ; it!=it_end; ++it ) {
+
+ //Fill Offset container
+ if( m_TRTHelper->straw_layer_hash_max() <= (*it)->identifyHash() )
+ ATH_MSG_ERROR("My assumption about the maximum size of the hash was wrong");
+ (*offsets)[ (*it)->identifyHash() ] = counter;
+
+ // skip empty collections
+ if( (*it)->empty() ) continue;
+
+ // loop over collection and convert to xAOD
+ for( const auto& prd : **it ){
+
+ ++counter;
+
+ // create and add xAOD object
+ xAOD::TrackMeasurementValidation* xprd = new xAOD::TrackMeasurementValidation();
+ xaod->push_back(xprd);
+
+ Identifier surfaceID = prd->identify();
+
+ xprd->setIdentifier(surfaceID.get_compact());
+
+
+ // Save ID info:
+ xprd->auxdata<int>("bec") = m_TRTHelper->barrel_ec( surfaceID ) ;
+ xprd->auxdata<int>("layer") = m_TRTHelper->layer_or_wheel(surfaceID ) ;
+ xprd->auxdata<int>("phi_module") = m_TRTHelper->phi_module( surfaceID ) ;
+ xprd->auxdata<int>("strawlayer") = m_TRTHelper->straw_layer( surfaceID ) ;
+ xprd->auxdata<int>("strawnumber") = m_TRTHelper->straw( surfaceID ) ;
+ int chip=0;
+ int board=-1;
+ m_neighbourSvc->getChip(surfaceID,chip);
+ if(abs(m_TRTHelper->barrel_ec(surfaceID))<2)
+ board=m_neighbourSvc->chipToBoardBarrel(chip,m_TRTHelper->layer_or_wheel(surfaceID));
+ else if (chip<12)
+ board=0;
+ else {
+ chip=chip-20;
+ board=1;
+ }
+ xprd->auxdata<int>("board") = board ;
+ xprd->auxdata<int>("chip") = chip ;
+
+
+ //Set Local Position
+ const Amg::Vector2D& locpos = prd->localPosition();
+ float locY(0.);
+ float locX = locpos.x();
+ if ( !(std::isinf(locpos.y()) || std::isnan(locpos.y())) ){
+ if (locpos.y()>=1e-07)
+ locY=locpos.y();
+ } else {
+ locY = -9999.;
+ }
+
+ // Set local error matrix
+ xprd->setLocalPosition(locX,locY);
+
+ const Amg::MatrixX& localCov = prd->localCovariance();
+ if(localCov.size() == 1){
+ xprd->setLocalPositionError( localCov(0,0), 0., 0. );
+ } else if(localCov.size() == 2){
+ xprd->setLocalPositionError( localCov(0,0), localCov(1,1), localCov(0,1) );
+ } else {
+ xprd->setLocalPositionError(0.,0.,0.);
+ }
+
+
+ //Set Global position
+ const Amg::Vector3D* gpos = prd->detectorElement()->surface(surfaceID).localToGlobal(prd->localPosition());
+ if(gpos){
+ xprd->setGlobalPosition(gpos->x(),gpos->y(),gpos->z());
+ delete gpos;
+ }
+
+
+ //TRT hit bit word
+ unsigned int word = prd->getWord();
+
+ //TRTCond::RtRelation const *rtr = m_trtcaldbSvc->getRtRelation(surfaceID);
+ double tot = prd->timeOverThreshold();
+ bool isvalid=false;
+ xprd->auxdata<float>("drifttime") = prd->driftTime(isvalid) ;
+ xprd->auxdata<int>("status") = isvalid;
+ xprd->auxdata<float>("tot") = tot ;
+ xprd->auxdata<char>("isHT") = prd->highLevel() ;
+ xprd->auxdata<float>("T0") = m_trtcaldbSvc->getT0(surfaceID) ;
+
+
+ // Save time info:
+ xprd->auxdata<float>("leadingEdge") = prd->rawDriftTime();
+
+ //Drift Time corrections
+ xprd->auxdata<float>("driftTimeToTCorrection") = m_driftFunctionTool->driftTimeToTCorrection(tot,surfaceID);
+ if(prd->highLevel())
+ xprd->auxdata<float>("driftTimeHTCorrection") = m_driftFunctionTool->driftTimeHTCorrection(surfaceID);
+ else
+ xprd->auxdata<float>("driftTimeHTCorrection") = 0;
+
+
+ //HT Bit patterens
+ unsigned int theWord = word & 0x04020100; //HT bits mask
+ char highThreshold = 0;
+ //this is the agreed-upon binning for HT bit patterns...
+ if (theWord == 0x04000000) //pattern 1
+ highThreshold = 1;
+ else if (theWord == 0x00020000) //pattern 2
+ highThreshold = 2;
+ else if (theWord == 0x00000100) //pattern 3
+ highThreshold = 3;
+ else if (theWord == 0x04020000) //pattern 4
+ highThreshold = 4;
+ else if (theWord == 0x00020100) //pattern 5
+ highThreshold = 5;
+ else if (theWord == 0x04000100) //pattern 6
+ highThreshold = 6;
+ else if (theWord == 0x04020100) //pattern 7
+ highThreshold = 7;
+
+ xprd->auxdata<char>("highThreshold") = highThreshold;
+
+ //Full bit patternword from the TRT
+ xprd->auxdata<unsigned int>("bitPattern") = word;
+
+ char isArgonStraw = 0;
+ if (!m_TRTStrawSummarySvc.empty()) {
+ if (m_TRTStrawSummarySvc->getStatusHT(surfaceID) != TRTCond::StrawStatus::Good) {
+ isArgonStraw = 1;
+ }
+ }
+ xprd->auxdata<char>("isArgon") = isArgonStraw ;
+
+
+ // Use the MultiTruth Collection to get a list of all true particle contributing to the DC
+ if(m_prdmtColl){
+ std::vector<int> barcodes;
+ auto range = m_prdmtColl->equal_range(surfaceID);
+ for (auto i = range.first; i != range.second; ++i) {
+ barcodes.push_back( i->second.barcode() );
+ }
+ xprd->auxdata< std::vector<int> >("truth_barcode") = barcodes;
+ }
+
+ if( m_sdoCollection ){
+ // find hit
+ auto pos = m_sdoCollection->find(surfaceID);
+ int sdo_word = -1000000;
+ if( pos != m_sdoCollection->end() ) {
+ sdo_word = pos->second.word();
+ }
+ xprd->auxdata<int>("sdo_word") = sdo_word;
+ }
+
+ }
+ }
+ ATH_MSG_DEBUG( " recorded TRT_PrepData obejcts: size " << xaod->size() );
+
+
+
+ // Code to test that something was added to SG:
+ StatusCode sc = StatusCode::SUCCESS;
+ if (msgLvl(MSG::DEBUG)){
+ const xAOD::TrackMeasurementValidationContainer* xaod2;
+ sc = evtStore()->retrieve(xaod2, m_driftcirclecontainer);
+ if (sc.isFailure()) {
+ ATH_MSG_DEBUG("Failed to retrieve " << m_driftcirclecontainer);
+ }
+ else ATH_MSG_DEBUG("xAOD info retreived " << m_driftcirclecontainer << "\t" << xaod2->size());
+
+ const xAOD::TrackMeasurementValidationAuxContainer* aux2;
+ sc = evtStore()->retrieve(aux2, (m_driftcirclecontainer+ "Aux."));
+ if (sc.isFailure()) {
+ ATH_MSG_DEBUG("Failed to retrieve " << m_driftcirclecontainer<< "Aux.");
+ }
+ else ATH_MSG_DEBUG("xAOD info retrieved " << m_driftcirclecontainer << "Aux. \t");
+ }
+
+ // --- end of event. Disable one-time warnings
+ m_firstEventWarnings = false;
+
+ return sc;
+}
+
+/////////////////////////////////////////////////////////////////////
+//
+// Finalize method:
+//
+/////////////////////////////////////////////////////////////////////
+StatusCode TRT_PrepDataToxAOD::finalize()
+{
+ return StatusCode::SUCCESS;
+}
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.h b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.h
new file mode 100644
index 00000000000..8a3fd43579b
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/TRT_PrepDataToxAOD.h
@@ -0,0 +1,66 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// TRT_PrepDataToxAOD.h
+// Header file for class TRT_PrepDataToxAOD
+///////////////////////////////////////////////////////////////////
+
+#ifndef TRT_PREPDATATOXAOD_H
+#define TRT_PREPDATATOXAOD_H
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "TRT_ConditionsServices/ITRT_StrawNeighbourSvc.h"
+#include "TRT_ConditionsServices/ITRT_CalDbSvc.h"
+#include "TRT_ConditionsServices/ITRT_StrawStatusSummarySvc.h"
+#include "TRT_DriftFunctionTool/ITRT_DriftFunctionTool.h"
+
+
+#include <string>
+
+class TRT_ID;
+class ITRT_CalDbSvc ;
+class ITRT_DriftFunctionTool;
+class ITRT_StrawSummarySvc;
+
+
+class TRT_PrepDataToxAOD : public AthAlgorithm {
+
+public:
+ // Constructor with parameters:
+ TRT_PrepDataToxAOD(const std::string &name,ISvcLocator *pSvcLocator);
+
+ // Basic algorithm methods:
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+
+private:
+
+ // --- Steering and configuration flags
+ bool m_useTruthInfo;
+ bool m_writeSDOs;
+
+ // --- Configuration keys
+ std::string m_driftcirclecontainer;
+ std::string m_SDOcontainer;
+ std::string m_multiTruth;
+
+ // --- Services and Tools
+ ToolHandle< ITRT_DriftFunctionTool > m_driftFunctionTool ; //!< DriftFunctionTool
+ ServiceHandle<ITRT_CalDbSvc> m_trtcaldbSvc ;
+ ServiceHandle<ITRT_StrawNeighbourSvc> m_neighbourSvc ;
+ ServiceHandle<ITRT_StrawStatusSummarySvc> m_TRTStrawSummarySvc;
+ const TRT_ID *m_TRTHelper;
+
+ // ---- Internal members
+ bool m_firstEventWarnings;
+
+};
+
+
+#endif
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_entries.cxx b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_entries.cxx
new file mode 100755
index 00000000000..5313b9af720
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_entries.cxx
@@ -0,0 +1,18 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "../TRT_PrepDataToxAOD.h"
+#include "../SCT_PrepDataToxAOD.h"
+#include "../PixelPrepDataToxAOD.h"
+
+
+//using namespace InDet;
+
+DECLARE_ALGORITHM_FACTORY( TRT_PrepDataToxAOD )
+DECLARE_ALGORITHM_FACTORY( SCT_PrepDataToxAOD )
+DECLARE_ALGORITHM_FACTORY( PixelPrepDataToxAOD )
+
+DECLARE_FACTORY_ENTRIES( InDetPrepRawDataToxAOD )
+{
+ DECLARE_ALGORITHM( TRT_PrepDataToxAOD )
+ DECLARE_ALGORITHM( SCT_PrepDataToxAOD )
+ DECLARE_ALGORITHM( PixelPrepDataToxAOD )
+}
diff --git a/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_load.cxx b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_load.cxx
new file mode 100755
index 00000000000..249ab411bec
--- /dev/null
+++ b/InnerDetector/InDetEventCnv/InDetPrepRawDataToxAOD/src/components/InDetPrepRawDataToxAOD_load.cxx
@@ -0,0 +1,4 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( InDetPrepRawDataToxAOD )
+
--
GitLab