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