BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ExtBesEmcGeometry.cxx
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2// BESIII Object_Oreiented Simulation and Reconstruction Tool //
3//---------------------------------------------------------------------------//
4// Descpirtion: Geometry of EMC detector
5// Author: Fu Chengdong
6// Created: Oct 23, 2003
7// Comment:
8//---------------------------------------------------------------------------//
9//
10#include "G4ThreeVector.hh"
11
12#include "ExtBesEmcGeometry.h"
13#include "ExtBesEmcParameter.h"
14#include "G4SystemOfUnits.hh"
15#include "G4ios.hh"
16
17ExtBesEmcGeometry::ExtBesEmcGeometry() { verboseLevel = 0; }
18
19ExtBesEmcGeometry::~ExtBesEmcGeometry() { ; }
20
21void ExtBesEmcGeometry::ReadEMCParameters() {
22 // Read EMC parameters from database
23 ExtBesEmcParameter& emcPara = ExtBesEmcParameter::GetInstance();
24 // BesEmcParameter emcPara;
25 // emcPara.ReadData();
26
27 TaperRingDz = emcPara.GetTaperRingDz();
28 TaperRingThickness1 = emcPara.GetTaperRingThickness1();
29 TaperRingThickness2 = emcPara.GetTaperRingThickness2();
30 TaperRingThickness3 = emcPara.GetTaperRingThickness3();
31 TaperRingTheta = emcPara.GetTaperRingTheta() * deg;
32 TaperRingInnerLength = emcPara.GetTaperRingInnerLength();
33 TaperRingOuterLength = emcPara.GetTaperRingOuterLength();
34
35 BSCRmin = emcPara.GetBSCRmin();
36 BSCDz = emcPara.GetBSCDz() - TaperRingThickness3;
37
38 BSCRmin1 = emcPara.GetBSCRmin1();
39 BSCRmax1 = emcPara.GetBSCRmax1();
40 BSCRmin2 = emcPara.GetBSCRmin2();
41 // BSCRmax2 = emcPara.GetBSCRmax2();
42 BSCDz1 = emcPara.GetBSCDz1() - TaperRingThickness3;
43
44 BSCAngleRotat = emcPara.GetBSCAngleRotat() * deg;
45 BSCNbPhi = emcPara.GetBSCNbPhi();
46 BSCNbTheta = emcPara.GetBSCNbTheta();
47
48 BSCCryLength = emcPara.GetCrystalLength();
49 BSCYFront0 = emcPara.GetBSCYFront0();
50 BSCYFront = emcPara.GetBSCYFront();
51 BSCYFront1 = emcPara.GetBSCYFront1();
52 BSCPosition0 = emcPara.GetBSCPosition0();
53 BSCPosition1 = emcPara.GetBSCPosition1();
54
55 fTyvekThickness = emcPara.GetTyvekThickness();
56 fAlThickness = emcPara.GetAlThickness();
57 fMylarThickness = emcPara.GetMylarThickness();
58
59 rearBoxLength = emcPara.GetRearBoxLength();
60 rearBoxDz = emcPara.GetRearBoxDz();
61
62 HangingPlateDz = emcPara.GetHangingPlateDz();
63 OCGirderAngle = emcPara.GetOCGirderAngle() * deg;
64
65 rearCasingThickness = emcPara.GetRearCasingThickness();
66
67 orgGlassLengthX = emcPara.GetOrgGlassLengthX();
68 orgGlassLengthY = emcPara.GetOrgGlassLengthY();
69 orgGlassLengthZ = emcPara.GetOrgGlassLengthZ();
70
71 PDLengthX = emcPara.GetPDLengthX();
72 PDLengthY = emcPara.GetPDLengthY();
73 PDLengthZ = emcPara.GetPDLengthZ();
74
75 AlPlateDz = emcPara.GetAlPlateDz();
76 PABoxDz = emcPara.GetPABoxDz();
77 PABoxThickness = emcPara.GetPABoxThickness();
78
79 cableDr = emcPara.GetCableDr();
80 waterPipeDr = emcPara.GetWaterPipeDr();
81 waterPipeThickness = emcPara.GetWaterPipeThickness();
82
83 SPBarThickness = emcPara.GetSPBarThickness();
84 SPBarThickness1 = emcPara.GetSPBarThickness1();
85 SPBarwidth = emcPara.GetSPBarwidth();
86
87 EndRingDz = emcPara.GetEndRingDz();
88 EndRingDr = emcPara.GetEndRingDr();
89 EndRingRmin = emcPara.GetEndRingRmin();
90
91 for ( G4int i = 0; i < BSCNbTheta * 6; i++ )
92 {
93 zHalfLength[i] = 0;
94 thetaAxis[i] = 0;
95 phiAxis[i] = 0;
96 yHalfLength1[i] = 0;
97 xHalfLength1[i] = 0;
98 xHalfLength2[i] = 0;
99 tanAlpha1[i] = 0;
100 yHalfLength2[i] = 0;
101 xHalfLength3[i] = 0;
102 xHalfLength4[i] = 0;
103 tanAlpha2[i] = 0;
104 thetaPosition[i] = 0;
105 xPosition[i] = 0;
106 yPosition[i] = 0;
107 zPosition[i] = 0;
108 }
109}
110
111void ExtBesEmcGeometry::PrintEMCParameters() {
112 // Print EMC parameters
113 ;
114}
115
116void ExtBesEmcGeometry::ComputeEMCParameters() {
117 ReadEMCParameters();
118
119 // Compute derived parameters of the calorimeter
120 G4double theta0;
121 G4int i;
122 // for debug
123 // G4Exception("ExtBesEmcGeometry::ComputeEMCParameters() starting.......");
124 //
125 // support frame
126 //
127 const G4double delta = 1. * mm; // for opening-cut girder
128 TaperRingRmin1 = BSCRmax1 - TaperRingThickness1;
129 TaperRingRmin2 = TaperRingRmin1 + TaperRingDz * tan( TaperRingTheta );
130 TaperRingDr = TaperRingThickness2 / cos( TaperRingTheta );
131 TaperRingOuterLength1 = TaperRingOuterLength + TaperRingThickness3 * tan( TaperRingTheta );
132
133 BSCRmax2 = TaperRingRmin2 + TaperRingDr - TaperRingThickness3 * tan( TaperRingTheta );
134 BSCRmax = BSCRmin + 33.8 * cm;
135 BSCPhiDphi = 360. * deg / BSCNbPhi;
136 BSCPhiSphi = -BSCPhiDphi / 2.;
137 BSCPhiDz = BSCDz;
138 BSCPhiRmax = BSCRmax - 0.5 * cm;
139 BSCPhiRmin = BSCRmin / sin( 90. * deg + BSCPhiDphi / 2. ) *
140 sin( 90. * deg + BSCPhiDphi / 2. - BSCAngleRotat ); // Rmin after rotate
141 SPBarDphi = 2 * atan( SPBarwidth / ( 2 * BSCRmax ) );
142
143 // crystal No.1(from z=0 to big)
144 zHalfLength[0] = BSCCryLength / 2.;
145 yHalfLength1[0] = BSCYFront0 / 2.;
146 xHalfLength1[0] = BSCPhiRmin * tan( BSCPhiDphi ) / 2.;
147 xHalfLength2[0] = xHalfLength1[0];
148 G4double rminprojected = BSCRmin * cos( BSCAngleRotat );
149 // rminprojected=BSCRmin;//according to hardware design
150 yHalfLength2[0] =
151 yHalfLength1[0] + ( BSCYFront0 - BSCPosition0 ) / rminprojected * BSCCryLength / 2.;
152 xHalfLength3[0] = xHalfLength1[0] * ( BSCPhiRmin + BSCCryLength ) / BSCPhiRmin;
153 xHalfLength4[0] = xHalfLength2[0] * ( BSCPhiRmin + BSCCryLength ) / BSCPhiRmin;
154
155 thetaPosition[0] = 90. * deg;
156 xPosition[0] = BSCPhiRmin + BSCCryLength / 2.;
157 yPosition[0] =
158 -( xHalfLength1[0] + xHalfLength3[0] + xHalfLength2[0] + xHalfLength4[0] ) / 4.;
159 zPosition[0] = ( yHalfLength1[0] + yHalfLength2[0] ) / 2.;
160
161 theta0 = 90. * deg - atan( ( BSCYFront0 - BSCPosition0 ) / rminprojected );
162 if ( verboseLevel > 1 )
163 {
164 G4cout << "------------------------------------>1" << G4endl
165 << "The direction of crystal is:" << theta0 / deg << "(deg)." << G4endl;
166 }
167
168 rearBoxPosX[0] = xPosition[0] + BSCCryLength / 2. + rearBoxDz / 2.;
169 rearBoxPosY[0] = -( xHalfLength3[0] + xHalfLength4[0] ) / 2.;
170 rearBoxPosZ[0] = yHalfLength2[0];
171
172 OCGirderRmin1[0] = rearBoxPosX[0] + rearBoxDz / 2. + delta;
173 OCGirderRmin2[0] = rearBoxPosX[0] + rearBoxDz / 2. + delta;
174 OCGirderDz[0] = rearBoxLength;
175 OCGirderPosZ[0] = rearBoxLength / 2;
176
177 cableLength[0] = BSCDz - rearBoxPosZ[0];
178 cablePosX[0] = BSCPhiRmax - cableDr - 2. * mm;
179 cablePosY[0] = rearBoxPosY[0] - 3 * cableDr;
180 cablePosZ[0] = BSCDz - cableLength[0] / 2;
181
182 // crystal No.2
183 zHalfLength[1] = BSCCryLength / 2.;
184 yHalfLength1[1] = BSCYFront0 / 2.;
185 G4double deltaR = BSCYFront0 * cos( theta0 );
186 xHalfLength1[1] = xHalfLength1[0];
187 xHalfLength2[1] = xHalfLength1[1] * ( BSCPhiRmin + deltaR ) / BSCPhiRmin;
188 yHalfLength2[1] =
189 yHalfLength1[1] +
190 tan( theta0 - atan( rminprojected /
191 ( BSCYFront0 * ( 1. + 1. / sin( theta0 ) ) - BSCPosition1 ) ) ) *
192 BSCCryLength / 2.;
193 deltaR = deltaR + BSCCryLength * sin( theta0 );
194 xHalfLength4[1] = xHalfLength1[1] * ( BSCPhiRmin + deltaR ) / BSCPhiRmin;
195 deltaR = deltaR - yHalfLength2[1] * 2. * cos( theta0 );
196 xHalfLength3[1] = xHalfLength1[1] * ( BSCPhiRmin + deltaR ) / BSCPhiRmin;
197
198 thetaPosition[1] = theta0;
199 xPosition[1] = BSCPhiRmin + BSCCryLength / 2. * sin( theta0 ) +
200 ( 3. * yHalfLength1[1] - yHalfLength2[1] ) / 2. * cos( theta0 );
201 yPosition[1] =
202 -( xHalfLength1[1] + xHalfLength3[1] + xHalfLength2[1] + xHalfLength4[1] ) / 4.;
203 zPosition[1] = yHalfLength1[0] * 2. +
204 ( yHalfLength1[1] * 2. / tan( theta0 ) + BSCCryLength / 2. ) * cos( theta0 ) +
205 ( yHalfLength1[1] + yHalfLength2[1] ) / 2. * sin( theta0 );
206
207 rearBoxPosX[1] = xPosition[1] + ( BSCCryLength / 2. + rearBoxDz / 2. ) * sin( theta0 );
208 rearBoxPosY[1] = -( xHalfLength3[1] + xHalfLength4[1] ) / 2.;
209 rearBoxPosZ[1] = zPosition[1] + ( BSCCryLength / 2. + rearBoxDz / 2. ) * cos( theta0 );
210
211 OCGirderRmin1[1] = rearBoxPosX[1] + rearBoxDz * sin( theta0 ) / 2. +
212 rearBoxLength * cos( theta0 ) / 2 + delta;
213 OCGirderRmin2[1] = rearBoxPosX[1] + rearBoxDz * sin( theta0 ) / 2. -
214 rearBoxLength * cos( theta0 ) / 2 + delta;
215 OCGirderDz[1] = rearBoxLength * sin( theta0 );
216 OCGirderPosZ[1] = rearBoxPosZ[1] + rearBoxDz * cos( theta0 ) / 2.;
217
218 cableLength[1] = BSCDz - rearBoxPosZ[1];
219 cablePosX[1] = cablePosX[0];
220 cablePosY[1] = cablePosY[0] + 2 * cableDr;
221 cablePosZ[1] = BSCDz - cableLength[1] / 2;
222
223 theta0 = theta0 - atan( ( yHalfLength2[1] - yHalfLength1[1] ) * 2. / BSCCryLength );
224
225 for ( i = 2; i < BSCNbTheta; i++ )
226 {
227 if ( verboseLevel > 1 )
228 {
229 G4cout << "------------------------------------>" << i << G4endl
230 << "The direction of crystal is:" << theta0 / deg << "(deg)." << G4endl;
231 }
232 // crystal No.i+1
233 zHalfLength[i] = zHalfLength[0];
234 yHalfLength1[i] = BSCYFront / 2.;
235 if ( i == BSCNbTheta - 1 ) { yHalfLength1[i] = BSCYFront1 / 2.; } // the rightest crystal
236 deltaR = yHalfLength1[i] * 2. * cos( theta0 );
237 xHalfLength1[i] = xHalfLength1[0];
238 xHalfLength2[i] = xHalfLength1[i] / BSCPhiRmin * ( BSCPhiRmin + deltaR );
239 yHalfLength2[i] =
240 yHalfLength1[i] * ( 1. + BSCCryLength / ( rminprojected / sin( theta0 ) +
241 yHalfLength1[i] * 2. / tan( theta0 ) ) );
242 deltaR = deltaR + BSCCryLength * sin( theta0 );
243 xHalfLength4[i] = xHalfLength1[i] / BSCPhiRmin * ( BSCPhiRmin + deltaR );
244 deltaR = deltaR - yHalfLength2[i] * 2. * cos( theta0 );
245 xHalfLength3[i] = xHalfLength1[i] / BSCPhiRmin * ( BSCPhiRmin + deltaR );
246
247 thetaPosition[i] = theta0;
248 xPosition[i] = BSCPhiRmin + zHalfLength[i] * sin( theta0 ) +
249 ( 3. * yHalfLength1[i] - yHalfLength2[i] ) / 2. * cos( theta0 );
250 yPosition[i] =
251 -( xHalfLength1[i] + xHalfLength3[i] + xHalfLength2[i] + xHalfLength4[i] ) / 4.;
252 zPosition[i] = BSCPosition1 + rminprojected / tan( theta0 ) +
253 ( 2. * yHalfLength1[i] / tan( theta0 ) + zHalfLength[i] ) * cos( theta0 ) +
254 ( yHalfLength1[i] + yHalfLength2[i] ) / 2. * sin( theta0 );
255
256 rearBoxPosX[i] = xPosition[i] + ( BSCCryLength / 2. + rearBoxDz / 2. ) * sin( theta0 );
257 rearBoxPosY[i] = -( xHalfLength3[i] + xHalfLength4[i] ) / 2.;
258 rearBoxPosZ[i] = zPosition[i] + ( BSCCryLength / 2. + rearBoxDz / 2. ) * cos( theta0 );
259
260 OCGirderRmin1[i] = rearBoxPosX[i] + rearBoxDz * sin( theta0 ) / 2. +
261 rearBoxLength * cos( theta0 ) / 2 + delta;
262 OCGirderRmin2[i] = rearBoxPosX[i] + rearBoxDz * sin( theta0 ) / 2. -
263 rearBoxLength * cos( theta0 ) / 2 + delta;
264 OCGirderDz[i] = rearBoxLength * sin( theta0 );
265 OCGirderPosZ[i] = rearBoxPosZ[i] + rearBoxDz * cos( theta0 ) / 2.;
266
267 G4int xCable, yCable;
268 xCable = i / 4;
269 yCable = i - 4 * (G4int)( i / 4 );
270 cableLength[i] = BSCDz - ( rearBoxPosZ[i] + rearBoxDz / 2. * cos( theta0 ) );
271 cablePosX[i] = cablePosX[0] - xCable * 2 * cableDr;
272 cablePosY[i] = cablePosY[0] + yCable * 2 * cableDr;
273 cablePosZ[i] = BSCDz - cableLength[i] / 2;
274
275 theta0 = theta0 -
276 atan( 2. * yHalfLength1[i] /
277 ( rminprojected / sin( theta0 ) + 2. * yHalfLength1[i] / tan( theta0 ) ) );
278 }
279 thetaPosition[BSCNbTheta] = theta0;
280 if ( verboseLevel > 1 )
281 {
282 G4cout << "------------------------------------>" << i << G4endl
283 << "The direction of crystal is:" << theta0 / deg << "(deg)." << G4endl;
284 }
285
286 G4double oop;
287 for ( i = 0; i < BSCNbTheta; i++ )
288 {
289 oop = sqrt( ( yHalfLength2[i] - yHalfLength1[i] ) * ( yHalfLength2[i] - yHalfLength1[i] ) +
290 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) *
291 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) /
292 4 );
293 thetaAxis[i] = atan( oop / BSCCryLength );
294 phiAxis[i] =
295 180. * deg +
296 atan( ( yHalfLength2[i] - yHalfLength1[i] ) /
297 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) * 2. );
298 tanAlpha2[i] = -( xHalfLength4[i] - xHalfLength3[i] ) / yHalfLength2[i] / 2.;
299 tanAlpha1[i] = -( xHalfLength2[i] - xHalfLength1[i] ) / yHalfLength1[i] / 2.;
300 // G4cout <<thetaAxis[i]<<", "
301 // <<phiAxis[i]<<", "
302 // <<tanAlpha1[i]<<", "
303 // <<tanAlpha2[i]<<G4endl;
304 }
305}
306
307void ExtBesEmcGeometry::ModifyForCasing() {
308 // Compute the sizes of the naked crystals and the casing
309 // Casing size
310 // BSCNbTheta---->2*BSCNbTheta-1 : before crystal
311 // 2*BSCNbTheta-->3*BSCNbTheta-1 : a side of crystal
312 // 3*BSCNbTheta-->4*BSCNbTheta-1 : b side of crystal
313 // 4*BSCNbTheta-->5*BSCNbTheta-1 : c side of crystal
314 // 5*BSCNbTheta-->6*BSCNbTheta-1 : d side of crystal
315 // d
316 // ----------
317 // | |
318 // | |
319 // c | | a
320 // | |
321 // | |
322 // | /
323 // | /
324 // | /
325 // b
326 G4double totalThickness = fTyvekThickness + fAlThickness + fMylarThickness; //
327 G4double delta = 0., angle1 = 0. * deg, angle2 = 0. * deg;
328 G4double oop; // the distance of the centers of the two parallel plane
329
330 G4double rminprojected = BSCRmin * cos( BSCAngleRotat );
331 rminprojected = BSCRmin; // accord with hardware design
332
333 G4int i;
334 for ( i = 0; i < BSCNbTheta; i++ )
335 {
336 zHalfLength[BSCNbTheta + i] = totalThickness / 2.;
337 yHalfLength1[BSCNbTheta + i] = yHalfLength1[i];
338 yHalfLength2[BSCNbTheta + i] = yHalfLength1[i] + ( yHalfLength2[i] - yHalfLength1[i] ) *
339 totalThickness / BSCCryLength;
340 xHalfLength1[BSCNbTheta + i] = xHalfLength1[i];
341 xHalfLength2[BSCNbTheta + i] = xHalfLength2[i];
342 xHalfLength1[BSCNbTheta * 2 + i] = xHalfLength3[BSCNbTheta + i] =
343 xHalfLength1[i] * ( BSCCryLength - totalThickness ) / BSCCryLength +
344 xHalfLength3[i] * totalThickness / BSCCryLength;
345 xHalfLength2[BSCNbTheta * 4 + i] = xHalfLength4[BSCNbTheta + i] =
346 xHalfLength2[i] * ( BSCCryLength - totalThickness ) / BSCCryLength +
347 xHalfLength4[i] * totalThickness / BSCCryLength;
348
349 zHalfLength[BSCNbTheta * 5 + i] = zHalfLength[BSCNbTheta * 4 + i] =
350 zHalfLength[BSCNbTheta * 3 + i] = zHalfLength[BSCNbTheta * 2 + i] =
351 zHalfLength[i] - totalThickness / 2.;
352
353 yHalfLength2[BSCNbTheta * 2 + i] = yHalfLength1[BSCNbTheta * 2 + i] =
354 totalThickness / cos( thetaPosition[i + 1] - thetaPosition[i] ) / 2.;
355 xHalfLength3[BSCNbTheta * 2 + i] = xHalfLength3[i];
356 xHalfLength4[BSCNbTheta * 2 + i] =
357 xHalfLength3[i] + ( xHalfLength4[i] - xHalfLength3[i] ) *
358 yHalfLength2[BSCNbTheta * 2 + i] / yHalfLength2[i];
359 xHalfLength2[BSCNbTheta * 2 + i] =
360 xHalfLength3[BSCNbTheta + i] +
361 ( xHalfLength4[BSCNbTheta + i] - xHalfLength3[BSCNbTheta + i] ) *
362 yHalfLength1[BSCNbTheta * 2 + i] / yHalfLength2[BSCNbTheta * 1 + i];
363
364 yHalfLength2[BSCNbTheta * 4 + i] = yHalfLength1[BSCNbTheta * 4 + i] = totalThickness / 2.;
365 xHalfLength4[BSCNbTheta * 4 + i] = xHalfLength4[i];
366 xHalfLength3[BSCNbTheta * 4 + i] =
367 xHalfLength4[i] - ( xHalfLength4[i] - xHalfLength3[i] ) *
368 yHalfLength2[BSCNbTheta * 4 + i] / yHalfLength2[i];
369 xHalfLength1[BSCNbTheta * 4 + i] =
370 xHalfLength4[BSCNbTheta + i] -
371 ( xHalfLength4[BSCNbTheta + i] - xHalfLength3[BSCNbTheta + i] ) *
372 yHalfLength1[BSCNbTheta * 4 + i] / yHalfLength2[BSCNbTheta * 1 + i];
373
374 delta = totalThickness / 2. + yHalfLength1[BSCNbTheta * 2 + i];
375 angle1 = atan( yHalfLength2[i] / ( xHalfLength4[i] - xHalfLength3[i] ) );
376 angle2 = atan( 2. * ( xHalfLength4[i] - xHalfLength2[i] ) * sin( angle1 ) / BSCCryLength );
377
378 yHalfLength1[BSCNbTheta * 5 + i] = yHalfLength1[BSCNbTheta * 3 + i] =
379 yHalfLength1[i] - delta;
380 yHalfLength2[BSCNbTheta * 5 + i] = yHalfLength2[BSCNbTheta * 3 + i] =
381 yHalfLength2[i] - delta;
382 xHalfLength4[BSCNbTheta * 3 + i] = xHalfLength3[BSCNbTheta * 3 + i] =
383 xHalfLength2[BSCNbTheta * 3 + i] = xHalfLength1[BSCNbTheta * 3 + i] =
384 totalThickness / cos( angle2 ) / sin( angle1 ) / 2.;
385 xHalfLength4[BSCNbTheta * 5 + i] = xHalfLength3[BSCNbTheta * 5 + i] =
386 xHalfLength2[BSCNbTheta * 5 + i] = xHalfLength1[BSCNbTheta * 5 + i] =
387 totalThickness / 2.;
388
389 zHalfLength[i] = zHalfLength[i] - totalThickness / 2.;
390 yHalfLength1[i] = yHalfLength1[i] - delta;
391 yHalfLength2[i] = yHalfLength2[i] - delta;
392 delta = totalThickness * ( 1. + 1. / cos( angle2 ) / sin( angle1 ) ) / 2.;
393 xHalfLength1[i] = xHalfLength1[i] - delta;
394 xHalfLength2[i] = xHalfLength2[i] - delta;
395 xHalfLength3[i] = xHalfLength3[i] - delta;
396 xHalfLength4[i] = xHalfLength4[i] - delta;
397
398 oop = sqrt( ( yHalfLength2[i] - yHalfLength1[i] ) * ( yHalfLength2[i] - yHalfLength1[i] ) +
399 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) *
400 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) /
401 4 );
402 thetaAxis[i] = atan( oop / BSCCryLength );
403 // -phi --->180+phi
404 phiAxis[i] =
405 180. * deg +
406 atan( ( yHalfLength2[i] - yHalfLength1[i] ) /
407 ( xHalfLength3[i] + xHalfLength4[i] - xHalfLength1[i] - xHalfLength2[i] ) * 2. );
408
409 oop = sqrt( ( yHalfLength2[BSCNbTheta + i] - yHalfLength1[BSCNbTheta + i] ) *
410 ( yHalfLength2[BSCNbTheta + i] - yHalfLength1[BSCNbTheta + i] ) +
411 ( xHalfLength3[BSCNbTheta + i] + xHalfLength4[BSCNbTheta + i] -
412 xHalfLength1[BSCNbTheta + i] - xHalfLength2[BSCNbTheta + i] ) *
413 ( xHalfLength3[BSCNbTheta + i] + xHalfLength4[BSCNbTheta + i] -
414 xHalfLength1[BSCNbTheta + i] - xHalfLength2[BSCNbTheta + i] ) /
415 4 );
416 thetaAxis[BSCNbTheta + i] = atan( oop / totalThickness );
417 phiAxis[BSCNbTheta + i] =
418 -atan( ( yHalfLength2[BSCNbTheta + i] - yHalfLength1[BSCNbTheta + i] ) /
419 ( xHalfLength3[BSCNbTheta + i] + xHalfLength4[BSCNbTheta + i] -
420 xHalfLength1[BSCNbTheta + i] - xHalfLength2[BSCNbTheta + i] ) *
421 2. );
422
423 oop = sqrt( ( yHalfLength2[i] - yHalfLength1[i] ) * ( yHalfLength2[i] - yHalfLength1[i] ) *
424 4 +
425 ( xHalfLength3[BSCNbTheta * 2 + i] + xHalfLength4[BSCNbTheta * 2 + i] -
426 xHalfLength1[BSCNbTheta * 2 + i] - xHalfLength2[BSCNbTheta * 2 + i] ) *
427 ( xHalfLength3[BSCNbTheta * 2 + i] + xHalfLength4[BSCNbTheta * 2 + i] -
428 xHalfLength1[BSCNbTheta * 2 + i] - xHalfLength2[BSCNbTheta * 2 + i] ) /
429 4 );
430 thetaAxis[BSCNbTheta * 2 + i] = atan( oop / ( zHalfLength[BSCNbTheta * 2 + i] * 2 ) );
431 phiAxis[BSCNbTheta * 2 + i] =
432 -atan( ( yHalfLength2[i] - yHalfLength1[i] ) /
433 ( xHalfLength3[BSCNbTheta * 2 + i] + xHalfLength4[BSCNbTheta * 2 + i] -
434 xHalfLength1[BSCNbTheta * 2 + i] - xHalfLength2[BSCNbTheta * 2 + i] ) *
435 4 );
436
437 oop = sqrt(
438 ( yHalfLength2[i] - yHalfLength1[i] ) * ( yHalfLength2[i] - yHalfLength1[i] ) * 4 +
439 ( xHalfLength4[i] - xHalfLength2[i] ) * ( xHalfLength4[i] - xHalfLength2[i] ) * 4 );
440 thetaAxis[BSCNbTheta * 3 + i] = atan( oop / ( zHalfLength[BSCNbTheta * 3 + i] * 2 ) );
441 phiAxis[BSCNbTheta * 3 + i] =
442 -atan( ( yHalfLength2[i] - yHalfLength1[i] ) / ( xHalfLength4[i] - xHalfLength2[i] ) );
443
444 thetaAxis[BSCNbTheta * 4 + i] =
445 atan( ( xHalfLength4[BSCNbTheta * 4 + i] + xHalfLength3[BSCNbTheta * 4 + i] -
446 xHalfLength2[BSCNbTheta * 4 + i] - xHalfLength1[BSCNbTheta * 4 + i] ) /
447 2. / ( zHalfLength[BSCNbTheta * 4 + i] * 2 ) );
448 phiAxis[BSCNbTheta * 4 + i] = 0;
449
450 thetaAxis[BSCNbTheta * 5 + i] =
451 atan( ( xHalfLength3[BSCNbTheta * 5 + i] - xHalfLength1[BSCNbTheta * 5 + i] ) /
452 ( zHalfLength[BSCNbTheta * 5 + i] * 2 ) );
453 phiAxis[BSCNbTheta * 5 + i] = -90. * deg;
454
455 tanAlpha2[BSCNbTheta + i] = tanAlpha1[BSCNbTheta + i] = tanAlpha1[i] =
456 -( xHalfLength2[i] - xHalfLength1[i] ) / yHalfLength1[i] / 2.;
457 tanAlpha1[BSCNbTheta * 2 + i] =
458 ( xHalfLength2[BSCNbTheta * 2 + i] - xHalfLength1[BSCNbTheta * 2 + i] ) /
459 yHalfLength1[BSCNbTheta * 2 + i] / 2.;
460 tanAlpha1[BSCNbTheta * 3 + i] = tanAlpha1[i] * 2.;
461 tanAlpha1[BSCNbTheta * 4 + i] =
462 ( xHalfLength2[BSCNbTheta * 4 + i] - xHalfLength1[BSCNbTheta * 4 + i] ) /
463 yHalfLength1[BSCNbTheta * 4 + i] / 2.;
464 tanAlpha1[BSCNbTheta * 5 + i] =
465 ( xHalfLength2[BSCNbTheta * 5 + i] - xHalfLength1[BSCNbTheta * 5 + i] ) /
466 yHalfLength1[BSCNbTheta * 5 + i] / 2.;
467
468 tanAlpha2[i] = -( xHalfLength4[i] - xHalfLength3[i] ) / yHalfLength2[i] / 2.;
469 // tanAlpha2[BSCNbTheta+i]=(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])/yHalfLength2[BSCNbTheta+i]/2.;
470 tanAlpha2[BSCNbTheta * 2 + i] =
471 ( xHalfLength4[BSCNbTheta * 2 + i] - xHalfLength3[BSCNbTheta * 2 + i] ) /
472 yHalfLength2[BSCNbTheta * 2 + i] / 2.;
473 tanAlpha2[BSCNbTheta * 3 + i] = tanAlpha2[i] * 2.;
474 tanAlpha2[BSCNbTheta * 4 + i] =
475 ( xHalfLength4[BSCNbTheta * 4 + i] - xHalfLength3[BSCNbTheta * 4 + i] ) /
476 yHalfLength2[BSCNbTheta * 4 + i] / 2.;
477 tanAlpha2[BSCNbTheta * 5 + i] =
478 ( xHalfLength4[BSCNbTheta * 5 + i] - xHalfLength3[BSCNbTheta * 5 + i] ) /
479 yHalfLength2[BSCNbTheta * 5 + i] / 2.;
480
481 zPosition[BSCNbTheta * 5 + i] = zPosition[BSCNbTheta * 3 + i] = zPosition[i] =
482 zPosition[i] + totalThickness / 2. * cos( thetaPosition[i] ) -
483 yHalfLength1[BSCNbTheta * 2 + i] * sin( thetaPosition[i] );
484 zPosition[i] = totalThickness / 2.; // for the newest method
485 xPosition[BSCNbTheta * 5 + i] = xPosition[BSCNbTheta * 3 + i] = xPosition[i] =
486 xPosition[i] + totalThickness / 2. * sin( thetaPosition[i] ) +
487 totalThickness * ( 1. / cos( thetaPosition[i + 1] - thetaPosition[i] ) - 1 ) / 2. *
488 cos( thetaPosition[i] );
489 xPosition[i] = totalThickness * ( 1. - 1. / cos( angle2 ) / sin( angle1 ) ) / 2.;
490 // for the newest method
491 yPosition[i] =
492 yPosition[i] + totalThickness * ( 1. - 1. / cos( angle2 ) / sin( angle1 ) ) / 2.;
493 yPosition[i] =
494 yHalfLength1[BSCNbTheta * 2 + i] - totalThickness / 2.; // for the newest method
495 yPosition[BSCNbTheta * 3 + i] = yPosition[i] * 2. + xHalfLength1[BSCNbTheta * 3 + i];
496 yPosition[BSCNbTheta * 5 + i] = xHalfLength1[BSCNbTheta * 5 + i];
497
498 xPosition[BSCNbTheta + i] =
499 BSCPhiRmin + zHalfLength[BSCNbTheta + i] * sin( thetaPosition[i] ) +
500 ( 3. * yHalfLength1[BSCNbTheta + i] - yHalfLength2[BSCNbTheta + i] ) / 2. *
501 cos( thetaPosition[i] );
502 yPosition[BSCNbTheta + i] =
503 ( xHalfLength1[BSCNbTheta + i] + xHalfLength3[BSCNbTheta + i] +
504 xHalfLength2[BSCNbTheta + i] + xHalfLength4[BSCNbTheta + i] ) /
505 4.;
506 zPosition[BSCNbTheta + i] =
507 BSCPosition1 + rminprojected / tan( thetaPosition[i] ) +
508 ( 2. * yHalfLength1[BSCNbTheta + i] / tan( thetaPosition[i] ) +
509 zHalfLength[BSCNbTheta + i] ) *
510 cos( thetaPosition[i] ) +
511 ( yHalfLength1[BSCNbTheta + i] + yHalfLength2[BSCNbTheta + i] ) / 2. *
512 sin( thetaPosition[i] );
513
514 xPosition[BSCNbTheta * 2 + i] =
515 xPosition[i] +
516 ( ( yHalfLength1[i] + yHalfLength2[i] ) / 2. + yHalfLength1[BSCNbTheta * 2 + i] ) *
517 cos( thetaPosition[i] );
518 zPosition[BSCNbTheta * 2 + i] =
519 zPosition[i] -
520 ( ( yHalfLength1[i] + yHalfLength2[i] ) / 2. + yHalfLength1[BSCNbTheta * 2 + i] ) *
521 sin( thetaPosition[i] );
522 yPosition[BSCNbTheta * 2 + i] =
523 ( xHalfLength1[BSCNbTheta * 2 + i] + xHalfLength3[BSCNbTheta * 2 + i] +
524 xHalfLength2[BSCNbTheta * 2 + i] + xHalfLength4[BSCNbTheta * 2 + i] ) /
525 4.;
526
527 xPosition[BSCNbTheta * 4 + i] =
528 xPosition[i] -
529 ( ( yHalfLength1[i] + yHalfLength2[i] ) / 2. + yHalfLength1[BSCNbTheta * 4 + i] ) *
530 cos( thetaPosition[i] );
531 zPosition[BSCNbTheta * 4 + i] =
532 zPosition[i] -
533 ( ( yHalfLength1[i] + yHalfLength2[i] ) / 2. + yHalfLength1[BSCNbTheta * 4 + i] ) *
534 sin( thetaPosition[i] );
535 yPosition[BSCNbTheta * 4 + i] =
536 ( xHalfLength1[BSCNbTheta * 4 + i] + xHalfLength3[BSCNbTheta * 4 + i] +
537 xHalfLength2[BSCNbTheta * 4 + i] + xHalfLength4[BSCNbTheta * 4 + i] ) /
538 4.;
539 }
540
541 if ( verboseLevel > 1 )
542 for ( i = 0; i < BSCNbTheta * 6; i++ )
543 {
544 G4cout << "The sizes of the " << i + 1 << " crystal are:" << G4endl
545 << "zHalfLength =" << zHalfLength[i] / cm << "(cm)," << G4endl
546 << "thetaAxis =" << thetaAxis[i] / deg << "(deg)," << G4endl
547 << "phiAxis =" << phiAxis[i] / deg << "(deg)," << G4endl
548 << "yHalfLength1=" << yHalfLength1[i] / cm << "(cm)," << G4endl
549 << "xHalfLength1=" << xHalfLength1[i] / cm << "(cm)," << G4endl
550 << "xHalfLength2=" << xHalfLength2[i] / cm << "(cm)," << G4endl
551 << "tanAlpha1 =" << tanAlpha1[i] << G4endl
552 << "yHalfLength2=" << yHalfLength2[i] / cm << "(cm)," << G4endl
553 << "xHalfLength3=" << xHalfLength3[i] / cm << "(cm)," << G4endl
554 << "xHalfLength4=" << xHalfLength4[i] / cm << "(cm)," << G4endl
555 << "tanAlpha2 =" << tanAlpha2[i] << "." << G4endl << "The position of the "
556 << i + 1 << " crystal is:" << G4endl << "(" << xPosition[i] / cm << ","
557 << yPosition[i] / cm << "," << zPosition[i] / cm << ")cm" << G4endl;
558 }
559}
560
561void ExtBesEmcGeometry::SetCasingThickness( G4ThreeVector val ) {
562 // change Gap thickness and recompute the calorimeter parameters
563 fTyvekThickness = val( 'X' );
564 fAlThickness = val( 'Y' );
565 fMylarThickness = val( 'Z' );
566}
567
568G4double ExtBesEmcGeometry::GetXPosition( G4int NbCrystal ) {
569 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return xPosition[NbCrystal]; }
570 else { return 0.; }
571}
572
573G4double ExtBesEmcGeometry::GetYPosition( G4int NbCrystal ) {
574 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return yPosition[NbCrystal]; }
575 else { return 0.; }
576}
577
578G4double ExtBesEmcGeometry::GetZPosition( G4int NbCrystal ) {
579 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return zPosition[NbCrystal]; }
580 else { return 0.; }
581}
582
583G4double ExtBesEmcGeometry::GetThetaPosition( G4int NbCrystal ) {
584 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return thetaPosition[NbCrystal]; }
585 else { return 0.; }
586}
587
588G4double ExtBesEmcGeometry::GetZHalfLength( G4int NbCrystal ) {
589 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return zHalfLength[NbCrystal]; }
590 else { return 0.; }
591}
592
593G4double ExtBesEmcGeometry::GetThetaAxis( G4int NbCrystal ) {
594 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return thetaAxis[NbCrystal]; }
595 else { return 0.; }
596}
597
598G4double ExtBesEmcGeometry::GetPhiAxis( G4int NbCrystal ) {
599 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return phiAxis[NbCrystal]; }
600 else { return 0.; }
601}
602
603G4double ExtBesEmcGeometry::GetYHalfLength1( G4int NbCrystal ) {
604 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return yHalfLength1[NbCrystal]; }
605 else { return 0.; }
606}
607
608G4double ExtBesEmcGeometry::GetXHalfLength1( G4int NbCrystal ) {
609 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return xHalfLength1[NbCrystal]; }
610 else { return 0.; }
611}
612
613G4double ExtBesEmcGeometry::GetXHalfLength2( G4int NbCrystal ) {
614 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return xHalfLength2[NbCrystal]; }
615 else { return 0.; }
616}
617
618G4double ExtBesEmcGeometry::GetTanAlpha1( G4int NbCrystal ) {
619 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return tanAlpha1[NbCrystal]; }
620 else { return 0.; }
621}
622
623G4double ExtBesEmcGeometry::GetYHalfLength2( G4int NbCrystal ) {
624 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return yHalfLength2[NbCrystal]; }
625 else { return 0.; }
626}
627
628G4double ExtBesEmcGeometry::GetXHalfLength3( G4int NbCrystal ) {
629 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return xHalfLength3[NbCrystal]; }
630 else { return 0.; }
631}
632
633G4double ExtBesEmcGeometry::GetXHalfLength4( G4int NbCrystal ) {
634 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return xHalfLength4[NbCrystal]; }
635 else { return 0.; }
636}
637
638G4double ExtBesEmcGeometry::GetTanAlpha2( G4int NbCrystal ) {
639 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return tanAlpha2[NbCrystal]; }
640 else { return 0.; }
641}
642
643G4VPhysicalVolume* ExtBesEmcGeometry::GetPhysiBSCCrystal( G4int NbCrystal ) {
644 if ( NbCrystal >= 0 && NbCrystal < BSCNbTheta * 6 ) { return physiBSCCrystal[NbCrystal]; }
645 else { return NULL; }
646}
tan
double tan(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:189
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
ExtBesEmcGeometry::GetXHalfLength1
G4double GetXHalfLength1(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:608
ExtBesEmcGeometry::GetThetaAxis
G4double GetThetaAxis(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:593
ExtBesEmcGeometry::GetTanAlpha2
G4double GetTanAlpha2(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:638
ExtBesEmcGeometry::GetThetaPosition
G4double GetThetaPosition(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:583
ExtBesEmcGeometry::GetPhiAxis
G4double GetPhiAxis(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:598
ExtBesEmcGeometry::GetTanAlpha1
G4double GetTanAlpha1(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:618
ExtBesEmcGeometry::GetYHalfLength2
G4double GetYHalfLength2(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:623
ExtBesEmcGeometry::GetYPosition
G4double GetYPosition(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:573
ExtBesEmcGeometry::GetZPosition
G4double GetZPosition(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:578
ExtBesEmcGeometry::GetXHalfLength3
G4double GetXHalfLength3(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:628
ExtBesEmcGeometry::SetCasingThickness
void SetCasingThickness(G4ThreeVector)
Definition ExtBesEmcGeometry.cxx:561
ExtBesEmcGeometry::GetPhysiBSCCrystal
G4VPhysicalVolume * GetPhysiBSCCrystal(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:643
ExtBesEmcGeometry::GetZHalfLength
G4double GetZHalfLength(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:588
ExtBesEmcGeometry::GetXHalfLength2
G4double GetXHalfLength2(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:613
ExtBesEmcGeometry::GetXHalfLength4
G4double GetXHalfLength4(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:633
ExtBesEmcGeometry::GetXPosition
G4double GetXPosition(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:568
ExtBesEmcGeometry::GetYHalfLength1
G4double GetYHalfLength1(G4int NbCrystal)
Definition ExtBesEmcGeometry.cxx:603
ExtBesEmcParameter::GetInstance
static ExtBesEmcParameter & GetInstance()
Definition ExtBesEmcParameter.cxx:29
ExtBesEmcParameter::GetTaperRingInnerLength
G4double GetTaperRingInnerLength()
Definition ExtBesEmcParameter.h:84
ExtBesEmcParameter::GetTaperRingThickness3
G4double GetTaperRingThickness3()
Definition ExtBesEmcParameter.h:82
ExtBesEmcParameter::GetRearCasingThickness
G4double GetRearCasingThickness()
Definition ExtBesEmcParameter.h:92
ExtBesEmcParameter::GetTaperRingThickness2
G4double GetTaperRingThickness2()
Definition ExtBesEmcParameter.h:81
ExtBesEmcParameter::GetTaperRingThickness1
G4double GetTaperRingThickness1()
Definition ExtBesEmcParameter.h:80
ExtBesEmcParameter::GetTaperRingOuterLength
G4double GetTaperRingOuterLength()
Definition ExtBesEmcParameter.h:85