Merge pull request #8 from brownd1978/libfix

Libfix

CMakeLists.txt

@@ -165,13 +165,13 @@ if (CMAKE_BUILD_TYPE MATCHES "^Debug$")
       # using Clang / LLVM
       link_libraries(-fprofile-instr-generate)
       add_custom_target(coverage
-        COMMAND gcovr --gcov-executable "llvm-cov gcov" -r ${CMAKE_SOURCE_DIR} ${PROJECT_BINARY_DIR} -s --exclude-directories Tests
+        COMMAND gcovr --gcov-executable "llvm-cov gcov" -r ${CMAKE_SOURCE_DIR} ${PROJECT_BINARY_DIR} -s --exclude-directories Mains
         WORKING_DIRECTORY ${PROJECT_BINARY_DIR}
       )
     elseif (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
       link_libraries(-lgcov) # gcov symbols
       add_custom_target(coverage
-        COMMAND gcovr -r ${CMAKE_SOURCE_DIR} ${PROJECT_BINARY_DIR} -s --exclude-directories Tests
+        COMMAND gcovr -r ${CMAKE_SOURCE_DIR} ${PROJECT_BINARY_DIR} -s --exclude-directories Mains
         WORKING_DIRECTORY ${PROJECT_BINARY_DIR}
       )
     endif()
@@ -182,11 +182,12 @@ endif()
 add_subdirectory(Detector)
 add_subdirectory(Spectra)
 add_subdirectory(General)
-add_subdirectory(Tests)
+add_subdirectory(Test)
+add_subdirectory(Mains)
 
-install(TARGETS Spectra Detector General
+install(TARGETS Spectra Detector General Test
     LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+  PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
 
 # install headers
 # Globbing here is fine because it does not influence build behaviour
@@ -211,14 +212,14 @@ fi
 export DEBUG_LEVEL=${CMAKE_BUILD_TYPE}
 
 # Linux:
-export LD_LIBRARY_PATH=\${PWD}/lib:\${LD_LIBRARY_PATH}
+export LD_LIBRARY_PATH=${CMAKE_CURRENT_BINARY_DIR}/lib:\${LD_LIBRARY_PATH}
 
 # MacOS:
-export DYLD_FALLBACK_LIBRARY_PATH=\${PWD}/lib:\${ROOT_LIB}:\${DYLD_FALLBACK_LIBRARY_PATH}
+export DYLD_FALLBACK_LIBRARY_PATH=${CMAKE_CURRENT_BINARY_DIR}/lib:\${ROOT_LIB}:\${DYLD_FALLBACK_LIBRARY_PATH}
 
 # Define the TrackToy package for others:
 export TRACKTOY_SOURCE_DIR=${CMAKE_SOURCE_DIR}
-export TRACKTOY_LIBRARY_DIR=\${PWD}/lib
+export TRACKTOY_LIBRARY_DIR=${CMAKE_CURRENT_BINARY_DIR}/lib
 
 ")
 

Data/Mu2eCalo.dat

@@ -0,0 +1,9 @@
+#
+#  Parameters describing the Mu2e tcalorimeter
+# disk1 information: rmin rmax zpos zhalf
+374.0 660.0 1667. 100.
+# disk2 information: rmin rmax zpos zhalf
+374.0 660.0 2367. 100.
+# crystal information: light propagation velocity (mm/ns), showermax depth, intrinsic time resolution (ns), spatial resolution (mm), minimum path to create a hit
+200.0 80.0 0.5 10.0 10
+

Data/Mu2eTracker.dat

@@ -4,6 +4,6 @@
 350.0 700.0 0.0 1500.0
 # Cell information: orientation (0=azimuthal, 1=axial), Ncells, cell radius, cell length, wall thickness, wire radius
 0 20736 2.5 800.0 0.015 0.0125
-# Hit information: drift velocity (mm/ns), propagation velocity (mm/ns), drift resolution (ns), propagation resolution (ns), min LR doca (mm) hit efficiency
-0.065 200 6.0 0.1 0.5 0.9
+# Hit information: drift velocity (mm/ns), propagation velocity (mm/ns), drift resolution (ns), propagation resolution (ns), min LR doca (mm), hit efficiency
+0.065 200 6.0 0.1 -10.0 0.95
 

Detector/CMakeLists.txt

@@ -1,9 +1,11 @@
 add_library(Detector SHARED
-   Tracker.cc
-   Target.cc
-   IPA.cc
-   EStar.cc
-)
+    Calorimeter.cc
+    Tracker.cc
+    Target.cc
+    HollowCylinder.cc
+    IPA.cc
+    EStar.cc
+    )
 
 # set top-level directory as include root
 target_include_directories(Detector PRIVATE ${PROJECT_SOURCE_DIR}/..)

Detector/Calorimeter.cc

@@ -0,0 +1,55 @@
+//
+//  Implementation of calorimeter
+//
+#include "TrackToy/Detector/Calorimeter.hh"
+#include "TrackToy/General/FileFinder.hh"
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <regex>
+#include <stdexcept>
+
+namespace TrackToy {
+  Calorimeter::Calorimeter(std::string const& calofile): tres_(-1.0), pres_(-1.0), minpath_(-1.0)
+  {
+    FileFinder filefinder;
+    std::string fullfile = filefinder.fullFile(calofile);
+    std::string line;
+    static std::string comment("#");
+    std::ifstream tgt_stream(fullfile,std::ios_base::in);
+    if(tgt_stream.fail()){
+      std::string errmsg = std::string("File doesn't exist" )+ fullfile;
+      throw std::invalid_argument(errmsg.c_str());
+    }
+    while (std::getline(tgt_stream, line)) {
+      // skip comments and blank lines
+      if (line.compare(0,1,comment) != 0 && line.size() > 0 ) {
+        // strip leading whitespace
+        line = std::regex_replace(line, std::regex("^ +"), "");
+        std::istringstream iss(line);
+        // first, disk 1 geometry
+        double rmin, rmax, zpos, zhalf, sprop;
+        if(disks_[0].rmin() < 0.0){
+          iss >> rmin >> rmax >> zpos >> zhalf;
+          disks_[0] = HollowCylinder(rmin,rmax,zpos,zhalf);
+          // disk 2 geometry
+        } else if (disks_[1].rmin() < 0.0){
+          iss >> rmin >> rmax >> zpos >> zhalf;
+          disks_[1] = HollowCylinder(rmin,rmax,zpos,zhalf);
+        } else {
+          iss >> sprop >> shmax_ >> tres_ >> pres_ >> minpath_;
+          // propagation velocity assumes crystals axis is along z
+          vprop_ = KinKal::VEC3(0.0,0.0,sprop);
+        }
+      }
+    }
+  }
+
+  void Calorimeter::print(std::ostream& os ) const {
+    std::cout << "Calorimeter  disk 1 " << disk(0) << " disk 2 " << disk(1) << std::endl;
+    std::cout << " Vprop " << vProp() << " showermax " << showerMax() << " time resolution " << timeResolution()
+    << " position resolution " << positionResolution() << " min hit path " << minPath() << std::endl;
+  }
+}
+
+

Detector/Calorimeter.hh

@@ -0,0 +1,98 @@
+//
+// twin disk calorimeter
+//
+#ifndef TrackToy_Detector_Calorimeter_hh
+#define TrackToy_Detector_Calorimeter_hh
+#include "TrackToy/Detector/HollowCylinder.hh"
+#include "TrackToy/General/TrajUtilities.hh"
+#include "KinKal/General/Vectors.hh"
+#include "KinKal/General/TimeRange.hh"
+#include "KinKal/Trajectory/Line.hh"
+#include "TRandom3.h"
+#include "Math/VectorUtil.h"
+#include "KinKal/Examples/ScintHit.hh"
+#include <array>
+#include <string>
+#include <vector>
+#include <iostream>
+namespace TrackToy {
+  class Calorimeter {
+    public:
+      Calorimeter(std::string const& tfile);
+      auto const& disk(size_t idisk) const { return disks_[idisk]; }
+      auto const& vProp() const { return vprop_; }
+      double showerMax() const { return shmax_; }
+      double timeResolution() const { return tres_; }
+      double positionResolution() const { return pres_; }
+      double minPath() const { return minpath_; }
+      void print(std::ostream& os) const;
+      template<class KTRAJ> unsigned simulateHits(KinKal::BFieldMap const& bfield,
+          KinKal::ParticleTrajectory<KTRAJ>& mctraj,
+          std::vector<std::shared_ptr<KinKal::Hit<KTRAJ>>>& hits, double tol) const;
+      template <class KTRAJ> void simulateHit(KinKal::ParticleTrajectory<KTRAJ> const& mctraj,
+          KinKal::TimeRange const& trange, std::vector<std::shared_ptr<KinKal::Hit<KTRAJ>>>& hits) const;
+    private:
+      std::array<HollowCylinder,2> disks_; // 2 calorimeter disks
+      KinKal::VEC3 vprop_; // light propagation velocity
+      double tres_, pres_; // time and (transverse) position resolution
+      double shmax_; // shower max depth
+      double minpath_; // minimum path to generate a hit
+      mutable TRandom3 tr_; // random number generator
+  };
+  template<class KTRAJ> unsigned Calorimeter::simulateHits(KinKal::BFieldMap const& bfield,
+      KinKal::ParticleTrajectory<KTRAJ>& mctraj,
+      std::vector<std::shared_ptr<KinKal::Hit<KTRAJ>>>& hits, double tol) const {
+    unsigned retval(0);
+   // extend through the first disk
+    extendZ(mctraj,bfield,disk(0).zmax(),tol);
+    double tstart = ztime(mctraj,mctraj.back().range().begin(),disk(0).zmin()-10.0);
+    double speed = mctraj.speed(tstart);
+    double tstep = 0.1*pres_/speed; // set step to fraction of transverse size
+    // find intersections with the first disk
+    std::vector<KinKal::TimeRange> tinters;
+    disk(0).intersect(mctraj,tinters,tstart,tstep);
+    if(tinters.size()== 0) {
+    // go to the second disk
+      extendZ(mctraj,bfield,disk(0).zmax(),tol);
+      tstart = ztime(mctraj,mctraj.back().range().begin(),disk(1).zmin()-10.0);
+      disk(1).intersect(mctraj,tinters,tstart,tstep);
+    }
+    if(tinters.size() > 0 && tinters.front().range()*speed > minpath_){
+      simulateHit(mctraj,tinters.front(),hits);
+      retval ++;
+    }
+    return retval;
+  }
+
+  template <class KTRAJ> void Calorimeter::simulateHit(KinKal::ParticleTrajectory<KTRAJ> const& mctraj,
+    KinKal::TimeRange const& trange, std::vector<std::shared_ptr<KinKal::Hit<KTRAJ>>>& hits) const {
+    using PTCA = KinKal::PiecewiseClosestApproach<KTRAJ,KinKal::Line>;
+    using SCINTHIT = KinKal::ScintHit<KTRAJ>;
+    using KinKal::Line;
+    using KinKal::CAHint;
+    // put the hit at showermax, or the end of the range (with a buffer)
+    double speed = mctraj.speed(trange.begin());
+    double tshmax = std::min(trange.begin() + shmax_/speed, trange.end()-0.1);
+    auto shmaxpos = mctraj.position3(tshmax);
+    auto hitpos = shmaxpos;
+    // figure out which disk
+    size_t idisk = (hitpos.Z() > disks_[0].zmax()) ? 1 : 0;
+    // smear the transverse position
+    hitpos.SetX(tr_.Gaus(hitpos.X(),pres_));
+    hitpos.SetY(tr_.Gaus(hitpos.Y(),pres_));
+    // set the z position to the sensor plane (back of the disk)
+    hitpos.SetZ(disks_[idisk].zmax());
+    double dz = disks_[idisk].zmax()-shmaxpos.Z();
+    // set the measurement time to correspond to the light propagation from showermax_, smeared by the resolution
+    double tmeas = tr_.Gaus(tshmax+dz/vprop_.Z(),tres_);
+    Line hline(hitpos,tmeas,vprop_,disks_[idisk].length());
+    // create the hit and add it; the hit has no material, since the original particle doesn't exist after showering so we don't need to simulate its interactions
+    CAHint tphint(tmeas,tmeas);
+    static double tprec(1e-8); // TPOCA precision
+    PTCA tp(mctraj,hline,tphint,tprec);
+    hits.push_back(std::make_shared<SCINTHIT>(tp, tres_*tres_, pres_*pres_,tprec));
+  }
+}
+
+#endif
+

Detector/HollowCylinder.cc

@@ -0,0 +1,7 @@
+#include "TrackToy/Detector/HollowCylinder.hh"
+namespace TrackToy {
+  std::ostream& operator <<(std::ostream& ost, HollowCylinder const& hcyl) {
+    ost << "Hollow Cylinder rmin " << hcyl.rmin() << " rmax " << hcyl.rmax() << " zmin " << hcyl.zmin() << " zmax " << hcyl.zmax();
+    return ost;
+  }
+}

Detector/HollowCylinder.hh

@@ -9,6 +9,7 @@
 #include "TrackToy/General/TrajUtilities.hh"
 #include <vector>
 #include <algorithm>
+#include <iostream>
 namespace TrackToy {
   using TimeRanges = std::vector<KinKal::TimeRange>;
   class HollowCylinder {
@@ -34,6 +35,7 @@ namespace TrackToy {
         return rho > rmin() && rho < rmax() && pos.Z() > zmin() && pos.Z() < zmax();
       }
   };
+  std::ostream& operator <<(std::ostream& ost, HollowCylinder const& hcyl);
 
   template<class PKTRAJ> void HollowCylinder::intersect(PKTRAJ const& pktraj, TimeRanges& tranges, double tstart, double tstep) const {
     tranges.clear();
@@ -91,6 +93,6 @@ namespace TrackToy {
     }
     // finish the last range
     if(inside)tranges.back().end() = std::min(ttest,pktraj.range().end());
-  }
+ }
 }
 #endif

Detector/Tracker.cc

@@ -10,15 +10,19 @@
 #include <stdexcept>
 
 namespace TrackToy {
-  Tracker::Tracker(MatEnv::MatDBInfo const& matdbinfo,std::string const& tgtfile):
+  Tracker::Tracker(MatEnv::MatDBInfo const& matdbinfo,std::string const& trkfile):
     ncells_(0), cellDensity_(-1.0), density_(-1.0), smat_(0),
     vdrift_(-1.0), vprop_(-1.0), sigt_(-1.0), sigl_(-1.0), lrdoca_(-1.0), hiteff_(-1.0)
   {
     FileFinder filefinder;
-    std::string fullfile = filefinder.fullFile(tgtfile);
+    std::string fullfile = filefinder.fullFile(trkfile);
     std::string line;
     static std::string comment("#");
-    std::ifstream tgt_stream(fullfile);
+    std::ifstream tgt_stream(fullfile,std::ios_base::in);
+    if(tgt_stream.fail()){
+      std::string errmsg = std::string("File doesn't exist" )+ fullfile;
+      throw std::invalid_argument(errmsg.c_str());
+    }
     double rmin(-1.0), rmax(-1.0), zpos(-1.0), zhalf(-1.0);
     double rcell(-1.0), lcell(-1.0), rwire(-1.0), wthick(-1.0);
     unsigned orient(0);
@@ -45,6 +49,7 @@ namespace TrackToy {
           double cmass = smat_->wallMaterial().density()*2.0*M_PI*rcell*lcell*wthick
             + smat_->gasMaterial().density()*M_PI*rcell*rcell*lcell
             + smat_->wireMaterial().density()*M_PI*rwire*rwire*lcell;
+          std::cout << "cell mass = " << cmass << std::endl;
           density_ = cmass*ncells_/cyl_.volume();
         } else if(vdrift_ < 0.0){
           // hit properties