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