CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemMilleAlign.cxx
Go to the documentation of this file.
4
5#include "GaudiKernel/IMessageSvc.h"
6#include "GaudiKernel/StatusCode.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/Bootstrap.h"
9
10#include "GaudiKernel/INTupleSvc.h"
11#include "GaudiKernel/NTuple.h"
12
13#include <iostream>
14#include <fstream>
15#include <iomanip>
16#include <string>
17#include <cstring>
18#include <cmath>
19
20#include "TMath.h"
21#include "TSpline.h"
22#include "TString.h"
23#include "TCanvas.h"
24#include "TGraph.h"
25
26#ifndef ENABLE_BACKWARDS_COMPATIBILITY
27// backwards compatblty wll be enabled ONLY n CLHEP 1.9
29#endif
30
31using namespace std;
32
34 cout << "CgemMilleAlign::CgemMilleAlign()" << endl;
35 for(int lay=0; lay<LAYERNMAX; lay++){
36 m_resiCut_phi[lay] = 7.0; //how to optimize this number guoaq-19-03-21
37 m_resiCut_v[lay] = 7.0; //how to optimize this number guoaq-19-03-21
38 m_docaCut[lay] = 0.5;
39 m_docaCut[lay] = 5.5;
40 }
41}
42
44 cout << "CgemMilleAlign::~CgemMilleAlign()" << endl;
45}
46
48 delete m_hresAll_phi;
49 delete m_hresAll_v;
50 delete m_hresAll_z;
51 for(int lay=0; lay<LAYERNMAX; lay++)
52 {
53 delete m_hresLay_phi[lay];
54 delete m_hresLay_v[lay];
55 delete m_hresLay_z[lay];
56 };
57 delete m_hddoca;
58 delete m_pMilleAlign;
59}
60
61
62void CgemMilleAlign::initialize(TObjArray* hlist, ICgemGeomSvc* cgemGeomSvc, ICgemCalibFunSvc* cgemFunSvc){
63 IMessageSvc* msgSvc;
64 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
65 MsgStream log(msgSvc, "MilleAlign");
66 log << MSG::INFO << "MilleAlign::initialize()" << endreq;
67
68 CgemAlignBase::initialize(hlist, cgemGeomSvc, cgemFunSvc);
69 m_hlist = hlist;
70 m_cgemGeomSvc = cgemGeomSvc;
71 m_cgemFunSvc = cgemFunSvc;
72
73 for(int il=0; il<6; il++)
74 {
75 int ilgeo = int(il/2);
76 m_r[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getMiddleROfGapD());
77 m_z[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getLengthOfCgemLayer()/2.0);
78 Ste_ang[il]=(m_cgemGeomSvc->getCgemLayer(ilgeo)->getAngleOfStereo()*TMath::Pi()/180.0);
79 }
80
81 int DataType = 1;
82
83 if(DataType==1)
84 {
85
86 sig_phi[0] = 0.13/m_r[0];
87 sig_phi[1] = 0.13/m_r[1];
88 sig_phi[2] = 0.13/m_r[2];
89 sig_phi[3] = 0.13/m_r[3];
90 sig_phi[4] = 0.13/m_r[4];
91 sig_phi[5] = 0.13/m_r[5];
92
93 sig_v[0] = 0.13;
94 sig_v[1] = 0.13;
95 sig_v[2] = 0.13;
96 sig_v[3] = 0.13;
97 sig_v[4] = 0.13;
98 sig_v[5] = 0.13;
99
100 }
101 else if(DataType==2)
102 {
103 // for DATA
104 sig_phi[0] = 0.03;
105 sig_phi[1] = 0.03;
106 sig_phi[2] = 0.03;
107 sig_phi[3] = 0.03;
108 sig_phi[4] = 0.03;
109 sig_phi[5] = 0.03;
110
111 sig_v[0] = 1.7;
112 sig_v[1] = 1.7;
113 sig_v[2] = 1.2;
114 sig_v[3] = 1.2;
115 sig_v[4] = 1.2;
116 sig_v[5] = 1.2;
117 }
118
119 for(int i=0; i<30; i++)
120 {
121 count[i] = 0 ;
122 }
123
124 // cuts of derivatives
125 cut_derLC_Phi[0][0] = -0.05;
126 cut_derLC_Phi[1][0] = 0.98;
127 cut_derLC_Phi[2][0] = -0.01;
128 cut_derLC_Phi[3][0] = -0.01;
129 cut_derLC_Phi[4][0] = -0.01;
130
131 cut_derLC_Phi[0][1] = 0.05;
132 cut_derLC_Phi[1][1] = 1.02;
133 cut_derLC_Phi[2][1] = 0.01;
134 cut_derLC_Phi[3][1] = 0.01;
135 cut_derLC_Phi[4][1] = 0.01;
136
137 cut_derLC_V[0][0] = -5.0;
138 cut_derLC_V[1][0] = 60.0;
139 cut_derLC_V[2][0] = -0.01;
140 cut_derLC_V[3][0] = 0.45;
141 cut_derLC_V[4][0] = -90.0;
142
143 cut_derLC_V[0][1] = 5.0;
144 cut_derLC_V[1][1] = 150;
145 cut_derLC_V[2][1] = 0.01;
146 cut_derLC_V[3][1] = 0.80;
147 cut_derLC_V[4][1] = 90.0;
148
149 cut_derGB_Phi[0][0] = -0.05;
150 cut_derGB_Phi[1][0] = -0.02;
151 cut_derGB_Phi[2][0] = -0.01;
152 cut_derGB_Phi[3][0] = -1.02;
153
154 cut_derGB_Phi[0][1] = 0.05;
155 cut_derGB_Phi[1][1] = 0.02;
156 cut_derGB_Phi[2][1] = 0.01;
157 cut_derGB_Phi[3][1] = -0.98;
158
159 cut_derGB_V[0][0] = -5.00;
160 cut_derGB_V[1][0] = -5.00;
161 cut_derGB_V[2][0] = -0.80;
162 cut_derGB_V[3][0] = -150;
163
164 cut_derGB_V[0][1] = 5.00;
165 cut_derGB_V[1][1] = 5.00;
166 cut_derGB_V[2][1] = -0.50;
167 cut_derGB_V[3][1] = -60;
168
169 // initialize hitograms
170 m_hresAll_phi = new TH1F("Resi_All_Phi", "", 200, -2.0, 2.0);
171 m_hlist->Add(m_hresAll_phi);
172
173 m_hresAll_v = new TH1F("Resi_All_v", "", 100, -30.0, 30.0);
174 m_hlist->Add(m_hresAll_v);
175
176 m_hresAll_z = new TH1F("Resi_All_z", "", 100, -30.0, 30.0);
177 m_hlist->Add(m_hresAll_z);
178
179 char hname_phi[200];
180 char hname_v[200];
181 char hname_z[200];
182
183 for(int lay=0; lay<LAYERNMAX; lay++){
184 sprintf(hname_phi, "Res_Phi_Layer%02d", lay);
185 m_hresLay_phi[lay] = new TH1F(hname_phi, "", 200, -2.0, 2.0);
186 m_hlist->Add(m_hresLay_phi[lay]);
187
188 sprintf(hname_v, "Res_V_Layer%02d", lay);
189 m_hresLay_v[lay] = new TH1F(hname_v, "", 100, -30.0, 30.0);
190 m_hlist->Add(m_hresLay_v[lay]);
191
192 sprintf(hname_z, "Res_Z_Layer%02d", lay);
193 m_hresLay_z[lay] = new TH1F(hname_z, "", 100, -30.0, 30.0);
194 m_hlist->Add(m_hresLay_z[lay]);
195 }
196
197 // for debug
198 m_hddoca = new TH1F("delt_doca", "", 200, -1.0, 1.0);
199 m_hlist->Add(m_hddoca);
200
201 char hname[200];
202 for(int lay=0; lay<LAYERNMAX; lay++){
203 sprintf(hname, "delt_docaLay%02d", lay);
204 m_hddocaLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
205 m_hlist->Add(m_hddocaLay[lay]);
206 }
207
208 m_GeoMethod = GEO_METHOD;
209 Plot_der = false;
210 debug = false;
211 re_3lay = false;
212 record_pos = true;
213 use_phi = true;
214 use_v = true;
215 mul_R = false;
216 mul_Ste_ang = false;
217
218 nTrack_tot = 0;
219
220 // initialize millepede
221
222 m_nlay = NLAYER;
223 m_nloc = NTRKPAR;
224 m_nglo_on_lay = NDOFALIGN; // NDOFALIGN : for each hit, 4 alignment parameters on each layer: dx, dy, dz, rz
225 m_npar = NDOFALIGN * m_nlay;// that is the global parameters, isn't it?
226
227 if(debug)cout<<"m_nlay, m_nloc, m_nglo_on_lay, m_npar = "<<m_nlay<<", "<<m_nloc<<", "<<m_nglo_on_lay<<", "<<m_npar<<endl;
228
229 m_middriftplane = m_cgemGeomSvc->getMidDriftPtr();
230 m_geoalign = m_cgemGeomSvc->getAlignPtr();
231
232 // get the initial alignment parameter and save it to a list
233 for(int i=0; i<NLAYER; i++)
234 {
235
236 ini_alig_par[i][0] = m_geoalign->getDx(i);
237 ini_alig_par[i][1] = m_geoalign->getDy(i);
238 ini_alig_par[i][2] = m_geoalign->getDz(i);
239 ini_alig_par[i][3] = m_geoalign->getRx(i);
240 ini_alig_par[i][4] = m_geoalign->getRy(i);
241 ini_alig_par[i][5] = m_geoalign->getRz(i);
242
243 }
244
245
246 int i;
247 for(i=0; i<NDOFALIGN; i++){
248 m_dofs[i] = g_dofs[i];
249 for (int j=i*m_nlay; j<(i+1)*m_nlay; j++)
250 {
251 m_sigm[j] = g_Sigm[j];
252 }
253 }
254
255 m_pMilleAlign = new Millepede();
256 m_pMilleAlign -> InitMille(&m_dofs[0], &m_sigm[0], m_nlay,m_nglo_on_lay, m_nloc,
258
259 m_derLC_Phi.resize(m_nloc);
260 m_derLC_V.resize(m_nloc);
261 m_derGB_Phi.resize(m_npar);
262 m_derGB_V.resize(m_npar);
263 m_derNonLin_Phi.resize(m_npar);
264 m_derNonLin_V.resize(m_npar);
265 m_par.resize(m_npar);
266 m_error.resize(m_npar);
267 m_pull.resize(m_npar);
268
269 log << MSG::INFO << "MilleAlign::initialize() finished!!" << endreq;
270
271 // contraints
272 std::vector<double> constDY1;
273 std::vector<double> constDY2;
274 std::vector<double> constDY3;
275 std::vector<double> constDY4;
276 std::vector<double> constDY5;
277 std::vector<double> constDY6;
278
279 std::vector<double> constDX5;
280 std::vector<double> constDZ5;
281 std::vector<double> constRX5;
282 std::vector<double> constRY5;
283 std::vector<double> constRZ5;
284
285 std::vector<double> constDX6;
286 std::vector<double> constDZ6;
287 std::vector<double> constRX6;
288 std::vector<double> constRY6;
289 std::vector<double> constRZ6;
290
291 std::vector<double> constDX1;
292 std::vector<double> constDZ1;
293 std::vector<double> constRX1;
294 std::vector<double> constRY1;
295 std::vector<double> constRZ1;
296
297 std::vector<double> constDX2;
298 std::vector<double> constDZ2;
299 std::vector<double> constRX2;
300 std::vector<double> constRY2;
301 std::vector<double> constRZ2;
302
303 std::vector<double> constDX3DX4;
304 std::vector<double> constDZ3DZ4;
305 std::vector<double> constRX3RX4;
306 std::vector<double> constRY3RY4;
307 std::vector<double> constRZ3RZ4;
308
309 std::vector<double> constRX3;
310 std::vector<double> constRY3;
311 std::vector<double> constRX4;
312 std::vector<double> constRY4;
313
314 // fix Dy of all layers
315 constDY1.resize(m_npar);
316 constDY2.resize(m_npar);
317 constDY3.resize(m_npar);
318 constDY4.resize(m_npar);
319 constDY5.resize(m_npar);
320 constDY6.resize(m_npar);
321
322 // fix Dx, Dz, Rz of the outmost layers
323 constDX5.resize(m_npar);
324 constDZ5.resize(m_npar);
325 constRX5.resize(m_npar);
326 constRY5.resize(m_npar);
327 constRZ5.resize(m_npar);
328
329 constDX6.resize(m_npar);
330 constDZ6.resize(m_npar);
331 constRX6.resize(m_npar);
332 constRY6.resize(m_npar);
333 constRZ6.resize(m_npar);
334
335 // fix Dx, Dz, Rz of the midle layers
336 constDX1.resize(m_npar);
337 constDZ1.resize(m_npar);
338 constRX1.resize(m_npar);
339 constRY1.resize(m_npar);
340 constRZ1.resize(m_npar);
341
342 constDX2.resize(m_npar);
343 constDZ2.resize(m_npar);
344 constRX2.resize(m_npar);
345 constRY2.resize(m_npar);
346 constRZ2.resize(m_npar);
347
348 constDX3DX4.resize(m_npar);
349 constDZ3DZ4.resize(m_npar);
350 constRX3RX4.resize(m_npar);
351 constRY3RY4.resize(m_npar);
352 constRZ3RZ4.resize(m_npar);
353
354 constRX3.resize(m_npar);
355 constRY3.resize(m_npar);
356 constRX4.resize(m_npar);
357 constRY4.resize(m_npar);
358
359 for(i=0; i<m_npar; i++){
360
361 constDY1[i] = 0.0;
362 constDY2[i] = 0.0;
363 constDY3[i] = 0.0;
364 constDY4[i] = 0.0;
365 constDY5[i] = 0.0;
366 constDY6[i] = 0.0;
367
368 constDX5[i] = 0.0;
369 constDZ5[i] = 0.0;
370 constRX5[i] = 0.0;
371 constRY5[i] = 0.0;
372 constRZ5[i] = 0.0;
373
374 constDX6[i] = 0.0;
375 constDZ6[i] = 0.0;
376 constRX6[i] = 0.0;
377 constRY6[i] = 0.0;
378 constRZ6[i] = 0.0;
379
380 constDX1[i] = 0.0;
381 constDZ1[i] = 0.0;
382 constRX1[i] = 0.0;
383 constRY1[i] = 0.0;
384 constRZ1[i] = 0.0;
385
386 constDX2[i] = 0.0;
387 constDZ2[i] = 0.0;
388 constRX2[i] = 0.0;
389 constRY2[i] = 0.0;
390 constRZ2[i] = 0.0;
391
392 constDX3DX4[i] = 0.0;
393 constDZ3DZ4[i] = 0.0;
394 constRX3RX4[i] = 0.0;
395 constRY3RY4[i] = 0.0;
396 constRZ3RZ4[i] = 0.0;
397
398 constRX3[i] = 0.0;
399 constRY3[i] = 0.0;
400 constRX4[i] = 0.0;
401 constRY4[i] = 0.0;
402
403 }
404
405 // Fix the position of 2nd and 3rd layer and dy in all 3 layers
406 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
407 // DX DY DZ RX RY RZ
408 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
409 // DX1 DX2 DX3 DX4 DX5 DX6 DY1 DY2 DY3 DY4 DY5 DY6 DZ1 DZ2 DZ3 DZ4 DZ5 DZ6 RX1 RX2 RX3 RX4 RX5 RX6 RY1 RY2 RY3 RY4 RY5 RY6 RZ1 RZ2 RZ3 RZ4 RZ5 RZ6
410 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
411 // index = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
412 // -------------------------- --------------------------- --------------------------- --------------------------- --------------------------- ---------------------------
413
414 // fix the DY on all layers
415 constDY1[6] = 1.0;
416 constDY2[7] = 1.0;
417 constDY3[8] = 1.0;
418 constDY4[9] = 1.0;
419 constDY5[10] = 1.0;
420 constDY6[11] = 1.0;
421
422 // fix the 3rd layer
423 constDX5[4] = 1.0;
424 constDZ5[16] = 1.0;
425 constRX5[22] = 1.0;
426 constRY5[28] = 1.0;
427 constRZ5[34] = 1.0;
428
429 constDX6[5] = 1.0;
430 constDZ6[17] = 1.0;
431 constRX6[23] = 1.0;
432 constRY6[29] = 1.0;
433 constRZ6[35] = 1.0;
434
435 // change to fix the 1st layer
436 constDX1[0] = 1.0;
437 constDZ1[12] = 1.0;
438 constRX1[18] = 1.0;
439 constRY1[24] = 1.0;
440 constRZ1[30] = 1.0;
441
442 constDX2[1] = 1.0;
443 constDZ2[13] = 1.0;
444 constRX2[19] = 1.0;
445 constRY2[25] = 1.0;
446 constRZ2[31] = 1.0;
447
448 // require the Dx3 = Dx4
449 constDX3DX4[2] = 1.0;
450 constDX3DX4[3] = -1.0;
451 // require the Dz3 = Dz4
452 constDZ3DZ4[14] = 1.0;
453 constDZ3DZ4[15] = -1.0;
454 // require the Rx3 = Rx4
455 constRX3RX4[20] = 1.0;
456 constRX3RX4[21] = -1.0;
457 constRX3[20] = 1.0;
458 constRX4[21] = 1.0;
459 // require the Ry3 = Ry4
460 constRY3RY4[26] = 1.0;
461 constRY3RY4[27] = -1.0;
462 constRY3[26] = 1.0;
463 constRY4[27] = 1.0;
464 // require the Rz3 = Rz4
465 constRZ3RZ4[32] = 1.0;
466 constRZ3RZ4[33] = -1.0;
467
468 // apply the constrains
469 m_pMilleAlign -> ConstF(&constDY1[0], 0.0);
470 m_pMilleAlign -> ConstF(&constDY2[0], 0.0);
471 m_pMilleAlign -> ConstF(&constDY3[0], 0.0);
472 m_pMilleAlign -> ConstF(&constDY4[0], 0.0);
473 m_pMilleAlign -> ConstF(&constDY5[0], 0.0);
474 m_pMilleAlign -> ConstF(&constDY6[0], 0.0);
475
476 m_pMilleAlign -> ConstF(&constDX5[0], 0.0);
477 m_pMilleAlign -> ConstF(&constDZ5[0], 0.0);
478 m_pMilleAlign -> ConstF(&constRX5[0], 0.0);
479 m_pMilleAlign -> ConstF(&constRY5[0], 0.0);
480 //m_pMilleAlign -> ConstF(&constRZ5[0], 0.0);
481
482 m_pMilleAlign -> ConstF(&constDX6[0], 0.0);
483 m_pMilleAlign -> ConstF(&constDZ6[0], 0.0);
484 m_pMilleAlign -> ConstF(&constRX6[0], 0.0);
485 m_pMilleAlign -> ConstF(&constRY6[0], 0.0);
486 //m_pMilleAlign -> ConstF(&constRZ6[0], 0.0);
487
488 m_pMilleAlign -> ConstF(&constDX1[0], 0.0);
489 m_pMilleAlign -> ConstF(&constDZ1[0], 0.0);
490 m_pMilleAlign -> ConstF(&constRX1[0], 0.0);
491 m_pMilleAlign -> ConstF(&constRY1[0], 0.0);
492 m_pMilleAlign -> ConstF(&constRZ1[0], 0.0);
493
494 m_pMilleAlign -> ConstF(&constDX2[0], 0.0);
495 m_pMilleAlign -> ConstF(&constDZ2[0], 0.0);
496 m_pMilleAlign -> ConstF(&constRX2[0], 0.0);
497 m_pMilleAlign -> ConstF(&constRY2[0], 0.0);
498 m_pMilleAlign -> ConstF(&constRZ2[0], 0.0);
499
500 m_pMilleAlign -> ConstF(&constDX3DX4[0], 0.0);
501 //m_pMilleAlign -> ConstF(&constDZ3DZ4[0], 0.0);
502
503 //m_pMilleAlign -> ConstF(&constRX3RX4[0], 0.0);
504 //m_pMilleAlign -> ConstF(&constRY3RY4[0], 0.0);
505
506 //m_pMilleAlign -> ConstF(&constRX3[0], 0.0);
507 //m_pMilleAlign -> ConstF(&constRX4[0], 0.0);
508 //m_pMilleAlign -> ConstF(&constRY3[0], 0.0);
509 //m_pMilleAlign -> ConstF(&constRY4[0], 0.0);
510
511 //m_pMilleAlign -> ConstF(&constRZ3RZ4[0], 0.0);
512
513}
514
516
517 IMessageSvc* msgSvc;
518 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
519 MsgStream log(msgSvc, "MilleAlign");
520 log << MSG::INFO << "MilleAlign::fillTree()" << endreq;
521
522 int recFlag;
523 int itrk;
524 int ihit;
525 int lay;
526 int Sheet_Meas;
527 int lay_virtual;
528 int Sheet_virtual;
529 double Phi_Meas;
530 double V_Meas;
531 double Z_Meas;
532 double resi_phi;
533 double resi_v;
534 double resi_phi_O;
535 double resi_v_O;
536 double resi_z;
537 double exp_phi;
538 double exp_v;
539 double exp_z;
540 double resi_comb;
541 double deri_Phi;
542 double deri_V;
543 int direction;
544 double hitSigm_Phi=0.0009;
545 double hitSigm_V=0.2239;
546 double lay_pos[6]; // 4 mis-alignment parameters for each layer
547
548 double trkpar[NTRKPAR];
549 double trkparms[NTRKPARALL]; // track parameters and errors
550
551 // numerical derivative
552 int ipar;
553 int iparGB;
554
555 CgemAliRecTrk* rectrk;
556 CgemAliRecHit* rechit;
557
558 HepPoint3D posUp, posDown;
559 double phiVUp[3], phiVDown[3];
560
561 HepPoint3D posUp_O, posDown_O;
562 double phiVUp_O[3], phiVDown_O[3];
563
564 int ntrk = event -> getNTrk();
565
566 log << MSG::INFO << "Number of tracks : "<<ntrk << endreq;
567
568 for(itrk=0; itrk<ntrk; itrk++){
569 rectrk = event->getRecTrk(itrk);
570 recFlag = rectrk->getStat();
571
572 trkpar[0] = rectrk -> getDr();
573 trkpar[1] = rectrk -> getPhi0();
574 trkpar[2] = rectrk -> getKappa();
575 trkpar[3] = rectrk -> getDz();
576 trkpar[4] = rectrk -> getTanLamda();
577 double chisq = rectrk -> getChisq();
578
579 int nHit = 0;
580 nHit = rectrk -> getNcluster(); //Why this number is two times larger than the inpur value ??? Guoaq -19June08
581
582 HepVector rechelix = rectrk->getHelix();
583 HepVector helix = rectrk->getHelix();
584 HepSymMatrix helixErr = rectrk->getHelixErr();
585
586 if(re_3lay&&nHit<12) continue; // only use the track can pass through all 3 layers of Cgem
587
588 int nHitUse = 0;
589 for(ihit=0; ihit<nHit/2; ihit++){
590 rechit = rectrk -> getRecHit(ihit);
591
592 lay = rechit -> getlayerid();
593 Phi_Meas = rechit -> getrecphi();
594 V_Meas = rechit -> getrecv();
595 Z_Meas = rechit -> getRecZ();
596 Sheet_Meas = rechit -> getsheetid();
597
598 if(lay==0) Sheet_virtual = (Phi_Meas<0)? 0:1;
599 else Sheet_virtual = Sheet_Meas;
600 lay_virtual = lay*2+Sheet_virtual;
601
602 lay_pos[0] = alignPar->getDx(lay_virtual);
603 lay_pos[1] = alignPar->getDy(lay_virtual);
604 lay_pos[2] = alignPar->getDz(lay_virtual);
605 lay_pos[3] = alignPar->getRx(lay_virtual);
606 lay_pos[4] = alignPar->getRy(lay_virtual);
607 lay_pos[5] = alignPar->getRz(lay_virtual);
608
609
610 if(m_GeoMethod==1)
611 {
612 // six or 3 in m_geoalign
613 m_geoalign->setDx(lay_virtual, lay_pos[0]);
614 m_geoalign->setDy(lay_virtual, lay_pos[1]);
615 m_geoalign->setDz(lay_virtual, lay_pos[2]);
616 m_geoalign->setRx(lay_virtual, lay_pos[3]);
617 m_geoalign->setRy(lay_virtual, lay_pos[4]);
618 m_geoalign->setRz(lay_virtual, lay_pos[5]);
619 m_middriftplane->setAlignment(m_geoalign);
620
621 StraightLine* trk = new StraightLine(trkpar[0], trkpar[1], trkpar[3], trkpar[4]);
622 m_middriftplane->getPointAligned(lay_virtual, trk, posUp, posDown, phiVUp, phiVDown);
623
624 phiVUp[0] = fmod((phiVUp[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
625 phiVDown[0] = fmod((phiVDown[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
626
627 phiVUp_O[0] = fmod((phiVUp_O[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
628 phiVDown_O[0] = fmod((phiVDown_O[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
629 Phi_Meas = fmod((Phi_Meas+2.0*TMath::Pi()),(2.0*TMath::Pi()));
630
631 resi_phi = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (Phi_Meas - phiVUp[0]) : (Phi_Meas - phiVDown[0]);
632 resi_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (V_Meas - phiVUp[1]) : (V_Meas - phiVDown[1]);
633
634 resi_phi_O = (fabs(Phi_Meas - phiVUp_O[0])<fabs(Phi_Meas - phiVDown_O[0]))? (Phi_Meas - phiVUp_O[0]) : (Phi_Meas - phiVDown_O[0]);
635 resi_v_O = (fabs(Phi_Meas - phiVUp_O[0])<fabs(Phi_Meas - phiVDown_O[0]))? (V_Meas - phiVUp_O[1]) : (V_Meas - phiVDown_O[1]);
636
637 resi_z = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (Z_Meas - posUp.z()) : (Z_Meas - posDown.z());
638
639 exp_phi = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[0]) : (phiVDown[0]);
640 exp_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[1]) : (phiVDown[1]);
641 exp_v = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (phiVUp[1]) : (phiVDown[1]);
642 exp_z = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (posUp.z()) : (posDown.z());
643
644 if(mul_R) resi_phi = resi_phi*m_r[lay_virtual];
645 if(mul_Ste_ang) resi_v = resi_v*sin(Ste_ang[lay_virtual]);
646
647 direction = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? -1 : 1;
648 delete trk;
649
650 }
651
652 // fill histograms
653 m_hresAll_phi->Fill(resi_phi);
654 m_hresAll_v->Fill(resi_v);
655 m_hresAll_z->Fill(resi_z);
656
657 m_hresLay_phi[lay_virtual]->Fill(resi_phi);
658 m_hresLay_v[lay_virtual]->Fill(resi_v);
659 m_hresLay_z[lay_virtual]->Fill(resi_z);
660
661 if (fabs(resi_phi) > m_resiCut_phi[lay_virtual]) continue;
662 if (fabs(resi_v) > m_resiCut_v[lay_virtual]) continue;
663
664 // reset the derivatives arrays
665 m_pMilleAlign -> ZerLoc(&m_derGB_Phi[0], &m_derLC_Phi[0], &m_derNonLin_Phi[0]);
666 m_pMilleAlign -> ZerLoc(&m_derGB_V[0], &m_derLC_V[0], &m_derNonLin_V[0]);
667
668 // derivatives of local parameters
669 for(ipar=0; ipar<m_nloc; ipar++){
670 if( ! getDeriLoc(nTrack_tot,ipar, lay_virtual, Sheet_Meas, direction, Phi_Meas, V_Meas, lay_pos, rechelix, helixErr, deri_Phi, deri_V) ){
671 cout << "getDeriLoc == false!" << setw(12) << itrk << setw(12) << ipar << endl;
672 return false;
673 }
674 if(debug)cout<<"deriLoc : ipar = "<<ipar<<", deri_Phi = "<<deri_Phi<<", deri_V = "<<deri_V<<endl;
675 m_derLC_Phi[ipar] = deri_Phi;
676 m_derLC_V[ipar] = deri_V;
677 }
678
679 // derivatives of global parameters
680 // ipar 0 - 3: dx, dy, dz, rz
681 for(ipar=0; ipar<m_nglo_on_lay; ipar++){
682 iparGB = getAlignParId(lay_virtual, ipar); // I need to use virtual layer here Guoaq/2020-09-08
683 if( ! getDeriGlo(nTrack_tot, ipar, iparGB, lay_virtual, Sheet_Meas, direction, Phi_Meas, V_Meas, lay_pos, helix, helixErr, deri_Phi, deri_V) )
684 {
685 cout << "getDeriGlo == false!" << setw(12) << itrk << setw(12) << ipar << endl;
686 return false;
687 }
688 m_derGB_Phi[iparGB] = deri_Phi;
689 m_derGB_V[iparGB] = deri_V;
690 if(debug)cout<<"deriGlo : ipar, iparGB, layer_virual = "<<ipar<<", "<<iparGB<<", "<<lay_virtual<<", deri_Phi = "<<deri_Phi<<", deri_V = "<<deri_V<<endl;
691 }
692 // Form one equation for each hit, what's that? How about this item? and this guy? guoaq-19-09-14
693
694 if(fabs(Phi_Meas-TMath::Pi())<0.025) continue;
695 if(fabs(Phi_Meas-0)<0.1) continue;
696 if(fabs(Phi_Meas-2.0*TMath::Pi())<0.025) continue;
697 if(fabs(fabs(trkpar[0])-m_r[lay_virtual])<10) continue;
698 if(fabs(fabs(Z_Meas)-m_z[lay_virtual])<5) continue;
699 if(m_derGB_V[getAlignParId(lay_virtual, 3)]<-1000||m_derGB_V[getAlignParId(lay_virtual, 3)]>1000) continue;
700 if(m_derGB_V[getAlignParId(lay_virtual, 4)]<-1000||m_derGB_V[getAlignParId(lay_virtual, 4)]>1000) continue;
701
702 double err_phi = 0;
703 double err_v = 0;
704
705 if(use_phi&&mul_R) err_phi = sig_phi[lay_virtual]*m_r[lay_virtual];
706 if(use_phi&&!mul_R) err_phi = sig_phi[lay_virtual];
707 if(use_v&&mul_Ste_ang) err_v = sig_v[lay_virtual]*sin(Ste_ang[lay_virtual]);
708 if(use_v&&!mul_Ste_ang) err_v = sig_v[lay_virtual];
709
710 m_pMilleAlign -> EquLoc(&m_derGB_Phi[0], &m_derLC_Phi[0], &m_derNonLin_Phi[0], resi_phi, err_phi);
711 m_pMilleAlign -> EquLoc(&m_derGB_V[0], &m_derLC_V[0], &m_derNonLin_V[0], resi_v, err_v);
712 nHitUse++;
713
714 } // loop of nhit
715
716 bool sc = m_pMilleAlign -> FitLoc(m_pMilleAlign->GetTrackNumber(), trkparms, 0);
717 if(sc) {m_pMilleAlign -> SetTrackNumber( m_pMilleAlign->GetTrackNumber()+1 );count[24]++;}
718
719 } // track loop
720
721 nTrack_tot++;
722 return true;
723}
724
725bool CgemMilleAlign::iniGeoAlignPar(){
726
727 for(int i=0; i<NLAYER; i++)
728 {
729
730 m_geoalign ->setDx(i,ini_alig_par[i][0]);
731 m_geoalign ->setDy(i,ini_alig_par[i][1]);
732 m_geoalign ->setDz(i,ini_alig_par[i][2]);
733 m_geoalign ->setRx(i,ini_alig_par[i][3]);
734 m_geoalign ->setRy(i,ini_alig_par[i][4]);
735 m_geoalign ->setRz(i,ini_alig_par[i][5]);
736
737 }
738
739
740}
742
743 IMessageSvc* msgSvc;
744 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
745 MsgStream log(msgSvc, "MilleAlign");
746 log << MSG::INFO << "MilleAlign::updateConst()" << endreq;
747
748 m_pMilleAlign -> MakeGlobalFit(&m_par[0], &m_error[0], &m_pull[0]);
749
750 int iLAYER;
751 int ipar;
752 double val;
753 double err;
754 for(int i=0; i<NDOFALIGN; i++){ //loop dx,dy,dz,rz,i
755 for(iLAYER=0; iLAYER<NLAYER; iLAYER++){ //loop 3 layer , iLAYER
756 ipar = i * NLAYER + iLAYER;
757
758 if(m_dofs[i]){
759 val = m_par[ipar];
760 err = m_error[ipar];
761 } else{
762 val = 0.0;
763 err = 0.0;
764 }
765 if(0 == i){
766 alignPar->setDelDx(iLAYER, val);
767 alignPar->setErrDx(iLAYER, err);
768 } else if(1 == i){
769 alignPar->setDelDy(iLAYER, val);
770 alignPar->setErrDy(iLAYER, err);
771 } else if(2 == i){
772 alignPar->setDelDz(iLAYER, val);
773 alignPar->setErrDz(iLAYER, err);
774 } else if(3 == i){
775 alignPar->setDelRx(iLAYER, val);
776 alignPar->setErrRx(iLAYER, err);
777 } else if(4 == i){
778 alignPar->setDelRy(iLAYER, val);
779 alignPar->setErrRy(iLAYER, err);
780 } else if(5 == i){
781 alignPar->setDelRz(iLAYER, val);
782 alignPar->setErrRz(iLAYER, err);
783 }
784
785 }
786 }
787 return true;
788}
789
790
791int CgemMilleAlign::getAlignParId(int lay, int iparHit){
792 int ipar = iparHit*6 + lay; // for 6 layer Guoaq/2020-09-08
793 return ipar;
794}
795
796bool CgemMilleAlign::getDeriLoc(int Intra, int ipar, int lay, int sheet, int dir, double mPhi, double mV, double laypos[], HepVector rechelix, HepSymMatrix &helixErr, double &deri_phi, double &deri_v){
797
798 // Defination of parameters
799 // ================================================
800 // Intra : index of track
801 // ipar : index of local paramters
802 // lay : id of geometry layer
803 // sheet: id of geometry sheet
804 // dir: direction of track: -1 using upper hit, 1 using bottom hit
805 // mPhi: measured phi
806 // mV: measured V
807 // laypos[]: alignment parameter of layer
808 // rechelix: reconstructed track helix parameters
809 // &heliErr: error of track parameters
810 // deri_phi: D[phi]/D[ipar] at measured phi
811 // deri_v: D[v]/D[ipar] at measured v
812 // ================================================
813 int layer_geo = int(lay/2);
814 int i;
815 double doca_all;
816 double radius;
817 double resi_phi;
818 double resi_v;
819 double exp_phi;
820 double exp_v;
821 double exp_z;
822 HepPoint3D posUp, posDown;
823 double phiVUp[3], phiVDown[3];
824 HepVector sampar(NTRKPAR, 0);
825 double xxLC_v[gNsamLC];
826 double xxLC[gNsamLC];
827 double yyLC_phi[gNsamLC];
828 double yyLC_v[gNsamLC];
829 double yyLC_z[gNsamLC];
830
831 CgemGeoReadoutPlane* m_readoutplane = m_cgemGeomSvc->getReadoutPlane(layer_geo, sheet);
832
833 // temp track para original track para
834 for(i=0; i<m_nloc; i++) sampar[i] = rechelix[i];
835 // get the start point of scan
836 double startpar = rechelix[ipar] - 0.5*gStepLC[ipar]*(double)gNsamLC;
837 for(i=0; i<gNsamLC; i++){
838 sampar[ipar] = startpar + (double)i * gStepLC[ipar];
839 xxLC[i] = sampar[ipar];
840
841 HepVector helix = sampar;
842
843 if(m_GeoMethod==1)
844 {
845 m_geoalign->setDx(lay, laypos[0]);
846 m_geoalign->setDy(lay, laypos[1]);
847 m_geoalign->setDz(lay, laypos[2]);
848 m_geoalign->setRx(lay, laypos[3]);
849 m_geoalign->setRy(lay, laypos[4]);
850 m_geoalign->setRz(lay, laypos[5]);
851
852 m_middriftplane->setAlignment(m_geoalign);
853
854 StraightLine* trk = new StraightLine(helix[0], helix[1],helix[3], helix[4]);
855 m_middriftplane->getPointAligned(lay, trk, posUp, posDown, phiVUp, phiVDown);
856
857 phiVUp[0] = fmod((phiVUp[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
858 phiVDown[0] = fmod((phiVDown[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
859
860 exp_phi = (dir==-1)? phiVUp[0]:phiVDown[0];
861 exp_v = (dir==-1)? phiVUp[1]:phiVDown[1];
862 exp_z = (dir==-1)? (posUp.z()) : (posDown.z());
863 resi_phi = (dir==-1)? phiVUp[0]:phiVDown[0] - mPhi;
864 resi_v = (dir==-1)? phiVUp[1]:phiVDown[1] - mV;
865 if(mul_R) resi_phi = resi_phi*m_r[lay];
866 if(mul_Ste_ang) resi_v = resi_v*sin(Ste_ang[lay]);
867
868 delete trk;
869
870 }
871 if(NULL == resi_phi||NULL==resi_v){
872 return false;
873 }
874 yyLC_phi[i] = exp_phi;
875 yyLC_v[i] = exp_v;
876 yyLC_z[i] = exp_z;
877 xxLC_v[i] = exp_v;
878 }
879
880 TSpline3* pSpline3_phi = new TSpline3("deri", xxLC, yyLC_phi, gNsamLC);
881 deri_phi = pSpline3_phi->Derivative(rechelix[ipar]);
882
883 TSpline3* pSpline3_v = new TSpline3("deri", xxLC, yyLC_v, gNsamLC);
884 deri_v = pSpline3_v->Derivative(rechelix[ipar]);
885
886 delete pSpline3_phi;
887 delete pSpline3_v;
888 return true;
889}
890
891bool CgemMilleAlign::getDeriGlo(int Intra, int iparHit, int iparGB, int lay, int sheet, int dir, double mPhi, double mV, double laypos[], HepVector helix,
892 HepSymMatrix &helixErr, double &deri_phi, double &deri_v){
893
894 // Defination of parameters
895 // ================================================
896 // Intra : index of track
897 // iparHit : index of alignment parameter for the layer contain the hit : [dx, dy, dz, rz]
898 // iparGB : index of global paramters, which need to be sorted properly !!!
899 // lay : id of geometry layer ? or virtual layer ?
900 // sheet: id of geometry sheet
901 // dir: direction of track: -1 using upper hit, 1 using bottom hit
902 // mPhi: measured phi
903 // mV: measured V
904 // laypos[]: alignment parameter of layer
905 // rechelix: reconstructed track helix parameters
906 // &heliErr: error of track parameters
907 // deri_phi: D[phi]/D[iparGB] at measured phi
908 // deri_v: D[v]/D[iparGB] at measured v
909 // ================================================
910
911 int layer_geo = int(lay/2);
912 int i;
913 double resi_phi;
914 double resi_v;
915 double exp_phi;
916 double exp_v;
917 double exp_z;
918 double radius;
919 HepPoint3D posUp, posDown;
920 double phiVUp[3], phiVDown[3];
921 double xxGB[gNsamGB];
922 double xxGB_v[gNsamGB];
923 double yyGB_phi[gNsamGB];
924 double yyGB_v[gNsamGB];
925 double yyGB_z[gNsamGB];
926 double dAlignPar;
927 double dAlignParini = laypos[iparHit];
928 // variation range interval step number
929 double startpar = dAlignParini - 0.5*gStepGB[iparGB]*(double)gNsamGB;
930 double layposSam[6];
931 CgemGeoReadoutPlane* m_readoutplane = m_cgemGeomSvc->getReadoutPlane(layer_geo, sheet);
932
933 // temp global para original global para.
934 for(i=0; i<6; i++) layposSam[i] = laypos[i]; // 0-2:dx,dy,dz; 3: theta_z
935
936 for(i=0; i<gNsamGB; i++){
937 dAlignPar = startpar + (double)i * gStepGB[iparGB];
938 xxGB[i] = dAlignPar;
939 if(0 == iparHit){ // dx
940 layposSam[0] = laypos[0] + dAlignPar;
941 } else if(1 == iparHit){ // dy
942 layposSam[1] = laypos[1] + dAlignPar;
943 } else if(2 == iparHit){ // dz
944 layposSam[2] = laypos[2] + dAlignPar;
945 } else if(3 == iparHit){ // rx
946 layposSam[3] = laypos[3] + dAlignPar;
947 } else if(4 == iparHit){ // ry
948 layposSam[4] = laypos[4] + dAlignPar;
949 } else if(5 == iparHit){ // rz
950 layposSam[5] = laypos[5] + dAlignPar;
951 }
952
953
954 if(m_GeoMethod==1)
955 {
956 m_geoalign->setDx(lay, layposSam[0]);
957 m_geoalign->setDy(lay, layposSam[1]);
958 m_geoalign->setDz(lay, layposSam[2]);
959 m_geoalign->setRx(lay, layposSam[3]);
960 m_geoalign->setRy(lay, layposSam[4]);
961 m_geoalign->setRz(lay, layposSam[5]);
962
963 m_middriftplane->setAlignment(m_geoalign);
964
965 StraightLine* trk = new StraightLine(helix[0], helix[1],helix[3], helix[4]);
966 m_middriftplane->getPointAligned(lay, trk, posUp, posDown, phiVUp, phiVDown);
967
968 phiVUp[0] = fmod((phiVUp[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
969 phiVDown[0] = fmod((phiVDown[0]+2.0*TMath::Pi()),(2.0*TMath::Pi()));
970
971 exp_phi = (dir==-1)? phiVUp[0]:phiVDown[0]; // global
972 exp_v = (dir==-1)? phiVUp[1]:phiVDown[1];
973 exp_z = (dir==-1)? (posUp.z()) : (posDown.z());
974
975 resi_phi = (dir==-1)? phiVUp[0]:phiVDown[0] - mPhi;
976 resi_v = (dir==-1)? phiVUp[1]:phiVDown[1] - mV;
977
978 if(mul_R) resi_phi = resi_phi*m_r[lay];
979 if(mul_Ste_ang) resi_v = resi_v*sin(Ste_ang[lay]);
980
981 delete trk;
982 }
983
984
985 if(NULL == resi_phi || NULL == resi_v) return false;
986
987 yyGB_phi[i] = exp_phi;
988 yyGB_v[i] = exp_v;
989 yyGB_z[i] = exp_z;
990 xxGB_v[i] = exp_v;
991 }
992
993 iniGeoAlignPar();
994
995 TSpline3* pSpline3_phi = new TSpline3("deri", xxGB, yyGB_phi, gNsamGB);
996 deri_phi = pSpline3_phi -> Derivative(dAlignParini);
997 TSpline3* pSpline3_v = new TSpline3("deri", xxGB, yyGB_v, gNsamGB);
998 deri_v = pSpline3_v -> Derivative(dAlignParini);
999
1000 delete pSpline3_v;
1001 delete pSpline3_phi;
1002 return true;
1003}
1004
1005
double sin(const BesAngle a)
Definition BesAngle.h:210
HepGeom::Point3D< double > HepPoint3D
IMessageSvc * msgSvc()
HepVector getHelix() const
int getStat() const
HepSymMatrix getHelixErr() const
virtual void initialize(TObjArray *hlist, ICgemGeomSvc *cgemGeomSvc, ICgemCalibFunSvc *cgemFunSvc)=0
void setErrRy(int iEP, double val)
void setDelDy(int iEP, double val)
void setErrDy(int iEP, double val)
void setErrRx(int iEP, double val)
double getDz(int iEP)
void setDelRy(int iEP, double val)
double getRy(int iEP)
void setDelDz(int iEP, double val)
void setDelRz(int iEP, double val)
void setErrDz(int iEP, double val)
double getDy(int iEP)
double getDx(int iEP)
void setErrDx(int iEP, double val)
void setErrRz(int iEP, double val)
void setDelDx(int iEP, double val)
double getRz(int iEP)
double getRx(int iEP)
void setDelRx(int iEP, double val)
double getDx(int layer_vir)
void setRy(int layer_vir, double v)
void setRz(int layer_vir, double v)
void setDz(int layer_vir, double v)
double getRx(int layer_vir)
double getRy(int layer_vir)
void setDy(int layer_vir, double v)
double getRz(int layer_vir)
double getDz(int layer_vir)
void setDx(int layer_vir, double v)
void setRx(int layer_vir, double v)
double getDy(int layer_vir)
void setAlignment(CgemGeoAlign *alignPtr)
bool getPointAligned(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
bool fillHist(CgemAliEvent *event, CgemAlignPar *alignPar)
void initialize(TObjArray *hlist, ICgemGeomSvc *cgemGeomSvc, ICgemCalibFunSvc *cgemFunSvc)
bool updateAlignPar(CgemAlignPar *alignPar)
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
virtual CgemGeoLayer * getCgemLayer(int i) const =0
virtual CgemGeoAlign * getAlignPtr() const =0
virtual CgemMidDriftPlane * getMidDriftPtr() const =0
const int NTRKPAR
const double g_res_cut
const int GEO_METHOD
const double gStepLC[5]
const int NLAYER
const double gStepGB[36]
const double g_start_chi_cut
const int NTRKPARALL
const double g_res_cut_init
const bool g_dofs[6]
const int NDOFALIGN
const double g_Sigm[36]
const int gNsamGB
const int gNsamLC
const int LAYERNMAX