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