BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEmcEndGeometry.cc
Go to the documentation of this file.
1#include "G4ThreeVector.hh"
2
3#include "BesEmcEndGeometry.hh"
4#include "BesEmcParameter.hh"
5#include "G4ios.hh"
6
7#include <iomanip>
8#include <fstream>
9#include <sstream>
10#include <iostream>
11#include <assert.h>
12
13
14using namespace std;
15
17{;}
18
20{;}
21
23{
25
26 WorldRmin1 = emcPara.GetWorldRmin1();
27 WorldRmax1 = emcPara.GetWorldRmax1();
28 WorldRmin2 = emcPara.GetWorldRmin2();//+5.*mm; //add 5mm to avoid warning
29 WorldRmax2 = emcPara.GetWorldRmax2();
30 WorldDz = emcPara.GetWorldDz();
31 WorldZPosition = emcPara.GetWorldZPosition();
32 CrystalLength = emcPara.GetCrystalLength();
33
34 SectorRmin1 = 489.27*mm;
35 SectorRmax1 = 880.*mm;
36 SectorRmin2 = 621.57*mm;
37 SectorRmax2 = 1113.53*mm;
38 SectorDz = WorldDz-44.*mm;
39 SectorZPosition = -18.*mm;
40
41 for(G4int i=0;i<6;i++)
42 cryNumInOneLayer[i] = emcPara.GetCryInOneLayer(i);
43 for(G4int i=0;i<5;i++)
44 pentaInOneSector[i] = emcPara.GetPentaInOneSector(i);
45
46 fTyvekThickness = emcPara.GetTyvekThickness();
47 fAlThickness = emcPara.GetAlThickness();
48 fMylarThickness = emcPara.GetMylarThickness();
49 totalThickness=fTyvekThickness+fAlThickness+fMylarThickness;
50
51 G4String ParaPath = getenv("EMCSIMROOT");
52 if(!ParaPath){
53 G4Exception("BOOST environment not set!");
54 }
55 ParaPath += "/dat/EmcEndGeometry.dat";
56
57 ifstream fin;
58 fin.open(ParaPath);
59 assert(fin);
60 for(G4int i=0;i<30;i++)
61 for (G4int j=0;j<24;j++)
62 fin>>param[i][j];
63 for(G4int i=0;i<5;i++)
64 for (G4int j=0;j<6;j++)
65 fin>>penta[i][j];
66 fin.close();
67}
68
69
71{
72
73 ////////////////////////////////////////////////////////////////////////
74 // emc endcap crystal //
75 //////////////////////////////////////////////////////////////////////////
76 // 1 //
77 // 0 __ -- //
78 // -- \ //
79 // | \ \ //
80 // | \ __ - 2 //
81 // | \ -- | //
82 // | | 3 | //
83 // | | | //
84 // | | | //
85 // | | | //
86 // | | 5 | //
87 // 4 \ | | //
88 // \ | __ - 6 //
89 // \| -- //
90 // 7 //
91 /////////////////////////////////////////////////////////////////////////
92
94
95 G4int pentaNb=0;
96 for(G4int i=0;i<30;i++)
97 {
98 for (G4int j=0;j<8;j++)
99 {
100 //use 24 parameters to construct 8 point of one crystal
101 G4ThreeVector* pPnt = new G4ThreeVector(param[i][j],param[i][j+8],param[i][j+16]-WorldZPosition-SectorZPosition);
102 fPnt[i][j] = *pPnt;
103 }
104
105 if(i==5||i==6||i==14||i==15||i==16)
106 //if(0) //no pentagonal crystal now
107 {
108 for(G4int j=0;j<8;j++)
109 fPnt[30+pentaNb][j] = fPnt[i][j];
110
111 //compute the 5th point of the pentagonal crystal
112 G4double y0,y1,y4,y5;
113 G4ThreeVector v0,v1,v4,v5;
114 y0 = penta[pentaNb][0];
115 y1 = penta[pentaNb][1];
116 y4 = penta[pentaNb][2];
117 y5 = penta[pentaNb][3];
118 v0 = (y0-fPnt[i][3].getY())*(fPnt[i][0]-fPnt[i][3])/(fPnt[i][0].getY()-fPnt[i][3].getY())
119 +fPnt[i][3];
120 v1 = (y1-fPnt[i][2].getY())*(fPnt[i][1]-fPnt[i][2])/(fPnt[i][1].getY()-fPnt[i][2].getY())
121 +fPnt[i][2];
122 v4 = (y4-fPnt[i][7].getY())*(fPnt[i][4]-fPnt[i][7])/(fPnt[i][4].getY()-fPnt[i][7].getY())
123 +fPnt[i][7];
124 v5 = (y5-fPnt[i][6].getY())*(fPnt[i][5]-fPnt[i][6])/(fPnt[i][5].getY()-fPnt[i][6].getY())
125 +fPnt[i][6];
126
127 G4double x1,x5;
128 x1 = penta[pentaNb][4];
129 x5 = penta[pentaNb][5];
130 v1 = (x1-v0.getX())*(v1-v0)/(v1.getX()-v0.getX())+v0; //v1', the fifth point
131 v5 = (x5-v4.getX())*(v5-v4)/(v5.getX()-v4.getX())+v4; //v5'
132
133 fPnt[i][0] = v0;
134 fPnt[i][1] = v1;
135 fPnt[i][4] = v4;
136 fPnt[i][5] = v5;
137
138 fPnt[30+pentaNb][2] = v1;
139 fPnt[30+pentaNb][3] = v0;
140 fPnt[30+pentaNb][6] = v5;
141 fPnt[30+pentaNb][7] = v4;
142
143 pentaNb++;
144 }
145 }
146
147 //reflect point in one sector according to new design
148 G4ThreeVector temp[35][8];
149 for(G4int i=0;i<35;i++)
150 {
151 for (G4int j=0;j<8;j++)
152 {
153 temp[i][j]=fPnt[i][j];
154 fPnt[i][j].rotateZ(157.5*deg);
155 fPnt[i][j].setX(-fPnt[i][j].getX());
156 }
157
158 // point 0<-->3, 1<-->2, 4<-->7, 5<-->6
159 for (G4int j=0;j<8;j++)
160 {
161 if(j<2)
162 {
163 G4ThreeVector v=fPnt[i][j];
164 fPnt[i][j]=fPnt[i][3-j];
165 fPnt[i][3-j]=v;
166 }
167 else if(j>=4&&j<6)
168 {
169 G4ThreeVector v=fPnt[i][j];
170 fPnt[i][j]=fPnt[i][11-j];
171 fPnt[i][11-j]=v;
172 }
173 }
174 }
175
176 //exchange sequence in the same layer
177 //Exchange(0,3); Exchange(1,2); Exchange(4,7); Exchange(5,6);
178 //Exchange(8,12); Exchange(9,11); Exchange(13,17); Exchange(14,16);
179 //Exchange(18,23); Exchange(19,22); Exchange(20,21); Exchange(24,29); Exchange(25,28); Exchange(26,27);
180
181 Exchange(0,3); Exchange(1,2); Exchange(4,7); Exchange(5,31); Exchange(6,30);
182 Exchange(8,12); Exchange(9,11); Exchange(13,17); Exchange(14,34); Exchange(15,33); Exchange(16,32);
183 Exchange(18,23); Exchange(19,22); Exchange(20,21); Exchange(24,29); Exchange(25,28); Exchange(26,27);
184
185 /*for(G4int i=0;i<35;i++)
186 {
187 G4cout<<"crystal number: "<<i<<G4endl;
188 for(G4int j=0;j<8;j++)
189 G4cout<<fPnt[i][j]<<G4endl;
190 }*/
191
192
193 //compute the 6 crystal beside the 20mm gap
194 for(G4int i=0;i<35;i++)
195 {
196 for (G4int j=0;j<8;j++)
197 {
198 G4ThreeVector pPnt1 = temp[i][j];
199 fPnt1[i][j] = pPnt1.rotateZ(67.5*deg);
200 }
201 if((i==3)||(i==7)||(i==12)||(i==17)||(i==23)||(i==29))
202 {
203 fPnt1[i][0].setX(10);
204 fPnt1[i][1].setX(10);
205 fPnt1[i][4].setX(10);
206 fPnt1[i][5].setX(10);
207
208 G4double y0 = fPnt1[i][0].getY()+10*(fPnt1[i][3].getY()-fPnt1[i][0].getY())/fPnt1[i][3].getX();
209 G4double y1 = fPnt1[i][1].getY()+10*(fPnt1[i][2].getY()-fPnt1[i][1].getY())/fPnt1[i][2].getX();
210 G4double y4 = fPnt1[i][4].getY()+10*(fPnt1[i][7].getY()-fPnt1[i][4].getY())/fPnt1[i][7].getX();
211 G4double y5 = fPnt1[i][5].getY()+10*(fPnt1[i][6].getY()-fPnt1[i][5].getY())/fPnt1[i][6].getX();
212
213 G4double z0 = fPnt1[i][0].getZ()+10*(fPnt1[i][3].getZ()-fPnt1[i][0].getZ())/fPnt1[i][3].getX();
214 G4double z1 = fPnt1[i][1].getZ()+10*(fPnt1[i][2].getZ()-fPnt1[i][1].getZ())/fPnt1[i][2].getX();
215 G4double z4 = fPnt1[i][4].getZ()+10*(fPnt1[i][7].getZ()-fPnt1[i][4].getZ())/fPnt1[i][7].getX();
216 G4double z5 = fPnt1[i][5].getZ()+10*(fPnt1[i][6].getZ()-fPnt1[i][5].getZ())/fPnt1[i][6].getX();
217
218 fPnt1[i][0].setY(y0);
219 fPnt1[i][1].setY(y1);
220 fPnt1[i][4].setY(y4);
221 fPnt1[i][5].setY(y5);
222
223 fPnt1[i][0].setZ(z0);
224 fPnt1[i][1].setZ(z1);
225 fPnt1[i][4].setZ(z4);
226 fPnt1[i][5].setZ(z5);
227 }
228 }
229}
230
231void BesEmcEndGeometry::Exchange(G4int cry1, G4int cry2)
232{
233 G4ThreeVector v;
234 for(G4int i=0;i<8;i++)
235 {
236 v = fPnt[cry1][i];
237 fPnt[cry1][i] = fPnt[cry2][i];
238 fPnt[cry2][i] = v;
239 }
240}
241
242void BesEmcEndGeometry::ExchangeSector7(G4int cry1, G4int cry2)
243{
244 G4ThreeVector v;
245 for(G4int i=0;i<8;i++)
246 {
247 v = fPnt1[cry1][i];
248 fPnt1[cry1][i] = fPnt1[cry2][i];
249 fPnt1[cry2][i] = v;
250 }
251}
252
253//reflect x axis to construct sectors on the left
255{
256 for(G4int i=0;i<35;i++)
257 {
258 for(G4int j=0;j<8;j++)
259 {
260 fPnt1[i][j].setX(-fPnt1[i][j].getX());
261 }
262
263 // point 0<-->3, 1<-->2, 4<-->7, 5<-->6
264 for (G4int j=0;j<8;j++)
265 {
266 if(j<2)
267 {
268 G4ThreeVector v=fPnt1[i][j];
269 fPnt1[i][j]=fPnt1[i][3-j];
270 fPnt1[i][3-j]=v;
271 }
272 else if(j>=4&&j<6)
273 {
274 G4ThreeVector v=fPnt1[i][j];
275 fPnt1[i][j]=fPnt1[i][11-j];
276 fPnt1[i][11-j]=v;
277 }
278 }
279
280 //change the sequence of eight points according to the requirment of IrregBox
281 // point 0<-->1, 2<-->3, 4<-->5, 6<-->7
282 /*for(G4int j=0;j<8;j++)
283 {
284 if((j%2)==0)
285 {
286 G4ThreeVector v2=fPnt1[i][j];
287 fPnt1[i][j]=fPnt1[i][j+1];
288 fPnt1[i][j+1]=v2;
289 }
290 }*/
291 }
294 ExchangeSector7(9,11); ExchangeSector7(13,17); ExchangeSector7(14,34);
295 ExchangeSector7(15,33); ExchangeSector7(16,32); ExchangeSector7(18,23);
296 ExchangeSector7(19,22); ExchangeSector7(20,21); ExchangeSector7(24,29);
297 ExchangeSector7(25,28); ExchangeSector7(26,27);
298}
299
300void BesEmcEndGeometry::Zoom(const G4ThreeVector pos[8], const G4double factor)
301{
302 G4ThreeVector center1(0,0,0);
303 G4ThreeVector center2(0,0,0);
304 for(G4int i=0;i<8;i++)
305 zoomPoint[i]=0;
306
307 for(G4int i=0;i<8;i++)
308 {
309 if(i<4) center1+=pos[i];
310 else center2+=pos[i];
311 }
312 center1/=4;
313 center2/=4;
314
315 for(G4int i=0;i<8;i++)
316 {
317 if(i<4) zoomPoint[i]=factor*pos[i]+(1-factor)*center1;
318 else zoomPoint[i]=factor*pos[i]+(1-factor)*center2;
319 }
320}
321
322//subtract the casing
323void BesEmcEndGeometry::ModifyForCasing(G4ThreeVector pos[8], G4int CryNb)
324{
325 G4ThreeVector center1=0; //the center of large surface
326 G4ThreeVector center2=0; //the center of small surface
327
328 const G4double dt=1e-5; //avoid overlap between crystal and casing
329
330 if(CryNb==5||CryNb==6||CryNb==14||CryNb==15||CryNb==16)
331 {
332 center1 = (pos[0]+pos[1])*(1-dt)/2+(pos[2]+pos[3])*dt/2;
333 center2 = (pos[4]+pos[5])*(1-dt)/2+(pos[6]+pos[7])*dt/2;
334 }
335 else if(CryNb>=30&&CryNb<35)
336 {
337 center1 = (pos[0]+pos[1])*dt/2+(pos[2]+pos[3])*(1-dt)/2;
338 center2 = (pos[4]+pos[5])*dt/2+(pos[6]+pos[7])*(1-dt)/2;
339 }
340 else
341 {
342 center1 = (pos[0]+pos[1]+pos[2]+pos[3])/4;
343 center2 = (pos[4]+pos[5]+pos[6]+pos[7])/4;
344 }
345
346 G4double r1=(pos[1]-center1).r();
347 G4double r2=(pos[2]-center1).r();
348 G4double r12=(pos[1]-pos[2]).r();
349 G4double theta=acos((r2*r2+r12*r12-r1*r1)/(2*r2*r12));
350 G4double h=r2*sin(theta); //distance from the center to the vertical side
351 G4double t1=totalThickness/h;
352
353 r1=(pos[5]-center2).r();
354 r2=(pos[6]-center2).r();
355 r12=(pos[5]-pos[6]).r();
356 theta=acos((r2*r2+r12*r12-r1*r1)/(2*r2*r12));
357 h=r2*sin(theta);
358 G4double t2=totalThickness/h;
359
360 for(G4int i=0;i<8;i++)
361 {
362 if(i<4)
363 {
364 cryPoint[i] = (1-t1)*pos[i]+t1*center1;
365 }
366 else
367 {
368 G4ThreeVector temp = (1-t2)*pos[i]+t2*center2;
369 cryPoint[i] = (1-totalThickness/CrystalLength)*temp+(totalThickness/CrystalLength)*pos[i-4];
370 }
371 //G4cout<<"casing="<<pos[i]<<"\t"<<"crystal="<<cryPoint[i]<<G4endl;
372 }
373}
374
376(const G4int partId, const G4int numTheta, const G4int numPhi)
377{
378 //ComputeParameters();
379 G4int sector=-1, cryNb=-1;
380 G4int leftFlag=0;
381 G4int downFlag=0;
382 G4int pentaFlag=0;
383 G4int flag=0;
384 G4double A1=0,a1=0,B1=0,b1=0,C1=0,c1=0,D1=0,d1=0,E1=0,e1=0; //dimension of pentagonal crystal
385 G4int m_numPhi=0;
386 G4ThreeVector position=0;
387 G4int cryInOneSector = cryNumInOneLayer[numTheta]/16; //number of crystal in one layer in one sector
388 G4ThreeVector pos[8];
389
390 if(partId==2) //west end
391 {
392 if(numPhi>=0&&numPhi<8*cryInOneSector)
393 m_numPhi=8*cryInOneSector-1-numPhi;
394 else if(numPhi>=8*cryInOneSector&&numPhi<16*cryInOneSector)
395 m_numPhi=16*cryInOneSector-1-numPhi;
396 }
397 else //east end
398 m_numPhi=numPhi;
399
400 if(numPhi>=4*cryInOneSector&&numPhi<5*cryInOneSector) //crystal in 4th sector
401 {
402 leftFlag = 1;
403 m_numPhi=8*cryInOneSector-1-numPhi;
404 }
405 else if(numPhi>=11*cryInOneSector&&numPhi<12*cryInOneSector) //crystal in 11th sector
406 {
407 leftFlag = 1;
408 downFlag = 1;
409 m_numPhi=numPhi-8*cryInOneSector;
410 }
411 if(numPhi>=12*cryInOneSector&&numPhi<13*cryInOneSector) //crystal in 12th sector
412 {
413 downFlag = 1;
414 m_numPhi=16*cryInOneSector-1-numPhi;
415 }
416
417 //translate numTheta and numPhi to sector and cryNb
418 G4int cryNbOffset = 0;
419 for(G4int i=0;i<numTheta;i++)
420 cryNbOffset += cryNumInOneLayer[i]/16;
421 if(cryInOneSector)
422 sector = m_numPhi/cryInOneSector;
423 cryNb = m_numPhi-cryInOneSector*sector+cryNbOffset;
424 sector += 3;
425 if(sector>15&&sector<=18)
426 sector -= 16;
427
428 //sector beside the gap
429 if(sector==6)
430 for(G4int i=0;i<8;i++)
431 pos[i]=fPnt1[cryNb][i];
432 else
433 for(G4int i=0;i<8;i++)
434 {
435 pos[i]=fPnt[cryNb][i];
436 pos[i].rotateZ(-67.5*deg+sector*22.5*deg);
437 }
438
439 //crystal dimension
440 Zoom(pos,0.999);
441 ModifyForCasing(zoomPoint,cryNb);
442 G4double A = (cryPoint[0]-cryPoint[3]).r();
443 G4double a = (cryPoint[4]-cryPoint[7]).r();
444 G4double B = (cryPoint[1]-cryPoint[2]).r();
445 G4double b = (cryPoint[5]-cryPoint[6]).r();
446 G4double C = (cryPoint[0]-cryPoint[1]).r();
447 G4double c = (cryPoint[4]-cryPoint[5]).r();
448 G4double D = (cryPoint[2]-cryPoint[3]).r();
449 G4double d = (cryPoint[6]-cryPoint[7]).r();
450
451 //reflect the axis according to the flag
452 for(G4int i=0;i<8;i++)
453 {
454 pos[i].setZ(pos[i].getZ()+WorldZPosition+SectorZPosition); //give the absolute z coordinate
455 if(leftFlag==1)
456 pos[i].setX(-pos[i].getX());
457 if(downFlag==1)
458 pos[i].setY(-pos[i].getY());
459 if(partId==2)
460 {
461 pos[i].setX(-pos[i].getX());
462 pos[i].setZ(-pos[i].getZ());
463 }
464 }
465
466 //compute the position
467 for(G4int j=4;j<8;j++)
468 position += pos[j];
469 position /= 4;
470
471 //compute pentagonal crystal
472 for(G4int i=0;i<5;i++)
473 {
474 if(cryNb==pentaInOneSector[i])
475 {
476 pentaFlag = 1;
477 G4ThreeVector penta[8];
478
479 //sector beside the gap
480 if(sector==6)
481 for(G4int j=0;j<8;j++)
482 penta[j]=fPnt1[30+i][j];
483 else
484 for(G4int j=0;j<8;j++)
485 {
486 penta[j]=fPnt[30+i][j];
487 penta[j].rotateZ(-67.5*deg+sector*22.5*deg);
488 }
489
490 //crystal dimension
491 ModifyForCasing(penta,30+i);
492 A1 = (cryPoint[1]-cryPoint[2]).r();
493 a1 = (cryPoint[5]-cryPoint[6]).r();
494 B1 = (cryPoint[1]-cryPoint[0]).r();
495 b1 = (cryPoint[5]-cryPoint[4]).r();
496 C1 = (cryPoint[0]-cryPoint[3]).r()+A;
497 c1 = (cryPoint[4]-cryPoint[7]).r()+a;
498 D1 = D;
499 d1 = d;
500 E1 = B;
501 e1 = b;
502
503 //reflect the axis according to the flag
504 for(G4int j=0;j<8;j++)
505 {
506 penta[j].setZ(penta[j].getZ()+WorldZPosition+SectorZPosition); //give the absolute z coordinate
507 if(leftFlag==1)
508 penta[j].setX(-penta[j].getX());
509 if(downFlag==1)
510 penta[j].setY(-penta[j].getY());
511 if(partId==2)
512 {
513 penta[j].setX(-penta[j].getX());
514 penta[j].setZ(-penta[j].getZ());
515 }
516 }
517
518 //compute the position
519 for(G4int j=4;j<8;j++)
520 {
521 if(j!=0&&j!=4)
522 position += pos[j];
523 if(j==0||j==1||j==4||j==5)
524 position += penta[j];
525 }
526 position /= 10;
527
528 flag = leftFlag+downFlag;
529 if(flag==1)
530 {
531 G4double temp1 = B1; B1 = D1; D1 = temp1;
532 temp1 = A1; A1 = E1; E1 = temp1;
533 temp1 = b1; b1 = d1; d1 = temp1;
534 temp1 = a1; a1 = e1; e1 = temp1;
535 }
536
537 break;
538 }
539 }
540
541 flag = leftFlag+downFlag+partId/2;
542 if(flag==1||flag==3)
543 {
544 G4double temp = C; C = D; D = temp;
545 temp = c; c = d; d = temp;
546 }
547
548 /*G4cout<<"##########################################################################"<<G4endl;
549 G4cout<<"###sector="<<sector<<",cryNb="<<cryNb<<",left flag="<<leftFlag<<",down flag="<<downFlag<<G4endl;
550 G4cout<<"partId="<<partId<<"\t"<<"theta number="<<numTheta<<"\t"<<"phi number="<<numPhi<<G4endl;
551 G4cout<<"crystal point:"<<G4endl;
552 for(G4int i=0;i<8;i++)
553 G4cout<<pos[i]<<G4endl;
554 G4cout<<"crystal position:"<<"\t"<<position<<"\t"<<"phi="<<position.phi()/deg<<G4endl;
555 G4cout<<"crystal dimension:"<<G4endl;
556 if(pentaFlag==1)
557 G4cout<<"A="<<A1<<",a="<<a1<<",B="<<B1<<",b="<<b1<<",C="<<C1<<",c="<<c1<<",D="<<D1<<",d="<<d1<<",E="<<E1<<",e="<<e1<<G4endl;
558 else
559 G4cout<<"A="<<A<<",a="<<a<<",B="<<B<<",b="<<b<<",C="<<C<<",c="<<c<<",D="<<D<<",d="<<d<<G4endl;
560 G4cout<<"##########################################################################"<<G4endl;*/
561
562 return position;
563}
Double_t e1
double sin(const BesAngle a)
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
TGraph2DErrors * dt
Definition: McCor.cxx:45
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition: RRes.h:29
void Exchange(G4int cry1, G4int cry2)
void ExchangeSector7(G4int cry1, G4int cry2)
void ModifyForCasing(G4ThreeVector pos[8], G4int CryNb)
void Zoom(const G4ThreeVector pos[8], const G4double factor)
G4ThreeVector ComputeDimAndPos(const G4int, const G4int, const G4int)
static BesEmcParameter & GetInstance()
const double b
Definition: slope.cxx:9