algorithms/calorimetry/HEXPLIT.cc

0001 // Copyright (C) 2023 Sebouh Paul
0002 // SPDX-License-Identifier: LGPL-3.0-or-later
0003
0004 // References:
0005 //   https://arxiv.org/abs/2308.06939
0006
0007 #include <DD4hep/Alignments.h>
0008 #include <DD4hep/DetElement.h>
0009 #include <DD4hep/Objects.h>
0010 #include <DD4hep/VolumeManager.h>
0011 #include <Evaluator/DD4hepUnits.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <stdlib.h>
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <gsl/pointers>                            // for not_null
0019 #include <vector>
0020
0021 #include "HEXPLIT.h"
0022 #include "algorithms/calorimetry/HEXPLITConfig.h"
0023
0024 namespace eicrecon {
0025
0026 //positions where the overlapping cells are relative to a given cell (in units of hexagon side length)
0027 const std::vector<double> HEXPLIT::neighbor_offsets_x =[]() {
0028   std::vector<double> x;
0029   double rs[2] ={1.5, sqrt(3)/2.};
0030   double offsets[2]={0, M_PI/2};
0031   for (int i = 0; i<2; i++){
0032     for (int j = 0; j < 6; j += 1)
0033       x.push_back(rs[i]*cos(j*M_PI/3+offsets[i]));
0034   }
0035   return x;
0036 }();
0037
0038 const std::vector<double> HEXPLIT::neighbor_offsets_y =[]() {
0039   std::vector<double> y;
0040   double rs[2] ={1.5, sqrt(3)/2.};
0041   double offsets[2]={0, M_PI/2};
0042   for (int i = 0; i<2; i++){
0043     for (int j = 0; j < 6; j += 1)
0044       y.push_back(rs[i]*sin(j*M_PI/3+offsets[i]));
0045   }
0046   return y;
0047 }();
0048
0049 //indices of the neighboring cells which overlap to produce a given subcell
0050 const int HEXPLIT::neighbor_indices[SUBCELLS][OVERLAP]={{0, 11,10}, {1, 6, 11},{2, 7, 6}, {3,8,7}, {4,9,8}, {5,10,9},
0051                          {6, 11, 7}, {7, 6, 8}, {8, 7, 9}, {9,8,10},{10,9,11},{11,10,6}};
0052
0053 //positions of the centers of subcells
0054 const std::vector<double> HEXPLIT::subcell_offsets_x =[]() {
0055   std::vector<double> x;
0056   double rs[2] ={0.75, sqrt(3)/4.};
0057   double offsets[2]={0, M_PI/2};
0058   for (int i = 0; i<2; i++){
0059       for (int j = 0; j < 6; j += 1)
0060         x.push_back(rs[i]*cos(j*M_PI/3+offsets[i]));
0061   }
0062   return x;
0063 }();
0064
0065 const std::vector<double> HEXPLIT::subcell_offsets_y =[]() {
0066   std::vector<double> y;
0067   double rs[2] ={0.75, sqrt(3)/4.};
0068   double offsets[2]={0, M_PI/2};
0069   for (int i = 0; i<2; i++){
0070     for (int j = 0; j < 6; j += 1)
0071       y.push_back(rs[i]*sin(j*M_PI/3+offsets[i]));
0072   }
0073   return y;
0074 }();
0075
0076 void HEXPLIT::init() { }
0077
0078 void HEXPLIT::process(const HEXPLIT::Input& input,
0079                       const HEXPLIT::Output& output) const {
0080
0081   const auto [hits] = input;
0082   auto [subcellHits] = output;
0083
0084   double MIP=m_cfg.MIP/dd4hep::GeV;
0085   double Emin=m_cfg.Emin_in_MIPs*MIP;
0086   double tmax=m_cfg.tmax/dd4hep::ns;
0087
0088   auto volman = m_detector->volumeManager();
0089
0090   for(const auto& hit : *hits){
0091     //skip hits that do not pass E and t cuts
0092     if (hit.getEnergy()<Emin || hit.getTime()>tmax)
0093       continue;
0094
0095     //keep track of the energy in each neighboring cell
0096     std::vector<double> Eneighbors(NEIGHBORS, 0.0);
0097
0098     double sl = hit.getDimension().x/2.;
0099     for (const auto& other_hit : *hits){
0100       // maximum distance between where the neighboring cell is and where it should be
0101       // based on an ideal geometry using the staggered tessellation pattern.
0102       // Deviations could arise from rounding errors or from detector misalignment.
0103       double tol=0.1; // in units of side lengths.
0104
0105       //only look at hits nearby within two layers of the current layer
0106       int dz=abs(hit.getLayer()-other_hit.getLayer());
0107       if (dz>2 || dz==0)
0108         continue;
0109       if (other_hit.getEnergy()<Emin || other_hit.getTime()>tmax)
0110         continue;
0111       //difference in transverse position (in units of side lengths)
0112       double dx=(other_hit.getLocal().x-hit.getLocal().x)/sl;
0113       double dy=(other_hit.getLocal().y-hit.getLocal().y)/sl;
0114       if (abs(dx)>2 || abs(dy)>sqrt(3))
0115         continue;
0116
0117       //loop over locations of the neighboring cells
0118       //and check if the jth hit matches this location
0119       for(int k=0;k<NEIGHBORS;k++){
0120         if(abs(dx-neighbor_offsets_x[k])<tol && abs(dy-neighbor_offsets_y[k])<tol){
0121           Eneighbors[k]+=other_hit.getEnergy();
0122           break;
0123         }
0124       }
0125     }
0126     double weights[SUBCELLS];
0127     for(int k=0; k<NEIGHBORS; k++){
0128       Eneighbors[k]=std::max(Eneighbors[k],MIP);
0129     }
0130     double sum_weights=0;
0131     for(int k=0; k<SUBCELLS; k++){
0132       weights[k]=Eneighbors[neighbor_indices[k][0]]*Eneighbors[neighbor_indices[k][1]]*Eneighbors[neighbor_indices[k][2]];
0133       sum_weights+=weights[k];
0134     }
0135     for(int k=0; k<SUBCELLS;k++){
0136
0137       //create the subcell hits.  First determine their positions in local coordinates.
0138       const decltype(edm4eic::CalorimeterHitData::local) local(hit.getLocal().x+subcell_offsets_x[k]*sl, hit.getLocal().y+subcell_offsets_y[k]*sl, hit.getLocal().z);
0139
0140       //convert this to a position object so that the global position can be determined
0141       dd4hep::Position local_position;
0142       local_position.SetX(local.x*dd4hep::mm);
0143       local_position.SetY(local.y*dd4hep::mm);
0144       local_position.SetZ(local.z*dd4hep::mm);
0145
0146       dd4hep::Position global_position;
0147       try {
0148
0149         //also convert this to the detector's global coordinates.  To do: check if this is correct
0150         auto alignment = volman.lookupDetElement(hit.getCellID()).nominal();
0151
0152         global_position = alignment.localToWorld(local_position);
0153
0154       }
0155       catch (...){
0156         // do this to prevent errors when running the test on the mock detector
0157         warning("Cannot find transformation from local to global coordinates.");
0158         global_position = local_position;
0159       }
0160
0161       //convert this from position object to a vector object
0162       const decltype(edm4eic::CalorimeterHitData::position) position = {static_cast<float>(global_position.X()/dd4hep::mm), static_cast<float>(global_position.Y()/dd4hep::mm), static_cast<float>(global_position.Z()/dd4hep::mm)};
0163
0164       //bounding box dimensions depend on the orientation of the rhombus
0165       int orientation = k%3==0;
0166       const decltype(edm4eic::CalorimeterHitData::dimension) dimension(sl*(orientation?1:1.5), sl*sqrt(3)/2.*(orientation?2:1),
0167                                                                        hit.getDimension()[2]);
0168
0169       subcellHits->create(
0170             hit.getCellID(),
0171             hit.getEnergy()*weights[k]/sum_weights,
0172             0,
0173             hit.getTime(),
0174             0,
0175             position,
0176             dimension,
0177             hit.getSector(),
0178             hit.getLayer(),
0179             local);
0180     }
0181   }
0182 }
0183
0184
0185 } // namespace eicrecon