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