diff --git a/include/ChargedParticle.h b/include/ChargedParticle.h index e34ea1c..c48d26f 100644 --- a/include/ChargedParticle.h +++ b/include/ChargedParticle.h @@ -72,7 +72,6 @@ class ChargedParticle: public TransientParticle { // Group steps by radiator of course; in the easiest case an entry and exit points; std::vector > m_RadiatorHistory; - bool m_StopTracing; //! #ifndef DISABLE_ROOT_IO diff --git a/include/CherenkovDetector.h b/include/CherenkovDetector.h index 964124f..80d4924 100644 --- a/include/CherenkovDetector.h +++ b/include/CherenkovDetector.h @@ -117,7 +117,7 @@ class CherenkovDetector: public TObject { // readout ID -> pixel position converter (for external usage) std::function m_ReadoutIDToPosition; //! - + std::function m_ReadoutIDToNormal; //! private: TString m_Name; // This is needed for dd4hep cell index decoding; diff --git a/include/OpticalPhoton.h b/include/OpticalPhoton.h index c8c0eb6..710d358 100644 --- a/include/OpticalPhoton.h +++ b/include/OpticalPhoton.h @@ -31,8 +31,9 @@ class OpticalPhoton: public TransientParticle { inline void SetVertexParentMomentum(const TVector3 &momentum){ m_VertexParentMomentum = momentum; }; inline void SetVolumeCopy(uint64_t copy) { m_VolumeCopy = copy; }; inline void SetDetectionPosition(const TVector3 &position) { m_DetectionPosition = position; }; + //inline void SetDetectionMomentum(const TVector3 &momentum) { m_DetectionMomentum = momentum;} /////////////////////////// inline void SetPhotonDetector(CherenkovPhotonDetector *pd) { m_PhotonDetector = pd; }; - + inline void SetNormalSurface(const TVector3 &nSurface) {m_NormalSurface = nSurface;};//////////////////////////////// inline void AddReflectionPoint(ReflectionPoint *reflection) { m_ReflectionPoints.push_back(reflection); }; @@ -54,7 +55,9 @@ class OpticalPhoton: public TransientParticle { inline const TVector3 &GetVertexMomentum( void ) const { return m_VertexMomentum; }; inline const TVector3 &GetVertexParentMomentum( void ) const { return m_VertexParentMomentum; }; inline const TVector3 &GetDetectionPosition( void ) const { return m_DetectionPosition; }; - + //inline const TVector3 &GetDetectionMomentum( void ) const { return m_DetectionMomentum; };/////////////////////////// + inline const TVector3 &GetNormalSurface( void ) const { return m_NormalSurface;};//////////////////////////////// + inline bool WasDetected( void ) const { return m_Detected; }; inline uint64_t GetVolumeCopy( void ) const { return m_VolumeCopy; }; @@ -78,6 +81,8 @@ class OpticalPhoton: public TransientParticle { TRef m_PhotonDetector; uint64_t m_VolumeCopy; TVector3 m_DetectionPosition; + //TVector3 m_DetectionMomentum; + TVector3 m_NormalSurface; double m_DetectionTime; // 'true' if the photon was actually detected; QE-based and geometric efficiency accounted; @@ -92,7 +97,10 @@ class OpticalPhoton: public TransientParticle { // Average estimated phi angle; no need to know it precisely (?); std::map m_Phi; //! - + + //Incident angle + std::map m_Angle; //! + #ifndef DISABLE_ROOT_IO ClassDef(OpticalPhoton, 3); #endif diff --git a/source/ChargedParticle.cc b/source/ChargedParticle.cc index fa194b9..cfd1721 100644 --- a/source/ChargedParticle.cc +++ b/source/ChargedParticle.cc @@ -35,7 +35,7 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) auto pd = photon->GetPhotonDetector(); const auto irt = pd->GetIRT(photon->GetVolumeCopy()); - if (!irt) { + if (!irt) { printf("No photosensor with this cellID found!\n"); continue; } //if @@ -46,7 +46,17 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) if (radiator->m_Locations.size() != zdim+1) continue; TVector3 phx = photon->GetDetectionPosition(); - + // ADding Stuff + const auto norm = (photon->GetNormalSurface()); + const auto mom = (1e9*photon->GetVertexMomentum()); + photon->m_Angle[radiator]=norm.Dot(mom); + //float angle = photon->m_Angle[radiator]; + //if((angle*(180/M_PI)) > 0) + //printf("-->IRT %lf %lf %lf %lf \n",norm.X(),norm.Y(),norm.Z(),180 - (angle*(180/M_PI))); + //if((angle*(180/M_PI)) < 0) + //printf("-->IRT %lf %lf %lf %lf \n",norm.X(),norm.Y(),norm.Z(),180 + (angle*(180/M_PI))); + // up to here + // Get effective attenuation length for this radiator, as well as the // parameterization of its rear surface in this particular sector; this is // not really a clean procedure for dRICH aerogel, but should be good enough