Bug Summary

File:ecl/geometry/src/ECLGeometryPar.cc
Warning:line 187, column 7
Use of memory after it is freed

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -O3 -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name ECLGeometryPar.cc -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/data/b2soft/buildbot/development/build -fcoverage-compilation-dir=/data/b2soft/buildbot/development/build -resource-dir /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++ -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/x86_64-redhat-linux -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/backward -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/python3.12 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/CLHEP -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/Geant4 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/root -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/belle_legacy -I include/ -D _PACKAGE_="ecl" -D G4UI_USE_TCSH -D RaveDllExport= -D HAS_SQLITE -D HAS_CALLGRIND -I include -I /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/libxml2 -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++ -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/x86_64-redhat-linux -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/backward -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../x86_64-redhat-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -Wno-missing-braces -Wno-unused-command-line-argument -std=c++20 -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fno-implicit-modules -fskip-odr-check-in-gmf -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /scan_build/2026-05-31-004316-385593-1 -x c++ ecl/geometry/src/ECLGeometryPar.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9/* Own header. */
10#include <ecl/geometry/ECLGeometryPar.h>
11
12/* ECL headers. */
13#include <ecl/dataobjects/ECLElementNumbers.h>
14
15/* Basf2 headers. */
16#include <framework/logging/Logger.h>
17
18/* Geant4 headers. */
19#include <G4NavigationHistory.hh>
20#include <G4PhysicalVolumeStore.hh>
21#include <G4Point3D.hh>
22#include <G4Transform3D.hh>
23#include <G4Vector3D.hh>
24#include <G4VTouchable.hh>
25
26/* C++ headers. */
27#include <iomanip>
28
29using namespace std;
30using namespace Belle2;
31using namespace ECL;
32
33ECLGeometryPar* ECLGeometryPar::m_B4ECLGeometryParDB = 0;
34
35/** Mapping class */
36class Mapping_t {
37public:
38 /** Retrieving theta and phi id of crystal */
39 static void Mapping(int id, int& ThetaId, int& PhiId)
40 {
41 ThetaId = m_Theta[((unsigned int)id) >> 4];
42 PhiId = id - m_dTheta[ThetaId] * 16 - ThetaId * 128;
43 }
44
45 /** Retrieving theta id, phi id, reciprocal shift and index */
46 static void Mapping(int id, int& ThetaId, int& PhiId, int& nrep, int& indx)
47 {
48 Mapping(id, ThetaId, PhiId);
49
50 int off = m_offsets[ThetaId];
51 int i = m_tbl[ThetaId];
52
53 int r = m_recip[i];
54 int d = m_denom[i];
55
56 nrep = (PhiId * r) >> m_RECIPROCAL_SHIFT;
57 indx = off + (PhiId - nrep * d);
58 }
59
60 /** return cell id as a function of theta id and phi id */
61 static int CellID(int ThetaId, int PhiId)
62 {
63 return PhiId + m_dTheta[ThetaId] * 16 + ThetaId * 128;
64 }
65
66 /** return offset based on theta id */
67 static int Offset(int ThetaId)
68 {
69 return m_dTheta[ThetaId] + ThetaId * 8;
70 }
71
72 /** getter for theta */
73 static int Indx2ThetaId(int indx)
74 {
75 return m_Theta[indx];
76 }
77
78 /** getter for number of crystals */
79 static int ThetaId2NCry(int ThetaId) // Theta Id to the number of crystals @ this Id
80 {
81 return m_denom[m_tbl[ThetaId]];
82 }
83
84private:
85 static const char m_dTheta[69];/**< array of theta offsets */
86 static const unsigned char m_Theta[546]; /**< array of theta */
87 static const unsigned char m_tbl[69]; /**< array of crystals per phi sector */
88 static const unsigned char m_offsets[69]; /**< array of offsets */
89
90 static const unsigned char m_RECIPROCAL_SHIFT = 16; /**< reciprocal shift */
91 static const unsigned int m_recip[5]; /**< array of reciprocal values */
92 static const unsigned int m_denom[5]; /**< array of denominator values */
93};
94
95const unsigned char Mapping_t::m_Theta[546] = {
96 0, 0, 0, 1, 1, 1,
97 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4,
98 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
99 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10,
100 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12,
101
102 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15,
103 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18,
104 19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 21, 21, 21, 21,
105 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24,
106 25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27,
107 28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30,
108 31, 31, 31, 31, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 33, 33, 33, 33,
109 34, 34, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36,
110 37, 37, 37, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39,
111 40, 40, 40, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 42, 42,
112 43, 43, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45,
113 46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 48, 48, 48,
114 49, 49, 49, 49, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 51,
115 52, 52, 52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54,
116 55, 55, 55, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 57, 57, 57,
117 58, 58, 58, 58, 58, 58, 58, 58, 58,
118
119 59, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 60, 60, 60, 60,
120 61, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 62,
121 63, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 64,
122 65, 65, 65, 65, 65, 65,
123 66, 66, 66, 66, 67, 67, 67, 67, 68, 68, 68, 68
124};
125
126const char Mapping_t::m_dTheta[69] = {
127 0, -5, -10, -14, -18, -22, -24, -26, -28, -30, -32, -34, -33, // forward
128
129 -32, -31, -30, -29, -28, -27, -26, -25, -24,
130 -23, -22, -21, -20, -19, -18, -17, -16, -15,
131 -14, -13, -12, -11, -10, -9, -8, -7, -6,
132 -5, -4, -3, -2, -1, 0, 1, 2, 3,
133 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
134
135 14, 15, 16, 14, 12, 10, 8, 6, 2, -2 // backward
136};
137
138const unsigned char Mapping_t::m_tbl[69] = { // pointer to the denominator/reciprocal arrays or how many crystals in a phi sector
139 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, // forward
140
141 0, 0, 0, 0, 0, 0, 0, 0, 0,
142 0, 0, 0, 0, 0, 0, 0, 0, 0,
143 0, 0, 0, 0, 0, 0, 0, 0, 0,
144 0, 0, 0, 0, 0, 0, 0, 0, 0,
145 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
146
147 4, 4, 3, 3, 3, 3, 3, 2, 2, 2 // backward
148};
149
150const unsigned char Mapping_t::m_offsets[69] = {
151 0, 3, 6, 10, 14, 18, 24, 30, 36, 42, 48, 54, 63, // forward
152
153 132, 134, 136, 138, 140, 142, 144, 146, 148,
154 150, 152, 154, 156, 158, 160, 162, 164, 166,
155 168, 170, 172, 174, 176, 178, 180, 182, 184,
156 186, 188, 190, 192, 194, 196, 198, 200, 202,
157 204, 206, 208, 210, 212, 214, 216, 218, 220, 222,
158
159 72, 81, 90, 96, 102, 108, 114, 120, 124, 128 // backward
160};
161
162#define pack(x) ((1<<Mapping_t::m_RECIPROCAL_SHIFT)/x+1)
163const unsigned int Mapping_t::m_recip[] = {pack(2), pack(3), pack(4), pack(6), pack(9)};
164const unsigned int Mapping_t::m_denom[] = { (2), (3), (4), (6), (9)};
165#undef pack
166
167
168ECLGeometryPar* ECLGeometryPar::Instance()
169{
170 if (!m_B4ECLGeometryParDB) m_B4ECLGeometryParDB = new ECLGeometryPar();
171 return m_B4ECLGeometryParDB;
172}
173
174ECLGeometryPar::ECLGeometryPar()
175{
176 clear();
177 // delay crystal positions fetching to a moment when it actually needed and geometry is already built
178 // read();
179}
180
181ECLGeometryPar::~ECLGeometryPar()
182{
183 if (m_B4ECLGeometryParDB
3.1
'm_B4ECLGeometryParDB' is non-null
) {
1
Assuming 'm_B4ECLGeometryParDB' is non-null
→
2
←
Taking true branch
→
4
←
Taking true branch
→
184 delete m_B4ECLGeometryParDB;
3
←
Calling '~ECLGeometryPar'
→
5
←
Memory is released
→
185 B2DEBUG(150, "m_B4ECLGeometryParDB deleted ")do { if (Belle2::LogSystem::debugEnabled()) do { if (Belle2::
LogSystem::Instance().isLevelEnabled(Belle2::LogConfig::c_Debug
, 150, "ecl")) { { LogVariableStream varStream; varStream <<
"m_B4ECLGeometryParDB deleted "; Belle2::LogSystem::Instance
().sendMessage(Belle2::LogMessage(Belle2::LogConfig::c_Debug,
std::move(varStream), "ecl", __PRETTY_FUNCTION__, "ecl/geometry/src/ECLGeometryPar.cc"
, 185, 150)); }; } } while(false); } while(false);
6
←
Assuming the condition is false
→
7
←
Taking false branch
→
8
←
Loop condition is false. Exiting loop
→
186 }
187 if (m_ECLForwardGlobalT) delete m_ECLForwardGlobalT;
9
←
Use of memory after it is freed
188 if (m_ECLBarrelGlobalT) delete m_ECLBarrelGlobalT;
189 if (m_ECLBackwardGlobalT) delete m_ECLBackwardGlobalT;
190}
191
192void ECLGeometryPar::clear()
193{
194 m_ini_cid = -1;
195 mPar_cellID = 0;
196 mPar_thetaID = 0;
197 mPar_phiID = 0;
198}
199
200// There is no way to get world coordinates of a local point of a physical volume in Geant.
201// The only way to check geometry is to trace particle and check volumes it crosses.
202// Here particle gun parameters are produced to check crystal positions with the center @ r0 and direction n
203void ParticleGunParameters(const G4String& comment, const G4ThreeVector& n, const G4ThreeVector& r0, double dphi)
204{
205 cout << comment << endl;
206 cout << "Center position = " << r0 << ", Direction = " << n << endl;
207 // closest point to z-axis
208 double t = -(n.z() * r0.z() - n * r0) / (n.z() * n.z() - n * n);
209 G4ThreeVector r = n * t + r0;
210 cout << "Closest point to z-axis = " << r << endl; // at the moment I do not see tilt in phi
211 const double r2d = 180 / M_PI3.14159265358979323846;
212 double th = r2d * n.theta();
213 double phi = r2d * n.phi();
214 double z = r.z();
215 dphi *= r2d;
216 cout << " 'thetaParams': [" << th << ", " << th << "]," << endl;
217 cout << " 'phiParams': [" << phi << "+0*" << dphi << ", " << phi << "+0*" << dphi << "]," << endl;
218 cout << " 'zVertexParams': [" << z << ", " << z << "]" << endl;
219}
220
221G4Transform3D getT(const G4VPhysicalVolume& v)
222{
223 return G4Transform3D(v.GetObjectRotationValue(), v.GetObjectTranslation());
224}
225
226G4Transform3D getT(const G4String& n)
227{
228 G4PhysicalVolumeStore* store = G4PhysicalVolumeStore::GetInstance();
229 // int N = store->size();
230 // cout<<N<<endl;
231 // for(int i=0;i<N;i++){
232 // G4VPhysicalVolume *vv = (*store)[i];
233 // if(vv->GetName()==n){
234 // cout<<i<<" "<<vv->GetName()<<" "<<vv->GetMultiplicity()<<" "<<vv->GetCopyNo()<<endl;
235 // G4PVReplica *r = dynamic_cast<G4PVReplica*>(vv);
236 // if(r) cout<<"Instance = "<<r->GetInstanceID()<<endl;
237 // }
238 // }
239
240 G4VPhysicalVolume* v = store->GetVolume(n);
241 return getT(*v);
242}
243
244void print_trans(const G4Transform3D& t)
245{
246 for (int i = 0; i < 4; i++) {
247 for (int j = 0; j < 3; j++)
248 cout << t(j, i) << " ";
249 cout << endl;
250 }
251}
252
253void ECLGeometryPar::read()
254{
255 m_crystals.clear();
256 m_crystals.reserve(46 * 2 + 132); // 10752 bytes
257
258 // Endcap sectors are rotated since behaviour of G4Replica class. It
259 // requires a physical volume during phi replication to be symmetric
260 // around phi=0. So we need to rotate it by -((2*pi)/16)/2 angle at the
261 // geometry description are return back here.
262
263 G4Point3D p0(0, 0, 0); G4Vector3D v0(0, 0, 1);
264
265 {
266 m_ECLForwardGlobalT = new G4Transform3D(getT("ECLForwardPhysical"));
267 G4Transform3D T = getT("ECLForwardCrystalSectorPhysical_1");
268 G4String tnamef("ECLForwardWrappedCrystal_Physical_");
269 for (int i = 0; i < 72; i++) {
270 G4String vname(tnamef); vname += to_string(i);
271 G4Transform3D cT = T * getT(vname);
272 CrystalGeom_t c = {cT * p0, cT * v0};
273 m_crystals.push_back(c);
274 // G4cout << i << " " << c.pos << " " << c.dir << G4endl;
275 }
276 }
277
278 {
279 m_ECLBackwardGlobalT = new G4Transform3D(getT("ECLBackwardPhysical"));
280 G4Transform3D T = getT("ECLBackwardCrystalSectorPhysical_0");
281 G4String tnamef("ECLBackwardWrappedCrystal_Physical_");
282 for (int i = 0; i < 60; i++) {
283 G4String vname(tnamef); vname += to_string(i);
284 G4Transform3D cT = T * getT(vname);
285 CrystalGeom_t c = {cT * p0, cT * v0};
286 m_crystals.push_back(c);
287 // G4cout << i << " " << c.pos << " " << c.dir << G4endl;
288 }
289 }
290
291 {
292 // get barrel sector (between two septums) transformation
293
294 // extract global barrel movements assuming we know transformation of the sector in axial symmetric situation
295 m_ECLBarrelGlobalT = new G4Transform3D(getT("ECLBarrelSectorPhysical_0") * G4RotateZ3D(M_PI3.14159265358979323846 / 2));
296 G4Transform3D Ts = G4RotateZ3D(-M_PI3.14159265358979323846 / 2);
297 // since barrel sector is symmetric around phi=0 we need to
298 // translate crystal with negative phi back to positive rotating
299 // crystal position by (2*M_PI/72) angle
300 G4Transform3D Ts1 = G4RotateZ3D(M_PI3.14159265358979323846 / 36) * Ts;
301
302 G4String tname("ECLBarrelWrappedCrystal_Physical_");
303 for (int i = 0; i < 2 * 46; i++) {
304 G4String vname(tname); vname += to_string(i);
305 G4Transform3D cT = ((i % 2) ? Ts1 : Ts) * getT(vname);
306 CrystalGeom_t c = {cT * p0, cT * v0};
307 m_crystals.push_back(c);
308 // G4cout << i << " " << c.pos << " " <<c.dir<<G4endl;
309 }
310 }
311
312 B2DEBUG(150, "ECLGeometryPar::read() initialized with " << m_crystals.size() << " crystals.")do { if (Belle2::LogSystem::debugEnabled()) do { if (Belle2::
LogSystem::Instance().isLevelEnabled(Belle2::LogConfig::c_Debug
, 150, "ecl")) { { LogVariableStream varStream; varStream <<
"ECLGeometryPar::read() initialized with " << m_crystals
.size() << " crystals."; Belle2::LogSystem::Instance().
sendMessage(Belle2::LogMessage(Belle2::LogConfig::c_Debug, std
::move(varStream), "ecl", __PRETTY_FUNCTION__, "ecl/geometry/src/ECLGeometryPar.cc"
, 312, 150)); }; } } while(false); } while(false);
313}
314
315template <int n>
316void sincos(const double* ss, int iphi, double& s, double& c)
317{
318 int n4 = n / 4;
319 int iq = iphi / n4; // quadrant
320 int is = iphi % n4;
321 int ic = n4 - is;
322 if (iq & 1) {int t = is; is = ic; ic = t;}
323 double ls = ss[is];
324 double lc = ss[ic];
325 // check quadrant and assign sign bit accordingly
326 if ((((0 << 0) + (0 << 1) + (1 << 2) + (1 << 3)) >> iq) & 1) ls = -ls;
327 if ((((0 << 0) + (1 << 1) + (1 << 2) + (0 << 3)) >> iq) & 1) lc = -lc;
328 s = ls;
329 c = lc;
330}
331
332const double ss72[] = {
333 0.0000000000000000000000,
334 0.0871557427476581735581,
335 0.1736481776669303488517,
336 0.2588190451025207623489,
337 0.3420201433256687330441,
338 0.4226182617406994361870,
339 0.5000000000000000000000,
340 0.5735764363510460961080,
341 0.6427876096865393263226,
342 0.7071067811865475244008,
343 0.7660444431189780352024,
344 0.8191520442889917896845,
345 0.8660254037844386467637,
346 0.9063077870366499632426,
347 0.9396926207859083840541,
348 0.9659258262890682867497,
349 0.9848077530122080593667,
350 0.9961946980917455322950,
351 1.0000000000000000000000
352};
353
354const double ss16[] = {
355 0.0000000000000000000000,
356 0.3826834323650897717285,
357 0.7071067811865475244008,
358 0.9238795325112867561282,
359 1.0000000000000000000000
360};
361
362void ECLGeometryPar::InitCrystal(int cid)
363{
364 if (m_crystals.size() == 0) read();
365 int thetaid, phiid, nreplica, indx;
366 Mapping_t::Mapping(cid, thetaid, phiid, nreplica, indx);
367 // cout<<cid<<" "<<thetaid<<" "<<phiid<<" "<<nreplica<<" "<<indx<<endl;
368 const CrystalGeom_t& t = m_crystals[indx];
369 double s, c;
370 if (indx > 131)
371 sincos<72>(ss72, nreplica, s, c);
372 else
373 sincos<16>(ss16, nreplica, s, c);
374
375 G4Transform3D* T;
376 if (ECLElementNumbers::isForward(cid + 1)) {
377 T = m_ECLForwardGlobalT;
378 } else if (ECLElementNumbers::isBarrel(cid + 1)) {
379 T = m_ECLBarrelGlobalT;
380 } else {
381 T = m_ECLBackwardGlobalT;
382 c = -c;
383 }
384 double xp = c * t.pos.x() - s * t.pos.y();
385 double yp = s * t.pos.x() + c * t.pos.y();
386 // m_current_crystal.pos.set(xp, yp, t.pos.z());
387
388 double xv = c * t.dir.x() - s * t.dir.y();
389 double yv = s * t.dir.x() + c * t.dir.y();
390 // m_current_crystal.dir.set(xv, yv, t.dir.z());
391
392 m_current_crystal.pos = (1 / CLHEP::cm) * (*T * G4Point3D(xp, yp, t.pos.z()));
393 m_current_crystal.dir = *T * G4Vector3D(xv, yv, t.dir.z());
394
395 // cout<<t.pos<<" "<<t.dir<<" "<<m_current_crystal.pos<<" "<<m_current_crystal.dir<<endl;
396 m_ini_cid = cid;
397}
398
399int ECLGeometryPar::GetCellID(int ThetaId, int PhiId)
400{
401 return Mapping_t::CellID(ThetaId, PhiId);
402}
403
404void ECLGeometryPar::Mapping(int cid)
405{
406 mPar_cellID = cid;
407 Mapping_t::Mapping(mPar_cellID, mPar_thetaID, mPar_phiID);
408}
409
410int ECLGeometryPar::TouchableDiodeToCellID(const G4VTouchable* touch)
411{
412 return TouchableToCellID(touch);
413}
414
415int ECLGeometryPar::TouchableToCellID(const G4VTouchable* touch)
416{
417 // touch->GetCopyNumber() is a virtual call, avoid it by using
418 // directly G4NavigationHistory so we will have only one virtual
419 // call instead of three here
420 const G4NavigationHistory* h = touch->GetHistory();
421 int hd = h->GetDepth();
422 int i1 = h->GetReplicaNo(hd - 1); // index of each volume is set at physical volume creation
423 int i2 = h->GetReplicaNo(hd - 2); // go up in volume hierarchy
424
425 int ThetaId = Mapping_t::Indx2ThetaId(i1);
426 int NCryst = Mapping_t::ThetaId2NCry(ThetaId); // the number of crystals in a sector at given ThetaId
427 int Offset = Mapping_t::Offset(ThetaId);
428
429 int ik = i1 - Offset;
430 int PhiId;
431 if (NCryst == 2) {
432 PhiId = (i2 - (ik % 2)) * NCryst + ik;
433 if (PhiId < 0) PhiId += 144;
434 } else {
435 int i3 = h->GetReplicaNo(hd - 3); // go up in volume hierarchy
436 if (ThetaId < 13)
437 PhiId = (i2 + 2 * i3) * NCryst + ik;
438 else {
439 // int tt[] = {3,2,1,0,7,6,5,4};
440 // if(i3r!=tt[i3]){
441 // cout<<i3<<" "<<i3<<" "<<tt[i3]<<endl;
442 // exit(0);
443 // }
444 int i3r = (3 - (i3 % 4)) + 4 * (i3 / 4);
445 PhiId = (2 * i3r + (1 - i2)) * NCryst + ik;
446 }
447 }
448
449 int cellID = Offset * 16 + PhiId;
450 // cout<<"ECLGeometryPar::TouchableToCellID "<<h->GetVolume(hd-1)->GetName()<<" "<<i1<<" "<<i2<<" "<<h->GetReplicaNo(hd - 3)<<" "<<ThetaId<<" "<<NCryst<<" "<<Offset<<" "<<PhiId<<endl;
451
452 // test of the position and direction of crystal
453 if (0) {
454 G4AffineTransform t = h->GetTopTransform();
455 G4ThreeVector o(0, 0, 0), n(0, 0, 1);
456 G4ThreeVector ro = t.Inverse().TransformPoint(o);
457 G4ThreeVector rn = t.Inverse().TransformAxis(n);
458
459 InitCrystal(cellID);
460 ro *= 1 / CLHEP::cm;
461
462 G4ThreeVector dr = m_current_crystal.pos - ro, dn = m_current_crystal.dir - rn;
463 if (dr.mag() > 1e-10 || dn.mag() > 1e-10) {
464 cout << "Mismatch " << h->GetVolume(hd - 1)->GetName() << " " << cellID << " " << ThetaId << " " << PhiId << " " << NCryst << " "
465 << hd << " " << i2 << " " << i1 << " " << m_current_crystal.pos << " " << ro << " " << rn << " " << dr << " " << dn << endl;
466
467 for (int i = 0; i < 144; i++) {
468 int ci = Mapping_t::CellID(ThetaId, i);
469 InitCrystal(ci);
470 dr = m_current_crystal.pos - ro;
471 if (dr.mag() < 1e-10) cout << "best PhiId = " << i << endl;
472 }
473 }
474 }
475 return cellID;
476}
477
478int ECLGeometryPar::ECLVolumeToCellID(const G4VTouchable* touch)
479{
480 int depth = touch->GetHistoryDepth();
481 if ((depth != 3) && (depth != 5)) {
482 B2WARNING("ECLGeometryPar::ECLVolumeToCellID: History depth = " << depth << " is out of range: should be 3 or 5.")do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Warning, 0, "ecl")) { { LogVariableStream varStream
; varStream << "ECLGeometryPar::ECLVolumeToCellID: History depth = "
<< depth << " is out of range: should be 3 or 5."
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Warning, std::move(varStream), "ecl", __PRETTY_FUNCTION__
, "ecl/geometry/src/ECLGeometryPar.cc", 482, 0)); }; } } while
(false);
483 return -1;
484 }
485 const G4String& vname = touch->GetVolume()->GetName();
486 std::size_t pos0 = vname.find("lv_forward_crystal_");
487 std::size_t pos1 = vname.find("lv_barrel_crystal_");
488 std::size_t pos2 = vname.find("lv_backward_crystal_");
489 if (pos0 == string::npos && pos1 == string::npos && pos2 == string::npos) {
490 B2WARNING("ECLGeometryPar::ECLVolumeToCellID: Volume name does not match pattern. NAME=" << vname)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Warning, 0, "ecl")) { { LogVariableStream varStream
; varStream << "ECLGeometryPar::ECLVolumeToCellID: Volume name does not match pattern. NAME="
<< vname; Belle2::LogSystem::Instance().sendMessage(Belle2
::LogMessage(Belle2::LogConfig::c_Warning, std::move(varStream
), "ecl", __PRETTY_FUNCTION__, "ecl/geometry/src/ECLGeometryPar.cc"
, 490, 0)); }; } } while(false);
491 return -1;
492 }
493 return TouchableToCellID(touch);
494}
495
496double ECLGeometryPar::time2sensor(int cid, const G4ThreeVector& hit_pos)
497{
498 if (cid != m_ini_cid) InitCrystal(cid);
499 double z = 15. - (hit_pos - m_current_crystal.pos) * m_current_crystal.dir; // position along the vector of crystal axis
500 double dt = 6.05 + z * (0.0749 - z * 0.00112); // flight time to diode sensor in nanoseconds
501 return dt;
502}
503
504EclNbr::EclNbr() :
505 m_nbrs(*new std::vector< Identifier >)
506{
507 mNbr_cellID = 0;
508 mNbr_thetaID = 0;
509 mNbr_phiID = 0;
510}
511
512EclNbr::EclNbr(const EclNbr& aNbr) :
513 m_nbrs(*new std::vector< Identifier > (aNbr.m_nbrs)),
514 m_nearSize(aNbr.m_nearSize)
515{
516 mNbr_cellID = 0;
517 mNbr_thetaID = 0;
518 mNbr_phiID = 0;
519}
520
521EclNbr::EclNbr(
522 const std::vector< Identifier >& aNbrs,
523 const std::vector< Identifier >::size_type aNearSize
524) :
525 m_nbrs(*new std::vector< Identifier > (aNbrs)),
526 m_nearSize(aNearSize)
527{
528 // sort vector separately for near, nxt-near nbrs
529 std::sort(m_nbrs.begin(), m_nbrs.begin() + aNearSize, std::less< Identifier >()) ;
530 std::sort(m_nbrs.begin() + aNearSize, m_nbrs.end(), std::less< Identifier >()) ;
531}
532
533EclNbr::~EclNbr()
534{
535 delete &m_nbrs ;
536}
537std::ostream& operator<<(std::ostream& os, const EclNbr& aNbr)
538{
539 os << "N(" ;
540 unsigned short i(0) ;
541 for (std::vector< EclNbr::Identifier >::const_iterator
542 iNbr(aNbr.nbrs().begin()) ;
543 iNbr != aNbr.nbrs().end() ; ++iNbr) {
544 ++i;
545 if (iNbr != aNbr.nbrs().begin() && i != aNbr.nearSize() + 1) os << "," ;
546 if (i == aNbr.nearSize() + 1) os << "|" ;
547 os << std::setw(4) << (*iNbr) ;
548 }
549 os << ")" ;
550 return os ;
551}
552
553
554void EclNbr::printNbr()
555{
556 unsigned short Nri(0) ;
557 cout << "(";
558 for (std::vector< EclNbr::Identifier >::const_iterator
559 iNbr(m_nbrs.begin()) ;
560 iNbr != m_nbrs.end() ; ++iNbr) {
561 ++Nri;
562 if (iNbr != m_nbrs.begin() && Nri != m_nearSize + 1) cout << "," ;
563 if (Nri == m_nearSize + 1) cout << "|" ;
564 cout << std::setw(4) << (*iNbr) ;
565 }
566 cout << ")" << endl;
567}
568//
569// assignment operators
570
571
572EclNbr& EclNbr::operator=(const EclNbr& aNbr)
573{
574 if (this != &aNbr) {
575 mNbr_cellID = aNbr.mNbr_cellID;
576 mNbr_thetaID = aNbr.mNbr_thetaID;
577 mNbr_phiID = aNbr.mNbr_phiID;
578
579 m_nbrs = aNbr.m_nbrs ;
580 m_nearSize = aNbr.m_nearSize ;
581 }
582 return *this ;
583}
584
585//
586// member functions
587//
588
589//
590// const member functions
591//
592
593const std::vector< EclNbr::Identifier >&
594EclNbr::nbrs() const
595{
596 return m_nbrs ;
597}
598
599const std::vector< EclNbr::Identifier >::const_iterator
600EclNbr::nearBegin() const
601{
602 return m_nbrs.begin() ;
603}
604
605const std::vector< EclNbr::Identifier >::const_iterator
606EclNbr::nearEnd() const
607{
608 return m_nbrs.begin() + m_nearSize ;
609}
610
611const std::vector< EclNbr::Identifier >::const_iterator
612EclNbr::nextBegin() const
613{
614 return m_nbrs.begin() + m_nearSize ;
615}
616
617const std::vector< EclNbr::Identifier >::const_iterator
618EclNbr::nextEnd() const
619{
620 return m_nbrs.end() ;
621}
622
623std::vector< EclNbr::Identifier >::size_type
624EclNbr::nearSize() const
625{
626 return m_nearSize ;
627}
628
629std::vector< EclNbr::Identifier >::size_type
630EclNbr::nextSize() const
631{
632 return (m_nbrs.size() - m_nearSize) ;
633}
634
635
636int EclNbr::GetCellID(int ThetaId, int PhiId)
637{
638 mNbr_cellID = Mapping_t::CellID(ThetaId, PhiId);
639 mNbr_thetaID = ThetaId;
640 mNbr_phiID = PhiId ;
641 return mNbr_cellID;
642}
643
644void EclNbr::Mapping(int cid)
645{
646 mNbr_cellID = cid;
647 Mapping_t::Mapping(mNbr_cellID, mNbr_thetaID, mNbr_phiID);
648}
649
650EclNbr
651EclNbr::getNbr(const Identifier aCellId)
652{
653 // generate nbr lists. always easier here to work with theta-phi
654
655 const int cellID = aCellId;
656 Mapping(cellID);
657 const int thetaId = GetThetaID();
658 const int phiId = GetPhiID();
659 std::vector< EclNbr::Identifier >::size_type nearSize(0);
660 std::vector< EclNbr::Identifier > vNbr;
661
662 vNbr.reserve(24) ; // except for extreme endcaps, always 24
663
664 int t00 = thetaId;
665 int tm1 = thetaId - 1;
666 int tm2 = thetaId - 2;
667 int tp1 = thetaId + 1;
668 int tp2 = thetaId + 2;
669
670 if (ECLElementNumbers::isBarrel(aCellId + 1)) {
671 // Barrel.
672 //
673 // 12 13 14 15 16 ^ theta
674 // 11 2 3 4 17 |
675 // 10 1 0 5 18 +--> phi X--+ view from inside
676 // 9 8 7 6 19 | (foot pointing e- dir)
677 // 24 23 22 21 20 Z
678 int f00 = phiId;
679 int fm1 = (phiId + 143) % 144;
680 int fp1 = (phiId + 1) % 144;
681 int fm2 = (phiId + 142) % 144;
682 int fp2 = (phiId + 2) % 144;
683
684 vNbr.push_back(GetCellID(t00, fm1));
685 vNbr.push_back(GetCellID(tp1, fm1));
686 vNbr.push_back(GetCellID(tp1, f00));
687 vNbr.push_back(GetCellID(tp1, fp1));
688 vNbr.push_back(GetCellID(t00, fp1));
689 vNbr.push_back(GetCellID(tm1, fp1));
690 vNbr.push_back(GetCellID(tm1, f00));
691 vNbr.push_back(GetCellID(tm1, fm1));
692
693 nearSize = vNbr.size();
694
695 vNbr.push_back(GetCellID(tm1, fm2));
696 vNbr.push_back(GetCellID(t00, fm2));
697 vNbr.push_back(GetCellID(tp1, fm2));
698 vNbr.push_back(GetCellID(tp2, fm2));
699 vNbr.push_back(GetCellID(tp2, fm1));
700 vNbr.push_back(GetCellID(tp2, f00));
701 vNbr.push_back(GetCellID(tp2, fp1));
702 vNbr.push_back(GetCellID(tp2, fp2));
703 vNbr.push_back(GetCellID(tp1, fp2));
704 vNbr.push_back(GetCellID(t00, fp2));
705 vNbr.push_back(GetCellID(tm1, fp2));
706 vNbr.push_back(GetCellID(tm2, fp2));
707 vNbr.push_back(GetCellID(tm2, fp1));
708 vNbr.push_back(GetCellID(tm2, f00));
709 vNbr.push_back(GetCellID(tm2, fm1));
710 vNbr.push_back(GetCellID(tm2, fm2));
711 } else {
712 // Forward or backward.
713 // endcap -- not always 24!
714 int n00 = 1000;
715 int np1 = 1000;
716 int np2 = 1000;
717 int nm1 = 1000;
718 int nm2 = 1000;
719 if (ECLElementNumbers::isForward(aCellId + 1)) {
720 // Forward.
721 const EclIdentifier mPerRingForward[]
722 = { 48, 48, 64, 64, 64, 96, 96, 96, 96, 96, 96, 144, 144, 144, 144 };
723 if (thetaId > 1) nm2 = mPerRingForward[ thetaId - 2 ];
724 if (thetaId > 0) nm1 = mPerRingForward[ thetaId - 1 ];
725 n00 = mPerRingForward[ thetaId ];
726 np1 = mPerRingForward[ thetaId + 1 ];
727 np2 = mPerRingForward[ thetaId + 2 ];
728 } else {
729 // Backward.
730 const EclIdentifier mPerRingBackward[]
731 = { 64, 64, 64, 96, 96, 96, 96, 96, 144, 144, 144, 144 };
732 if (thetaId < 67) np2 = mPerRingBackward[ 66 - thetaId ];
733 if (thetaId < 68) np1 = mPerRingBackward[ 67 - thetaId ];
734 n00 = mPerRingBackward[ 68 - thetaId ];
735 nm1 = mPerRingBackward[ 69 - thetaId ];
736 nm2 = mPerRingBackward[ 70 - thetaId ];
737 }
738 // f-- are phi's, t-- are thetas
739 // all calculations should be integer arith - pcs
740 // f(th,phi)
741 // criterion: center -> next bin
742 int f0000 = phiId;
743 int fp100 = (f0000 * np1 + np1 / 2) / n00;
744 int fp200 = (f0000 * np2 + np2 / 2) / n00;
745 int fm100 = (f0000 * nm1 + nm1 / 2) / n00;
746 int fm200 = (f0000 * nm2 + nm2 / 2) / n00;
747
748 int f00m1 = (f0000 + n00 - 1) % n00; // should be exact
749 int f00m2 = (f0000 + n00 - 2) % n00;
750 int f00p1 = (f0000 + 1) % n00;
751 int f00p2 = (f0000 + 2) % n00;
752
753 int fp1m1 = (fp100 + np1 - 1) % np1;
754 int fp1m2 = (fp100 + np1 - 2) % np1;
755 int fp1p1 = (fp100 + 1) % np1;
756 int fp1p2 = (fp100 + 2) % np1;
757
758 int fm1m1 = (fm100 + nm1 - 1) % nm1;
759 int fm1m2 = (fm100 + nm1 - 2) % nm1;
760 int fm1p1 = (fm100 + 1) % nm1;
761 int fm1p2 = (fm100 + 2) % nm1;
762
763 int fp2m1 = (fp200 + np2 - 1) % np2;
764 int fp2m2 = (fp200 + np2 - 2) % np2;
765 int fp2p1 = (fp200 + 1) % np2;
766 int fp2p2 = (fp200 + 2) % np2;
767
768 int fm2m1 = (fm200 + nm2 - 1) % nm2;
769 int fm2m2 = (fm200 + nm2 - 2) % nm2;
770 int fm2p1 = (fm200 + 1) % nm2;
771 int fm2p2 = (fm200 + 2) % nm2;
772 int delta = n00 / 16;
773// int sector = phiId/delta; // 0..15
774 int nth = phiId % delta;
775
776 switch (thetaId) {
777 case 0:
778 if (nth == 1)
779 fp2p2 = 1000;
780 break;
781 case 1:
782 if (nth == 1) {
783 fp2p2 = 1000;
784 fp1p2 = fp1p1;
785 fp1p1 = 1000;
786 }
787 break;
788 case 2:
789 if ((nth == 0) || (nth == 1)) {
790 fm2p2 = 1000;
791 fm1p2 = fm1p1;
792 fm1p1 = 1000;
793 } else if ((nth == 2) || (nth == 3)) {
794 fm2m2 = 1000;
795 fm1m2 = fm1m1;
796 fm1m1 = 1000;
797 }
798 break;
799 case 3:
800 if ((nth == 0) || (nth == 3)) {
801 fm2p2 = fm2m2 = 1000;
802 } else if (nth == 1) {
803 fm2p2 = 1000;
804 } else if (nth == 2) {
805 fm2m2 = 1000;
806 }
807 break;
808 case 5:
809 if ((nth == 2) || (nth == 5)) {
810 fm2m2 = 1000;
811 fm1m2 = fm1m1;
812 fm1m1 = 1000;
813 } else {
814 fm2p2 = 1000;
815 fm1p2 = fm1p1;
816 fm1p1 = 1000;
817 }
818 break;
819 case 6:
820 fm2p2 = 1000;
821 if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5)) {
822 fm2m2 = 1000;
823 }
824 break;
825 case 11:
826 if ((nth == 2) || (nth == 5) || (nth == 8)) {
827 fm2m2 = 1000;
828 fm1m2 = fm1m1;
829 fm1m1 = 1000;
830 } else {
831 fm2p2 = 1000;
832 fm1p2 = fm1p1;
833 fm1p1 = 1000;
834 }
835 break;
836 case 12:
837 fm2p2 = 1000;
838 if ((nth == 0) || (nth == 2) || (nth == 3)
839 || (nth == 5) || (nth == 6) || (nth == 8))
840 fm2m2 = 1000;
841 break;
842 case 65:
843 if ((nth == 2) || (nth == 5)) {
844 fp2m2 = 1000;
845 fp1m2 = fp1m1;
846 fp1m1 = 1000;
847 } else {
848 fp2p2 = 1000;
849 fp1p2 = fp1p1;
850 fp1p1 = 1000;
851 }
852 break;
853 case 64:
854 fp2p2 = 1000;
855 if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5))
856 fp2m2 = 1000;
857 break;
858 case 60:
859 if ((nth == 2) || (nth == 5) || (nth == 8)) {
860 fp2m2 = 1000;
861 fp1m2 = fp1m1;
862 fp1m1 = 1000;
863 } else {
864 fp2p2 = 1000;
865 fp1p2 = fp1p1;
866 fp1p1 = 1000;
867 }
868 break;
869 case 59:
870 fp2p2 = 1000;
871 if ((nth == 0) || (nth == 2) || (nth == 3) || (nth == 5)
872 || (nth == 6) || (nth == 8))
873 fp2m2 = 1000;
874 break;
875 }//switch
876
877 // insert near-nbrs
878 vNbr.push_back(GetCellID(t00, f00m1));
879 vNbr.push_back(GetCellID(t00, f00p1));
880 if (nm1 < 999) {
881 vNbr.push_back(GetCellID(tm1, fm100));
882 if (fm1m1 < 999)
883 vNbr.push_back(GetCellID(tm1, fm1m1));
884 if (fm1p1 < 999)
885 vNbr.push_back(GetCellID(tm1, fm1p1));
886 }
887 if (np1 < 999) {
888 vNbr.push_back(GetCellID(tp1, fp100));
889 if (fp1m1 < 999)
890 vNbr.push_back(GetCellID(tp1, fp1m1));
891 if (fp1p1 < 999)
892 vNbr.push_back(GetCellID(tp1, fp1p1));
893 }
894 nearSize = vNbr.size() ;
895
896 // now on to next-near neighbors
897 if (nm2 < 999) {
898 vNbr.push_back(GetCellID(tm2, fm200));
899 if (fm2m1 < 999)
900 vNbr.push_back(GetCellID(tm2, fm2m1));
901 if (fm2p1 < 999)
902 vNbr.push_back(GetCellID(tm2, fm2p1));
903 if (fm2m2 < 999)
904 vNbr.push_back(GetCellID(tm2, fm2m2));
905 if (fm2p2 < 999)
906 vNbr.push_back(GetCellID(tm2, fm2p2));
907 }
908 if (nm1 < 999) {
909 if (fm1m2 < 999)
910 vNbr.push_back(GetCellID(tm1, fm1m2));
911 if (fm1p2 < 999)
912 vNbr.push_back(GetCellID(tm1, fm1p2));
913 }
914 vNbr.push_back(GetCellID(t00, f00m2));
915 vNbr.push_back(GetCellID(t00, f00p2));
916 if (np1 < 999) {
917 if (fp1m2 < 999)
918 vNbr.push_back(GetCellID(tp1, fp1m2));
919 if (fp1p2 < 999)
920 vNbr.push_back(GetCellID(tp1, fp1p2));
921 }
922 if (np2 < 999) {
923 vNbr.push_back(GetCellID(tp2, fp200));
924 if (fp2m1 < 999)
925 vNbr.push_back(GetCellID(tp2, fp2m1));
926 if (fp2p1 < 999)
927 vNbr.push_back(GetCellID(tp2, fp2p1));
928 if (fp2m2 < 999)
929 vNbr.push_back(GetCellID(tp2, fp2m2));
930 if (fp2p2 < 999)
931 vNbr.push_back(GetCellID(tp2, fp2p2));
932 }
933 }
934 return EclNbr(vNbr, nearSize);
935}