File indexing completed on 2024-09-27 07:02:57
0001
0002
0003
0004
0005
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
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
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
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
0092 if (hit.getEnergy()<Emin || hit.getTime()>tmax)
0093 continue;
0094
0095
0096 std::vector<double> Eneighbors(NEIGHBORS, 0.0);
0097
0098 double sl = hit.getDimension().x/2.;
0099 for (const auto& other_hit : *hits){
0100
0101
0102
0103 double tol=0.1;
0104
0105
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
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
0118
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
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
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
0150 auto alignment = volman.lookupDetElement(hit.getCellID()).nominal();
0151
0152 global_position = alignment.localToWorld(local_position);
0153
0154 }
0155 catch (...){
0156
0157 warning("Cannot find transformation from local to global coordinates.");
0158 global_position = local_position;
0159 }
0160
0161
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
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 }