Changed version

CMakeLists.txt

@@ -11,7 +11,7 @@ project(uLib)
 
 # The version number.
 set(PROJECT_VERSION_MAJOR 0)
-set(PROJECT_VERSION_MINOR 4)
+set(PROJECT_VERSION_MINOR 5)
 set(PROJECT_VERSION "${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}")
 set(PROJECT_SOVERSION "${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}")
 

src/Math/VoxRaytracer.cpp

@@ -25,11 +25,12 @@
 
 
 #include <iostream>
-#include <unordered_map>
 
 #include "VoxRaytracer.h"
 #include "Utils.h"
 
+#define unlikely(expr) __builtin_expect(!!(expr), 0)
+
 inline float fast_sign(float f) { return 1 - 2 * (f < 0); }
 
 namespace uLib {
@@ -48,17 +49,16 @@ void VoxRaytracer::RayData::AddElement(Id_t id, float L)
 
 void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in)
 {
-    if (!in.m_Data.size())
+    if (unlikely(!in.m_Data.size())) {
-    {
         std::cout << "Warinig: PoCA on exit border!\n";
+        return;
     }
-    else if (!m_Data.size())
+    else if (unlikely(!m_Data.size())) {
-    {
         m_Data = in.m_Data;
         std::cout << "Warinig: PoCA on entrance border!\n";
+        return;
     }
-    else
+    else {
-    {
         // Opzione 1) un voxel in piu' //
         m_Data.reserve(m_Data.size() + in.m_Data.size());
         m_Data.insert(m_Data.end(), in.m_Data.begin(), in.m_Data.end());
@@ -82,22 +82,16 @@ void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in)
 void VoxRaytracer::RayData::PrintSelf(std::ostream &o)
 {
     o << "Ray: total lenght " << m_TotalLength << "\n";
-    for(std::vector<Element>::iterator it = m_Data.begin(); it < m_Data.end(); ++it)
+    std::vector<Element>::iterator it;
-    {
+    for(it = m_Data.begin(); it < m_Data.end(); ++it)
         o << "[ " << (*it).vox_id << ", " << (*it).L << "] \n";
-    }
 }
 
+
 ////////////////////////////////////////////////////////////////////////////////
 //// RAY TRACER ////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 
-VoxRaytracer::VoxRaytracer(StructuredGrid &image) : m_Image(&image)
-{
-    m_scale << (m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
-               (m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
-               (m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
-}
 
 bool VoxRaytracer::GetEntryPoint(const HLine3f &line, HPoint3f &pt)
 {
@@ -156,7 +150,9 @@ bool VoxRaytracer::GetExitPoint(const HLine3f &line, HPoint3f &pt)
 }
 
 
-VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const
+VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in,
+                                                       const HPoint3f &out)
+const
 {
     RayData ray;
     Vector4f pt1 = m_Image->GetLocalPoint(in);
@@ -166,138 +162,53 @@ VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in, const
     float l = s.head(3).norm();
     Vector3f L(l/s(0), l/s(1), l/s(2));
 
+    // Vector3f scale; // FIXXX
+    // scale << (m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
+    //          (m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
+    //          (m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
+
     Vector3f offset;
-    for(int i=0;i<3;++i)
+    for(int i=0;i<3;++i) offset(i) = (s(i)>=0) - (pt1(i)-floor(pt1(i))) ;
-    {
-        offset(i) = (s(i)>=0) - (pt1(i)-floor(pt1(i)));
-    }
     offset = offset.cwiseProduct(L).cwiseAbs();
     L = L.cwiseAbs();
 
     //---- Check if the ray only crosses one voxel
     Vector3i vid = m_Image->Find(in);
-    if (vid == m_Image->Find(out))
+    if(vid == m_Image->Find(out)){
-    {
+      ray.AddElement(m_Image->Map(vid),s.norm());
-        ray.AddElement(m_Image->Map(vid),s.norm());
+      return ray;
-        return ray;
     }
 
     //---- Otherwise, loop until ray is finished
     int id; float d;
-    while (l>0)
+    while(l>0){
-    {
-        d = offset.minCoeff(&id);
 
-        if (m_Image->IsInsideGrid(vid))
+      d = offset.minCoeff(&id);
-        {
-            ray.AddElement(m_Image->Map(vid), d * m_scale(id));
-        }
 
-        vid(id) += (int)fast_sign(s(id));
+      if(m_Image->IsInsideGrid(vid)){
+	ray.AddElement(m_Image->Map(vid), d * m_scale(id) );
+      }
 
-        l -= d;
+      // nan check //
-        offset.array() -= d;
+      //        if(unlikely(!isFinite(d * scale(id)))) {
-        offset(id) = fmin(L(id),l);
+      //            std:: cout << "NAN in raytracer\n";
+      //            exit(1);
+      //        }
+
+      vid(id) += (int)fast_sign(s(id));
+
+      l -= d;
+      offset.array() -= d;
+      offset(id) = fmin(L(id),l);
     }
     return ray;
 }
 
-static int encode_v(Vector3i in)
-{
-    return ((in[0] + 1) << 4) + ((in[1] + 1) << 2) + in[2] + 1;
-}
-
-static Vector3i decode_v(int in)
-{
-    Vector3i result {
-     ((in & 48) >> 4) - 1,
-     ((in & 12) >> 2) - 1,
-     (in & 3) - 1
-    };
-    return result;
-}
-
-VoxRaytracer::RayData VoxRaytracer::BeamBetweenPoints(const HPoint3f &in, const HPoint3f &out, Vector3i thickness) const
-{
-    if (thickness[0] < 0) thickness[0] = 0;
-    if (thickness[1] < 0) thickness[1] = 0;
-    if (thickness[2] < 0) thickness[2] = 0;
-
-    Vector3i zero_v { 0, 0, 0 };
-
-    RayData ray = TraceBetweenPoints(in, out);
-    if (thickness == zero_v || ray.Data().size() == 0) return ray;
-
-    /*
-     * Calculate the forbidden relocations
-     */
-
-    std::unordered_map<int, int> ban_points(26);
-
-    Vector3i prevPos = m_Image->UnMap(ray.Data()[0].vox_id);
-    Vector3i currDir = zero_v;
-    int currLen = 1;
-
-    for (int k = 1; k < ray.Data().size(); k++)
-    {
-        Vector3i currPos = m_Image->UnMap(ray.Data()[k].vox_id);
-        Vector3i offset = currPos - prevPos;
-        prevPos = currPos;
-
-        if (k == 1) currDir = offset;
-
-        if (offset == currDir)
-        {
-            currLen++;
-            continue;
-        }
-
-        int enc_v = encode_v(currDir);
-        if (ban_points.find(enc_v) == ban_points.end())
-        {
-            ban_points.emplace(enc_v, currLen);
-        }
-        else if (currLen > ban_points[enc_v])
-        {
-            ban_points[enc_v] = currLen;
-        }
-
-        currDir = offset;
-        currLen = 2;
-    }
-
-    /*
-     * Calculate the beam section
-     */
-    std::vector<Vector3i> relocs;
-    relocs.push_back(zero_v);
-
-    /*
-     * Compose the beam
-     */
-    RayData beam;
-    for (auto iter : ray.Data())
-    {
-        Vector3i rPos = m_Image->UnMap(iter.vox_id);
-
-        for (Vector3i reloc : relocs)
-        {
-            Vector3i cPos = rPos + reloc;
-            if (!m_Image->IsInsideGrid(cPos)) continue;
-
-            beam.AddElement(m_Image->Map(cPos), iter.L);
-        }
-    }
-
-    return beam;
-}
-
-
-
 // 20150528 SV for absorbed muons
 VoxRaytracer::RayData VoxRaytracer::TraceLine(const HLine3f &line) const
 {
     RayData ray;
+
     Vector4f pt = m_Image->GetLocalPoint(line.origin);
     Vector4f s = m_Image->GetLocalPoint(line.direction);
 
@@ -339,4 +250,4 @@ VoxRaytracer::RayData VoxRaytracer::TraceLine(const HLine3f &line) const
     return ray;
 }
 
-} //end of namespace uLib
+}

src/Math/VoxRaytracer.h

@@ -38,8 +38,7 @@ namespace uLib {
 class VoxRaytracer {
 
 public:
-    class RayData
+    class RayData {
-    {
     public:
         RayData() : m_TotalLength(0) {}
 
@@ -63,8 +62,14 @@ public:
         Scalarf              m_TotalLength;
     };
 
-public:
+
-    VoxRaytracer(StructuredGrid &image);
+ public:
+ VoxRaytracer(StructuredGrid &image) : m_Image(&image) {
+      m_scale <<
+	(m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
+	(m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
+	(m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
+    }
 
     bool GetEntryPoint(const HLine3f &line, HPoint3f &pt);
 
@@ -72,8 +77,6 @@ public:
 
     RayData TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const;
 
-    RayData BeamBetweenPoints(const HPoint3f &in, const HPoint3f &out, Vector3i thickness) const;
-
     RayData TraceLine(const HLine3f &line) const;
 
     inline StructuredGrid* GetImage() const { return this->m_Image; }

src/Root/CMakeLists.txt

@@ -1,12 +1,40 @@
 set(HEADERS RootMathDense.h
             RootMuonScatter.h
-            RootHitRaw.h)
+            RootHitRaw.h
+            muCastorMCTrack.h
+            muCastorHit.h
+            muCastorInfo.h
+            muCastorSkinHit.h)
 
-set(SOURCES ${HEADERS} RootMuonScatter.cpp)
+set(SOURCES ${HEADERS} RootMuonScatter.cpp
+                       muCastorMCTrack.cpp
+                       muCastorHit.cpp
+                       muCastorInfo.cpp
+                       muCastorSkinHit.cpp)
+
+set(DICTIONARY_HEADERS muCastorMCTrack.h
+                       muCastorHit.h
+                       muCastorInfo.h
+                       muCastorSkinHit.h)
 
 set(LIBRARIES ${ROOT_LIBRARIES}
               ${PACKAGE_LIBPREFIX}Math)
 
+set(rDictName ${PACKAGE_LIBPREFIX}RootDict)
+root_generate_dictionary(${rDictName} ${DICTIONARY_HEADERS}
+                         LINKDEF Linkdef.h)
+set_source_files_properties(${rDictName}.cxx
+                            PROPERTIES GENERATED TRUE)
+set_source_files_properties(${rDictName}.h
+                            PROPERTIES GENERATED TRUE)
+list(APPEND SOURCES ${rDictName}.cxx)
+
+# TODO use a custom target linked to root_generate_dictionary
+set(R_ARTIFACTS ${CMAKE_CURRENT_BINARY_DIR}/lib${rDictName}_rdict.pcm
+                ${CMAKE_CURRENT_BINARY_DIR}/lib${rDictName}.rootmap)
+install(FILES ${R_ARTIFACTS}
+        DESTINATION ${PACKAGE_INSTALL_LIB_DIR})
+
 set(libname ${PACKAGE_LIBPREFIX}Root)
 set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)
 set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Root PARENT_SCOPE)

src/Root/Linkdef.h

@@ -79,6 +79,7 @@ using namespace ROOT::Mutom;
 #pragma link C++ class  muCastorMCTrack+;
 #pragma link C++ class  muCastorHit+;
 #pragma link C++ class  muCastorInfo+;
+#pragma link C++ class  muCastorSkinHit+;
 
 #endif // __CINT__
 

src/Root/muCastorSkinHit.cpp

@@ -0,0 +1,43 @@
+//----------------------------------------------------------
+// Class : CastorSkinHit 
+// Date: October 2020
+// Author: Germano Bonomi   germano.bonomi@unibs.it
+//----------------------------------------------------------
+
+#include <iostream>
+
+#include "muCastorSkinHit.h"
+
+/// \cond CLASSIMP
+ClassImp(muCastorSkinHit)
+/// \endcond
+
+using namespace std;
+
+//_____________________________________________________________________________
+muCastorSkinHit::muCastorSkinHit() :
+    fDetID(-1),
+    fPdgCode(-1),
+    fMotherID(-1),
+    fMomX(0.),
+    fMomY(0.),
+    fMomZ(0.),
+    fPosX(0.),
+    fPosY(0.),
+    fPosZ(0.)
+{}
+
+//_____________________________________________________________________________
+muCastorSkinHit::~muCastorSkinHit() 
+{}
+//_____________________________________________________________________________
+void muCastorSkinHit::Print(const Option_t* /*opt*/) const
+{
+    cout << "  detID: " << fDetID 
+         << "  position (cm): (" 
+         << fPosX << ", " << fPosY << ", " << fPosZ << ")"
+         << "  momentum (MeV/c): (" 
+         << fMomX << ", " << fMomY << ", " << fMomZ << ")"
+         << endl;
+}
+

src/Root/muCastorSkinHit.h

@@ -0,0 +1,46 @@
+//----------------------------------------------------------
+// Class : CastorSkinHit 
+// Date: October 2020
+// Author: Germano Bonomi   germano.bonomi@unibs.it
+//----------------------------------------------------------
+
+#ifndef muCastor_SKINHIT_H
+#define muCastor_SKINHIT_H
+
+
+#include <TObject.h>
+#include <TVector3.h>
+
+class muCastorSkinHit : public TObject
+{
+public:
+    muCastorSkinHit();
+    virtual ~muCastorSkinHit();
+
+    // methods
+    virtual void Print(const Option_t* option = "") const;
+
+    // set methods
+    void SetDetID(Int_t id)         { fDetID = id; };  
+    void SetPdgCode(Int_t pdg)      { fPdgCode = pdg; };
+    void SetMotherID(Int_t mid)     { fMotherID = mid; };
+    void SetMom(TVector3 xyz)       { fMomX = xyz.X(); fMomY = xyz.Y(); fMomZ = xyz.Z(); };
+    void SetPos(TVector3 xyz)       { fPosX = xyz.X(); fPosY = xyz.Y(); fPosZ = xyz.Z(); };
+      
+private:
+    Int_t      fDetID;      // Detector module ID 
+    Int_t      fPdgCode;    // Particle PDG Code
+    Int_t      fMotherID;   // Particle mother ID (-1 = primary, 0 = secondary, etc..)
+    Double_t   fMomX;       // Track momentum when releasing the hit (X)
+    Double_t   fMomY;       // Track momentum when releasing the hit (Y)
+    Double_t   fMomZ;       // Track momentum when releasing the hit (Z)
+    Double_t   fPosX;       // Hit coordinates (at the entrance of the detector) (X)
+    Double_t   fPosY;       // Hit coordinates (at the entrance of the detector) (Y)
+    Double_t   fPosZ;       // Hit coordinates (at the entrance of the detector) (Z)
+    
+    ClassDef(muCastorSkinHit,1) //muCastorSkinHit  
+};
+
+#endif //muCastort_SKINHIT_H
+
+