20 update sandgeomanager for the ecal endcap modules

CMakeLists.txt

@@ -7,7 +7,7 @@
 cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
 project(SandReco LANGUAGES CXX VERSION 01.00.00)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_PLATFORM_INDEPENDENT_CODE ON)
@@ -37,9 +37,9 @@ target_link_libraries(Struct PUBLIC)
 ROOT_GENERATE_DICTIONARY(StructDict struct.h MODULE Struct LINKDEF include/StructLinkDef.h)
 
 # Create SANDGeoManager lib
-add_library(SANDGeoManager SHARED src/SANDGeoManager.cpp src/SANDSTTTubeInfo.cpp src/SANDECALCellInfo.cpp)
+add_library(SANDGeoManager SHARED src/SANDGeoManager.cpp src/SANDSTTTubeInfo.cpp src/SANDECALCellInfo.cpp src/SANDENDCAPModInfo.cpp)
 target_link_libraries(SANDGeoManager PUBLIC EDepSim::edepsim_io)
-ROOT_GENERATE_DICTIONARY(SANDGeoManagerDict SANDGeoManager.h SANDSTTTubeInfo.h SANDECALCellInfo.h MODULE SANDGeoManager LINKDEF include/SANDGeoManagerLinkDef.h)
+ROOT_GENERATE_DICTIONARY(SANDGeoManagerDict SANDGeoManager.h SANDSTTTubeInfo.h SANDECALCellInfo.h SANDENDCAPModInfo.h MODULE SANDGeoManager LINKDEF include/SANDGeoManagerLinkDef.h)
 
 # Creates a libUtils shared library
 add_library(Utils SHARED src/utils.cpp src/transf.cpp)

include/SANDClustering.h

@@ -30,7 +30,9 @@ cluster Calc_variables(std::vector<reco_cell>);
 cluster Create_cluster(std::vector<dg_cell>);
 
 bool RepetitionCheck(std::vector<int>, int);
-bool isNeighbour(int, int);
+//bool isNeighbour(int, int);
+bool isNeighbour(const dg_cell&, const dg_cell&);
+bool isNeighbour_dg_reco(const dg_cell&, const reco_cell&);
 
 double EfromADCsingle(double adc, double f);
 double DfromTDC(double, double);

include/SANDENDCAPModInfo.h

@@ -0,0 +1,96 @@
+#include <TGeoBBox.h>
+#include <TGeoManager.h>
+#include <TGeoNode.h>
+#include <TGeoTube.h>
+#include <TMath.h>
+#include <TObject.h>
+
+#include <iostream>
+#include <map>
+
+#ifndef SANDENDCAPModInfo_H
+#define SANDENDCAPModInfo_H
+
+// class for storing geometric info of the SAND ECAL cells
+class SANDENDCAPModInfo : public TObject
+{
+ public:
+  enum class Orient { kHorizontal, kVertical };
+
+ private:
+  int id_;    // id of the cell
+  double x_;  // x position of the center of the module
+  double y_;  // y position of the center of the module
+  double z_;  // z position of the center of the module
+  double l_hor_;
+  double l_vert_;
+  double r_max_;
+  double r_min_;
+  int n_sec_;     // number of sections in the module (either 4 or 5)
+  double lmin_;   // min (outer) cell path length
+  double lmax_;   // max (inner) cell path length
+  double width_;  // width of the module (locally along the x direction)
+  double mod_dz_; // total depth of the module (including Al thickness)
+  double Al_dz_;  // half_thickness of the Al plate
+  TString path_;
+  TGeoNode* mod_node_;       // corresponding node of the module
+  TGeoHMatrix mod_hmatrix_;  // h_matrix of the module
+  // std::map<int, TGeoNode> cell_sections_;  // map of the sections within the
+  //                                          // module cell (indicized by the
+  //                                          // section id) (could be something
+  //                                          // else apart from TGeoNode)
+
+  // double length_;       // length of the cell?
+  // Orient orientation_;  // orientation of the cell?
+  // private Al_thick setter from mod_node_
+  void get_Al_dz();
+  void compute_min_max_l();
+
+ public:
+  SANDENDCAPModInfo();  // Default constructor
+  // SANDENDCAPModInfo(int id, double x, double y, double z, double length,
+  //                   Orient orientation);  // parametric constructor
+  SANDENDCAPModInfo(int arg_id, TGeoNode* arg_mod_node,
+                    const TGeoHMatrix& arg_hmatrix);
+  // Setter methods for the attributes
+  void id(int arg_id);
+  void x(double arg_x);
+  void y(double arg_y);
+  void z(double arg_z);
+  void width(double arg_width);
+  void Al_dz(double arg_Al_thick);
+  // void orientation(Orient arg_orientation);
+  // Getter methods for the attributes
+  int id() const;
+  int n_sections() const;
+  double x() const;
+  double y() const;
+  double z() const;
+  double width() const;
+  double mod_dz() const;
+  double l_hor() const;
+  double l_vert() const;
+  double rmin() const;
+  double rmax() const;
+  double Al_dz() const;
+  TString path() const;
+  TGeoNode* mod_node() const;
+  TGeoHMatrix mod_hmatrix() const;
+  // PMT pos. and path length computation
+  // void get_ecal_endcap_cell_local_id(double x, double y, double z,
+  //                                    int& cell_id) const;+
+
+  // compute the arc length along the curved sections at a given depth
+  double get_curv_arc_len(double depth) const;
+  // compute the total cell path length given the depth along the module (w.r.t. the outer layer)
+  double get_cell_tot_len(double depth) const;
+
+  ClassDef(SANDENDCAPModInfo, 1);
+};
+
+#ifdef __MAKECINT__
+#pragma link C++ class SANDENDCAPModInfo + ;
+#endif
+
+#endif
+

include/SANDGeoManager.h

@@ -1,4 +1,5 @@
 #include "SANDECALCellInfo.h"
+#include "SANDENDCAPModInfo.h"
 #include "SANDSTTTubeInfo.h"
 
 #include <TGeoManager.h>
@@ -47,16 +48,26 @@ const char* const path_endcapL_template =
 const char* const path_endcapR_template =
     "volWorld_PV_1/rockBox_lv_PV_0/volDetEnclosure_PV_0/volSAND_PV_0/"
     "MagIntVol_volume_PV_0/kloe_calo_volume_PV_0/ECAL_end_lv_PV_1";
+const char* const endcap_mod_regex_string =
+    "ECAL_ec_mod_([0-9]+)_lv_PV_([0-9]+)(/|)";
+const char* const endcap_mod_path_regex_string =
+    "ECAL_endcap_lv_PV_([0-9]+)/ECAL_ec_mod_([0-9]+)_lv_PV_([0-9]+)(/|)";
+
 const char* const barrel_module_name = "ECAL_lv_PV";
 const char* const endcap_module_name = "ECAL_end_lv_PV";
+const char* const barrel_last_passive_slab_name = "volECALPassiveSlab_208_PV";
 
+const double endcap_cell_width = 44.4;
 const int number_of_layers = 5;
 const int number_of_cells_per_barrel_layer = 12;
 const int number_of_barrel_modules = 24;
-const int number_of_cells_per_endcap_layer = 90;
+const int number_of_cells_per_endcap_layer = 6;
 
-// thickness of the layers in mm
+// thickness of the layers in mm (barrel)
 const double layer_thickness[number_of_layers] = {44., 44., 44., 44., 54.};
+// thickness of the cell layers in mm (endcap)
+const double ec_layer_thickness[number_of_layers] = {44.4, 44.4, 44.4, 44.4,
+                                                     52.4};
 
 // endcap module id
 const int endcap_module_ids[2] = {30, 40};
@@ -83,6 +94,11 @@ class SANDGeoManager : public TObject
                                              // value: info on cell)
   std::map<int, SANDSTTTubeInfo> sttmap_;    // map of stt tube (key: id, value:
                                              // info on tube)
+
+  std::map<int, SANDENDCAPModInfo> endcapmap_;  // map of the endcap modules
+                                                // (key: mod id, value: mod
+                                                // info)
+
   mutable TPRegexp stt_tube_regex_{
       sand_geometry::stt::stt_tube_regex_string};  // regular expression
                                                    // to match relevant
@@ -102,23 +118,45 @@ class SANDGeoManager : public TObject
                                                      // match relevant info
                                                      // about module from volume
                                                      // path
+  mutable TPRegexp endcap_mod_regex_{
+      sand_geometry::ecal::endcap_mod_regex_string};  // regular expression to
+                                                      // match relevant info
+                                                      // about endcap module
+                                                      // from volume path
+  mutable TPRegexp endcap_mod_path_regex_{
+      sand_geometry::ecal::endcap_mod_path_regex_string};  // regular expression
+                                                           // to match relevant
+                                                           // info about endcap
+                                                           // module from volume
+                                                           // path
   std::map<int, std::map<double, int> >
       stt_tube_tranverse_position_map_;  // map (key: plane id, value: map (key:
                                          // tube id, value: 2D position [i.e. x
                                          // = z, y = transversal coord]))
 
   // ECAL
-  std::vector<double> get_levels_z(double half_module_height) const;
+  std::vector<double> get_levels_z(double half_module_height,
+                                   const double (&layers_thickness)[5]) const;
   int encode_ecal_barrel_cell_local_id(int layer, int cell) const;
   int encode_ecal_endcap_cell_local_id(int layer, int cell) const;
+  static int encode_endcap_mod_id(int module_id, int module_replica_id,
+                                  int endcap_side_id);
+  static void decode_endcap_mod_id(int endcap_mod_global_id, int& module_id,
+                                   int& module_replica_id, int& endcap_side_id);
   std::pair<int, int> decode_ecal_barrel_cell_local_id(int id) const;
   std::pair<int, int> decode_ecal_endcap_cell_local_id(int id) const;
   std::map<int, TVector3> get_ecal_barrel_cell_center_local_position(
       const std::vector<double>& zlevels, double m, double q) const;
   std::map<int, TVector3> get_ecal_endcap_cell_center_local_position(
       const std::vector<double>& zlevels, double rmin, double rmax) const;
+  // new (alternative version)
+  std::map<int, TVector3> get_ec_cell_center_local_position(
+      const std::vector<double>& zlevels,
+      const SANDENDCAPModInfo& module) const;
+
   bool is_ecal_barrel(const TString& volume_name) const;
   bool is_ecal_endcap(const TString& volume_name) const;
+  bool is_endcap_mod(const TString& volume_name) const;
   bool check_and_process_ecal_path(TString& volume_path) const;
   void get_ecal_barrel_module_and_layer(const TString& volume_name,
                                         const TString& volume_path,
@@ -131,10 +169,30 @@ class SANDGeoManager : public TObject
   void get_ecal_barrel_cell_local_id(double x, double y, double z,
                                      const TGeoNode* const node,
                                      int& cell_local_id) const;
+  //   void get_ecal_endcap_cell_local_id(double x, double y, double z,
+  //                                      const TGeoNode* const node,
+  //                                      int& cell_local_id) const;
   void get_ecal_endcap_cell_local_id(double x, double y, double z,
-                                     const TGeoNode* const node,
+                                     const int& endcap_mod_id,
                                      int& cell_local_id) const;
+  int get_barrel_path_len(const double& hx, const double& hy, const double& hz,
+                          double& d1, double& d2) const;
+  int get_barrel_hit_pos(const double& d1, const int& global_cellID,
+                         double& reco_x, double& reco_y, double& reco_z) const;
+  int get_endcap_path_len(const double& hx, const double& hy, const double& hz,
+                          const int& endcap_mod_id, double& d1,
+                          double& d2) const;
+  int get_endcap_hit_pos(const double& d1, const int& global_cellID,
+                         const int& modID, double& reco_x, double& reco_y,
+                         double& reco_z) const;
+  double compute_cell_d1(const double& cell_l, const double& tdc_1,
+                         const double& tdc_2) const;
+
+  // mod id for the new endcap modules
+  int get_endcap_mod_id(const TString& volume_path) const;
   void set_ecal_info();
+  void set_ecal_endcap_info(const TGeoHMatrix& matrix);
+  void set_ecal_endcap_info();
 
   // STT
   bool is_stt_tube(const TString& volume_name) const;
@@ -187,6 +245,12 @@ class SANDGeoManager : public TObject
   static void decode_ecal_cell_id(int cell_global_id, int& detector_id,
                                   int& module_id, int& layer_id,
                                   int& cell_local_id);
+  int get_hit_path_len(const double& hx, const double& hy, const double& hz,
+                       const int& global_cell_id, double& d1, double& d2) const;
+
+  int get_reco_hit_pos(const int& cellID, const double& cell_l,
+                       const double& tdc_1, const double& tdc_2, double& reco_x,
+                       double& reco_y, double& reco_z) const;
 
   // STT
   static int encode_stt_tube_id(int stt_plane_global_id, int stt_tube_local_id);

include/StructLinkDef.h

@@ -5,21 +5,28 @@
 //#pragma link C++ class std::map < int, std::vector < int>> + ;
 //#pragma link C++ class std::map < int, double> + ;
 #pragma link C++ class std::vector < pe> + ;
+#pragma link C++ class std::vector < truecluster> + ;
+#pragma link C++ class std::vector < cluster_generator> + ;
 #pragma link C++ class std::vector < dg_ps> + ;
 #pragma link C++ class std::vector < dg_cell> + ;
 #pragma link C++ class std::vector < reco_cell > +;
+#pragma link C++ class std::vector < incomplete_cell> + ;
 //#pragma link C++ class std::map < std::string, std::vector < hit>> + ;
 #pragma link C++ class std::vector < dg_tube> + ;
 #pragma link C++ class std::vector < track> + ;
 #pragma link C++ class std::vector < cluster> + ;
 //#pragma link C++ class std::vector < particle> + ;
+#pragma link C++ class truecluster + ;
 #pragma link C++ class pe + ;
+#pragma link C++ class cluster_generator + ;
+#pragma link C++ class truecluster + ;
 #pragma link C++ class dg_ps + ;
 #pragma link C++ class dg_tube + ;
 #pragma link C++ class dg_cell + ;
 #pragma link C++ class reco_cell + ;
+#pragma link C++ class incomplete_cell + ;
 #pragma link C++ class cluster + ;
 #pragma link C++ class track + ;
-//#pragma link C++ class particle + ;
+// #pragma link C++ class particle + ;
 #pragma link C++ class event + ;
 #endif

include/struct.h

@@ -10,6 +10,50 @@ struct pe {
   double time;
   int h_index;
 };
+struct cluster_generator{
+  int pdg_code; 
+  int parent_pdg_code;
+  int track_id;
+  int parent_track_id; 
+  double dep_energy;
+  double initial_energy;
+  double initial_momentum; 
+  double initial_x; 
+  double initial_y;
+  double initial_z;  
+};
+
+struct truecluster{
+    int tid;
+    double x;
+    double y;
+    double z;
+    double t; 
+    double e;
+    double sx;
+    double sy;
+    double sz;
+    int ntot_cell;
+    int cell_l0;
+    int cell_l1; 
+    int cell_l2;
+    int cell_l3;
+    int cell_l4;
+    double energy_l0;
+    double energy_l1;
+    double energy_l2;
+    double energy_l3;
+    double energy_l4;
+    double lay0_maxE;
+    double lay1_maxE;
+    double lay2_maxE;
+    double lay3_maxE;
+    double lay4_maxE;
+    double asymmetry; 
+    double Eoverp;
+    bool moregens= false;
+    std::vector<cluster_generator> vec_generator;
+};
 
 struct hit {
   std::string det;
@@ -63,6 +107,21 @@ struct reco_cell {
   dg_ps ps2;
 };
 
+struct incomplete_cell {
+  int id;
+  bool isbarrel;
+  int endcap; 
+  double z; // better to use a custom Vector2D
+  double y;
+  double x;
+  double l;
+  int mod;
+  int lay;
+  double e;
+  int fired_pmt; 
+  dg_ps ps;
+};
+
 struct dg_tube {
   std::string det;
   int did;
@@ -93,8 +152,8 @@ struct cluster {
   double varx;
   double vary;
   double varz;
-  //std::vector<dg_cell> cells;
   std::vector<reco_cell> reco_cells; 
+  std::vector<incomplete_cell> incomplete_cells; 
 };
 
 struct track {