BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRecEndCapGeo.cxx
Go to the documentation of this file.
1//
2// EmcRecEndCapGeo
3//
4// Dec 18, 2003, Created by Wang.Zhe
5//
6// unit: mm, radian
7//
8#include "EmcRecGeoSvc/EmcRecEndCapGeo.h"
9
10EmcRecEndCapGeo::EmcRecEndCapGeo() {
11 ParameterInitialize();
12 CalculateStandardSector1();
13 CalculateStandardSector2();
14 FillCCenterVector();
15}
16
17EmcRecEndCapGeo::~EmcRecEndCapGeo() {}
18
19void EmcRecEndCapGeo::ParameterInitialize() {
20 // EndCap is too complex. More careful and detailed initialization
21 // will be done to make the later calculation easier.
22 int i;
23
24 // Read constant for Database
25 // general
26 flength = 280;
27 fzshift = 100;
28 fz = 1285;
29 fr[0] = 880;
30 fr[1] = 813.60278843;
31 fr[2] = 748.39248406;
32 fr[3] = 684.54540341;
33 fr[4] = 621.94216056;
34 fr[5] = 560.46548641;
35 fr[6] = 500;
36
37 // for Ring5 and 4
38 fphi5[0] = 0;
39 fphi5[1] = 4.4;
40 fphi5[2] = 8.02;
41 fphi5[3] = 11.64;
42 fphi5[4] = 15.26;
43 fphi5[5] = 18.88;
44 fphi5[6] = 22.50;
45 for ( i = 0; i < 7; ++i ) { fphi5[i] = fphi5[i] * 3.14159265358979 / 180; }
46
47 // for Ring3 and 2
48 fphi3[0] = 0;
49 fphi3[1] = 5.124;
50 fphi3[2] = 9.468;
51 fphi3[3] = 13.812;
52 fphi3[4] = 18.156;
53 fphi3[5] = 22.50;
54 for ( i = 0; i < 6; ++i ) { fphi3[i] = fphi3[i] * 3.14159265358979 / 180; }
55
56 /// for Ring1 and 0
57 fphi1[0] = 0;
58 fphi1[1] = 6.21;
59 fphi1[2] = 11.64;
60 fphi1[3] = 17.07;
61 fphi1[4] = 22.50;
62 for ( i = 0; i < 5; ++i ) { fphi1[i] = fphi1[i] * 3.14159265358979 / 180; }
63}
64
65void EmcRecEndCapGeo::CalculateStandardSector1() {
66 int i, j;
67 EmcRecGeoPlane pl1, pl2, pl3;
68 HepPoint3D po0, po1, po2, po3, po4, po5, po6, po7, po8, po9;
69
70 double tantheta1, costheta1, sintheta1;
71
72 HepPoint3D O, n1; // very important
73 O.setX( 0 );
74 O.setY( 0 );
75 O.setZ( 0 );
76 n1.setX( 0 );
77 n1.setY( 0 );
78 n1.setZ( fzshift );
79
80 // ******** Ring 5 ********
81 HepPoint3D t, n2, n;
82 HepPoint3D w, m;
83 double dphi1 = fphi5[1] - fphi5[0];
84 double dphi2 = fphi5[2] - fphi5[1];
85 double dphi;
86 double l, ll, lp;
87
88 t.setX( fr[0] );
89 t.setY( 0 );
90 t.setZ( fzshift + fz );
91
92 tantheta1 = fr[0] / fz;
93 costheta1 = 1 / sqrt( 1 + tantheta1 * tantheta1 );
94 sintheta1 = costheta1 * tantheta1;
95
96 for ( i = 0; i < 6; ++i ) { fRing5[i].Type( EmcRecCrystal::SixPlane() ); }
97
98 // No. 0,1'
99 for ( i = 0; i < 2; ++i )
100 {
101 if ( i == 0 ) { dphi = dphi1; }
102 else { dphi = dphi2; }
103 // first, point 3
104 po3 = t;
105 fRing5[i].Set( 3, po3 );
106 // then, point 0
107 po0 = t;
108 fRing5[i].Set( 0, po0.rotateZ( dphi ) );
109 // point 7 and point 4
110 l = fRing5[i].Get( 3 ).distance( fRing5[i].Get( 0 ) ) / 2;
111 ll = sqrt( fz * fz + fr[0] * fr[0] );
112 lp = ll * flength / sqrt( ll * ll - l * l );
113 po7 = t;
114 po7.setX( po7.x() + lp * sintheta1 );
115 po7.setZ( po7.z() + lp * costheta1 );
116 fRing5[i].Set( 7, po7 );
117 po4 = po7;
118 fRing5[i].Set( 4, po4.rotateZ( dphi ) );
119 // point 1,2
120 pl1.Build( 0, 1, 0, 0 );
121 n2 = ( fRing5[i].Get( 0 ) + fRing5[i].Get( 3 ) ) / 2;
122 n = n2 - n1;
123 pl2.Build( n, fRing5[i].Get( 3 ) );
124 m.setX( fr[1] );
125 m.setY( 0 );
126 m.setZ( fzshift + fz );
127 w = m;
128 w.rotateZ( dphi );
129 pl3.Build( n1, m, w );
130 FindInt( pl1, pl2, pl3, po2 );
131
132 fRing5[i].Set( 2, po2 );
133 po1 = po2;
134 fRing5[i].Set( 1, po1.rotateZ( dphi ) );
135 // point 6,5
136 pl2.Build( n, fRing5[i].Get( 7 ) );
137 FindInt( pl1, pl2, pl3, po6 );
138 fRing5[i].Set( 6, po6 );
139 po5 = po6;
140 fRing5[i].Set( 5, po5.rotateZ( dphi ) );
141 }
142
143 // No. 1
144 for ( i = 0; i < 8; ++i ) { fRing5[1].Set( i, fRing5[1].Get( i ).rotateZ( dphi1 ) ); }
145
146 //===== No. 2--5 =====
147 for ( i = 2; i < 6; ++i )
148 {
149 for ( j = 0; j < 8; ++j )
150 { fRing5[i].Set( j, fRing5[1].Get( j ).rotateZ( fphi5[i] - fphi5[1] ) ); }
151 }
152 // Finally, ring 5 is done.
153 // Check for ring 5
154 // for(i=0;i<6;++i) {
155 // cout<<fRing5[i]<<endl;
156 // fRing5[i].EndCapCheckout();
157 // }
158
159 // ******** Ring 4 ********
160 EmcRecGeoPlane pl4, pl5, pl6;
161 double dphip;
162 HepPoint3D ttt;
163
164 t.setX( fr[1] );
165 t.setY( 0 );
166 t.setZ( fzshift + fz );
167
168 dphi1 = fphi5[1] - fphi5[0];
169 dphi2 = fphi5[2] - fphi5[1];
170
171 tantheta1 = fr[1] / fz;
172 costheta1 = 1 / sqrt( 1 + tantheta1 * tantheta1 );
173 sintheta1 = costheta1 * tantheta1;
174
175 for ( i = 0; i < 6; ++i ) { fRing4[i].Type( EmcRecCrystal::SixPlane() ); }
176
177 // It's too complicated. Boring!
178 // up and down for crystal 0,1 and 4,5
179 EmcRecCrystal up[2], down[2];
180 for ( i = 0; i < 2; ++i )
181 {
182 if ( i == 0 ) { dphi = dphi1; }
183 else { dphi = dphi2; }
184 // first point 3,0, up still needs a rotation
185 po3 = t;
186 down[i].Set( 3, po3 );
187 po0 = t;
188 down[i].Set( 0, po0.rotateZ( dphi ) );
189 //
190 po3 = t;
191 up[i].Set( 3, po3 );
192 po0 = t;
193 up[i].Set( 0, po0.rotateZ( dphi2 ) );
194 // then point 7,4, up still needs a rotation
195 l = down[i].Get( 3 ).distance( down[i].Get( 0 ) ) / 2;
196 ll = sqrt( fz * fz + fr[1] * fr[1] );
197 lp = ll * flength / sqrt( ll * ll - l * l );
198 po7 = t;
199 po7.setX( po7.x() + lp * sintheta1 );
200 po7.setZ( po7.z() + lp * costheta1 );
201 down[i].Set( 7, po7 );
202 po4 = po7;
203 down[i].Set( 4, po4.rotateZ( dphi ) );
204
205 l = up[i].Get( 3 ).distance( up[i].Get( 0 ) ) / 2;
206 ll = sqrt( fz * fz + fr[1] * fr[1] );
207 lp = ll * flength / sqrt( ll * ll - l * l );
208 po7 = t;
209 po7.setX( po7.x() + lp * sintheta1 );
210 po7.setZ( po7.z() + lp * costheta1 );
211 up[i].Set( 7, po7 );
212 po4 = po7;
213 up[i].Set( 4, po4.rotateZ( dphi2 ) );
214
215 up[i].Set( 0, up[i].Get( 0 ).rotateZ( dphi ) );
216 up[i].Set( 3, up[i].Get( 3 ).rotateZ( dphi ) );
217 up[i].Set( 4, up[i].Get( 4 ).rotateZ( dphi ) );
218 up[i].Set( 7, up[i].Get( 7 ).rotateZ( dphi ) );
219 // 0,3,4,7 is done.
220 // switch to 1,2,5,6
221 // for down
222 pl1.Build( 0, 1, 0, 0 );
223 n2 = ( down[i].Get( 0 ) + down[i].Get( 3 ) ) / 2;
224 n = n2 - n1;
225 pl2.Build( n, down[i].Get( 3 ) );
226 m.setX( fr[2] );
227 m.setY( 0 );
228 m.setZ( fzshift + fz );
229 w = m;
230 if ( i == 0 ) { dphip = fphi5[2] - fphi5[0]; }
231 else { dphip = fphi5[6] - fphi5[4]; }
232 w.rotateZ( dphip );
233 pl3.Build( n1, m, w ); // very important
234 FindInt( pl1, pl2, pl3, po2 );
235 down[i].Set( 2, po2 );
236 //
237 pl4.Build( O, n1, down[i].Get( 0 ) );
238 FindInt( pl4, pl2, pl3, po1 );
239 down[i].Set( 1, po1 );
240 //
241 // point 6,5
242 pl2.Build( n, down[i].Get( 7 ) );
243 FindInt( pl1, pl2, pl3, po6 );
244 down[i].Set( 6, po6 );
245 //
246 FindInt( pl4, pl2, pl3, po5 );
247 down[i].Set( 5, po5 );
248 // for up
249 // point 2, 1
250 n2 = ( up[i].Get( 0 ) + up[i].Get( 3 ) ) / 2;
251 n = n2 - n1;
252 pl1.Build( O, n1, up[i].Get( 3 ) );
253 pl2.Build( n, up[i].Get( 3 ) );
254 FindInt( pl1, pl2, pl3, po2 );
255 up[i].Set( 2, po2 );
256 //
257 pl4.Build( O, n1, up[i].Get( 0 ) );
258 FindInt( pl4, pl2, pl3, po1 );
259 up[i].Set( 1, po1 );
260 // point 5, 6
261 pl2.Build( n, up[i].Get( 7 ) );
262 FindInt( pl1, pl2, pl3, po6 );
263 up[i].Set( 6, po6 );
264 //
265 FindInt( pl4, pl2, pl3, po5 );
266 up[i].Set( 5, po5 );
267 }
268
269 for ( i = 0; i < 8; ++i )
270 {
271 fRing4[0].Set( i, down[0].Get( i ) );
272 fRing4[1].Set( i, up[0].Get( i ) );
273 fRing4[4].Set( i, down[1].Get( i ).rotateZ( fphi5[4] ) );
274 fRing4[5].Set( i, up[1].Get( i ).rotateZ( fphi5[4] ) );
275 }
276
277 // crystal 2 and 3
278 dphi = fphi5[3] - fphi5[2];
279 // first, point 3
280 po3 = t;
281 fRing4[2].Set( 3, po3 );
282 // then, point 0
283 po0 = t;
284 fRing4[2].Set( 0, po0.rotateZ( dphi ) );
285 // point 7 and point 4
286 l = fRing4[2].Get( 3 ).distance( fRing4[2].Get( 0 ) ) / 2;
287 ll = sqrt( fz * fz + fr[1] * fr[1] );
288 lp = ll * flength / sqrt( ll * ll - l * l );
289 po7 = t;
290 po7.setX( po7.x() + lp * sintheta1 );
291 po7.setZ( po7.z() + lp * costheta1 );
292 fRing4[2].Set( 7, po7 );
293 po4 = po7;
294 fRing4[2].Set( 4, po4.rotateZ( dphi ) );
295 // point 2, 1
296 pl1.Build( 0, 1, 0, 0 );
297 n2 = ( fRing4[2].Get( 0 ) + fRing4[2].Get( 3 ) ) / 2;
298 n = n2 - n1;
299 pl2.Build( n, fRing4[2].Get( 3 ) );
300 m.setX( fr[2] );
301 m.setY( 0 );
302 m.setZ( fzshift + fz );
303 w = m;
304 w.rotateZ( dphi );
305 pl3.Build( n1, m, w );
306 FindInt( pl1, pl2, pl3, po2 );
307 fRing4[2].Set( 2, po2 );
308 po1 = po2;
309 fRing4[2].Set( 1, po1.rotateZ( dphi ) );
310 // point 6,5
311 pl2.Build( n, fRing4[2].Get( 7 ) );
312 FindInt( pl1, pl2, pl3, po6 );
313 fRing4[2].Set( 6, po6 );
314 po5 = po6;
315 fRing4[2].Set( 5, po5.rotateZ( dphi ) );
316
317 // Crystal 2, 3 still need a rotation.
318 for ( i = 0; i < 8; ++i ) { fRing4[2].Set( i, fRing4[2].Get( i ).rotateZ( fphi5[2] ) ); }
319 for ( i = 0; i < 8; ++i )
320 { fRing4[3].Set( i, fRing4[2].Get( i ).rotateZ( fphi5[3] - fphi5[2] ) ); }
321
322 // Finally done. Check it out.
323 // for(i=0;i<6;++i) {
324 // cout<<fRing4[i]<<endl;
325 // fRing4[i].EndCapCheckout();
326 // }
327
328 // ***************** Ring3 ********************
329 // Here I changed my way of calculation.
330 // Don't care for it. Still BORING!
331 // do some preparation
332 HepPoint3D base3[5];
333 base3[0].setX( fr[2] );
334 base3[0].setY( 0 );
335 base3[0].setZ( fz + fzshift );
336 for ( i = 1; i < 5; ++i ) { base3[i] = base3[0]; }
337 base3[1].rotateZ( fphi5[2] );
338 base3[2].rotateZ( fphi5[3] );
339 base3[3].rotateZ( fphi5[4] );
340 base3[4].rotateZ( fphi5[6] );
341 // for(i=0;i<5;++i) {
342 // cout<<base3[i]<<endl;
343 // }
344
345 HepPoint3D base2[6];
346 base2[0].setX( fr[3] );
347 base2[0].setY( 0 );
348 base2[0].setZ( fz + fzshift );
349 for ( i = 1; i < 6; ++i )
350 {
351 base2[i] = base2[0];
352 base2[i].rotateZ( fphi3[i] );
353 }
354 // for(i=0;i<6;++i) {
355 // cout<<base2[i]<<endl;
356 // }
357 // cout<<endl;
358
359 EmcRecGeoPlane pphi[6];
360 for ( i = 0; i < 6; ++i ) { pphi[i].Build( O, n1, base2[i] ); }
361 EmcRecGeoPlane ptheta[4];
362 for ( i = 0; i < 4; ++i ) { ptheta[i].Build( n1, base3[i], base3[i + 1] ); }
363 EmcRecGeoPlane pthetap[5];
364 for ( i = 0; i < 5; ++i ) { pthetap[i].Build( n1, base2[i], base2[i + 1] ); }
365
366 // Once an error occor here. I just declare HepPoint3D nn[4];
367 // Finally, the operation of nn[4] gets out of range.
368 // It has overlaped another varible.
369 HepPoint3D nn[5];
370 nn[0] = ( base3[0] + base3[1] ) / 2 - n1;
371 nn[1] = base3[1] - n1;
372 nn[2] = base3[2] - n1;
373 nn[3] = base3[3] - n1;
374 nn[4] = ( base3[3] + base3[4] ) / 2 - n1;
375
376 EmcRecGeoPlane psection[5];
377 for ( i = 0; i < 5; ++i ) { psection[i].Build( nn[i], base3[i] ); }
378
379 EmcRecGeoPlane psection2[5];
380 HepPoint3D bp[5];
381 for ( i = 0; i < 5; ++i )
382 {
383 bp[i] = base3[i];
384 bp[i].setX( bp[i].x() + flength * nn[i].x() / nn[i].mag() );
385 bp[i].setY( bp[i].y() + flength * nn[i].y() / nn[i].mag() );
386 bp[i].setZ( bp[i].z() + flength * nn[i].z() / nn[i].mag() );
387 psection2[i].Build( nn[i], bp[i] );
388 }
389
390 EmcRecGeoPlane pthetatmp;
391 for ( i = 0; i < 5; ++i )
392 {
393 if ( i == 0 || i == 4 )
394 {
395 fRing3[i].Type( EmcRecCrystal::SixPlane() );
396 if ( i == 0 ) { pthetatmp = ptheta[0]; }
397 if ( i == 4 ) { pthetatmp = ptheta[3]; }
398 FindInt( pphi[i], pthetatmp, psection[i], po3 );
399 FindInt( pphi[i], pthetap[i], psection[i], po2 );
400 FindInt( pphi[i + 1], pthetatmp, psection[i], po0 );
401 FindInt( pphi[i + 1], pthetap[i], psection[i], po1 );
402 fRing3[i].Set( 0, po0 );
403 fRing3[i].Set( 1, po1 );
404 fRing3[i].Set( 2, po2 );
405 fRing3[i].Set( 3, po3 );
406 FindInt( pphi[i], pthetatmp, psection2[i], po7 );
407 FindInt( pphi[i], pthetap[i], psection2[i], po6 );
408 FindInt( pphi[i + 1], pthetatmp, psection2[i], po4 );
409 FindInt( pphi[i + 1], pthetap[i], psection2[i], po5 );
410 fRing3[i].Set( 4, po4 );
411 fRing3[i].Set( 5, po5 );
412 fRing3[i].Set( 6, po6 );
413 fRing3[i].Set( 7, po7 );
414 }
415 else
416 {
417 fRing3[i].Type( EmcRecCrystal::SevenPlane() );
418 po0 = base3[i];
419 FindInt( pphi[i], ptheta[i - 1], psection[i], po4 );
420 FindInt( pphi[i], pthetap[i], psection[i], po3 );
421 FindInt( pphi[i + 1], ptheta[i], psection[i], po1 );
422 FindInt( pphi[i + 1], pthetap[i], psection[i], po2 );
423 fRing3[i].Set( 0, po0 );
424 fRing3[i].Set( 1, po1 );
425 fRing3[i].Set( 2, po2 );
426 fRing3[i].Set( 3, po3 );
427 fRing3[i].Set( 4, po4 );
428 po5 = bp[i];
429 FindInt( pphi[i], ptheta[i - 1], psection2[i], po9 );
430 FindInt( pphi[i], pthetap[i], psection2[i], po8 );
431 FindInt( pphi[i + 1], ptheta[i], psection2[i], po6 );
432 FindInt( pphi[i + 1], pthetap[i], psection2[i], po7 );
433 fRing3[i].Set( 5, po5 );
434 fRing3[i].Set( 6, po6 );
435 fRing3[i].Set( 7, po7 );
436 fRing3[i].Set( 8, po8 );
437 fRing3[i].Set( 9, po9 );
438 }
439 }
440 /// Finally, Ring 3 is done.
441 /// The maximum descripancy is less than 1mm.
442 /// However, it's hard to know where they from.
443 /// There are too many approximation in structure design.
444 /// It's already less than the tolerance of production.
445 /// Then, omit them!
446 // for(i=0;i<5;++i) {
447 // cout<<fRing3[i]<<endl;
448 // fRing3[i].EndCapCheckout();
449 // }
450
451 /// Get into Ring2
452 HepPoint3D base1[4];
453 base1[0].setX( fr[4] );
454 base1[0].setY( 0 );
455 base1[0].setZ( fz + fzshift );
456 for ( i = 1; i < 4; ++i ) { base1[i] = base1[0]; }
457 base1[1].rotateZ( fphi3[2] );
458 base1[2].rotateZ( fphi3[3] );
459 base1[3].rotateZ( fphi3[5] );
460 EmcRecGeoPlane ptheta1[3];
461 for ( i = 0; i < 3; ++i ) { ptheta1[i].Build( n1, base1[i], base1[i + 1] ); }
462
463 HepPoint3D nn2[5];
464 for ( i = 0; i < 5; ++i ) { nn2[i] = ( base2[i] + base2[i + 1] ) / 2 - n1; }
465 EmcRecGeoPlane psec2[5];
466 for ( i = 0; i < 5; ++i ) { psec2[i].Build( nn2[i], base2[i] ); }
467
468 EmcRecGeoPlane psec2p[5];
469 HepPoint3D bpp[5];
470 for ( i = 0; i < 5; ++i )
471 {
472 bpp[i] = base2[i];
473 bpp[i].setX( bpp[i].x() + flength * nn2[i].x() / nn2[i].mag() );
474 bpp[i].setY( bpp[i].y() + flength * nn2[i].y() / nn2[i].mag() );
475 bpp[i].setZ( bpp[i].z() + flength * nn2[i].z() / nn2[i].mag() );
476 psec2p[i].Build( nn2[i], bpp[i] );
477 }
478
479 EmcRecGeoPlane ptheta1tmp;
480 for ( i = 0; i < 5; ++i )
481 {
482 fRing2[i].Type( EmcRecCrystal::SixPlane() );
483 if ( i < 2 ) { ptheta1tmp = ptheta1[0]; }
484 if ( i == 2 ) { ptheta1tmp = ptheta1[1]; }
485 if ( i > 2 ) { ptheta1tmp = ptheta1[2]; }
486 po3 = base2[i];
487 FindInt( pphi[i], ptheta1tmp, psec2[i], po2 );
488 po0 = base2[i + 1];
489 FindInt( pphi[i + 1], ptheta1tmp, psec2[i], po1 );
490 FindInt( pphi[i], pthetap[i], psec2p[i], po7 );
491 FindInt( pphi[i], ptheta1tmp, psec2p[i], po6 );
492 FindInt( pphi[i + 1], pthetap[i], psec2p[i], po4 );
493 FindInt( pphi[i + 1], ptheta1tmp, psec2p[i], po5 );
494 fRing2[i].Set( 0, po0 );
495 fRing2[i].Set( 1, po1 );
496 fRing2[i].Set( 2, po2 );
497 fRing2[i].Set( 3, po3 );
498 fRing2[i].Set( 4, po4 );
499 fRing2[i].Set( 5, po5 );
500 fRing2[i].Set( 6, po6 );
501 fRing2[i].Set( 7, po7 );
502 }
503
504 /// Ring2 is done. It seems a little easier.
505 // for(i=0;i<5;++i) {
506 // cout<<fRing2[i]<<endl;
507 // fRing2[i].EndCapCheckout();
508 // }
509
510 /// ************* Ring1 *************
511 HepPoint3D base0[5];
512 base0[0].setX( fr[5] );
513 base0[0].setY( 0 );
514 base0[0].setZ( fz + fzshift );
515 for ( i = 1; i < 5; ++i )
516 {
517 base0[i] = base0[0];
518 base0[i].rotateZ( fphi1[i] );
519 }
520
521 EmcRecGeoPlane pphi1[5];
522 for ( i = 0; i < 5; ++i ) { pphi1[i].Build( O, n1, base0[i] ); }
523
524 EmcRecGeoPlane ptheta0[4];
525 for ( i = 0; i < 4; ++i ) { ptheta0[i].Build( n1, base0[i], base0[i + 1] ); }
526
527 HepPoint3D nn1[4];
528 nn1[0] = ( base1[0] + base1[1] ) / 2 - n1;
529 nn1[1] = base1[1] - n1;
530 nn1[2] = base1[2] - n1;
531 nn1[3] = ( base1[2] + base1[3] ) / 2 - n1;
532
533 EmcRecGeoPlane psec1[4];
534 for ( i = 0; i < 4; ++i ) { psec1[i].Build( nn1[i], base1[i] ); }
535
536 EmcRecGeoPlane psec1p[4];
537 HepPoint3D qq[4];
538 for ( i = 0; i < 4; ++i )
539 {
540 qq[i] = base1[i];
541 qq[i].setX( qq[i].x() + flength * nn1[i].x() / nn1[i].mag() );
542 qq[i].setY( qq[i].y() + flength * nn1[i].y() / nn1[i].mag() );
543 qq[i].setZ( qq[i].z() + flength * nn1[i].z() / nn1[i].mag() );
544 psec1p[i].Build( nn1[i], qq[i] );
545 }
546
547 EmcRecGeoPlane pt1tmp;
548 for ( i = 0; i < 4; ++i )
549 {
550 if ( i == 0 || i == 3 )
551 {
552 fRing1[i].Type( EmcRecCrystal::SixPlane() );
553 if ( i == 0 ) { pt1tmp = ptheta1[0]; }
554 else { pt1tmp = ptheta1[2]; }
555 FindInt( pphi1[i], pt1tmp, psec1[i], po3 );
556 FindInt( pphi1[i], ptheta0[i], psec1[i], po2 );
557 FindInt( pphi1[i + 1], pt1tmp, psec1[i], po0 );
558 FindInt( pphi1[i + 1], ptheta0[i], psec1[i], po1 );
559 FindInt( pphi1[i], pt1tmp, psec1p[i], po7 );
560 FindInt( pphi1[i], ptheta0[i], psec1p[i], po6 );
561 FindInt( pphi1[i + 1], pt1tmp, psec1p[i], po4 );
562 FindInt( pphi1[i + 1], ptheta0[i], psec1p[i], po5 );
563 fRing1[i].Set( 0, po0 );
564 fRing1[i].Set( 1, po1 );
565 fRing1[i].Set( 2, po2 );
566 fRing1[i].Set( 3, po3 );
567 fRing1[i].Set( 4, po4 );
568 fRing1[i].Set( 5, po5 );
569 fRing1[i].Set( 6, po6 );
570 fRing1[i].Set( 7, po7 );
571 }
572 else
573 {
574 fRing1[i].Type( EmcRecCrystal::SevenPlane() );
575 po0 = base1[i];
576 FindInt( pphi1[i], ptheta1[i - 1], psec1[i], po4 );
577 FindInt( pphi1[i], ptheta0[i], psec1[i], po3 );
578 FindInt( pphi1[i + 1], ptheta1[i], psec1[i], po1 );
579 FindInt( pphi1[i + 1], ptheta0[i], psec1[i], po2 );
580 fRing1[i].Set( 0, po0 );
581 fRing1[i].Set( 1, po1 );
582 fRing1[i].Set( 2, po2 );
583 fRing1[i].Set( 3, po3 );
584 fRing1[i].Set( 4, po4 );
585 po5 = qq[i];
586 FindInt( pphi1[i], ptheta1[i - 1], psec1p[i], po9 );
587 FindInt( pphi1[i], ptheta0[i], psec1p[i], po8 );
588 FindInt( pphi1[i + 1], ptheta1[i], psec1p[i], po6 );
589 FindInt( pphi1[i + 1], ptheta0[i], psec1p[i], po7 );
590 fRing1[i].Set( 5, po5 );
591 fRing1[i].Set( 6, po6 );
592 fRing1[i].Set( 7, po7 );
593 fRing1[i].Set( 8, po8 );
594 fRing1[i].Set( 9, po9 );
595 }
596 }
597
598 /// Check it out!
599 // for(i=0;i<4;++i) {
600 // cout<<fRing1[i]<<endl;
601 // fRing1[i].EndCapCheckout();
602 // }
603 /// Last ring! Ring0
604 HepPoint3D basem1[5];
605 basem1[0].setX( fr[6] );
606 basem1[0].setY( 0 );
607 basem1[0].setZ( fz + fzshift );
608 for ( i = 1; i < 5; ++i )
609 {
610 basem1[i] = basem1[0];
611 basem1[i].rotateZ( fphi1[i] );
612 }
613
614 EmcRecGeoPlane pthetam1[4];
615 for ( i = 0; i < 4; ++i ) { pthetam1[i].Build( n1, basem1[i], basem1[i + 1] ); }
616
617 HepPoint3D nn0[4];
618 for ( i = 0; i < 4; ++i ) { nn0[i] = ( base0[i] + base0[i + 1] ) / 2 - n1; }
619
620 EmcRecGeoPlane psec0[4];
621 for ( i = 0; i < 4; ++i ) { psec0[i].Build( nn0[i], base0[i] ); }
622
623 EmcRecGeoPlane psec0p[4];
624 HepPoint3D qq0[4];
625 for ( i = 0; i < 4; ++i )
626 {
627 qq0[i] = base0[i];
628 qq0[i].setX( qq0[i].x() + flength * nn0[i].x() / nn0[i].mag() );
629 qq0[i].setY( qq0[i].y() + flength * nn0[i].y() / nn0[i].mag() );
630 qq0[i].setZ( qq0[i].z() + flength * nn0[i].z() / nn0[i].mag() );
631 psec0p[i].Build( nn0[i], qq0[i] );
632 }
633
634 for ( i = 0; i < 4; ++i )
635 {
636 fRing0[i].Type( EmcRecCrystal::SixPlane() );
637 po3 = base0[i];
638 FindInt( pphi1[i], pthetam1[i], psec0[i], po2 );
639 po0 = base0[i + 1];
640 FindInt( pphi1[i + 1], pthetam1[i], psec0[i], po1 );
641 FindInt( pphi1[i], ptheta0[i], psec0p[i], po7 );
642 FindInt( pphi1[i], pthetam1[i], psec0p[i], po6 );
643 FindInt( pphi1[i + 1], ptheta0[i], psec0p[i], po4 );
644 FindInt( pphi1[i + 1], pthetam1[i], psec0p[i], po5 );
645 fRing0[i].Set( 0, po0 );
646 fRing0[i].Set( 1, po1 );
647 fRing0[i].Set( 2, po2 );
648 fRing0[i].Set( 3, po3 );
649 fRing0[i].Set( 4, po4 );
650 fRing0[i].Set( 5, po5 );
651 fRing0[i].Set( 6, po6 );
652 fRing0[i].Set( 7, po7 );
653 }
654
655 /// Sector1 is finished. !!!
656 // for(i=0;i<4;++i) {
657 // cout<<fRing0[i]<<endl;
658 // fRing0[i].EndCapCheckout();
659 // }
660 for ( i = 0; i < 6; ++i )
661 {
662 if ( i < 4 )
663 {
664 fCrystal[0][0][i] = fRing0[i];
665 fCrystal[0][1][i] = fRing1[i];
666 }
667 if ( i < 5 )
668 {
669 fCrystal[0][2][i] = fRing2[i];
670 fCrystal[0][3][i] = fRing3[i];
671 }
672 fCrystal[0][4][i] = fRing4[i];
673 fCrystal[0][5][i] = fRing5[i];
674 }
675}
676
677void EmcRecEndCapGeo::CalculateStandardSector2() {
678 EmcRecCrystal edge[6];
679
680 int i, j;
681
682 for ( i = 1; i < 6; ++i )
683 {
684 fRing5p[i] = fRing5[i];
685 fRing4p[i] = fRing4[i];
686 }
687 for ( i = 1; i < 5; ++i )
688 {
689 fRing3p[i] = fRing3[i];
690 fRing2p[i] = fRing2[i];
691 }
692 for ( i = 1; i < 4; ++i )
693 {
694 fRing1p[i] = fRing1[i];
695 fRing0p[i] = fRing0[i];
696 }
697
698 edge[5] = fRing5[0];
699 edge[4] = fRing4[0];
700 edge[3] = fRing3[0];
701 edge[2] = fRing2[0];
702 edge[1] = fRing1[0];
703 edge[0] = fRing0[0];
704
705 EmcRecGeoPlane p10mm;
706 p10mm.Build( 0, 1, 0, -10 );
707
708 EmcRecGeoPlane sec1, sec2;
709 EmcRecGeoPlane theta1, theta2;
710
711 HepPoint3D po3, po2, po7, po6;
712
713 for ( i = 0; i < 6; ++i )
714 {
715 sec1.Build( edge[i].Get( 0 ), edge[i].Get( 1 ), edge[i].Get( 2 ) );
716 sec2.Build( edge[i].Get( 4 ), edge[i].Get( 5 ), edge[i].Get( 6 ) );
717 theta1.Build( edge[i].Get( 2 ), edge[i].Get( 5 ), edge[i].Get( 6 ) );
718 theta2.Build( edge[i].Get( 3 ), edge[i].Get( 4 ), edge[i].Get( 7 ) );
719
720 FindInt( sec1, theta1, p10mm, po2 );
721 FindInt( sec1, theta2, p10mm, po3 );
722 FindInt( sec2, theta1, p10mm, po6 );
723 FindInt( sec2, theta2, p10mm, po7 );
724
725 edge[i].Set( 2, po2 );
726 edge[i].Set( 3, po3 );
727 edge[i].Set( 6, po6 );
728 edge[i].Set( 7, po7 );
729 }
730
731 fRing5p[0] = edge[5];
732 fRing4p[0] = edge[4];
733 fRing3p[0] = edge[3];
734 fRing2p[0] = edge[2];
735 fRing1p[0] = edge[1];
736 fRing0p[0] = edge[0];
737
738 double pio2 = 3.14159265358979 / 2;
739
740 /// rotate 90 degrees
741 for ( i = 0; i < 6; ++i )
742 {
743 for ( j = 0; j < 10; ++j )
744 {
745 fRing5p[i].Set( j, fRing5p[i].Get( j ).rotateZ( pio2 ) );
746 fRing4p[i].Set( j, fRing4p[i].Get( j ).rotateZ( pio2 ) );
747 }
748 }
749 for ( i = 0; i < 5; ++i )
750 {
751 for ( j = 0; j < 10; ++j )
752 {
753 fRing3p[i].Set( j, fRing3p[i].Get( j ).rotateZ( pio2 ) );
754 fRing2p[i].Set( j, fRing2p[i].Get( j ).rotateZ( pio2 ) );
755 }
756 }
757 for ( i = 0; i < 4; ++i )
758 {
759 for ( j = 0; j < 10; ++j )
760 {
761 fRing1p[i].Set( j, fRing1p[i].Get( j ).rotateZ( pio2 ) );
762 fRing0p[i].Set( j, fRing0p[i].Get( j ).rotateZ( pio2 ) );
763 }
764 }
765
766 // check it!!!
767 // for(i=0;i<4;++i) {
768 // cout<<fRing1p[i]<<endl;
769 // }
770 // for(i=0;i<5;++i) {
771 // cout<<fRing3p[i]<<endl;
772 // }
773 // for(i=0;i<6;++i) {
774 // cout<<fRing5p[i]<<endl;
775 // }
776 for ( i = 0; i < 6; ++i )
777 {
778 if ( i < 4 )
779 {
780 fCrystal[1][0][i] = fRing0p[i];
781 fCrystal[1][1][i] = fRing1p[i];
782 }
783 if ( i < 5 )
784 {
785 fCrystal[1][2][i] = fRing2p[i];
786 fCrystal[1][3][i] = fRing3p[i];
787 }
788 fCrystal[1][4][i] = fRing4p[i];
789 fCrystal[1][5][i] = fRing5p[i];
790 }
791}
792
793void EmcRecEndCapGeo::FillCCenterVector() {
794 unsigned int module;
795 unsigned int phi, phimax, theta;
796 Identifier id;
797 EmcRecCrystal aCry;
798 HepPoint3D aCenter;
799 HepPoint3D aFrontCenter;
800
801 module = EmcID::getENDCAP_EAST();
802 for ( theta = 0; theta <= 5; ++theta )
803 {
804 phimax = EmcID::getPHI_ENDCAP_MAX( theta );
805 for ( phi = 0; phi <= phimax; ++phi )
806 {
807 id = EmcID::crystal_id( module, theta, phi );
808 aCry = GetCrystal( id );
809 aCenter = aCry.Center();
810 aFrontCenter = aCry.FrontCenter();
811 fCCenter.push_back( aCenter );
812 fCFrontCenter.push_back( aFrontCenter );
813 }
814 }
815
816 module = EmcID::getENDCAP_WEST();
817 for ( theta = 0; theta <= 5; ++theta )
818 {
819 phimax = EmcID::getPHI_ENDCAP_MAX( theta );
820 for ( phi = 0; phi <= phimax; ++phi )
821 {
822 id = EmcID::crystal_id( module, theta, phi );
823 aCry = GetCrystal( id );
824 aCenter = aCry.Center();
825 aFrontCenter = aCry.FrontCenter();
826 fCCenter.push_back( aCenter );
827 fCFrontCenter.push_back( aFrontCenter );
828 }
829 }
830}
831
832EmcRecCrystal EmcRecEndCapGeo::GetCrystal( const Identifier& id ) const {
833 int i;
834 EmcRecCrystal cry;
835 unsigned int module = EmcID::barrel_ec( id );
836 unsigned int theta = EmcID::theta_module( id );
837 unsigned int phi = EmcID::phi_module( id );
838
839 unsigned int phiMax, phiMax16;
840 unsigned int phiQuotient, phiRemainder;
841
842 phiMax = EmcID::getPHI_ENDCAP_MAX( theta );
843 phiMax16 = ( phiMax + 1 ) / 16;
844 phiQuotient = (unsigned int)( phi / phiMax16 );
845 phiRemainder = phi % phiMax16;
846
847 // cout<<phiQuotient<<" "<<phiRemainder<<endl;
848
849 if ( module == EmcID::getENDCAP_EAST() )
850 {
851 if ( phiQuotient != 3 && phiQuotient != 4 && phiQuotient != 11 && phiQuotient != 12 )
852 {
853 cry = fCrystal[0][theta][phiRemainder];
854 for ( i = 0; i < 10; ++i )
855 { cry.Set( i, cry.Get( i ).rotateZ( phiQuotient * fphi5[6] ) ); }
856 }
857 else
858 {
859 if ( phiQuotient == 4 ) { cry = fCrystal[1][theta][phiRemainder]; }
860 if ( phiQuotient == 3 )
861 {
862 cry = fCrystal[1][theta][phiMax16 - 1 - phiRemainder];
863 for ( i = 0; i < 10; ++i ) { cry.SetX( i, -cry.Get( i ).x() ); }
864 }
865 if ( phiQuotient == 11 )
866 {
867 cry = fCrystal[1][theta][phiMax16 - 1 - phiRemainder];
868 for ( i = 0; i < 10; ++i ) { cry.SetY( i, -cry.Get( i ).y() ); }
869 }
870 if ( phiQuotient == 12 )
871 {
872 cry = fCrystal[1][theta][phiRemainder];
873 for ( i = 0; i < 10; ++i )
874 {
875 cry.SetX( i, -cry.Get( i ).x() );
876 cry.SetY( i, -cry.Get( i ).y() );
877 }
878 }
879 }
880 }
881
882 if ( module == EmcID::getENDCAP_WEST() )
883 {
884 unsigned int phipp;
885 unsigned int phiMax2 = ( phiMax + 1 ) / 2;
886 if ( phi < phiMax2 ) { phipp = phiMax2 - 1 - phi; }
887 else { phipp = phiMax + phiMax2 - phi; }
888 Identifier idd = EmcID::crystal_id( EmcID::getENDCAP_EAST(), theta, phipp );
889 cry = GetCrystal( idd );
890 for ( i = 0; i < 10; ++i )
891 {
892 cry.SetX( i, -cry.Get( i ).x() );
893 cry.SetZ( i, -cry.Get( i ).z() );
894 }
895 }
896
897 return cry;
898}
899
900HepPoint3D EmcRecEndCapGeo::GetCCenter( const Identifier& id ) const {
901 unsigned int module, theta, phi;
902 unsigned int i, j;
903
904 module = EmcID::barrel_ec( id );
905 theta = EmcID::theta_module( id );
906 phi = EmcID::phi_module( id );
907
908 i = 0;
909 if ( module == EmcID::getENDCAP_EAST() )
910 {
911 for ( j = 0; j < theta; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
912 i += ( phi + 1 );
913 }
914 else
915 {
916 for ( j = 0; j < 6; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
917 for ( j = 0; j < theta; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
918 i += ( phi + 1 );
919 }
920
921 return fCCenter[i - 1];
922}
923
924HepPoint3D EmcRecEndCapGeo::GetCFrontCenter( const Identifier& id ) const {
925 unsigned int module, theta, phi;
926 unsigned int i, j;
927
928 module = EmcID::barrel_ec( id );
929 theta = EmcID::theta_module( id );
930 phi = EmcID::phi_module( id );
931
932 i = 0;
933 if ( module == EmcID::getENDCAP_EAST() )
934 {
935 for ( j = 0; j < theta; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
936 i += ( phi + 1 );
937 }
938 else
939 {
940 for ( j = 0; j < 6; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
941 for ( j = 0; j < theta; ++j ) { i += ( EmcID::getPHI_ENDCAP_MAX( j ) + 1 ); }
942 i += ( phi + 1 );
943 }
944
945 return fCFrontCenter[i - 1];
946}
947
948void EmcRecEndCapGeo::FindInt( const EmcRecGeoPlane& p1, const EmcRecGeoPlane& p2,
949 const EmcRecGeoPlane& p3, HepPoint3D& p ) {
950 // solve a system of linear equation
951 // The form is AX=B
952 HepMatrix A( 3, 3 );
953 HepVector B( 3 );
954 HepVector X( 3 );
955
956 A( 1, 1 ) = p1.a();
957 A( 1, 2 ) = p1.b();
958 A( 1, 3 ) = p1.c();
959 B( 1 ) = -p1.d();
960 A( 2, 1 ) = p2.a();
961 A( 2, 2 ) = p2.b();
962 A( 2, 3 ) = p2.c();
963 B( 2 ) = -p2.d();
964 A( 3, 1 ) = p3.a();
965 A( 3, 2 ) = p3.b();
966 A( 3, 3 ) = p3.c();
967 B( 3 ) = -p3.d();
968
969 X = solve( A, B );
970 // cout<<A;
971 // cout<<B;
972 // cout<<X;
973
974 p.setX( X( 1 ) );
975 p.setY( X( 2 ) );
976 p.setZ( X( 3 ) );
977 // cout<<p;
978}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
n
const Int_t n
Definition DataBase/tau_mode.c:57
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
w
double w
Definition EvtPhokharaDef.hh:28
theta2
@ theta2
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkKalDeriv.h:24
theta1
@ theta1
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkKalDeriv.h:24
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition InstallArea/x86_64-el9-gcc13-opt/include/EmcRecGeoSvc/EmcRecEndCapGeo.h:21
n2
int n2
Definition SD0Tag.cxx:59
n1
int n1
Definition SD0Tag.cxx:58
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition EmcID.cxx:63
EmcID::getENDCAP_WEST
static unsigned int getENDCAP_WEST()
Definition EmcID.cxx:99
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
EmcID::getENDCAP_EAST
static unsigned int getENDCAP_EAST()
Definition EmcID.cxx:95
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:41
EmcID::getPHI_ENDCAP_MAX
static unsigned int getPHI_ENDCAP_MAX(const unsigned int theta)
Definition EmcID.cxx:85
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:46
EmcRecCrystal::Set
void Set(int index, const HepPoint3D &aPoint)
Definition Emc/EmcRecGeoSvc/include/EmcRecGeoSvc/EmcRecCrystal.h:96
EmcRecEndCapGeo::GetCrystal
EmcRecCrystal GetCrystal(const Identifier &id) const
Definition EmcRecEndCapGeo.cxx:832
EmcRecEndCapGeo::GetCCenter
HepPoint3D GetCCenter(const Identifier &id) const
Definition EmcRecEndCapGeo.cxx:900
EmcRecEndCapGeo::GetCFrontCenter
HepPoint3D GetCFrontCenter(const Identifier &id) const
Definition EmcRecEndCapGeo.cxx:924
EmcRecGeoPlane::Build
void Build(double a=0, double b=0, double c=0, double d=0)
Definition Emc/EmcRecGeoSvc/include/EmcRecGeoSvc/EmcRecEndCapGeo.h:162
t
int t()
Definition t.c:1