CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemGeoAlign.cxx
Go to the documentation of this file.
1#include "CgemGeomSvc/CgemGeoAlign.h"
2
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/IMessageSvc.h"
5#include "GaudiKernel/StatusCode.h"
6#include "GaudiKernel/ISvcLocator.h"
7#include "GaudiKernel/Bootstrap.h"
8#include "GaudiKernel/SmartDataPtr.h"
9#include "GaudiKernel/IDataProviderSvc.h"
10#include "GaudiKernel/PropertyMgr.h"
11
12#include "TMath.h"
13
14#include <iostream>
15#include <fstream>
16#include <iomanip>
17
18using namespace std;
19
20CgemGeoAlign::CgemGeoAlign(){
21 m_r[0] = 87.915;
22 m_r[1] = 130.415;
23 m_r[2] = 172.915;
24}
25
26void CgemGeoAlign::initAlignPar(string alignFile){
27 ifstream fin(alignFile.c_str());
28 char strtmp[200];
29 fin.getline(strtmp,1000);
30 string str[6];
31 // It will be better to use a variable to instead of a specific number
32 for (int layer=0; layer<6; layer++){
33 fin >> str[layer] >> m_dx[layer] >> m_dy[layer] >> m_dz[layer]
34 >> m_rx[layer] >> m_ry[layer] >> m_rz[layer];
35
36 m_dxOrig[layer] = m_dx[layer];
37 m_dyOrig[layer] = m_dy[layer];
38 m_dzOrig[layer] = m_dz[layer];
39 m_rxOrig[layer] = m_rx[layer];
40 m_ryOrig[layer] = m_ry[layer];
41 m_rzOrig[layer] = m_rz[layer];
42 }
43 fin.close();
44}
45
46void CgemGeoAlign::resetAlignPar(){
47 for (int layer=0; layer<6; layer++){
48 m_dx[layer] = m_dxOrig[layer];
49 m_dy[layer] = m_dyOrig[layer];
50 m_dz[layer] = m_dzOrig[layer];
51 m_rx[layer] = m_rxOrig[layer];
52 m_ry[layer] = m_ryOrig[layer];
53 m_rz[layer] = m_rzOrig[layer];
54 }
55}
56
57
58
59StraightLine CgemGeoAlign::StraightLineConversion(int layer_vir, StraightLine lineOriginal){
60
61 int layer_geo = int(layer_vir/2);
62 // 1. Get 2 points on straight line:
63 HepPoint3D p1 = lineOriginal.xAtR(m_r[layer_geo], -1);
64 HepPoint3D p2 = lineOriginal.xAtR(m_r[layer_geo], 1);
65
66 double dist = p1.distance2(p2);
67 if(fabs(dist) < 0.001){
68 double s1 = lineOriginal.sAtR(m_r[layer_geo], -1);
69 double s2 = lineOriginal.sAtR(m_r[layer_geo], 1);
70
71 p1 = lineOriginal.x(s1);
72 p2 = lineOriginal.x(s2);
73 }
74
75 // 2. Transformation
76 // cout << "alignment par: " << setw(15) << m_dx[layer] << setw(15) << m_dy[layer]
77 // << setw(15) << m_dz[layer] << setw(15) << m_rz[layer] << endl;
78
79 // cout << "2 points before transform: "
80 // << setw(15) << p1.x() << setw(15) << p1.y() << setw(15) << p1.z()
81 // << setw(15) << p2.x() << setw(15) << p2.y() << setw(15) << p2.z() << endl;
82
83 HepPoint3D newP1 = point_transform(layer_vir, p1);
84 HepPoint3D newP2 = point_transform(layer_vir, p2);
85
86 // cout << "2 points after transform: "
87 // << setw(15) << newP1.x() << setw(15) << newP1.y() << setw(15) << newP1.z()
88 // << setw(15) << newP2.x() << setw(15) << newP2.y() << setw(15) << newP2.z() << endl;
89
90 // 3. New line parameters
91 StraightLine newLine(newP1, newP2);
92
93 // cout << "use StraightLine" << endl;
94 // cout << "track before transform: " << setw(15) << lineOriginal.dr() << setw(15) << lineOriginal.phi0()
95 // << setw(15) << lineOriginal.dz() << setw(15) << lineOriginal.tanl() << endl;
96 // cout << "track after transform: " << setw(15) << newLine.dr() << setw(15) << newLine.phi0()
97 // << setw(15) << newLine.dz() << setw(15) << newLine.tanl() << endl << endl;
98
99 return newLine;
100}
101
102void CgemGeoAlign::StraightLineConversion_v1(int layer_vir, double lineOriginal[], double lineConverted[]){
103
104 // m_dx[0] = 1;
105 // m_dy[0] = 1;
106 // m_dz[0] = 1;
107 // m_rz[0] = 0; //TMath::Pi()/2.;
108
109 // 0. Get original straight line parameters:
110
111 double drho = lineOriginal[0];
112 double phi0 = lineOriginal[1];
113 double dz = lineOriginal[2];
114 double tgl = lineOriginal[3];
115
116 // 1. Get 2 points on straight line:
117
118 //--- Line on x-y plane: y = ax + b
119 //--- Line on s-z plane: z = -1.*tgl*s + dz;
120 // s = sqrt((x-x0)^2+(y-y0)^2)
121
122 int flg_parallel_x = 0;
123 int flg_parallel_y = 0;
124
125 if(phi0 == TMath::Pi()/2. || phi0 == -1.*TMath::Pi()/2.) flg_parallel_x = 1;
126 if(phi0 == 0 || phi0 == TMath::Pi()) flg_parallel_y = 1;
127
128 double a, b, x0, x1, x2, y0, y1, y2, s1, s2, z1, z2;
129 x0 = drho*cos(phi0);
130 y0 = drho*sin(phi0);
131
132 if(flg_parallel_x == 0 && flg_parallel_y == 0){
133 a = tan(TMath::Pi()/2.+phi0);
134 b = drho / cos(phi0);
135
136 x1 = 0;
137 y1 = drho / cos(phi0);
138 s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
139 z1 = -1.*tgl*s1 + dz;
140
141 x2 = -1.*b/a;
142 y2 = 0;
143 s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
144 z2 = -1.*tgl*s2 + dz;
145 if(x1==x2 || y1==y2 || z1==z2){
146 x2 = (1.-b)/a;
147 y2 = 1;
148 s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
149 z2 = -1.*tgl*s2 + dz;
150 }
151 }
152 else if(flg_parallel_x == 1){
153 x1 = drho; y1 = 0; s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)); z1 = -1.*tgl*s1 + dz;
154 x2 = drho; y2 = 1; s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0)); z2 = -1.*tgl*s2 + dz;
155 }
156 else if(flg_parallel_y == 1){
157 x1 = 0; y1 = drho; s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)); z1 = -1.*tgl*s1 + dz;
158 x2 = 1; y2 = drho; s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0)); z2 = -1.*tgl*s2 + dz;
159 }
160
161 // 2. Transformation
162
163 double shift_x = m_dx[layer_vir];
164 double shift_y = m_dy[layer_vir];
165 double shift_z = m_dz[layer_vir];
166 double rotation_z = m_rz[layer_vir];
167
168 // cout << "alignment par: " << setw(15) << shift_x << setw(15) << shift_y
169 // << setw(15) << shift_z << setw(15) << rotation_z << endl;
170
171 // cout << "2 points before transform: "
172 // << setw(15) << x1 << setw(15) << y1 << setw(15) << z1
173 // << setw(15) << x2 << setw(15) << y2 << setw(15) << z2 << endl;
174
175 point_transform(x1, y1, z1, shift_x, shift_y, shift_z, rotation_z);
176 point_transform(x2, y2, z2, shift_x, shift_y, shift_z, rotation_z);
177
178 // cout << "2 points after transform: "
179 // << setw(15) << x1 << setw(15) << y1 << setw(15) << z1
180 // << setw(15) << x2 << setw(15) << y2 << setw(15) << z2 << endl;
181
182 // 3. New line parameters
183
184 //--- Line on x-y plane: y = ax + b
185 // y1 = ax1 + b
186 // y2 = ax2 + b
187 // a = (y1-y2)/(x1-x2) = -tan(TMath::Pi()/2.-phi0) = tan(phi0-TMath::Pi()/2.);
188 // b = y1-a*x1 = drho / cos(phi0);
189 // phi0 = atan(a)+TMath::Pi()/2.;
190 // drho = b*cos(phi0);
191 //--- Line on s-z plane: z = -1.*tgl*s + dz;
192 // s = sqrt((x-x0)^2+(y-y0)^2)
193 // z1 = -1.*tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))+dz
194 // z2 = -1.*tgl*sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))+dz
195 // tgl = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
196 // dz = z1 + tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
197
198 double new_a;
199 double new_b;
200 double new_phi0;
201 double new_drho;
202 double new_tgl;
203 double new_dz;
204
205 if(x1!=x2 && y1!=y2){
206 new_a = (y1-y2)/(x1-x2);
207 new_b = y1-new_a*x1;
208 new_phi0 = atan(new_a)+TMath::Pi()/2.;
209 new_drho = new_b*cos(phi0);
210 new_tgl = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
211 new_dz = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
212 }
213 else if(x1==x2){
214 new_phi0 = 0;
215 new_drho = x1;
216 new_tgl = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
217 new_dz = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
218 }
219 else if(y1==y2){
220 new_phi0 = TMath::Pi()/2.;
221 new_drho = y1;
222 new_tgl = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
223 new_dz = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
224 }
225
226 lineConverted[0] = new_drho;
227 lineConverted[1] = new_phi0;
228 lineConverted[2] = new_dz;
229 lineConverted[3] = new_tgl;
230
231 // cout << "track before transform: " << setw(15) << lineOriginal[0] << setw(15) << lineOriginal[1]
232 // << setw(15) << lineOriginal[2] << setw(15) << lineOriginal[3] << endl;
233 // cout << "track after transform: " << setw(15) << lineConverted[0] << setw(15) << lineConverted[1]
234 // << setw(15) << lineConverted[2] << setw(15) << lineConverted[3] << endl << endl;
235
236}
237
238void CgemGeoAlign::HelixConversion(int layer_vir, double helixOriginal[], double helixConverted[]){
239}
240
241HepPoint3D CgemGeoAlign::point_transform(int layer_vir, HepPoint3D pos){
242 // this function is aimed to transform global coordinate (with alignment) to local coordinate (without alignment).
243
244 double x_shift = pos.x() - m_dx[layer_vir];
245 double y_shift = pos.y() - m_dy[layer_vir];
246 double z_shift = pos.z() - m_dz[layer_vir];
247
248 Hep3Vector InputVector(x_shift,y_shift,z_shift);
249
250 InputVector.rotateZ(-m_rz[layer_vir]);
251 InputVector.rotateY(-m_ry[layer_vir]);
252 InputVector.rotateX(-m_rx[layer_vir]);
253
254 double x = InputVector.x();
255 double y = InputVector.y();
256 double z = InputVector.z();
257
258 return HepPoint3D(x,y,z);
259}
260
261
262HepPoint3D CgemGeoAlign::point_transform(int layer_geo, int sheet_geo, HepPoint3D pos){
263 // this function is aimed to transform global coordinate (with alignment) to local coordinate (without alignment).
264
265 int layer_vir = layer_geo*2+sheet_geo;
266
267 double x_shift = pos.x() - m_dx[layer_vir];
268 double y_shift = pos.y() - m_dy[layer_vir];
269 double z_shift = pos.z() - m_dz[layer_vir];
270
271 Hep3Vector InputVector(x_shift,y_shift,z_shift);
272
273 InputVector.rotateZ(-m_rz[layer_vir]);
274 InputVector.rotateY(-m_ry[layer_vir]);
275 InputVector.rotateX(-m_rx[layer_vir]);
276
277 double x = InputVector.x();
278 double y = InputVector.y();
279 double z = InputVector.z();
280
281 return HepPoint3D(x,y,z);
282}
283
284
285
286void CgemGeoAlign::point_transform(double& x, double& y, double& z, double shift_x,
287 double shift_y, double shift_z, double rotation_z){
288
289 double x_shift = x - shift_x;
290 double y_shift = y - shift_y;
291 double z_shift = z - shift_z;
292
293 x = x_shift*cos(rotation_z) + y_shift*sin(rotation_z);
294 y = -1.*x_shift*sin(rotation_z) + y_shift*cos(rotation_z);
295 z = z_shift;
296}
297
298
299
300// HepPoint3D CgemGeoAlign::point_invTransform(int layer_vir, HepPoint3D pos){
301// // this function is aimed to transform local coordinate (without alignment) to global coordinate (with alignment).
302// // this version is aimed only for 4 alignment parameters on each layer: dx,dy,dz,rz
303// double xp = pos.x();
304// double yp = pos.y();
305// double zp = pos.z();
306//
307// double rz = -1.*m_rz[layer_vir];
308// double xrot = xp*cos(rz) + yp*sin(rz);
309// double yrot = -1.*xp*sin(rz) + yp*cos(rz);
310//
311// double x = xrot + m_dx[layer_vir];
312// double y = yrot + m_dy[layer_vir];
313// double z = zp + m_dz[layer_vir];
314//
315// // cout << "before invTransform: " << setw(15) << xp<< setw(15) << yp << setw(15) << zp << endl;
316// // cout << " after invTransform: " << setw(15) << xp<< setw(15) << yp << setw(15) << zp << endl;
317// return HepPoint3D(x, y, z);
318// }
319
320
321
322HepPoint3D CgemGeoAlign::point_invTransform(int layer_vir, HepPoint3D pos){
323
324 // this function is aimed to transform local coordinate (without alignment) to global coordinate (with alignment).
325 double xp = pos.x();
326 double yp = pos.y();
327 double zp = pos.z();
328
329 Hep3Vector InputVector(xp,yp,zp);
330
331 InputVector.rotateX(m_rx[layer_vir]);
332 InputVector.rotateY(m_ry[layer_vir]);
333 InputVector.rotateZ(m_rz[layer_vir]);
334
335 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
336 Hep3Vector InputVector_t1 = InputVector+OO;
337
338 double x = InputVector_t1.x();
339 double y = InputVector_t1.y();
340 double z = InputVector_t1.z();
341
342 return HepPoint3D(x, y, z);
343}
344
345HepPoint3D CgemGeoAlign::point_invTransform(int layer_geo, int sheet_geo, HepPoint3D pos){
346
347 // this function is aimed to transform local coordinate (without alignment) to global coordinate (with alignment).
348 int layer_vir = layer_geo*2+sheet_geo;
349
350 double xp = pos.x();
351 double yp = pos.y();
352 double zp = pos.z();
353
354
355 Hep3Vector InputVector(xp,yp,zp);
356
357 InputVector.rotateX(m_rx[layer_vir]);
358 InputVector.rotateY(m_ry[layer_vir]);
359 InputVector.rotateZ(m_rz[layer_vir]);
360
361 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
362 Hep3Vector InputVector_t1 = InputVector+OO;
363
364 double x = InputVector_t1.x();
365 double y = InputVector_t1.y();
366 double z = InputVector_t1.z();
367
368 return HepPoint3D(x, y, z);
369}
370
371// ===========================================================================================================
372//
373// ===========================================================================================================
374//
375// The code below is used to calculate the intersections between a given line and a given Cylinder with
376// arbitrary orientation !!
377//
378// ===========================================================================================================
379
380
381void CgemGeoAlign::GetLocPos(int layer_vir, Hep3Vector GposUp, Hep3Vector GposDown, HepPoint3D& LposUp, HepPoint3D& LposDown)
382{
383 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
384 Hep3Vector LVposUp = GposUp - OO;
385 Hep3Vector LVposDown = GposDown - OO;
386
387 // the sign dependent on the defination of rotation direction
388 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
389 LVposUp.rotateZ(-m_rz[layer_vir]);
390 LVposUp.rotateY(-m_ry[layer_vir]);
391 LVposUp.rotateX(-m_rx[layer_vir]);
392
393 LVposDown.rotateZ(-m_rz[layer_vir]);
394 // Why should I do these two rotations
395 LVposDown.rotateY(-m_ry[layer_vir]);
396 LVposDown.rotateX(-m_rx[layer_vir]);
397
398 LposUp.setX(LVposUp.x());
399 LposUp.setY(LVposUp.y());
400 LposUp.setZ(LVposUp.z());
401
402 LposDown.setX(LVposDown.x());
403 LposDown.setY(LVposDown.y());
404 LposDown.setZ(LVposDown.z());
405
406}
407
408
409
410void CgemGeoAlign::GetLocPos(int layer_vir, Hep3Vector Gpos, HepPoint3D& Lpos)
411{
412 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
413 Hep3Vector LVpos = Gpos - OO;
414
415 // the sign dependent on the defination of rotation direction
416 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
417 LVpos.rotateZ(-m_rz[layer_vir]);
418 LVpos.rotateY(-m_ry[layer_vir]);
419 LVpos.rotateX(-m_rx[layer_vir]);
420
421 Lpos.setX(LVpos.x());
422 Lpos.setY(LVpos.y());
423 Lpos.setZ(LVpos.z());
424
425}
426
427
428
429void CgemGeoAlign::GetLocPos(int layer_geo, int sheet_geo, Hep3Vector GposUp, Hep3Vector GposDown, HepPoint3D& LposUp, HepPoint3D& LposDown)
430{
431
432 int layer_vir = int(layer_geo*2+sheet_geo);
433 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
434 Hep3Vector LVposUp = GposUp - OO;
435 Hep3Vector LVposDown = GposDown - OO;
436
437 // the sign dependent on the defination of rotation direction
438 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
439 LVposUp.rotateZ(-m_rz[layer_vir]);
440 LVposUp.rotateY(-m_ry[layer_vir]);
441 LVposUp.rotateX(-m_rx[layer_vir]);
442
443 LVposDown.rotateZ(-m_rz[layer_vir]);
444 // Why should I do these two rotations
445 LVposDown.rotateY(-m_ry[layer_vir]);
446 LVposDown.rotateX(-m_rx[layer_vir]);
447
448 LposUp.setX(LVposUp.x());
449 LposUp.setY(LVposUp.y());
450 LposUp.setZ(LVposUp.z());
451
452 LposDown.setX(LVposDown.x());
453 LposDown.setY(LVposDown.y());
454 LposDown.setZ(LVposDown.z());
455
456}
457
458
459
460void CgemGeoAlign::GetLocPos(int layer_geo, int sheet_geo, Hep3Vector Gpos, HepPoint3D& Lpos)
461{
462
463 int layer_vir = int(layer_geo*2+sheet_geo);
464 Hep3Vector OO(m_dx[layer_vir],m_dy[layer_vir],m_dz[layer_vir]);
465 Hep3Vector LVpos = Gpos - OO;
466
467 // the sign dependent on the defination of rotation direction
468 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
469 LVpos.rotateZ(-m_rz[layer_vir]);
470 LVpos.rotateY(-m_ry[layer_vir]);
471 LVpos.rotateX(-m_rx[layer_vir]);
472
473 Lpos.setX(LVpos.x());
474 Lpos.setY(LVpos.y());
475 Lpos.setZ(LVpos.z());
476
477}
478
479
480bool CgemGeoAlign::getinter(ToyRay ARay, int layer_vir, Hep3Vector& posUp, Hep3Vector& posDown)
481{
482 int layer_geo = int(layer_vir/2);
483
484 ToyCgem ACgem(layer_vir,m_r[layer_geo], m_dx[layer_vir], m_dy[layer_vir], m_dz[layer_vir], m_rx[layer_vir], m_ry[layer_vir], m_rz[layer_vir]);
485
486 Hep3Vector LV(cos(ARay.GetPhi()-TMath::Pi()/2.),sin(ARay.GetPhi()-TMath::Pi()/2.),ARay.GetTanL());
487 Hep3Vector N = LV.unit();
488 Hep3Vector O(ARay.GetDrho()*cos(ARay.GetPhi()), ARay.GetDrho()*sin(ARay.GetPhi()), ARay.GetDz());
489 Hep3Vector A = ACgem.getdir();
490 Hep3Vector B(ACgem.GetX0(), ACgem.GetY0(),ACgem.GetZ0());
491 double R = ACgem.GetR();
492
493 if (N.cross(A).mag()==0)
494 {
495 return false;
496 }
497
498
499 if (((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))<0)
500 {
501 return false;
502 }
503
504 // Calculate the sulution
505 double D1 = (((N.cross(A)).dot((B-O).cross(A)))+sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
506 double D2 = (((N.cross(A)).dot((B-O).cross(A)))-sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
507
508 // Calculate the Z positons of the intersections
509 double z1 = A.dot(N*D1-B+O);
510 double z2 = A.dot(N*D2-B+O);
511
512 // Get the final solution
513 Hep3Vector Inter1 = D1*N + O;
514 Hep3Vector Inter2 = D2*N + O;
515
516 // the sign dependent on the defination of ratation direction
517 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
518
519 posUp = (Inter1.phi()>Inter2.phi())? Inter1:Inter2;
520 posDown = (Inter1.phi()>Inter2.phi())? Inter2:Inter1;
521
522 return true;
523
524}
525
526
527
528
529bool CgemGeoAlign::getinter(ToyRay ARay, int layer_vir, Hep3Vector& pos)
530{
531 int layer_geo = int(layer_vir/2);
532
533 ToyCgem ACgem(layer_vir,m_r[layer_geo], m_dx[layer_vir], m_dy[layer_vir], m_dz[layer_vir], m_rx[layer_vir], m_ry[layer_vir], m_rz[layer_vir]);
534
535 Hep3Vector LV(cos(ARay.GetPhi()-TMath::Pi()/2.),sin(ARay.GetPhi()-TMath::Pi()/2.),ARay.GetTanL());
536 Hep3Vector N = LV.unit();
537 Hep3Vector O(ARay.GetDrho()*cos(ARay.GetPhi()), ARay.GetDrho()*sin(ARay.GetPhi()), ARay.GetDz());
538 Hep3Vector A = ACgem.getdir();
539 Hep3Vector B(ACgem.GetX0(), ACgem.GetY0(),ACgem.GetZ0());
540 double R = ACgem.GetR();
541
542 if (N.cross(A).mag()==0)
543 {
544 return false;
545 }
546
547
548 if (((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))<0)
549 {
550 return false;
551 }
552
553 // Calculate the sulution
554 double D1 = (((N.cross(A)).dot((B-O).cross(A)))+sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
555 double D2 = (((N.cross(A)).dot((B-O).cross(A)))-sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
556
557 // Calculate the Z positons of the intersections
558 double z1 = A.dot(N*D1-B+O);
559 double z2 = A.dot(N*D2-B+O);
560
561 // Get the final solution
562 Hep3Vector Inter1 = D1*N + O;
563 Hep3Vector Inter2 = D2*N + O;
564
565 // the sign dependent on the defination of ratation direction
566 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
567
568 Hep3Vector posUp = (Inter1.phi()>Inter2.phi())? Inter1:Inter2;
569 Hep3Vector posDown = (Inter1.phi()>Inter2.phi())? Inter2:Inter1;
570
571 pos = (layer_vir%2==0)? posDown:posUp;
572 return true;
573
574}
575
576
577
578
579
580
581
582
583bool CgemGeoAlign::getinter(ToyRay ARay, int layer_geo, int sheet_geo, Hep3Vector& posUp, Hep3Vector& posDown)
584{
585 int layer_vir = int(layer_geo*2+sheet_geo);
586
587 ToyCgem ACgem(layer_vir,m_r[layer_geo], m_dx[layer_vir], m_dy[layer_vir], m_dz[layer_vir], m_rx[layer_vir], m_ry[layer_vir], m_rz[layer_vir]);
588
589 Hep3Vector LV(cos(ARay.GetPhi()-TMath::Pi()/2.),sin(ARay.GetPhi()-TMath::Pi()/2.),ARay.GetTanL());
590 Hep3Vector N = LV.unit();
591 Hep3Vector O(ARay.GetDrho()*cos(ARay.GetPhi()), ARay.GetDrho()*sin(ARay.GetPhi()), ARay.GetDz());
592 Hep3Vector A = ACgem.getdir();
593 Hep3Vector B(ACgem.GetX0(), ACgem.GetY0(),ACgem.GetZ0());
594 double R = ACgem.GetR();
595
596 if (N.cross(A).mag()==0)
597 {
598 return false;
599 }
600
601 if (((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))<0)
602 {
603 return false;
604 }
605
606 // Calculate the sulution
607 double D1 = (((N.cross(A)).dot((B-O).cross(A)))+sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
608 double D2 = (((N.cross(A)).dot((B-O).cross(A)))-sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
609
610 // Calculate the Z positons of the intersections
611 double z1 = A.dot(N*D1-B+O);
612 double z2 = A.dot(N*D2-B+O);
613
614 // Get the final solution
615 Hep3Vector Inter1 = D1*N + O;
616 Hep3Vector Inter2 = D2*N + O;
617
618 // the sign dependent on the defination of ratation direction
619 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
620
621 posUp = (Inter1.phi()>Inter2.phi())? Inter1:Inter2;
622 posDown = (Inter1.phi()>Inter2.phi())? Inter2:Inter1;
623
624 return true;
625
626}
627
628
629
630
631bool CgemGeoAlign::getinter(ToyRay ARay, int layer_geo, int sheet_geo, Hep3Vector& pos)
632{
633 int layer_vir = int(layer_geo*2+sheet_geo);
634
635 ToyCgem ACgem(layer_vir,m_r[layer_geo], m_dx[layer_vir], m_dy[layer_vir], m_dz[layer_vir], m_rx[layer_vir], m_ry[layer_vir], m_rz[layer_vir]);
636
637 Hep3Vector LV(cos(ARay.GetPhi()-TMath::Pi()/2.),sin(ARay.GetPhi()-TMath::Pi()/2.),ARay.GetTanL());
638 Hep3Vector N = LV.unit();
639 Hep3Vector O(ARay.GetDrho()*cos(ARay.GetPhi()), ARay.GetDrho()*sin(ARay.GetPhi()), ARay.GetDz());
640 Hep3Vector A = ACgem.getdir();
641 Hep3Vector B(ACgem.GetX0(), ACgem.GetY0(),ACgem.GetZ0());
642 double R = ACgem.GetR();
643
644 if (N.cross(A).mag()==0)
645 {
646 return false;
647 }
648
649 if (((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))<0)
650 {
651 return false;
652 }
653
654 // Calculate the sulution
655 double D1 = (((N.cross(A)).dot((B-O).cross(A)))+sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
656 double D2 = (((N.cross(A)).dot((B-O).cross(A)))-sqrt(((N.cross(A)).dot(N.cross(A)))*R*R-(A.dot(A)*((B-O).dot(N.cross(A)))*((B-O).dot(N.cross(A))))))/((N.cross(A)).dot(N.cross(A)));
657
658 // Calculate the Z positons of the intersections
659 double z1 = A.dot(N*D1-B+O);
660 double z2 = A.dot(N*D2-B+O);
661
662 // Get the final solution
663 Hep3Vector Inter1 = D1*N + O;
664 Hep3Vector Inter2 = D2*N + O;
665
666 // the sign dependent on the defination of ratation direction
667 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
668
669 Hep3Vector posUp = (Inter1.phi()>Inter2.phi())? Inter1:Inter2;
670 Hep3Vector posDown = (Inter1.phi()>Inter2.phi())? Inter2:Inter1;
671
672 pos = (layer_vir%2==0)? posDown:posUp;
673
674 return true;
675
676}
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693ToyCgem::ToyCgem(int layer_vir, double r, double x0, double y0, double z0, double rx, double ry, double rz)
694{
695
696 int layer_geo = int(layer_vir/2.);
697 R = r;
698 Layer = layer_geo;
699 X0 = x0;
700 Y0 = y0;
701 Z0 = z0;
702 RX = rx;
703 RY = ry;
704 RZ = rz;
705
706}
707
708Hep3Vector ToyCgem::getdir()
709{
710 Hep3Vector orin(0,0,1);
711 // the sign dependent on the defination of ratation direction. Here we assume a clock wise rotation
712 // change the rotation to anti-clock wise direction By - Guoaq-Sep-28-2020
713 orin.rotateX(RX);
714 orin.rotateY(RY);
715 orin.rotateZ(RZ);
716 return orin;
717}
718
719ToyRay::ToyRay(double dr, double phi, double dz, double tanl)
720{
721 Drho = dr;
722 Phi = phi;
723 Dz = dz;
724 TanL = tanl;
725}
726
727Hep3Vector ToyRay::getdir()
728{
729 Hep3Vector LV(cos(Phi-TMath::Pi()/2.),sin(Phi-TMath::Pi()/2.),TanL);
730 return LV.unit();
731}
732
733
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
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
cross
EvtVector3R cross(const EvtVector3R &p1, const EvtVector3R &p2)
Definition EvtVector3R.cc:84
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37
CgemGeoAlign::StraightLineConversion_v1
void StraightLineConversion_v1(int layer_vir, double lineOriginal[], double lineConverted[])
Definition CgemGeoAlign.cxx:102
CgemGeoAlign::StraightLineConversion
StraightLine StraightLineConversion(int layer_vir, StraightLine lineOriginal)
Definition CgemGeoAlign.cxx:59
CgemGeoAlign::GetLocPos
void GetLocPos(int layer_vir, Hep3Vector GposUp, Hep3Vector GposDown, HepPoint3D &LposUp, HepPoint3D &LposDown)
Definition CgemGeoAlign.cxx:381
CgemGeoAlign::point_invTransform
HepPoint3D point_invTransform(int layer_vir, HepPoint3D pos)
Definition CgemGeoAlign.cxx:322
CgemGeoAlign::initAlignPar
void initAlignPar(std::string alignFile)
Definition CgemGeoAlign.cxx:26
CgemGeoAlign::HelixConversion
void HelixConversion(int layer_vir, double helixOriginal[], double helixConverted[])
Definition CgemGeoAlign.cxx:238
CgemGeoAlign::resetAlignPar
void resetAlignPar()
Definition CgemGeoAlign.cxx:46
CgemGeoAlign::getinter
bool getinter(ToyRay ARay, int layer_vir, Hep3Vector &posUp, Hep3Vector &posDown)
Definition CgemGeoAlign.cxx:480
CgemGeoAlign::point_transform
HepPoint3D point_transform(int layer_vir, HepPoint3D pos)
Definition CgemGeoAlign.cxx:241
StraightLine::x
HepPoint3D x(double s=0.) const
returns position after moving s in downwoards
Definition StraightLine.cxx:81
StraightLine::xAtR
HepPoint3D xAtR(double R, int direction=1) const
Definition StraightLine.cxx:100
StraightLine::sAtR
double sAtR(double R, int direction=1) const
Definition StraightLine.cxx:89
ToyCgem::getdir
Hep3Vector getdir()
Definition CgemGeoAlign.cxx:708
ToyCgem::GetZ0
double GetZ0()
Definition CgemGeoAlign.h:21
ToyCgem::GetX0
double GetX0()
Definition CgemGeoAlign.h:19
ToyCgem::GetR
double GetR()
Definition CgemGeoAlign.h:25
ToyCgem::GetY0
double GetY0()
Definition CgemGeoAlign.h:20
ToyCgem::ToyCgem
ToyCgem(int layer_vir, double r, double x0, double y0, double z0, double rx, double ry, double rz)
Definition CgemGeoAlign.cxx:693
ToyRay::GetDz
double GetDz()
Definition CgemGeoAlign.h:47
ToyRay::GetTanL
double GetTanL()
Definition CgemGeoAlign.h:48
ToyRay::getdir
Hep3Vector getdir()
Definition CgemGeoAlign.cxx:727
ToyRay::ToyRay
ToyRay(double dr, double phi, double dz, double tanl)
Definition CgemGeoAlign.cxx:719
ToyRay::GetDrho
double GetDrho()
Definition CgemGeoAlign.h:45
ToyRay::GetPhi
double GetPhi()
Definition CgemGeoAlign.h:46