diff --git a/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt b/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..8f1e720b867398acbdbdcba483abcf61ed84fefc --- /dev/null +++ b/Control/CalypsoExample/BarcodeChecker/CMakeLists.txt @@ -0,0 +1,11 @@ +################################################################################ +# Package: BarcodeChecker +################################################################################ + +# Declare the package name: +atlas_subdir( BarcodeChecker ) + +# Install files from the package: +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) + +atlas_install_joboptions( share/*.py ) \ No newline at end of file diff --git a/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py new file mode 100644 index 0000000000000000000000000000000000000000..cc92bd8ac574dbaccfe522366666462b6219d988 --- /dev/null +++ b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py @@ -0,0 +1,38 @@ +# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + +import AthenaPython.PyAthena as PyAthena +from AthenaPython.PyAthena import StatusCode, McEventCollection, HepMC, CLHEP +import McParticleEvent.Pythonizations + +__author__ = "Dave Caser <dcasper@uci.edu>" + +class BarcodeCheckerAlg(PyAthena.Alg): + def __init__(self, name="BarcodeChecker", MCEventKey="TruthEvent"): + super(BarcodeCheckerAlg,self).__init__(name=name) + self.MCEventKey = MCEventKey + return + + def initialize(self): + self.maxLow = 0 + self.maxMid = 0 + self.maxHi = 0 + return StatusCode.Success + + + def execute(self): + evtCollection = self.evtStore[self.MCEventKey] + for mcEvt in evtCollection: + for mcParticle in mcEvt.particles: + barCode = mcParticle.barcode() + self.maxLow = max(self.maxLow, barCode%200000) + if barCode%1000000 > 200000: + self.maxMid = max(self.maxMid, barCode%1000000 - 200000) + self.maxHi = max(self.maxHi, barCode//1000000) + return StatusCode.Success + + def finalize(self): + print("Low part: ", self.maxLow, " out of 200000 (",100*self.maxLow/200000,"% of overflow)") + print("Mid part: ", self.maxMid, " out of ", 1000000 - 200000, " (",100*self.maxMid/(1000000-200000),"% of overflow") + print("Hi part: ", self.maxHi, " out of ", (1<<31)//1000000, " (", 100*self.maxHi/((1<<31)//1000000),"% of overflow") + + return StatusCode.Success \ No newline at end of file diff --git a/Control/CalypsoExample/BarcodeChecker/python/__init__.py b/Control/CalypsoExample/BarcodeChecker/python/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..d13ae164caa4e93bf2899cea4e67d0ec515784a2 --- /dev/null +++ b/Control/CalypsoExample/BarcodeChecker/python/__init__.py @@ -0,0 +1 @@ +# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration diff --git a/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py b/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py new file mode 100644 index 0000000000000000000000000000000000000000..1e739beb0c3865e924721b21fe5683b9f18fdaec --- /dev/null +++ b/Control/CalypsoExample/BarcodeChecker/share/BarcodeChecker_jobOptions.py @@ -0,0 +1,29 @@ +# +# Usage: athena.py -c'INPUT=["faser.1.gfaser.root","faser.2.gfaser.root"]' BarcodeChecker/BarcodeChecker_jobOptions.py >& BarcodeChecker.log +# +# INPUT is mandatory (INPUT can be a list, as shown above) +# + +if not 'INPUT' in dir(): + print("Missing INPUT parameter") + exit() + +if not isinstance(INPUT, (list,tuple)): + INPUT = [INPUT] + pass + +from AthenaCommon.GlobalFlags import globalflags + +globalflags.InputFormat.set_Value_and_Lock('pool') + +import AthenaPoolCnvSvc.ReadAthenaPool + +svcMgr.EventSelector.InputCollections = INPUT + +from BarcodeChecker.BarcodeCheckerAlg import BarcodeCheckerAlg + +from AthenaCommon.AlgSequence import AlgSequence +job = AlgSequence() +job += BarcodeCheckerAlg() + +theApp.EvtMax = -1 diff --git a/Control/CalypsoExample/Digitization/CMakeLists.txt b/Control/CalypsoExample/Digitization/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..f8caebb598f6cea848edafea449a0bd19e657c6d --- /dev/null +++ b/Control/CalypsoExample/Digitization/CMakeLists.txt @@ -0,0 +1,26 @@ +################################################################################ +# Package: Digitization +################################################################################ + +# Declare the package name: +atlas_subdir( Digitization ) + +# Component(s) in the package: +#atlas_add_component( GeoModelTest +# src/GeoModelTestAlg.cxx +# src/components/GeoModelTest_entries.cxx +# INCLUDE_DIRS ${GEOMODEL_INCLUDE_DIRS} +# LINK_LIBRARIES ${GEOMODEL_LIBRARIES} AthenaBaseComps GeoModelFaserUtilities ScintReadoutGeometry TrackerReadoutGeometry MagFieldInterfaces MagFieldElements MagFieldConditions ) + +# Install files from the package: +#atlas_install_joboptions( share/*.py ) +#atlas_install_python_modules( python/*.py ) +atlas_install_scripts( scripts/*.sh scripts/*.py ) + +#atlas_add_test( ProdRecoTI12 +# SCRIPT scripts/faser_reco.py ${CMAKE_CURRENT_SOURCE_DIR}/../rawdata/Faser-Physics-001920-filtered.raw TI12Data +# PROPERTIES TIMEOUT 300 ) +#atlas_add_test( ProdRecoTestBeam +# SCRIPT scripts/faser_reco.py ${CMAKE_CURRENT_SOURCE_DIR}/../RAWDATA/Faser-Physics-003613-filtered.raw TestBeamData +# PROPERTIES TIMEOUT 300 ) + diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi.py b/Control/CalypsoExample/Digitization/scripts/faser_digi.py new file mode 100755 index 0000000000000000000000000000000000000000..b9d0708be4e7d8b5c443397c6df0826b3f95d2df --- /dev/null +++ b/Control/CalypsoExample/Digitization/scripts/faser_digi.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python +# +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Run with: +# ./faser_digi.py filepath runtype +# +# filepath - fully qualified path, including url if needed, to the input HITS file +# example: "root://eospublic.cern.ch//eos/experiment/faser/sim/GeniePilot/HITS/1/faser.150fbInv.1.001.HITS.pool.root" +# +# runtype - MANDATORY flag to specify the data type (TI12OldMC or TI12MC or TestBeamMC). +# Not extracted (yet) from file path for MC data +# +import sys +import argparse + +parser = argparse.ArgumentParser(description="Run FASER reconstruction") + +parser.add_argument("file_path", + help="Fully qualified path of the raw input file") +parser.add_argument("run_type", nargs="?", default="", + help="Specify run type (if it can't be parsed from path)") +parser.add_argument("-r", "--reco", default="", + help="Specify reco tag (to append to output filename)") +parser.add_argument("-n", "--nevents", type=int, default=-1, + help="Specify number of events to process (default: all)") +parser.add_argument("-v", "--verbose", action='store_true', + help="Turn on DEBUG output") + +args = parser.parse_args() + +from pathlib import Path + +filepath=Path(args.file_path) + +# runtype has been provided +if len(args.run_type) > 0: + runtype=args.run_type + +# Extract runtype from path +# Should be directory above run +# i.e.: TestBeamData/Run-004150/Faser-Physics-004150-00000.raw" +else: + if True or len(filepath.parts) < 3: + print("Can't determine run type from path - specify on command line ") + sys.exit(-1) + +# runtype = filepath.parts[-3] + +print(f"Starting digitization of {filepath.name} with type {runtype}") +if args.nevents > 0: + print(f"Reconstructing {args.nevents} events by command-line option") + +# Start digitization + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory +from AthenaCommon.Constants import VERBOSE, INFO + +from AthenaCommon.Configurable import Configurable +from CalypsoConfiguration.AllConfigFlags import ConfigFlags + +Configurable.configurableRun3Behavior = True + +# Flags for this job +ConfigFlags.Input.isMC = False # Needed to bypass autoconfig +ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + +ConfigFlags.Input.ProjectName = "mc20" +ConfigFlags.GeoModel.Align.Dynamic = False +ConfigFlags.Beam.NumberOfCollisions = 0. +ConfigFlags.Digitization.TruthOutput = True + +# TI12 old geometry +if runtype == "TI12OldMC": + ConfigFlags.GeoModel.FaserVersion = "FASER-01" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" + +# Testbeam setup +elif runtype == "TestBeamMC" : + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00" + +# New TI12 geometry (ugh) +elif runtype == "TI12MC": + ConfigFlags.GeoModel.FaserVersion = "FASERNU-02" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" + +else: + print("Invalid run type found:", runtype) + print("Specify correct type or update list") + sys.exit(-1) + + +# Must use original input string here, as pathlib mangles double // in path names +ConfigFlags.Input.Files = [ args.file_path ] + +filestem = filepath.stem +if len(args.reco) > 0: + filestem += f"-{args.reco}" + +# ConfigFlags.addFlag("Output.xAODFileName", f"{filestem}-xAOD.root") +ConfigFlags.Output.RDOFileName = f"{filestem}-RDO.root" + +# +# Play around with this? +# ConfigFlags.Concurrency.NumThreads = 2 +# ConfigFlags.Concurrency.NumConcurrentEvents = 2 +ConfigFlags.lock() + +# +# Configure components +from CalypsoConfiguration.MainServicesConfig import MainServicesCfg +from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg +from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + +acc = MainServicesCfg(ConfigFlags) +acc.merge(PoolReadCfg(ConfigFlags)) +acc.merge(PoolWriteCfg(ConfigFlags)) + +# +# Needed, or move to MainServicesCfg? +from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg +acc.merge(FaserGeometryCfg(ConfigFlags)) + +# Set up algorithms +from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_DigitizationCfg +acc.merge(FaserSCT_DigitizationCfg(ConfigFlags)) + +# from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg +# acc.merge(WaveformReconstructionCfg(ConfigFlags)) + +# # Not ready for primetime +# # from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg +# # acc.merge(CalorimeterReconstructionCfg(ConfigFlags)) + +# # Tracker clusters +# from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg +# acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags)) + +# # SpacePoints +# from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg +# acc.merge(TrackerSpacePointFinderCfg(ConfigFlags)) + +# # Try Dave's fitter +# from TrackerClusterFit.TrackerClusterFitConfig import ClusterFitAlgCfg +# acc.merge(ClusterFitAlgCfg(ConfigFlags)) + +# +# Configure output +# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg +# itemList = [ "xAOD::EventInfo#*" +# , "xAOD::EventAuxInfo#*" +# , "xAOD::FaserTriggerData#*" +# , "xAOD::FaserTriggerDataAux#*" +# , "FaserSCT_RDO_Container#*" +# # , "Tracker::FaserSCT_ClusterContainer#*" +# # , "FaserSCT_SpacePointContainer#*" +# # , "FaserSCT_SpacePointOverlapCollection#*" +# # , "TrackCollection#*" +# ] +# acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList)) + +# Waveform reconstruction output +# from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg +# acc.merge(WaveformReconstructionOutputCfg(ConfigFlags)) + +# Calorimeter reconstruction output +# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg +# acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags)) + +# Check what we have +# print( "Writing out xAOD objects:" ) +# print( acc.getEventAlgo("OutputStreamxAOD").ItemList ) + +# Hack to avoid problem with our use of MC databases when isMC = False +replicaSvc = acc.getService("DBReplicaSvc") +replicaSvc.COOLSQLiteVetoPattern = "" +replicaSvc.UseCOOLSQLite = True +replicaSvc.UseCOOLFrontier = False +replicaSvc.UseGeomSQLite = True + +# Configure verbosity +# ConfigFlags.dump() +if args.verbose: + acc.foreach_component("*").OutputLevel = VERBOSE + + #acc.getService("FaserByteStreamInputSvc").DumpFlag = True + #acc.getService("FaserEventSelector").OutputLevel = VERBOSE + #acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE + #acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE + #acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE + +else: + acc.foreach_component("*").OutputLevel = INFO + +acc.foreach_component("*ClassID*").OutputLevel = INFO + +acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M" + +# Execute and finish +sys.exit(int(acc.run(maxEvents=args.nevents).isFailure())) diff --git a/Generators/GenieReader/python/GenieReaderAlg.py b/Generators/GenieReader/python/GenieReaderAlg.py index c930ae66ae2aa8ef4eddf7f31873acc598d31a5c..c97ea494b089a25ff88e4c2a1f2abdd4bf472c6b 100644 --- a/Generators/GenieReader/python/GenieReaderAlg.py +++ b/Generators/GenieReader/python/GenieReaderAlg.py @@ -2,9 +2,8 @@ from AthenaCommon.AppMgr import ServiceMgr as svcMgr from GeneratorModules.EvgenAlg import EvgenAlg -from AthenaPython.PyAthena import StatusCode +from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType from AthenaCommon.SystemOfUnits import GeV, m - import ROOT __author__ = "Dave Caser <dcasper@uci.edu>" @@ -22,6 +21,22 @@ class GenieReader(EvgenAlg): from AthenaPython.PyAthena import HepMC as HepMC evt.weights().push_back(1.0) + treeEventInfo = self.evtStore["TTreeEventInfo"] + treeEventId = treeEventInfo.event_ID() + runNumber = treeEventId.run_number() + eventNumber = treeEventId.event_number() + # print("found run/event number ", runNumber, "/", eventNumber) + + mcEventType = EventType() + mcEventType.add_type(EventType.IS_SIMULATION) + + mcEventId = EventID(run_number = runNumber, event_number = eventNumber) + mcEventInfo = EventInfo(id = mcEventId, type = mcEventType) + self.evtStore.record(mcEventInfo, "McEventInfo", True, False) + ROOT.SetOwnership(mcEventType, False) + ROOT.SetOwnership(mcEventId, False) + ROOT.SetOwnership(mcEventInfo, False) + pos = HepMC.FourVector(self.evtStore["vx"]*m, self.evtStore["vy"]*m, self.evtStore["vz"]*m, 0) gv = HepMC.GenVertex(pos) ROOT.SetOwnership(gv, False) @@ -34,7 +49,7 @@ class GenieReader(EvgenAlg): py = list(self.evtStore["py"]) pz = list(self.evtStore["pz"]) E = list(self.evtStore["E"]) - M = list(self.evtStore["M"]) + M = list(self.evtStore["M"]) for i in range(nParticles): gp = HepMC.GenParticle() diff --git a/Generators/GenieReader/share/GenieReader_jobOptions.py b/Generators/GenieReader/share/GenieReader_jobOptions.py index 91ad7c1318b772fd8cb015eda03ea274b0a6815a..649e705b66eb05f7b5ac2b331890a914794bec60 100644 --- a/Generators/GenieReader/share/GenieReader_jobOptions.py +++ b/Generators/GenieReader/share/GenieReader_jobOptions.py @@ -30,9 +30,10 @@ job += GenieReader() from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream ostream = AthenaPoolOutputStream( "StreamEVGEN" , OUTPUT, noTag=True ) -ostream.ItemList += [ "xAOD::EventInfo#*", "xAOD::EventAuxInfo#*", - "EventInfo#*", +ostream.ItemList.remove("EventInfo#*") +ostream.ItemList += [ "EventInfo#McEventInfo", "McEventCollection#*" ] +print(ostream.ItemList) if not 'EVTMAX' in dir(): EVTMAX = -1 diff --git a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx index 88592d298b2eb974007a73c15ad2413beb193648..f9c937aa22c8a69153e7c32174d4c1009e9b5845 100644 --- a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx +++ b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.cxx @@ -3,6 +3,7 @@ */ #include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h" #include "ScintSimEventTPCnv/ScintHits/ScintHit_p1.h" #include "ScintHitCollectionCnv.h" @@ -18,11 +19,19 @@ ScintHitCollection_PERS* ScintHitCollectionCnv::createPersistent(ScintHitCollect ScintHitCollection* ScintHitCollectionCnv::createTransient() { MsgStream mlog(msgSvc(), "ScintHitCollectionConverter" ); ScintHitCollectionCnv_p1 converter_p1; + ScintHitCollectionCnv_p1a converter_p1a; static const pool::Guid p1_guid("B2573A16-4B46-4E1E-98E3-F93421680779"); + static const pool::Guid p1a_guid("0DFC461F-9FF5-4447-B4F0-7CC7157191D1"); ScintHitCollection *trans_cont(0); - if( this->compareClassGuid(p1_guid)) { + if( this->compareClassGuid(p1a_guid)) + { + std::unique_ptr< ScintHitCollection_p1a > col_vect( this->poolReadObject< ScintHitCollection_p1a >() ); + trans_cont = converter_p1a.createTransient( col_vect.get(), mlog ); + } + else if( this->compareClassGuid(p1_guid)) + { std::unique_ptr< ScintHitCollection_p1 > col_vect( this->poolReadObject< ScintHitCollection_p1 >() ); trans_cont = converter_p1.createTransient( col_vect.get(), mlog ); } else { diff --git a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h index ca9c329ddad032c35766ad6a271fd7637734b1c0..4cdaf2046a595f763c424baae5f6942ab91a33fd 100644 --- a/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h +++ b/Scintillator/ScintEventCnv/ScintSimEventAthenaPool/src/ScintHitCollectionCnv.h @@ -8,12 +8,14 @@ #include "ScintSimEvent/ScintHitCollection.h" #include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1.h" #include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h" #include "AthenaPoolCnvSvc/T_AthenaPoolCustomCnv.h" // Gaudi #include "GaudiKernel/MsgStream.h" // typedef to the latest persistent version -typedef ScintHitCollection_p1 ScintHitCollection_PERS; -typedef ScintHitCollectionCnv_p1 ScintHitCollectionCnv_PERS; +typedef ScintHitCollection_p1a ScintHitCollection_PERS; +typedef ScintHitCollectionCnv_p1a ScintHitCollectionCnv_PERS; class ScintHitCollectionCnv : public T_AthenaPoolCustomCnv<ScintHitCollection, ScintHitCollection_PERS > { friend class CnvFactory<ScintHitCollectionCnv>; diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h new file mode 100644 index 0000000000000000000000000000000000000000..4f8de2b4c704240402bdf5ea23a36cf3a12acefe --- /dev/null +++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h @@ -0,0 +1,41 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SCINTHITCOLLECTIONCNV_P1A_H +#define SCINTHITCOLLECTIONCNV_P1A_H + +// ScintHitCollectionCnv_p1, T/P separation of Scint Hits +// author D.Costanzo <davide.costanzo@cern.ch> +// author O.Arnaez <olivier.arnaez@cern.ch> + +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" +#include "ScintSimEvent/ScintHitCollection.h" +#include "ScintHitCollection_p1a.h" + + +class ScintHitCollectionCnv_p1a : public T_AthenaPoolTPCnvBase<ScintHitCollection, ScintHitCollection_p1a> +{ + public: + + ScintHitCollectionCnv_p1a() {}; + + virtual ScintHitCollection* createTransient(const ScintHitCollection_p1a* persObj, MsgStream &log); + + virtual void persToTrans(const ScintHitCollection_p1a* persCont, + ScintHitCollection* transCont, + MsgStream &log) ; + virtual void transToPers(const ScintHitCollection* transCont, + ScintHitCollection_p1a* persCont, + MsgStream &log) ; + + private: + + static const double m_persEneUnit; + static const double m_persLenUnit; + static const double m_persAngUnit; + static const double m_2bHalfMaximum; + static const int m_2bMaximum; +}; + +#endif diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h new file mode 100644 index 0000000000000000000000000000000000000000..8e714d6e571c746ae6376536ed62c5681476232e --- /dev/null +++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h @@ -0,0 +1,57 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SCINTHITCOLLECTION_P1A_H +#define SCINTHITCOLLECTION_P1A_H + +/* + +Authors: Davide Costanzo Rob Duxfield + +*/ + +#include <vector> +#include <string> + +class ScintHitCollection_p1a +{ + public: +/// Default constructor + ScintHitCollection_p1a (); + //private: + + std::vector<float> m_hit1_meanTime; // 1 element per string + std::vector<float> m_hit1_x0; // + std::vector<float> m_hit1_y0; // + std::vector<float> m_hit1_z0; // + std::vector<float> m_hit1_theta; // + std::vector<float> m_hit1_phi; // + std::vector<unsigned long> m_nHits; // + + std::vector<unsigned short> m_hitEne_2b; // 1 element per hit + std::vector<unsigned short> m_hitLength_2b; // + + std::vector<unsigned short> m_dTheta; // 1 element per hit except for first hit in string + std::vector<unsigned short> m_dPhi; // + + std::vector<float> m_hitEne_4b; // 1 element per hit with m_hitEne_2b[i] == 2**16 + + std::vector<float> m_hitLength_4b; // 1 element per hit with m_hitLength_2b[i] == 2**16 + + std::vector<unsigned long> m_barcode; + std::vector<unsigned long> m_mcEvtIndex; + std::vector<char> m_evtColl; + std::vector<unsigned long> m_nBC; + + std::vector<unsigned long> m_id; + std::vector<unsigned long> m_nId; +}; + + +// inlines + +inline +ScintHitCollection_p1a::ScintHitCollection_p1a () {} + +#endif diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h index f0042a97ec671d78b8fb783af5f516116570952a..1606077dd0efcd14371c1c7445d7a553b0b0579d 100644 --- a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h +++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/ScintSimEventTPCnvDict.h @@ -15,6 +15,8 @@ #include "ScintSimEventTPCnv/ScintHits/ScintHitCnv_p1.h" #include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1.h" #include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h" #include "ScintSimEventTPCnv/ScintHits/ScintHit_p1.h" #endif // SCINTEVENTTPCNV_SCINTSIMEVENTTPCNVDICT_H diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml index 6b0d915eac9a1e2fedb5f8254bbbf419dec16dac..345d08165aeb11723e94f19e7195f42fcb5a88af 100644 --- a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml +++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/ScintSimEventTPCnv/selection.xml @@ -4,4 +4,5 @@ <class name="ScintHit_p1" /> <class name="std::vector<ScintHit_p1>" /> <class name="ScintHitCollection_p1" id="B2573A16-4B46-4E1E-98E3-F93421680779" /> + <class name="ScintHitCollection_p1a" id="0DFC461F-9FF5-4447-B4F0-7CC7157191D1" /> </lcgdict> diff --git a/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cd162d1ad3c8064966a59c3517b19b233bc29c80 --- /dev/null +++ b/Scintillator/ScintEventCnv/ScintSimEventTPCnv/src/ScintHits/ScintHitCollectionCnv_p1a.cxx @@ -0,0 +1,334 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ScintSimEvent/ScintHit.h" +#include "ScintSimEvent/ScintHitCollection.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollection_p1a.h" +#include "ScintSimEventTPCnv/ScintHits/ScintHitCollectionCnv_p1a.h" +#include "GeneratorObjects/HepMcParticleLink.h" + +#include <cmath> + +//CLHEP +#include "CLHEP/Geometry/Point3D.h" +// Gaudi +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/ThreadLocalContext.h" + +// Athena +#include "AthenaKernel/ExtendedEventContext.h" +#include "StoreGate/StoreGateSvc.h" + +// * * * stolen from eflowRec * * * // +inline double phicorr(double a) +{ + if (a <= -M_PI) + { + return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI))); + } + else if (a > M_PI) + { + return a-(2*M_PI*floor((a+M_PI)/(2*M_PI))); + } + else + { + return a; + } +} + +// * * * stolen from eflowRec * * * // +inline double cycle(double a, double b) +{ + double del = b-a; + if (del > M_PI) + { + return a+2.0*M_PI; + } + else if (del < -M_PI) + { + return a-2.0*M_PI; + } + else + { + return a; + } +} + + +const double ScintHitCollectionCnv_p1a::m_persEneUnit = 1.0e-5; +const double ScintHitCollectionCnv_p1a::m_persLenUnit = 1.0e-5; +const double ScintHitCollectionCnv_p1a::m_persAngUnit = 1.0e-5; +const double ScintHitCollectionCnv_p1a::m_2bHalfMaximum = pow(2.0, 15.0); +const int ScintHitCollectionCnv_p1a::m_2bMaximum = (unsigned short)(-1); + + +void ScintHitCollectionCnv_p1a::transToPers(const ScintHitCollection* transCont, ScintHitCollection_p1a* persCont, MsgStream &log) +{ + // Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and + // persistifies the end point of each hit plus the start point of the first hit in each string. + // + // Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of + // first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and + // used to specify the hit direction. All hit lengths are stored as 2 byte numbers. + // + // Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean + // time of the first hit per string. + // + // See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info. + + static const double dRcut = 1.0e-7; + static const double dTcut = 1.0; + + const EventContext& ctx = Gaudi::Hive::currentContext(); + const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy(); + const HepMcParticleLink * lastLink=nullptr; + int lastId = -1; + double stringFirstTheta = 0.0; + double stringFirstPhi = 0.0; + double lastT = 0.0; + double persSumE = 0.0; + double transSumE = 0.0; + unsigned int idx = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + unsigned int endHit = 0; + HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0); + HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0); + + for (ScintHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) { + + ScintHitCollection::const_iterator siHit = it; + + + if ( !lastLink || (siHit->particleLink() != *lastLink) ) { + + // store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode + + lastLink = &(siHit->particleLink()); + persCont->m_barcode.push_back(lastLink->barcode()); + unsigned short index{0}; + const HepMcParticleLink::index_type position = + HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(), + lastLink->getEventCollection(), + proxy).at(0); + if (position!=0) { + index = lastLink->eventIndex(); + if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) { + log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg; + } + } + persCont->m_mcEvtIndex.push_back(index); + persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar()); + + if (idx > 0) { + persCont->m_nBC.push_back(idx - endBC); + endBC = idx; + } + } + + if ( (int)siHit->identify() != lastId ) { + + // store id once for set of consecutive hits with same barcode + + lastId = siHit->identify(); + persCont->m_id.push_back(lastId); + + if (idx > 0) { + persCont->m_nId.push_back(idx - endId); + endId = idx; + } + } + + HepGeom::Point3D<double> st = siHit->localStartPosition(); + HepGeom::Point3D<double> en = siHit->localEndPosition(); + + const double dx = st.x() - lastTransEnd.x(); + const double dy = st.y() - lastTransEnd.y(); + const double dz = st.z() - lastTransEnd.z(); + const double t = siHit->meanTime(); + + const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one + const double dTLast = fabs(t - lastT); + + CLHEP::Hep3Vector direction(0.0, 0.0, 0.0); + double theta = 0.0; + double phi = 0.0; + bool startNewString = false; + + if (dRLast < dRcut && dTLast < dTcut) { + + // hit is part of existing string + + direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() ); + + theta = direction.theta(); + phi = phicorr( direction.phi() ); + + const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 ); + const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 ); + + if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) { + persCont->m_dTheta.push_back(dTheta_2b); + persCont->m_dPhi.push_back(dPhi_2b); + theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit; + phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit; + phi = phicorr(phi); + } + else { + startNewString = true; + } + } + + if (startNewString || dRLast >= dRcut || dTLast >= dTcut) { + + // begin new hit string + + direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() ); + + theta = direction.theta(); + phi = phicorr( direction.phi() ); + + persCont->m_hit1_meanTime.push_back(t); + persCont->m_hit1_x0.push_back(st.x()); + persCont->m_hit1_y0.push_back(st.y()); + persCont->m_hit1_z0.push_back(st.z()); + persCont->m_hit1_theta.push_back(theta); + persCont->m_hit1_phi.push_back(phi); + + lastPersEnd = HepGeom::Point3D<double>(st.x(), st.y(), st.z()); + + stringFirstTheta = theta; + stringFirstPhi = phi; + + if (idx > 0) { + persCont->m_nHits.push_back(idx - endHit); + endHit = idx; + } + } + + lastTransEnd = HepGeom::Point3D<double>(en.x(), en.y(), en.z()); + transSumE += siHit->energyLoss(); + + const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over + // whole hit string to chosen precision + + const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to + // the chosen precision, NOT the length of the + // hit (small difference in practice). + double eneLoss = 0.0; + + if (eneLoss_2b >= m_2bMaximum) { + eneLoss = siHit->energyLoss(); + persCont->m_hitEne_2b.push_back(m_2bMaximum); + persCont->m_hitEne_4b.push_back(eneLoss); + } + else { + eneLoss = eneLoss_2b * m_persEneUnit; + persCont->m_hitEne_2b.push_back(eneLoss_2b); + } + + double length = 0.0; + + if (hitLength_2b >= m_2bMaximum) { + length = direction.mag(); + persCont->m_hitLength_2b.push_back(m_2bMaximum); + persCont->m_hitLength_4b.push_back(direction.mag()); + } + else { + length = hitLength_2b * m_persLenUnit; + persCont->m_hitLength_2b.push_back(hitLength_2b); + } + + CLHEP::Hep3Vector persDir(length, 0.0, 0.0); + persDir.setTheta(theta); + persDir.setPhi(phi); + + lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir; + persSumE += eneLoss; + lastT = t; + + ++idx; + } + + persCont->m_nBC.push_back(idx - endBC); + persCont->m_nId.push_back(idx - endId); + persCont->m_nHits.push_back(idx - endHit); +} + + +ScintHitCollection* ScintHitCollectionCnv_p1a::createTransient(const ScintHitCollection_p1a* persObj, MsgStream &log) { + std::unique_ptr<ScintHitCollection> trans(std::make_unique<ScintHitCollection>("DefaultCollectionName",persObj->m_nHits.size())); + persToTrans(persObj, trans.get(), log); + return(trans.release()); +} + + +void ScintHitCollectionCnv_p1a::persToTrans(const ScintHitCollection_p1a* persCont, ScintHitCollection* transCont, MsgStream &/*log*/) +{ + const EventContext& ctx = Gaudi::Hive::currentContext(); + + unsigned int hitCount = 0; + unsigned int angleCount = 0; + unsigned int idxBC = 0; + unsigned int idxId = 0; + unsigned int idxEne4b = 0; + unsigned int idxLen4b = 0; + unsigned int endHit = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + + for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) { + + if (persCont->m_nHits[i]) { + + const unsigned int start = endHit; + endHit += persCont->m_nHits[i]; + + const double t0 = persCont->m_hit1_meanTime[i]; + const double theta0 = persCont->m_hit1_theta[i]; + const double phi0 = persCont->m_hit1_phi[i]; + HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]); + + for (unsigned int j = start; j < endHit; j++) { + + if (j >= endBC + persCont->m_nBC[idxBC]) + endBC += persCont->m_nBC[idxBC++]; + + if (j >= endId + persCont->m_nId[idxId]) + endId += persCont->m_nId[idxId++]; + + const double eneLoss_2b = persCont->m_hitEne_2b[hitCount]; + const double hitLength_2b = persCont->m_hitLength_2b[hitCount]; + + const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++]; + const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++]; + + const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0; + const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0; + + const double meanTime = t0; + const double theta = theta0 + dTheta; + const double phi = phicorr(phi0 + dPhi); + + CLHEP::Hep3Vector r(length, 0.0, 0.0); + r.setTheta(theta); + r.setPhi(phi); + + HepGeom::Point3D<double> endThis( endLast + r ); + + HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX; + if (persCont->m_mcEvtIndex[idxBC] == 0) { + flag = HepMcParticleLink::IS_POSITION; + } + HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx ); + transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]); + + endLast = endThis; + + ++hitCount; + if (j > start) ++angleCount; + } + } + } +} diff --git a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx index 51514fb5b599630e9c32d04c406439c5ea891991..2149bdb33ed2ef55ea57a0c1bf55b025dd69dc16 100644 --- a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx +++ b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.cxx @@ -3,7 +3,6 @@ */ #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h" -//#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHit_p2.h" #include "FaserSiHitCollectionCnv.h" @@ -18,14 +17,23 @@ FaserSiHitCollection_PERS* FaserSiHitCollectionCnv::createPersistent(FaserSiHitC FaserSiHitCollection* FaserSiHitCollectionCnv::createTransient() { MsgStream mlog(msgSvc(), "FaserSiHitCollectionConverter" ); FaserSiHitCollectionCnv_p3 converter_p3; + FaserSiHitCollectionCnv_p3a converter_p3a; static const pool::Guid p3_guid("FF9508DE-3E25-425D-9556-16D319DCE0E1"); + static const pool::Guid p3a_guid("72FD9F51-AB1B-4DF7-B430-6CCAE0A994DB"); FaserSiHitCollection *trans_cont(0); - if( this->compareClassGuid(p3_guid)) { + if( this->compareClassGuid(p3a_guid)) + { + std::unique_ptr< FaserSiHitCollection_p3a > col_vect( this->poolReadObject< FaserSiHitCollection_p3a >() ); + trans_cont = converter_p3a.createTransient( col_vect.get(), mlog ); + } + else if( this->compareClassGuid(p3_guid)) + { std::unique_ptr< FaserSiHitCollection_p3 > col_vect( this->poolReadObject< FaserSiHitCollection_p3 >() ); trans_cont = converter_p3.createTransient( col_vect.get(), mlog ); - } else { + } + else { throw std::runtime_error("Unsupported persistent version of Data container"); } return trans_cont; diff --git a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h index 3a81e19fa0671a028948f7217aa256c06e94c9a1..47892fe8a5cd1738428aa90858c95e83d46fe232 100644 --- a/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h +++ b/Tracker/TrackerEventCnv/TrackerSimEventAthenaPool/src/FaserSiHitCollectionCnv.h @@ -8,12 +8,14 @@ #include "TrackerSimEvent/FaserSiHitCollection.h" #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3.h" #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h" #include "AthenaPoolCnvSvc/T_AthenaPoolCustomCnv.h" // Gaudi #include "GaudiKernel/MsgStream.h" // typedef to the latest persistent version -typedef FaserSiHitCollection_p3 FaserSiHitCollection_PERS; -typedef FaserSiHitCollectionCnv_p3 FaserSiHitCollectionCnv_PERS; +typedef FaserSiHitCollection_p3a FaserSiHitCollection_PERS; +typedef FaserSiHitCollectionCnv_p3a FaserSiHitCollectionCnv_PERS; class FaserSiHitCollectionCnv : public T_AthenaPoolCustomCnv<FaserSiHitCollection, FaserSiHitCollection_PERS > { friend class CnvFactory<FaserSiHitCollectionCnv>; diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h new file mode 100644 index 0000000000000000000000000000000000000000..8f453cf57ab52efb41fbe739d6543d75ba77c72a --- /dev/null +++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h @@ -0,0 +1,41 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef FASERSIHITCOLLECTIONCNV_P3A_H +#define FASERSIHITCOLLECTIONCNV_P3A_H + +// FaserSiHitCollectionCnv_p3, T/P separation of FaserSi Hits +// author D.Costanzo <davide.costanzo@cern.ch> +// author O.Arnaez <olivier.arnaez@cern.ch> + +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" +#include "TrackerSimEvent/FaserSiHitCollection.h" +#include "FaserSiHitCollection_p3a.h" + + +class FaserSiHitCollectionCnv_p3a : public T_AthenaPoolTPCnvBase<FaserSiHitCollection, FaserSiHitCollection_p3a> +{ + public: + + FaserSiHitCollectionCnv_p3a() {}; + + virtual FaserSiHitCollection* createTransient(const FaserSiHitCollection_p3a* persObj, MsgStream &log); + + virtual void persToTrans(const FaserSiHitCollection_p3a* persCont, + FaserSiHitCollection* transCont, + MsgStream &log) ; + virtual void transToPers(const FaserSiHitCollection* transCont, + FaserSiHitCollection_p3a* persCont, + MsgStream &log) ; + + private: + + static const double m_persEneUnit; + static const double m_persLenUnit; + static const double m_persAngUnit; + static const double m_2bHalfMaximum; + static const int m_2bMaximum; +}; + +#endif diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h new file mode 100644 index 0000000000000000000000000000000000000000..137ebd68b84c0ea6e6643132023d90c49b23df0f --- /dev/null +++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h @@ -0,0 +1,57 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef FASERSIHITCOLLECTION_P3A_H +#define FASERSIHITCOLLECTION_P3A_H + +/* + +Authors: Davide Costanzo Rob Duxfield + +*/ + +#include <vector> +#include <string> + +class FaserSiHitCollection_p3a +{ + public: +/// Default constructor + FaserSiHitCollection_p3a (); + //private: + + std::vector<float> m_hit1_meanTime; // 1 element per string + std::vector<float> m_hit1_x0; // + std::vector<float> m_hit1_y0; // + std::vector<float> m_hit1_z0; // + std::vector<float> m_hit1_theta; // + std::vector<float> m_hit1_phi; // + std::vector<unsigned long> m_nHits; // + + std::vector<unsigned short> m_hitEne_2b; // 1 element per hit + std::vector<unsigned short> m_hitLength_2b; // + + std::vector<unsigned short> m_dTheta; // 1 element per hit except for first hit in string + std::vector<unsigned short> m_dPhi; // + + std::vector<float> m_hitEne_4b; // 1 element per hit with m_hitEne_2b[i] == 2**16 + + std::vector<float> m_hitLength_4b; // 1 element per hit with m_hitLength_2b[i] == 2**16 + + std::vector<unsigned long> m_barcode; + std::vector<unsigned long> m_mcEvtIndex; + std::vector<char> m_evtColl; + std::vector<unsigned long> m_nBC; + + std::vector<unsigned long> m_id; + std::vector<unsigned long> m_nId; +}; + + +// inlines + +inline +FaserSiHitCollection_p3a::FaserSiHitCollection_p3a () {} + +#endif diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h index a730fca90074a9da71420e41651384173aa298cf..a898e3095bee7e4b25d2090b0e7365cbdcc4d00e 100644 --- a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h +++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnvDict.h @@ -15,6 +15,8 @@ #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCnv_p2.h" #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3.h" #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h" #include "TrackerSimEventTPCnv/TrackerHits/FaserSiHit_p2.h" #endif // TRACKERSIMEVENTTPCNV_TRACKERSIMEVENTTPCNVDICT_H diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml index 7684404421c37267df94656b57ff9067dcce9503..6b34f88fcb38eb134efb2a283ac6635860a76928 100644 --- a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml +++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/TrackerSimEventTPCnv/selection.xml @@ -3,4 +3,5 @@ <class name="FaserSiHit_p2" /> <class name="std::vector<FaserSiHit_p2>" /> <class name="FaserSiHitCollection_p3" id="FF9508DE-3E25-425D-9556-16D319DCE0E1" /> + <class name="FaserSiHitCollection_p3a" id="72FD9F51-AB1B-4DF7-B430-6CCAE0A994DB" /> </lcgdict> diff --git a/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0aebca47dc6ae6eb0d6662d4a2b4abe1a52dc1ae --- /dev/null +++ b/Tracker/TrackerEventCnv/TrackerSimEventTPCnv/src/TrackerHits/FaserSiHitCollectionCnv_p3a.cxx @@ -0,0 +1,351 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TrackerSimEvent/FaserSiHit.h" +#include "TrackerSimEvent/FaserSiHitCollection.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollection_p3a.h" +#include "TrackerSimEventTPCnv/TrackerHits/FaserSiHitCollectionCnv_p3a.h" +#include "GeneratorObjects/HepMcParticleLink.h" + +#include <cmath> + +//CLHEP +#include "CLHEP/Geometry/Point3D.h" +// Gaudi +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/ThreadLocalContext.h" + +// Athena +#include "StoreGate/StoreGateSvc.h" +#include "AthenaKernel/ExtendedEventContext.h" + + +// * * * stolen from eflowRec * * * // +inline double phicorr(double a) +{ + if (a <= -M_PI) + { + return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI))); + } + else if (a > M_PI) + { + return a-(2*M_PI*floor((a+M_PI)/(2*M_PI))); + } + else + { + return a; + } +} + +// * * * stolen from eflowRec * * * // +inline double cycle(double a, double b) +{ + double del = b-a; + if (del > M_PI) + { + return a+2.0*M_PI; + } + else if (del < -M_PI) + { + return a-2.0*M_PI; + } + else + { + return a; + } +} + + +const double FaserSiHitCollectionCnv_p3a::m_persEneUnit = 1.0e-5; +const double FaserSiHitCollectionCnv_p3a::m_persLenUnit = 1.0e-5; +const double FaserSiHitCollectionCnv_p3a::m_persAngUnit = 1.0e-5; +const double FaserSiHitCollectionCnv_p3a::m_2bHalfMaximum = pow(2.0, 15.0); +const int FaserSiHitCollectionCnv_p3a::m_2bMaximum = (unsigned short)(-1); + + +void FaserSiHitCollectionCnv_p3a::transToPers(const FaserSiHitCollection* transCont, FaserSiHitCollection_p3a* persCont, MsgStream &log) +{ + // Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and + // persistifies the end point of each hit plus the start point of the first hit in each string. + // + // Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of + // first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and + // used to specify the hit direction. All hit lengths are stored as 2 byte numbers. + // + // Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean + // time of the first hit per string. + // + // See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info. + + static const double dRcut = 1.0e-7; + static const double dTcut = 1.0; + + const EventContext& ctx = Gaudi::Hive::currentContext(); + const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy(); + const HepMcParticleLink * lastLink=nullptr; + int lastId = -1; + double stringFirstTheta = 0.0; + double stringFirstPhi = 0.0; + double lastT = 0.0; + double persSumE = 0.0; + double transSumE = 0.0; + unsigned int idx = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + unsigned int endHit = 0; + HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0); + HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0); + + for (FaserSiHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) { + + FaserSiHitCollection::const_iterator siHit = it; + + + if ( !lastLink || (siHit->particleLink() != *lastLink) ) { + + // store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode + + lastLink = &(siHit->particleLink()); + persCont->m_barcode.push_back(lastLink->barcode()); + // log << MSG::ALWAYS << "TP: Event collection from link: " << lastLink->getEventCollection() << " ; char value: " << lastLink->getEventCollectionAsChar() << " barcode: " << lastLink->barcode() << endmsg; + unsigned short index{0}; + const HepMcParticleLink::index_type position = + HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(), + lastLink->getEventCollection(), + proxy).at(0); + if (position!=0) { + index = lastLink->eventIndex(); + if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) { + log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg; + } + } + persCont->m_mcEvtIndex.push_back(index); + persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar()); + + if (idx > 0) { + persCont->m_nBC.push_back(idx - endBC); + endBC = idx; + } + } + + if ( (int)siHit->identify() != lastId ) { + + // store id once for set of consecutive hits with same barcode + + lastId = siHit->identify(); + persCont->m_id.push_back(lastId); + + if (idx > 0) { + persCont->m_nId.push_back(idx - endId); + endId = idx; + } + } + + HepGeom::Point3D<double> st = siHit->localStartPosition(); + HepGeom::Point3D<double> en = siHit->localEndPosition(); + + const double dx = st.x() - lastTransEnd.x(); + const double dy = st.y() - lastTransEnd.y(); + const double dz = st.z() - lastTransEnd.z(); + const double t = siHit->meanTime(); + + const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one + const double dTLast = fabs(t - lastT); + + CLHEP::Hep3Vector direction(0.0, 0.0, 0.0); + double theta = 0.0; + double phi = 0.0; + bool startNewString = false; + + if (dRLast < dRcut && dTLast < dTcut) { + + // hit is part of existing string + + direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() ); + + theta = direction.theta(); + phi = phicorr( direction.phi() ); + + const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 ); + const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 ); + + if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) { + persCont->m_dTheta.push_back(dTheta_2b); + persCont->m_dPhi.push_back(dPhi_2b); + theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit; + phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit; + phi = phicorr(phi); + } + else { + startNewString = true; + } + } + + if (startNewString || dRLast >= dRcut || dTLast >= dTcut) { + + // begin new hit string + + direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() ); + + theta = direction.theta(); + phi = phicorr( direction.phi() ); + + persCont->m_hit1_meanTime.push_back(t); + persCont->m_hit1_x0.push_back(st.x()); + persCont->m_hit1_y0.push_back(st.y()); + persCont->m_hit1_z0.push_back(st.z()); + persCont->m_hit1_theta.push_back(theta); + persCont->m_hit1_phi.push_back(phi); + + lastPersEnd = HepGeom::Point3D<double>(st.x(), st.y(), st.z()); + + stringFirstTheta = theta; + stringFirstPhi = phi; + + if (idx > 0) { + persCont->m_nHits.push_back(idx - endHit); + endHit = idx; + } + } + + lastTransEnd = HepGeom::Point3D<double>(en.x(), en.y(), en.z()); + transSumE += siHit->energyLoss(); + + const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over + // whole hit string to chosen precision + + const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to + // the chosen precision, NOT the length of the + // hit (small difference in practice). + double eneLoss = 0.0; + + if (eneLoss_2b >= m_2bMaximum) { + eneLoss = siHit->energyLoss(); + persCont->m_hitEne_2b.push_back(m_2bMaximum); + persCont->m_hitEne_4b.push_back(eneLoss); + } + else { + eneLoss = eneLoss_2b * m_persEneUnit; + persCont->m_hitEne_2b.push_back(eneLoss_2b); + } + + double length = 0.0; + + if (hitLength_2b >= m_2bMaximum) { + length = direction.mag(); + persCont->m_hitLength_2b.push_back(m_2bMaximum); + persCont->m_hitLength_4b.push_back(direction.mag()); + } + else { + length = hitLength_2b * m_persLenUnit; + persCont->m_hitLength_2b.push_back(hitLength_2b); + } + + CLHEP::Hep3Vector persDir(length, 0.0, 0.0); + persDir.setTheta(theta); + persDir.setPhi(phi); + + lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir; + persSumE += eneLoss; + lastT = t; + + ++idx; + } + + persCont->m_nBC.push_back(idx - endBC); + persCont->m_nId.push_back(idx - endId); + persCont->m_nHits.push_back(idx - endHit); + + // log << MSG::ALWAYS << "nBC: " << persCont->m_nBC << endmsg; + // log << MSG::ALWAYS << "nId: " << persCont->m_nId << endmsg; + // log << MSG::ALWAYS << "nHits:" << persCont->m_nHits << endmsg; + +} + + +FaserSiHitCollection* FaserSiHitCollectionCnv_p3a::createTransient(const FaserSiHitCollection_p3a* persObj, MsgStream &log) { + // for (size_t i = 0; i < persObj->m_evtColl.size(); i++) + // { + // if (persObj->m_evtColl[i] != 'a') + // log << MSG::ALWAYS << "Corrupted in createTransient: " << persObj->m_evtColl[i] << endmsg; + // } + std::unique_ptr<FaserSiHitCollection> trans(std::make_unique<FaserSiHitCollection>("DefaultCollectionName",persObj->m_nHits.size())); + persToTrans(persObj, trans.get(), log); + return(trans.release()); +} + + +void FaserSiHitCollectionCnv_p3a::persToTrans(const FaserSiHitCollection_p3a* persCont, FaserSiHitCollection* transCont, MsgStream &/*log*/) +{ + const EventContext& ctx = Gaudi::Hive::currentContext(); + + unsigned int hitCount = 0; + unsigned int angleCount = 0; + unsigned int idxBC = 0; + unsigned int idxId = 0; + unsigned int idxEne4b = 0; + unsigned int idxLen4b = 0; + unsigned int endHit = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + + // log << MSG::ALWAYS << "nHits.size(): " << persCont->m_nHits.size() << " nBC.size(): " << persCont->m_nBC.size() << endmsg; + // log << MSG::ALWAYS << " nBC: " << persCont->m_nBC << endmsg; + + for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) { + if (persCont->m_nHits[i]) { + const unsigned int start = endHit; + endHit += persCont->m_nHits[i]; + // log << MSG::ALWAYS << "i: " << i << " persCont->m_nHits[i]: " << persCont->m_nHits[i] << " start: " << start << " endHit: " << endHit << endmsg; + + const double t0 = persCont->m_hit1_meanTime[i]; + const double theta0 = persCont->m_hit1_theta[i]; + const double phi0 = persCont->m_hit1_phi[i]; + HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]); + + for (unsigned int j = start; j < endHit; j++) { + + if (j >= endBC + persCont->m_nBC[idxBC]) + endBC += persCont->m_nBC[idxBC++]; + + if (j >= endId + persCont->m_nId[idxId]) + endId += persCont->m_nId[idxId++]; + + const double eneLoss_2b = persCont->m_hitEne_2b[hitCount]; + const double hitLength_2b = persCont->m_hitLength_2b[hitCount]; + + const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++]; + const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++]; + + const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0; + const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0; + + const double meanTime = t0; + const double theta = theta0 + dTheta; + const double phi = phicorr(phi0 + dPhi); + + CLHEP::Hep3Vector r(length, 0.0, 0.0); + r.setTheta(theta); + r.setPhi(phi); + + HepGeom::Point3D<double> endThis( endLast + r ); + + HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX; + if (persCont->m_mcEvtIndex[idxBC] == 0) { + flag = HepMcParticleLink::IS_POSITION; + } + // log << MSG::ALWAYS << "PT: Event collection from persCont: " << persCont->m_evtColl[idxBC] << " ; idxBC: " << idxBC << " evtColl.size(): " << persCont->m_evtColl.size() << " evtIndex.size(): " << persCont->m_mcEvtIndex.size() << endmsg; + + HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx ); + transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]); + + endLast = endThis; + + ++hitCount; + if (j > start) ++angleCount; + } + } + } + // log << MSG::ALWAYS << "hitCount: " << hitCount << endmsg; +}