1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
/// @file CoaxialConesStruct.h
/// @author Raman Sehgal (raman.sehgal@cern.ch)
#ifndef VECGEOM_VOLUMES_COAXIALCONESSTRUCT_H_
#define VECGEOM_VOLUMES_COAXIALCONESSTRUCT_H_
#include "VecGeom/base/Global.h"
#include "VecGeom/base/Vector3D.h"
#include "VecGeom/base/Vector.h"
#include "VecGeom/volumes/ConeStruct.h"
namespace vecgeom {
inline namespace VECGEOM_IMPL_NAMESPACE {
template <typename T = double>
struct CoaxialConesStruct {
// CoaxialCone parameters
unsigned int fNumOfCones;
T fDz;
T fSPhi;
T fDPhi;
Vector<T> fRmin1Vect;
Vector<T> fRmax1Vect;
Vector<T> fRmin2Vect;
Vector<T> fRmax2Vect;
/*
** A function to Print the sequence of Coaxial cones,
** Useful for debugging the Coaxial cones
*/
void Print() const
{
for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
std::cout << std::endl << "====== Index of cone struct : " << i << " ======" << std::endl;
fConeStructVector[i]->Print();
std::cout << std::endl;
}
}
VECCORE_ATT_HOST_DEVICE
CoaxialConesStruct() {}
VECCORE_ATT_HOST_DEVICE
CoaxialConesStruct(unsigned int numOfCones, T *rmin1Vect, T *rmax1Vect, T *rmin2Vect, T *rmax2Vect, T dz, T sphi,
T dphi)
: fNumOfCones(numOfCones), fDz(dz), fSPhi(sphi), fDPhi(dphi)
{
for (unsigned int i = 0; i < fNumOfCones; i++) {
fRmin1Vect.push_back(rmin1Vect[i]);
fRmax1Vect.push_back(rmax1Vect[i]);
fRmin2Vect.push_back(rmin2Vect[i]);
fRmax2Vect.push_back(rmax2Vect[i]);
}
for (unsigned int i = 0; i < fNumOfCones; i++) {
fConeStructVector.push_back(
new ConeStruct<T>(fRmin1Vect[i], fRmax1Vect[i], fRmin2Vect[i], fRmax2Vect[i], dz, sphi, dphi));
if (i == 0) {
fMinR = fRmin1Vect[i] < fRmin2Vect[i] ? fRmin1Vect[i] : fRmin2Vect[i];
}
if (i == fNumOfCones - 1) {
fMaxR = fRmax1Vect[i] > fRmax2Vect[i] ? fRmax1Vect[i] : fRmax2Vect[i];
}
}
}
VECCORE_ATT_HOST_DEVICE
CoaxialConesStruct(Vector<Precision> rmin1Vect, Vector<Precision> rmax1Vect, Vector<Precision> rmin2Vect,
Vector<Precision> rmax2Vect, T dz, T sphi, T dphi)
: fNumOfCones(rmin1Vect.size()), fDz(dz), fSPhi(sphi), fDPhi(dphi), fRmin1Vect(rmin1Vect), fRmax1Vect(rmax1Vect),
fRmin2Vect(rmin2Vect), fRmax2Vect(rmax2Vect)
{
for (unsigned int i = 0; i < fNumOfCones; i++) {
fConeStructVector.push_back(
new ConeStruct<T>(fRmin1Vect[i], fRmax1Vect[i], fRmin2Vect[i], fRmax2Vect[i], dz, sphi, dphi));
if (i == 0) {
fMinR = fRmin1Vect[i] < fRmin2Vect[i] ? fRmin1Vect[i] : fRmin2Vect[i];
}
if (i == fNumOfCones - 1) {
fMaxR = fRmax1Vect[i] > fRmax2Vect[i] ? fRmax1Vect[i] : fRmax2Vect[i];
}
}
}
// Vector of Cones
Vector<ConeStruct<T> *> fConeStructVector;
T fSurfaceArea; // area of the surface
T fCubicVolume; // volume
T fMaxR;
T fMinR;
// Precalculated cached values
VECCORE_ATT_HOST_DEVICE
Precision Capacity()
{
Precision volume = 0.;
for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
volume += fConeStructVector[i]->Capacity();
}
return volume;
}
VECCORE_ATT_HOST_DEVICE
Precision ConicalSurfaceArea()
{
Precision conicalSurfaceArea = 0.;
Precision mmin, mmax, dmin, dmax;
for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
mmin = (fConeStructVector[i]->fRmin1 + fConeStructVector[i]->fRmin2) * 0.5;
mmax = (fConeStructVector[i]->fRmax1 + fConeStructVector[i]->fRmax2) * 0.5;
dmin = (fConeStructVector[i]->fRmin2 - fConeStructVector[i]->fRmin1);
dmax = (fConeStructVector[i]->fRmax2 - fConeStructVector[i]->fRmax1);
conicalSurfaceArea +=
fConeStructVector[i]->fDPhi *
(mmin * vecCore::math::Sqrt(dmin * dmin + 4 * fConeStructVector[i]->fDz * fConeStructVector[i]->fDz) +
mmax * vecCore::math::Sqrt(dmax * dmax + 4 * fConeStructVector[i]->fDz * fConeStructVector[i]->fDz));
}
return conicalSurfaceArea;
}
VECCORE_ATT_HOST_DEVICE
Precision SurfaceAreaLowerZPlanes(int index)
{
return fConeStructVector[index]->fDPhi * 0.5 *
(fConeStructVector[index]->fRmax1 * fConeStructVector[index]->fRmax1 -
fConeStructVector[index]->fRmin1 * fConeStructVector[index]->fRmin1);
}
VECCORE_ATT_HOST_DEVICE
Precision SurfaceAreaUpperZPlanes(int index)
{
return fConeStructVector[index]->fDPhi * 0.5 *
(fConeStructVector[index]->fRmax2 * fConeStructVector[index]->fRmax2 -
fConeStructVector[index]->fRmin2 * fConeStructVector[index]->fRmin2);
}
VECCORE_ATT_HOST_DEVICE
Precision SurfaceAreaOfZPlanes()
{
Precision surfaceAreaOfZPlanes = 0.;
for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
surfaceAreaOfZPlanes = SurfaceAreaLowerZPlanes(i) + SurfaceAreaUpperZPlanes(i);
}
return surfaceAreaOfZPlanes;
}
VECCORE_ATT_HOST_DEVICE
Precision SurfaceArea() { return (ConicalSurfaceArea() + SurfaceAreaOfZPlanes()); }
VECCORE_ATT_HOST_DEVICE
bool Normal(Vector3D<Precision> const &p, Vector3D<Precision> &norm) const
{
norm.Set(0.);
bool valid = false;
for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
bool validNormal = false;
Vector3D<Precision> normal(0., 0., 0.);
validNormal = fConeStructVector[i]->Normal(p, normal);
if (validNormal) {
norm += normal;
}
valid |= validNormal;
}
if (valid) {
norm.Normalize();
} else {
norm.Set(0., 0., 1.);
}
return valid;
}
};
} // namespace VECGEOM_IMPL_NAMESPACE
} // namespace vecgeom
#endif
|
