diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCnv_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCnv_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..27da6d00141b7dc2c44a1bd3b9c91dce4dfb23c2 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCnv_p1.h @@ -0,0 +1,34 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SIHITCNV_P1_H +#define SIHITCNV_P1_H + +/* +Transient/Persistent converter for SiHit class +Author: Davide Costanzo +*/ + +#include "InDetSimEvent/SiHit.h" +#include "SiHit_p1.h" + +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" + +class MsgStream; + + +class SiHitCnv_p1 : public T_AthenaPoolTPCnvBase<SiHit, SiHit_p1> +{ +public: + + SiHitCnv_p1() {} + + virtual void persToTrans(const SiHit_p1* persObj, SiHit* +transObj, MsgStream &log); + virtual void transToPers(const SiHit* transObj, SiHit_p1* +persObj, MsgStream &log); +}; + + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..69dafafadeb2f83f9e692059e0d307dd748ab9f8 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p1.h @@ -0,0 +1,12 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#include "SiHitCollection_p1.h" +#undef private +#include "InDetSimEvent/SiHitCollection.h" +#include "SiHitCnv_p1.h" +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" + +typedef T_AtlasHitsVectorCnv< SiHitCollection, SiHitCollection_p1, SiHitCnv_p1 > SiHitCollectionCnv_p1; diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p2.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p2.h new file mode 100755 index 0000000000000000000000000000000000000000..99119aaeede552e50f53d934bc43d5b90c4fefc4 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p2.h @@ -0,0 +1,40 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SIHITCOLLECTIONCNV_P2_H +#define SIHITCOLLECTIONCNV_P2_H + +// SiHitCollectionCnv_p2, T/P separation of Si Hits +// author D.Costanzo <davide.costanzo@cern.ch> + +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" +#include "InDetSimEvent/SiHitCollection.h" +#include "SiHitCollection_p2.h" + + +class SiHitCollectionCnv_p2 : public T_AthenaPoolTPCnvBase<SiHitCollection, SiHitCollection_p2> +{ + public: + + SiHitCollectionCnv_p2() {}; + + virtual SiHitCollection* createTransient(const SiHitCollection_p2* persObj, MsgStream &log); + + virtual void persToTrans(const SiHitCollection_p2* persCont, + SiHitCollection* transCont, + MsgStream &log) ; + virtual void transToPers(const SiHitCollection* transCont, + SiHitCollection_p2* 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/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..72d4498ce18bce0694d7b08e0696100528e3589c --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p1.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SIHITCOLLECTION_P1_H +#define SIHITCOLLECTION_P1_H + +/* + +Persistent represenation of a SiHitContainer, +Author: Davide Costanzo + +*/ + +#include <vector> +#include <string> +#include "SiHit_p1.h" + + +class SiHitCollection_p1 +{ +public: + /// typedefs + typedef std::vector<SiHit_p1> HitVector; + typedef HitVector::const_iterator const_iterator; + typedef HitVector::iterator iterator; + + + /// Default constructor + SiHitCollection_p1 (); + + // Accessors + const std::string& name() const; + const HitVector& getVector() const; +private: + std::vector<SiHit_p1> m_cont; + std::string m_name; +}; + + +// inlines + +inline +SiHitCollection_p1::SiHitCollection_p1 () {} + +inline +const std::string& +SiHitCollection_p1::name() const +{return m_name;} + +inline +const std::vector<SiHit_p1>& +SiHitCollection_p1::getVector() const +{return m_cont;} + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p2.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p2.h new file mode 100755 index 0000000000000000000000000000000000000000..f682a768584f571543e7320a92aa9aa061ee9beb --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHitCollection_p2.h @@ -0,0 +1,55 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SIHITCOLLECTION_P2_H +#define SIHITCOLLECTION_P2_H + +/* + +Authors: Davide Costanzo Rob Duxfield + +*/ + +#include <vector> +#include <string> + +class SiHitCollection_p2 +{ + public: +/// Default constructor + SiHitCollection_p2 (); + 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 short> 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 short> m_nBC; + + std::vector<unsigned long> m_id; + std::vector<unsigned short> m_nId; +}; + + +// inlines + +inline +SiHitCollection_p2::SiHitCollection_p2 () {} + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHit_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHit_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..aff1dbaa4d72fc7d778182998bfbcce9bde896f1 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/SiHit_p1.h @@ -0,0 +1,23 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef SIHIT_P1_H +#define SIHIT_P1_H + +#include "GeneratorObjectsTPCnv/HepMcParticleLink_p1.h" + +class SiHit_p1 { + public: + SiHit_p1() {}; + friend class SiHitCnv_p1; + + private: + float m_stX, m_stY, m_stZ; + float m_enX, m_enY, m_enZ; + float m_energyLoss; // deposited energy + float m_meanTime; // time of energy deposition + HepMcParticleLink_p1 m_partLink; + unsigned int m_ID; +}; +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCnv_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCnv_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..f944e914b10ac66d782725bd6b1cdbf2dff9b881 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCnv_p1.h @@ -0,0 +1,34 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRT_HITCNV_P1_H +#define TRT_HITCNV_P1_H + +/* +Transient/Persistent converter for TRT_Hit class +Author: Davide Costanzo +*/ + +#include "InDetSimEvent/TRTUncompressedHit.h" +#include "TRT_Hit_p1.h" + +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" + +class MsgStream; + + +class TRT_HitCnv_p1 : public T_AthenaPoolTPCnvBase<TRTUncompressedHit, TRT_Hit_p1> +{ +public: + + TRT_HitCnv_p1() {} + + virtual void persToTrans(const TRT_Hit_p1* persObj, TRTUncompressedHit* +transObj, MsgStream &log); + virtual void transToPers(const TRTUncompressedHit* transObj, TRT_Hit_p1* +persObj, MsgStream &log); +}; + + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..422565ef68751c8b9807031b91ad500374bd6bc8 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p1.h @@ -0,0 +1,12 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#include "TRT_HitCollection_p1.h" +#undef private +#include "InDetSimEvent/TRTUncompressedHitCollection.h" +#include "TRT_HitCnv_p1.h" +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" + +typedef T_AtlasHitsVectorCnv< TRTUncompressedHitCollection, TRT_HitCollection_p1, TRT_HitCnv_p1 > TRT_HitCollectionCnv_p1; diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p2.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p2.h new file mode 100644 index 0000000000000000000000000000000000000000..86eefd209f633e5330c7a06e7ababf73b0e5574e --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p2.h @@ -0,0 +1,32 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRTHITCOLLECTIONCNV_P2_H +#define TRTHITCOLLECTIONCNV_P2_H + +// SiHitCollectionCnv_p2, T/P separation of Si Hits +// author Robert Duxfield <r.duxfield@sheffield.ac.uk> + +#include "InDetSimEvent/TRTUncompressedHitCollection.h" +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" +#include "TRT_HitCollection_p2.h" + + +class TRT_HitCollectionCnv_p2 : public T_AthenaPoolTPCnvBase<TRTUncompressedHitCollection, TRT_HitCollection_p2> +{ + public: + + TRT_HitCollectionCnv_p2() {}; + + virtual TRTUncompressedHitCollection* createTransient(const TRT_HitCollection_p2* persObj, MsgStream &log); + + virtual void persToTrans(const TRT_HitCollection_p2* persCont, + TRTUncompressedHitCollection* transCont, + MsgStream &log) ; + virtual void transToPers(const TRTUncompressedHitCollection* transCont, + TRT_HitCollection_p2* persCont, + MsgStream &log) ; +}; + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p3.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p3.h new file mode 100644 index 0000000000000000000000000000000000000000..b3ae5d4b9d248d842a39ea106316e59c91c35451 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p3.h @@ -0,0 +1,32 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRTHITCOLLECTIONCNV_P3_H +#define TRTHITCOLLECTIONCNV_P3_H + +// SiHitCollectionCnv_p2, T/P separation of Si Hits +// author Robert Duxfield <r.duxfield@sheffield.ac.uk> + +#include "InDetSimEvent/TRTUncompressedHitCollection.h" +#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h" +#include "TRT_HitCollection_p3.h" + + +class TRT_HitCollectionCnv_p3 : public T_AthenaPoolTPCnvBase<TRTUncompressedHitCollection, TRT_HitCollection_p3> +{ + public: + + TRT_HitCollectionCnv_p3() {}; + + virtual TRTUncompressedHitCollection* createTransient(const TRT_HitCollection_p3* persObj, MsgStream &log); + + virtual void persToTrans(const TRT_HitCollection_p3* persCont, + TRTUncompressedHitCollection* transCont, + MsgStream &log) ; + virtual void transToPers(const TRTUncompressedHitCollection* transCont, + TRT_HitCollection_p3* persCont, + MsgStream &log) ; +}; + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..328a9950af27d59b2e0534dcb47565f1c00e4efa --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p1.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRT_HITCOLLECTION_P1_H +#define TRT_HITCOLLECTION_P1_H + +/* + +Persistent represenation of a TRT_HitContainer, +Author: Davide Costanzo + +*/ + +#include <vector> +#include <string> +#include "TRT_Hit_p1.h" + + +class TRT_HitCollection_p1 +{ +public: + /// typedefs + typedef std::vector<TRT_Hit_p1> HitVector; + typedef HitVector::const_iterator const_iterator; + typedef HitVector::iterator iterator; + + + /// Default constructor + TRT_HitCollection_p1 (); + + // Accessors + const std::string& name() const; + const HitVector& getVector() const; +private: + std::vector<TRT_Hit_p1> m_cont; + std::string m_name; +}; + + +// inlines + +inline +TRT_HitCollection_p1::TRT_HitCollection_p1 () {} + +inline +const std::string& +TRT_HitCollection_p1::name() const +{return m_name;} + +inline +const std::vector<TRT_Hit_p1>& +TRT_HitCollection_p1::getVector() const +{return m_cont;} + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p2.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p2.h new file mode 100644 index 0000000000000000000000000000000000000000..be9ba614f50187ca5a213ecd76ea9f930807783c --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p2.h @@ -0,0 +1,50 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRTHITCOLLECTION_P2_H +#define TRTHITCOLLECTION_P2_H + +/* + +Author: Rob Duxfield <r.duxfield@sheffield.ac.uk> + +*/ + +#include <vector> +#include <string> + +class TRT_HitCollection_p2 +{ + public: + /// Default constructor + TRT_HitCollection_p2 (); + private: + + std::vector<float> m_hit1_startX; // 1 element per string + std::vector<float> m_hit1_startY; // + std::vector<float> m_hit1_startZ; // + std::vector<unsigned short> m_nHits; // + + std::vector<float> m_hitEne; // 1 element per hit + std::vector<float> m_meanTime; // + std::vector<float> m_kinEne; // + std::vector<float> m_endX; // + std::vector<float> m_endY; // + std::vector<float> m_endZ; // + std::vector<unsigned long> m_hitId; // + + std::vector<unsigned long> m_barcode; + std::vector<unsigned short> m_nBC; + + std::vector<long> m_id; + std::vector<unsigned short> m_nId; +}; + + +// inlines + +inline +TRT_HitCollection_p2::TRT_HitCollection_p2 () {} + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p3.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p3.h new file mode 100644 index 0000000000000000000000000000000000000000..688c2ec1a08375e05849b28810dee01b16ee5262 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p3.h @@ -0,0 +1,64 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRTHITCOLLECTION_P3_H +#define TRTHITCOLLECTION_P3_H + +/* + Author: Rob Duxfield <r.duxfield@sheffield.ac.uk> Spring 2008 + _p3 Integer compression: <Andrew.Beddall@cern.ch> Spring 2009 + See http://cern.ch/beddall/TRThitCompression/ +*/ + +#include <vector> +#include <string> + +class TRT_HitCollection_p3 +{ + public: + /// Default constructor + TRT_HitCollection_p3 (); + + private: + + // + // 1 element per string (a string resides in one straw; there may be more than one string in a straw) + // + + std::vector<unsigned short> m_nHits; // number of hits in the string (0,1,2 ... ,hundreds). + std::vector<unsigned short> m_strawId2b; // straw id | 24-bit + std::vector<unsigned char> m_strawId1b; // straw id | integer. + std::vector<unsigned char> m_startR; // hit start radius (0, 2 mm) [not always stored]. + std::vector<unsigned char> m_startPhi; // hit start phi (-pi, pi). + std::vector<unsigned char> m_startZ; // hit start z (-365, +365 mm), and 1-bit startRflag. + + // + // 1 element per hit, there are typically 1 or 2 hits per string, but can be hundreds! + // + + std::vector<unsigned short> m_kinEne; // short float, kinematic energy of the particle causing the hit. + std::vector<unsigned short> m_steplength; // short float, g4 step length; endZ is derived from this. + std::vector<unsigned char> m_endR; // hit end radius (0, 2 mm) [Not always stored]. + std::vector<unsigned char> m_endPhi; // hit end phi (-pi, pi). + std::vector<unsigned short> m_meanTime; // time to center of the hit, and 1-bit idZsign and 1-bit endRflag. + std::vector<float> m_meanTimeof; // t >= 75 ns overflow to a float. + + // + // much less frequent + // + + std::vector<float> m_hitEne; // energy deposited; *only stored for photons* (m_id=22) + std::vector<unsigned short> m_nId; + std::vector<unsigned int> m_barcode; + std::vector<unsigned short> m_nBC; + std::vector<int> m_id; // particle code. + +}; + +// inlines + +inline +TRT_HitCollection_p3::TRT_HitCollection_p3 () {} + +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_Hit_p1.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_Hit_p1.h new file mode 100755 index 0000000000000000000000000000000000000000..b1171be05e648f346bf1a6d456b059128c94144f --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetHits/TRT_Hit_p1.h @@ -0,0 +1,29 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRT_HIT_P1_H +#define TRT_HIT_P1_H + +#include "GeneratorObjectsTPCnv/HepMcParticleLink_p1.h" + +class TRT_Hit_p1 { + public: + TRT_Hit_p1() {}; + friend class TRT_HitCnv_p1; + + private: + int hitID; // To identify the hit + HepMcParticleLink_p1 m_partLink; // link to the particle generating the hit + int particleEncoding; // PDG id + float kineticEnergy; // kin energy of the particle + float energyDeposit; // energy deposit by the hit + float preStepX; + float preStepY; + float preStepZ; + float postStepX; + float postStepY; + float postStepZ; + float globalTime; +}; +#endif diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetSimEventTPCnvDict.h b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetSimEventTPCnvDict.h new file mode 100644 index 0000000000000000000000000000000000000000..c75cae12da1d0e740b88be487a7a8d81c4d98ce4 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/InDetSimEventTPCnvDict.h @@ -0,0 +1,30 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef INDETEVENTTPCNV_INDETSIMEVENTTPCNVDICT_H +#define INDETEVENTTPCNV_INDETSIMEVENTTPCNVDICT_H + +//----------------------------------------------------------------------------- +// +// file: InDetSimEventTPCnvDict_p1.h +// +//----------------------------------------------------------------------------- + + +#include "InDetSimEventTPCnv/InDetHits/SiHitCnv_p1.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p1.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p2.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCollection_p1.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCollection_p2.h" +#include "InDetSimEventTPCnv/InDetHits/SiHit_p1.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCnv_p1.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p1.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p2.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p3.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p1.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p2.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p3.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_Hit_p1.h" + +#endif // INDETEVENTTPCNV_INDETSIMEVENTTPCNVDICT_H diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/selection.xml b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/selection.xml new file mode 100755 index 0000000000000000000000000000000000000000..92b636b81e4aa10fa880c38ce0167be3b1ea67fb --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/InDetSimEventTPCnv/selection.xml @@ -0,0 +1,14 @@ +<lcgdict> + + <!-- InDetHits --> + <class name="SiHit_p1" /> + <class name="std::vector<SiHit_p1>" /> + <class name="SiHitCollection_p1" id="36D1FF8E-5734-4A93-A133-F286CF47DB72" /> + <class name="SiHitCollection_p2" id="BD1469C5-C904-40B8-82B9-43D25888D884" /> + <class name="TRT_Hit_p1" /> + <class name="std::vector<TRT_Hit_p1>" /> + <class name="TRT_HitCollection_p1" id="6688E934-157E-421A-B6D1-A35FC8BD651C" /> + <class name="TRT_HitCollection_p2" id="473FF621-3466-4D87-9469-4780A6A77023" /> + <class name="TRT_HitCollection_p3" id="FB5F5BFC-43E5-44E1-B79C-C330C1480E2E" /> + +</lcgdict> diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/cmt/requirements b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/cmt/requirements new file mode 100755 index 0000000000000000000000000000000000000000..a1ce407981024866c1c5e453d9125295dfa47f05 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/cmt/requirements @@ -0,0 +1,35 @@ +package InDetSimEventTPCnv + +author Andreas Wildauer <Andreas.Wildauer@cern.ch> +author Edward Moyse <edward.moyse@cern.ch> + +public +use AtlasPolicy AtlasPolicy-* +use AthenaPoolCnvSvc AthenaPoolCnvSvc-* Database/AthenaPOOL +use GaudiInterface GaudiInterface-* External +use GeneratorObjectsTPCnv GeneratorObjectsTPCnv-* Generators +use InDetSimEvent InDetSimEvent-* InnerDetector + +private +use AtlasCLHEP AtlasCLHEP-* External +use AtlasReflex AtlasReflex-* External -no_auto_imports +use Identifier Identifier-* DetectorDescription +use StoreGate StoreGate-* Control +end_private + +public + +apply_pattern tpcnv_library +library InDetSimEventTPCnv *.cxx \ + InDetHits/*.cxx + + +# The following use is to get the lcgdict pattern. +# This is "hidden" behind "private" and "no_auto_imports" to keep +# clients of EventInfo from seeing excess dependencies +private +use AtlasReflex AtlasReflex-* External -no_auto_imports + +# Pattern to build the dict lib. User should create a single heade +# file: <package>Dict.h which includes all other .h files. See EventInfoDict +apply_pattern lcgdict dict=InDetSimEventTPCnv selectionfile=selection.xml headerfiles="../InDetSimEventTPCnv/InDetSimEventTPCnvDict.h" \ No newline at end of file diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCnv_p1.cxx b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCnv_p1.cxx new file mode 100755 index 0000000000000000000000000000000000000000..c9052504924e0dfd57dba30b55bdf282eee23e72 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCnv_p1.cxx @@ -0,0 +1,58 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#define protected public +#include "InDetSimEvent/SiHit.h" +#undef private +#undef protected +#include "Identifier/Identifier.h" +#include "GeneratorObjectsTPCnv/HepMcParticleLinkCnv_p1.h" + +#include "InDetSimEventTPCnv/InDetHits/SiHit_p1.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCnv_p1.h" + + +void +SiHitCnv_p1::persToTrans(const SiHit_p1* persObj, SiHit* transObj, +MsgStream &log) +{ +// if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "SiHitCnv_p1::persToTrans called " << endreq; + HepMcParticleLinkCnv_p1 HepMcPLCnv; + + transObj->m_stX = persObj->m_stX; + transObj->m_stY = persObj->m_stY; + transObj->m_stZ = persObj->m_stZ; + + transObj->m_enX = persObj->m_enX; + transObj->m_enY = persObj->m_enY; + transObj->m_enZ = persObj->m_enZ; + + transObj->m_energyLoss = persObj->m_energyLoss; + transObj->m_meanTime = persObj->m_meanTime; + transObj->m_ID = persObj->m_ID; + HepMcPLCnv.persToTrans(&(persObj->m_partLink),&(transObj->m_partLink), log); +} + + +void +SiHitCnv_p1::transToPers(const SiHit* transObj, SiHit_p1* persObj, +MsgStream &log) +{ +// if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "SiHitCnv_p1::transToPers called " << endreq; + HepMcParticleLinkCnv_p1 HepMcPLCnv; + + persObj->m_stX = transObj->m_stX; + persObj->m_stY = transObj->m_stY; + persObj->m_stZ = transObj->m_stZ; + + persObj->m_enX = transObj->m_enX; + persObj->m_enY = transObj->m_enY; + persObj->m_enZ = transObj->m_enZ; + + persObj->m_energyLoss = transObj->m_energyLoss; + persObj->m_meanTime = transObj->m_meanTime; + persObj->m_ID = transObj->m_ID; + HepMcPLCnv.transToPers(&(transObj->m_partLink),&(persObj->m_partLink), log); +} diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCollectionCnv_p2.cxx b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCollectionCnv_p2.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9557dcd0e91adac8d6bd6628d18783ba3019bb90 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/SiHitCollectionCnv_p2.cxx @@ -0,0 +1,311 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#define protected public +#include "InDetSimEvent/SiHit.h" +#include "InDetSimEvent/SiHitCollection.h" +#include "InDetSimEventTPCnv/InDetHits/SiHitCollection_p2.h" +#undef private +#undef protected +#include "InDetSimEventTPCnv/InDetHits/SiHitCollectionCnv_p2.h" + +#include <cmath> + +//CLHEP +#include "CLHEP/Geometry/Point3D.h" +// Gaudi +#include "GaudiKernel/MsgStream.h" +// Athena +#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 SiHitCollectionCnv_p2::m_persEneUnit = 1.0e-5; +const double SiHitCollectionCnv_p2::m_persLenUnit = 1.0e-5; +const double SiHitCollectionCnv_p2::m_persAngUnit = 1.0e-5; +const double SiHitCollectionCnv_p2::m_2bHalfMaximum = pow(2.0, 15.0); +const int SiHitCollectionCnv_p2::m_2bMaximum = (unsigned short)(-1); + + +void SiHitCollectionCnv_p2::transToPers(const SiHitCollection* transCont, SiHitCollection_p2* 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; + + int lastBarcode = -1; + 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 (SiHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) { + + // const SiHit* siHit = *it; + + SiHitCollection::const_iterator siHit = it; + + + if ( siHit->m_partLink.barcode() != lastBarcode ) { + + // store barcode once for set of consecutive hits with same barcode + + lastBarcode = siHit->m_partLink.barcode(); + persCont->m_barcode.push_back(lastBarcode); + + if (idx > 0) { + persCont->m_nBC.push_back(idx - endBC); + endBC = idx; + } + } + + if ( (int)siHit->m_ID != lastId ) { + + // store id once for set of consecutive hits with same barcode + + lastId = siHit->m_ID; + persCont->m_id.push_back(lastId); + + if (idx > 0) { + persCont->m_nId.push_back(idx - endId); + endId = idx; + } + } + + const double dx = siHit->m_stX - lastTransEnd.x(); + const double dy = siHit->m_stY - lastTransEnd.y(); + const double dz = siHit->m_stZ - lastTransEnd.z(); + const double t = siHit->m_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( siHit->m_enX - lastPersEnd.x(), siHit->m_enY - lastPersEnd.y(), siHit->m_enZ - 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( siHit->m_enX - siHit->m_stX, siHit->m_enY - siHit->m_stY, siHit->m_enZ - siHit->m_stZ ); + + theta = direction.theta(); + phi = phicorr( direction.phi() ); + + persCont->m_hit1_meanTime.push_back(t); + persCont->m_hit1_x0.push_back(siHit->m_stX); + persCont->m_hit1_y0.push_back(siHit->m_stY); + persCont->m_hit1_z0.push_back(siHit->m_stZ); + persCont->m_hit1_theta.push_back(theta); + persCont->m_hit1_phi.push_back(phi); + + lastPersEnd = HepGeom::Point3D<double>(siHit->m_stX, siHit->m_stY, siHit->m_stZ); + + stringFirstTheta = theta; + stringFirstPhi = phi; + + if (idx > 0) { + persCont->m_nHits.push_back(idx - endHit); + endHit = idx; + } + } + + lastTransEnd = HepGeom::Point3D<double>(siHit->m_enX, siHit->m_enY, siHit->m_enZ); + transSumE += siHit->m_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->m_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); +} + + +SiHitCollection* SiHitCollectionCnv_p2::createTransient(const SiHitCollection_p2* persObj, MsgStream &log) { + std::auto_ptr<SiHitCollection> trans(new SiHitCollection("DefaultCollectionName",persObj->m_nHits.size())); + persToTrans(persObj, trans.get(), log); + return(trans.release()); +} + + +void SiHitCollectionCnv_p2::persToTrans(const SiHitCollection_p2* persCont, SiHitCollection* transCont, MsgStream &/*log*/) +{ + 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 ); + + transCont->Emplace( endLast, endThis, eneLoss, meanTime, persCont->m_barcode[idxBC], persCont->m_id[idxId]); + + endLast = endThis; + + ++hitCount; + if (j > start) ++angleCount; + } + } + } +} diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCnv_p1.cxx b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCnv_p1.cxx new file mode 100755 index 0000000000000000000000000000000000000000..f27cf61f971454f18824c7ca3239547081bc6627 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCnv_p1.cxx @@ -0,0 +1,62 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#define protected public +#include "InDetSimEvent/TRTUncompressedHit.h" +#undef private +#undef protected +#include "Identifier/Identifier.h" +#include "GeneratorObjectsTPCnv/HepMcParticleLinkCnv_p1.h" + +#include "InDetSimEventTPCnv/InDetHits/TRT_Hit_p1.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCnv_p1.h" + + +void +TRT_HitCnv_p1::persToTrans(const TRT_Hit_p1* persObj, TRTUncompressedHit* transObj, +MsgStream &log) +{ +// if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "TRT_HitCnv_p1::persToTrans called " << endreq; + HepMcParticleLinkCnv_p1 HepMcPLCnv; + transObj->hitID = persObj-> hitID; + HepMcPLCnv.persToTrans(&(persObj->m_partLink),&(transObj->m_partLink), log); + transObj->particleEncoding = persObj->particleEncoding; + transObj->kineticEnergy = persObj->kineticEnergy; + transObj->energyDeposit = persObj->energyDeposit; + + transObj->preStepX = persObj->preStepX; + transObj->preStepY = persObj->preStepY; + transObj->preStepZ = persObj->preStepZ; + + transObj->postStepX = persObj->postStepX; + transObj->postStepY = persObj->postStepY; + transObj->postStepZ = persObj->postStepZ; + transObj->globalTime = persObj->globalTime; + +} + + +void +TRT_HitCnv_p1::transToPers(const TRTUncompressedHit* transObj, TRT_Hit_p1* persObj, +MsgStream &log) +{ +// if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "TRT_HitCnv_p1::transToPers called " << endreq; + HepMcParticleLinkCnv_p1 HepMcPLCnv; + persObj->hitID = transObj-> hitID; + HepMcPLCnv.transToPers(&(transObj->m_partLink),&(persObj->m_partLink), log); + persObj->particleEncoding = transObj->particleEncoding; + persObj->kineticEnergy = transObj->kineticEnergy; + persObj->energyDeposit = transObj->energyDeposit; + + persObj->preStepX = transObj->preStepX; + persObj->preStepY = transObj->preStepY; + persObj->preStepZ = transObj->preStepZ; + + persObj->postStepX = transObj->postStepX; + persObj->postStepY = transObj->postStepY; + persObj->postStepZ = transObj->postStepZ; + persObj->globalTime = transObj->globalTime; + +} diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p2.cxx b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p2.cxx new file mode 100644 index 0000000000000000000000000000000000000000..34c2cd85df87b835010de0aeec95d7682eba6572 --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p2.cxx @@ -0,0 +1,174 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#define protected public +#include "InDetSimEvent/TRTUncompressedHit.h" +#include "InDetSimEvent/TRTUncompressedHitCollection.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p2.h" +#undef private +#undef protected +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p2.h" + +#include <cmath> + +// CLHEP +#include "CLHEP/Geometry/Point3D.h" + +// Gaudi +#include "GaudiKernel/MsgStream.h" + +// Athena +#include "StoreGate/StoreGateSvc.h" + + + + + +void TRT_HitCollectionCnv_p2::transToPers(const TRTUncompressedHitCollection* transCont, TRT_HitCollection_p2* 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; + + // if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "In TRT_HitCollectionCnv_p2::transToPers()" << endreq; + + int lastBarcode = -1; + int lastId = -1; + double lastT = 0.0; + unsigned int idx = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + unsigned int endHit = 0; + HepGeom::Point3D<double> lastEnd(0.0, 0.0, 0.0); + + for (TRTUncompressedHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) { + + + if ( it->m_partLink.barcode() != lastBarcode ) { + + // store barcode once for set of consecutive hits with same barcode + + lastBarcode = it->m_partLink.barcode(); + persCont->m_barcode.push_back(lastBarcode); + + if (idx > 0) { + persCont->m_nBC.push_back(idx - endBC); + endBC = idx; + } + } + + if ( (int)it->particleEncoding != lastId ) { + + // store id once for set of consecutive hits with same id + + lastId = it->particleEncoding; + persCont->m_id.push_back(lastId); + + if (idx > 0) { + persCont->m_nId.push_back(idx - endId); + endId = idx; + } + } + + const HepGeom::Point3D<double> hitStart(it->preStepX, it->preStepY, it->preStepZ); + const double t = it->globalTime; + + const double dRLast = lastEnd.distance(hitStart); // dR between end of previous hit and start of current one + const double dTLast = fabs(t - lastT); + + if (dRLast >= dRcut || dTLast >= dTcut) { + + // begin new hit string + + persCont->m_hit1_startX.push_back( hitStart.x() ); + persCont->m_hit1_startY.push_back( hitStart.y() ); + persCont->m_hit1_startZ.push_back( hitStart.z() ); + + if (idx > 0) { + persCont->m_nHits.push_back(idx - endHit); + endHit = idx; + } + } + + persCont->m_hitId.push_back( it->hitID ); + persCont->m_kinEne.push_back( it->kineticEnergy ); + persCont->m_meanTime.push_back( t ); + + persCont->m_hitEne.push_back( it->energyDeposit ); + persCont->m_endX.push_back( it->postStepX ); + persCont->m_endY.push_back( it->postStepY ); + persCont->m_endZ.push_back( it->postStepZ ); + + lastEnd = HepGeom::Point3D<double>(it->postStepX, it->postStepY, it->postStepZ); + lastT = t; + ++idx; + } + + persCont->m_nBC.push_back(idx - endBC); + persCont->m_nId.push_back(idx - endId); + persCont->m_nHits.push_back(idx - endHit); +} + + +TRTUncompressedHitCollection* TRT_HitCollectionCnv_p2::createTransient(const TRT_HitCollection_p2* persObj, MsgStream &log) { + std::auto_ptr<TRTUncompressedHitCollection> trans(new TRTUncompressedHitCollection("DefaultCollectionName",persObj->m_nHits.size())); + persToTrans(persObj, trans.get(), log); + return(trans.release()); +} + + +void TRT_HitCollectionCnv_p2::persToTrans(const TRT_HitCollection_p2* persCont, TRTUncompressedHitCollection* transCont, MsgStream& /*log*/) +{ + // if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "In TRT_HitCollectionCnv_p2::persToTrans()" << endreq; + + unsigned int hitCount = 0; + unsigned int idxBC = 0; + unsigned int idxId = 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]; + + HepGeom::Point3D<double> endLast(persCont->m_hit1_startX[i], persCont->m_hit1_startY[i], persCont->m_hit1_startZ[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 = persCont->m_hitEne[hitCount]; + + HepGeom::Point3D<double> endThis(persCont->m_endX[j], persCont->m_endY[j], persCont->m_endZ[j]); + + transCont->Emplace( persCont->m_hitId[hitCount], persCont->m_barcode[idxBC], persCont->m_id[idxId], persCont->m_kinEne[hitCount], + eneLoss, endLast.x(), endLast.y(), endLast.z(), endThis.x(), endThis.y(), endThis.z(), + persCont->m_meanTime[hitCount] ); + + endLast = endThis; + ++hitCount; + } + } + } +} diff --git a/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p3.cxx b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p3.cxx new file mode 100644 index 0000000000000000000000000000000000000000..23e8ce7818b9ad3369a398be52bb8310efa7236d --- /dev/null +++ b/InnerDetector/InDetEventCnv/InDetSimEventTPCnv/src/InDetHits/TRT_HitCollectionCnv_p3.cxx @@ -0,0 +1,498 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#define private public +#define protected public +#include "InDetSimEvent/TRTUncompressedHit.h" +#include "InDetSimEvent/TRTUncompressedHitCollection.h" +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollection_p3.h" +#undef private +#undef protected +#include "InDetSimEventTPCnv/InDetHits/TRT_HitCollectionCnv_p3.h" + +#include <cmath> + +// CLHEP +#include "CLHEP/Geometry/Point3D.h" +#include "CLHEP/Units/SystemOfUnits.h" + +// Gaudi +#include "GaudiKernel/MsgStream.h" + +// Athena +#include "StoreGate/StoreGateSvc.h" + +// Transient(Geant) to Persistent(Disk) +void TRT_HitCollectionCnv_p3::transToPers(const TRTUncompressedHitCollection* transCont, TRT_HitCollection_p3* persCont, MsgStream &log) +{ + + /* + Spring 2009 + Andrew Beddall - lossy TRT G4hit compression [p3] + + In p1, p2 versions, GEANT hits are persistified on disk as floats. + In this p3 version, floats are compressed to "integers"/"short-floats" before persistifying. + The saving is about 75%; see http://cern.ch/beddall/TRThitCompression/ + + Spring 2008 + Rob Duxfield - lossless TRT G4hit compression [p2] + + 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. + */ + + // The original units from the hit simulation are indicated in comments; + // they are all in CLHEP units except for hitEne which is in keV. + // I sometimes make use of CLHEP scales *CLHEP::mm and *CLHEP::ns (both=1) for clarity (I hope!). + // See also https://twiki.cern.ch/twiki/bin/view/Atlas/TrtSoftware#Production_of_Hits + + static const double dRcut = 1.0e-7*CLHEP::mm; + static const double dTcut = 1.0*CLHEP::ns; // redundant? + + // if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "In TRT_HitCollectionCnv_p3::transToPers()" << endreq; + + int lastBarcode = -1; + int lastId = -1; + double lastT = 0.0*CLHEP::ns; + unsigned int idx = 0; + unsigned int endBC = 0; + unsigned int endId = 0; + unsigned int endHit = 0; + HepGeom::Point3D<double> lastEnd(0.0, 0.0, 0.0); // mm + + for (TRTUncompressedHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) { + + // const TRTUncompressedHit* trtHit = *it; + TRTUncompressedHitCollection::const_iterator trtHit = it; + + if ( trtHit->m_partLink.barcode() != lastBarcode || idx - endBC > 65500) { // max unsigned short = 65535; + // store barcode once for set of consecutive hits with same barcode + lastBarcode = trtHit->m_partLink.barcode(); + persCont->m_barcode.push_back(lastBarcode); + if ( idx > 0 ) { + persCont->m_nBC.push_back(idx - endBC); + endBC = idx; + } + } + + if ( (int)trtHit->particleEncoding != lastId || idx - endId > 65500) { // max unsigned short = 65535; + // store id once for set of consecutive hits with same id + lastId = trtHit->particleEncoding; + persCont->m_id.push_back(lastId); + if ( idx > 0 ) { + persCont->m_nId.push_back(idx - endId); + endId = idx; + } + } + + const HepGeom::Point3D<double> hitStart(trtHit->preStepX, trtHit->preStepY, trtHit->preStepZ); // mm + + const double meanTime = trtHit->globalTime; // ns // Time of flight from the I.P. to the center of the hit. + const double dTLast = fabs(meanTime - lastT); // |d(meantime)| between the previous hit and the current one. + const double dRLast = lastEnd.distance(hitStart); // Distance between end of previous hit and start of current one; + // this is zero if the hit is a continuation of the same particle in the same straw. + + // Begin a new string if the current and previous hits are disconnected; + // it looks like dTcut is redundant (but not sure about this). + if ( dRLast >= dRcut || dTLast >= dTcut ) { + + // if ( dRLast < dRcut) std::cout << "AJBdTLastTriggeredNewString " << dRLast << " " << dTLast << std::endl; + + //////////////////// + // new hit string // + //////////////////// + + // + // Persistify string *strawId* using 24 bits. + // Assumes 0 <= strawId <= 16,777,215 (strawId appears to be < 4,000,000) + // + const unsigned int strawId = trtHit->hitID; + persCont->m_strawId1b.push_back( (unsigned char)(strawId % 256) ); // 8 bits + persCont->m_strawId2b.push_back( (unsigned short)(strawId / 256) ); // 16 bits + if ( strawId>16777215 ) + log << MSG::WARNING << "TRT_HitCollectionCnv: strawId > 2^24-1 cannot be persistified correctly! " << endreq; + + // + // Persistify string start radius using 1 bit (istartRflag) or 8 bits (startR) + // Note that the smallest value of R is the wire radius (0.0155 mm) + // + // R will be flagged as 2 mm if it is within 0.1 um of the straw wall => max error = 0.1 um, + // otherwise compress with 8 bits => max error = 3.9 um (0.078 ns), RMS error = 1.1 um (0.022 ns) + // + const double startR = sqrt( hitStart.x()*hitStart.x() + hitStart.y()*hitStart.y() ); // mm + unsigned short istartRflag; + if ( startR > 1.9999*CLHEP::mm ) { + istartRflag=1; // persistify as a 1-bit flag + } + else { + istartRflag=0; // compress to 8 bits with a span of 2 mm + persCont->m_startR.push_back( (unsigned char)(startR/(2.0*CLHEP::mm)*256.0) ); + } + + // + // Persistify string *startPhi* using 8 bits (min=-pi, max=+pi) + // Max. error = 12 mrad (< 24 um, 0.48 ns); RMS error = 7 mrad (< 14 um, 0.28 ns) + // + const double startPhi = atan2( hitStart.y(), hitStart.x() ); // returns range -pi to +pi rad + persCont->m_startPhi.push_back( (unsigned char)( (startPhi+M_PI)/(2.0*M_PI)*256.0 ) ); + + // + // Persistify *startZ* using a 4 bits (min = -365 mm, max= +365 mm) + // Max. error = 25 mm (25e-3/(0.75c) = 0.111 ns * 2 reflect = 0.222 ns) + // RMS error = 14 mm (14e-3/(0.75c) = 0.062 ns * 2 reflect = 0.124 ns) + // Also the 1-bit *istartRflag* is packed into this variable. + // + // Note: + // In the digi code we need to allow for something like 22.5 mm outside straw. + // Also because we have short straws, + // short straws are about < +-180 mm, long straws are about < +-350 mm + // The following compressions can give a large "out of straw" value; + // *don't* use these: (2.0), 32.0, 128.0, 256.0. + + unsigned char istartZ = (unsigned char)( (hitStart.z()+365.0*CLHEP::mm)/(730.0*CLHEP::mm)*16.0 ); + istartZ = (istartZ << 1) | istartRflag; + persCont->m_startZ.push_back( istartZ ); + + if ( idx > 0 ) { + persCont->m_nHits.push_back( idx - endHit ); + endHit = idx; + } + /* + // Validation output + std::cout.precision(15); + std::cout << "AJBTtoPstrawId " << strawId << std::endl; + std::cout << "AJBTtoPstartR " << startR << std::endl; + std::cout << "AJBTtoPstartPhi " << startPhi << std::endl; + std::cout << "AJBTtoPstartX " << hitStart.x() << std::endl; + std::cout << "AJBTtoPstartY " << hitStart.y() << std::endl; + std::cout << "AJBTtoPstartZ " << hitStart.z() << std::endl; + */ + } // end of "begin new hit string" + + ////////////////////////// + // Now for the end hits // + ////////////////////////// + + const HepGeom::Point3D<double> hitEnd(trtHit->postStepX, trtHit->postStepY, trtHit->postStepZ); // mm + const HepGeom::Point3D<double> hitLength = (hitEnd - hitStart); + + // + // Here both *kinEne* (kinetic energy of the particle causing the hit) and + // *steplength* (g4hit length) are persistified using a 15-bit "short float" + // (9 bit unsigned mantissa, 6 bit unsigned exponent). + // This stores values in the range 0.51*2^0 = 0.51 to 1.00*2^63 = 9.2e18. + // I enforce the limits 1.0 and 9.0e18; see below. + // Max relative error = 0.0010, RMS = 0.0004 + // + // Notes: + // + // - G4 gives kinEne in MeV; I sometimes see values ~ 1e-7 MeV (100 meV) [float round-off?] + // So I multiply by 1e9 and store in units of meV => range 1.0 meV to 9.0e18 meV (9000 TeV!) + // - About 1 in 10000 hits have steplength ~ 1e-7 mm [float round-off?] + // so again I multiply by 1e9 and store in units of pm => range 1.0 pm to 9.0e18 pm (9000 km) + // - The mantissa has maximum 9 bits, the exponent has maximum 6 bits, + // Note: a rare condition causes an 10-bit mantissa (mantissa=512). + // + double kinEne = trtHit->kineticEnergy * 1.0e9; // nano Mev = meV. + double steplength = hitLength.distance() * 1.0e9; // nano mm = pm. + if ( kinEne < 1.0 ) kinEne=1.0; // Keep the value + if ( steplength < 1.0 ) steplength=1.0; // well within the + if ( kinEne > 9.0e18 ) kinEne=9.0e18; // range of the + if ( steplength > 9.0e18 ) steplength=9.0e18; // short float. + const unsigned int kexponent = (unsigned int)ceil(log10(kinEne)/0.30102999566398); + const unsigned int sexponent = (unsigned int)ceil(log10(steplength)/0.30102999566398); + const unsigned int kmantissa = (unsigned int)(kinEne/pow(2.0,kexponent)*1024) - 512; + const unsigned int smantissa = (unsigned int)(steplength/pow(2.0,sexponent)*1024) - 512; + persCont->m_kinEne.push_back( (kmantissa << 6) | kexponent ); + persCont->m_steplength.push_back( (smantissa << 6) | sexponent ); + + // + // Persistify hit end radius using 1 bit (iendRflag) or 8 bits (endR). + // Note that the smallest value of R is the wire radius (0.0155 mm) + // + // R will be flagged as 2 mm if it is within 0.1 um of the straw wall => max error = 0.1 um, + // otherwise compress with 8 bits. The errors are as for startR, but can increased greatly + // after steplength preservation in PtoT. + // + const double endR = sqrt( hitEnd.x()*hitEnd.x() + hitEnd.y()*hitEnd.y() ); // mm + unsigned short iendRflag; + if ( endR > 1.9999*CLHEP::mm ) { + iendRflag=1; // persistify as a 1-bit flag + } + else { + iendRflag=0; // compress to 8 bits with a span of 2 mm + persCont->m_endR.push_back( (unsigned char)(endR/(2.0*CLHEP::mm)*256.0) ); + } + + // + // Persistify string *endPhi* using 8 bits (min=-pi, max=+pi) + // The errors are as for startPhi, but are very different after steplength + // preservation in PtoT. + // + const double endPhi = atan2( hitEnd.y(), hitEnd.x() ); // returns range -pi to +pi rad + persCont->m_endPhi.push_back( (unsigned char)( (endPhi+M_PI)/(2.0*M_PI)*256.0 ) ); + + // + // Persistify hit *meanTime* using 10 bits (min=0.,span=75 ns) + // with float overflow for meanTime >= 75ns (the tail of the distribution). + // Max. error = 0.037 ns; RMS error = 0.021 ns. + // Also the 1-bit *iendRflag* and 1-bit *idZsign* are packed into this variable. + // + unsigned short idZsign = (hitLength.z()>0.0) ? 1 : 0; // flag the sign of dZ + unsigned short imeanTime = ( meanTime < 75.0*CLHEP::ns ) ? (unsigned short)(meanTime/(75.0*CLHEP::ns)*1024.0) : 1023; + if ( imeanTime == 1023 ) persCont->m_meanTimeof.push_back( (float)meanTime ); // "overflow flag" + imeanTime = (imeanTime << 2) | (idZsign << 1) | iendRflag; + persCont->m_meanTime.push_back( imeanTime ); + + // + // Persistify hit *hitEne* (the energy deposited by the hit in keV) using a float but only for photons + // (relatively very few of these). Digitisation does not use hitEne for charged particles. + // + if ( lastId == 22 || + (int)(abs(lastId)/100000) == 41 || + (int)(abs(lastId)/10000000) == 1 + ) persCont->m_hitEne.push_back( (float)(trtHit->energyDeposit) ); // keV + + lastEnd = hitEnd; + lastT = meanTime; + ++idx; + /* + // Validation output + std::cout.precision(15); + std::cout << "AJBTtoPendR " << endR << std::endl; + std::cout << "AJBTtoPendPhi " << endPhi << std::endl; + std::cout << "AJBTtoPendX " << hitEnd.x() << std::endl; + std::cout << "AJBTtoPendY " << hitEnd.y() << std::endl; + std::cout << "AJBTtoPendZ " << hitEnd.z() << std::endl; + std::cout << "AJBTtoPmeanTime " << meanTime << std::endl; + std::cout << "AJBTtoPkinEne " << trtHit->kineticEnergy << std::endl; + std::cout << "AJBTtoPhitEne " << trtHit->energyDeposit << std::endl; + std::cout << "AJBTtoPsteplength " << hitLength.distance() << std::endl; + */ + } + + persCont->m_nBC.push_back(idx - endBC); + persCont->m_nId.push_back(idx - endId); + persCont->m_nHits.push_back( idx - endHit ); + +} // transToPers + + +// Create Transient +TRTUncompressedHitCollection* TRT_HitCollectionCnv_p3::createTransient(const TRT_HitCollection_p3* persObj, MsgStream &log) { + std::auto_ptr<TRTUncompressedHitCollection> trans(new TRTUncompressedHitCollection("DefaultCollectionName",persObj->m_nHits.size())); + persToTrans(persObj, trans.get(), log); + return(trans.release()); +} //createTransient + + +// Persistent(Disk) to Transient +void TRT_HitCollectionCnv_p3::persToTrans(const TRT_HitCollection_p3* persCont, TRTUncompressedHitCollection* transCont, MsgStream& /*log*/) +{ + + // if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "In TRT_HitCollectionCnv_p3::persToTrans()" << endreq; + + // some values are read less than once per hit, these need counters. + unsigned int meanTimeofCount=0, startRCount=0, endRCount=0, hitEneCount=0; + unsigned int idxBC=0, idxId=0, endHit=0, endBC=0, endId=0; + + // + // loop over strings - index [i] + // + + for ( unsigned int i = 0; i < persCont->m_nHits.size(); i++ ) { + + if ( persCont->m_nHits[i] ) { // at least one hit in the string + + const unsigned int startHit = endHit; + endHit += persCont->m_nHits[i]; + + // + // string strawId + // + const unsigned int i1 = persCont->m_strawId1b[i]; // 8 bits + const unsigned int i2 = persCont->m_strawId2b[i]; // 16 bits + const unsigned int strawId = i2*256+i1; // => 24 bits (0 to 16,777,215) + + // + // string startPhi + // + const unsigned int istartPhi = persCont->m_startPhi[i]; // 8 bits + const double startPhi = -M_PI + (istartPhi+0.5)*2.0*M_PI/256.0; // rad (min = -pi, max = +pi) + + // + // string startZ + // + const unsigned int istartZ = persCont->m_startZ[i] >> 1; // 4 bits + double startZ = -365.0*CLHEP::mm + (istartZ+0.5)*730.0*CLHEP::mm/16.0; // (min = -365 mm, max = +365 mm) + + // + // start Rflag + // + const unsigned int istartRflag = persCont->m_startZ[i] & 1; // 1 bit + + // + // string startR + // + double startR; + if ( istartRflag == 1 ) { + startR = 2.0*CLHEP::mm; // 1 bit + } + else { + const unsigned int istartR = persCont->m_startR[startRCount++]; // 8 bits + startR = (istartR+0.5)*2.0*CLHEP::mm/256.0; // (range 0 - 2 mm) + if ( startR < 0.0155*CLHEP::mm ) startR = 0.0155*CLHEP::mm; // The wire radius + } + + // + // string startX, startY (derived from R,Phi) + // + double startX = startR*cos(startPhi); + double startY = startR*sin(startPhi); + /* + // Validation output + std::cout.precision(15); + std::cout << "AJBPtoTstrawId " << strawId << std::endl; + std::cout << "AJBPtoTstartR " << startR << std::endl; + std::cout << "AJBPtoTstartPhi " << startPhi << std::endl; + std::cout << "AJBPtoTstartX " << startX << std::endl; + std::cout << "AJBPtoTstartY " << startY << std::endl; + std::cout << "AJBPtoTstartZ " << startZ << std::endl; + std::cout << "AJBPtoTnHits " << persCont->m_nHits[i] << std::endl; + */ + // + // loop over end hits in the string - index [j] + // + + for ( unsigned int j = startHit; 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++]; + + // + // hit meanTime + // + const unsigned int imeanTime = persCont->m_meanTime[j] >> 2; // 10 bits + double meanTime = (imeanTime+0.5)*75.0*CLHEP::ns/1024.0; // (min = 0.0 ns, max = 75.0 ns) + if ( imeanTime == 1023 ) meanTime = (double)persCont->m_meanTimeof[meanTimeofCount++]; // ns, 32-bit float overflow + + // + // dZ sign + // + const unsigned int idZsign = (persCont->m_meanTime[j] >> 1 ) & 1; // 1 bit + + // + // endR flag + // + const unsigned int iendRflag = persCont->m_meanTime[j] & 1; // 1 bit + + // + // hit energy deposited in keV (only relevant for photons) 32-bit float + // + const double hitEne = ( persCont->m_id[idxId] == 22 || + (int)(abs(persCont->m_id[idxId])/100000) == 41 || + (int)(abs(persCont->m_id[idxId])/10000000) == 1 + ) ? (double)persCont->m_hitEne[hitEneCount++] : 0.0; + + // + // hit endPhi (can be modified later during "steplength preservation") + // + const unsigned int iendPhi = persCont->m_endPhi[j]; // 8 bits + double endPhi = -M_PI + (iendPhi+0.5)*2.0*M_PI/256.0; // rad (min = -pi, max = +pi) + + // + // string endR (can be modified later during "steplength preservation") + // + double endR; + if ( iendRflag==1 ) { + endR = 2.0*CLHEP::mm; // 1 bit + } + else { + const unsigned int iendR = persCont->m_endR[endRCount++]; + endR = (iendR+0.5)*2.0*CLHEP::mm/256.0; // 8 bits + if ( endR < 0.0155*CLHEP::mm ) endR = 0.0155*CLHEP::mm; // the wire radius + } + + // + // hit endX, endY (derived from R,Phi) + // + double endX = endR*cos(endPhi); // can be modified later during "steplength preservation" + double endY = endR*sin(endPhi); // can be modified later during "steplength preservation" + + // Save the (o)riginal endX, endY values for the next hit start because + // they might get shrunk to fit the g4 steplength of the current hit. + double endXo = endX; + double endYo = endY; + + // + // g4 step length of the hit, m_steplength, and + // kinetic energy of the hit, m_kinEne, are both 15-bit short floats. + // Note: a rare condition causes a 16-bit short float (mantissa=512). + // + const int kmantissa = persCont->m_kinEne[j] >> 6; // 9 bits (expected) + const int smantissa = persCont->m_steplength[j] >> 6; + const int kexponent = persCont->m_kinEne[j] & 0x3F; // 6 bits + const int sexponent = persCont->m_steplength[j] & 0x3F; + const double kinEne = (kmantissa+512.5)/1024 * pow(2.0,kexponent) / 1.0e9; // MeV + double g4steplength = (smantissa+512.5)/1024 * pow(2.0,sexponent) / 1.0e9; // mm + if ( idZsign==0 ) g4steplength = -g4steplength; + + // + // Preserving the steplength of the hit by setting endZ or shrinking dX,dY. + // + double dX = endX-startX; + double dY = endY-startY; + double dZ; + double dXY2 = dX*dX+dY*dY; + double dL2 = g4steplength*g4steplength; + if ( dL2 > dXY2 ) { // define dZ such that steplength = g4steplength + dZ = sqrt(dL2-dXY2); + if (g4steplength<0.0) dZ=-dZ; + } + else { // dL2 < dXY2 // shrink dX,dY such that dXY = g4steplength + dX = dX * sqrt(dL2/dXY2); // this includes the cases where dL2=0! + dY = dY * sqrt(dL2/dXY2); + dZ = 0.0*CLHEP::mm; + endX = startX + dX; + endY = startY + dY; + //endR = sqrt( endX*endX + endY*endY ); // for validation information + //endPhi = atan2(endY,endX); // for validation information + } + double endZ = startZ + dZ; + //dX = endX-startX; // for validation information + //dY = endY-startY; // for validation information + /* + // Validation output + std::cout.precision(15); + std::cout << "AJBPtoTendR " << endR << std::endl; + std::cout << "AJBPtoTendPhi " << endPhi << std::endl; + std::cout << "AJBPtoTendX " << endX << std::endl; + std::cout << "AJBPtoTendY " << endY << std::endl; + std::cout << "AJBPtoTendZ " << endZ << std::endl; + std::cout << "AJBPtoTmeanTime " << meanTime << std::endl; + std::cout << "AJBPtoTkinEne " << kinEne << std::endl; + std::cout << "AJBPtoThitEne " << hitEne << std::endl; + std::cout << "AJBPtoTsteplength " << sqrt(dX*dX+dY*dY+dZ*dZ) << std::endl; + */ + // + // Notes: + // - All units are CLHEP, except hitEne which is in keV. + // - For charged particles kinEne is *zero*! + // + + transCont->Emplace( strawId, persCont->m_barcode[idxBC], persCont->m_id[idxId], + kinEne, hitEne, startX, startY, startZ, + endX, endY, endZ, meanTime ); + // + // End of this hit becomes the start of the next; + // use the original (uncorrected) values for X,Y + // but the derived value for Z. + // + startX = endXo; startY = endYo; startZ = endZ; + + } + } // nhits>0 + } // straw loop +} // persToTrans