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