64 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
65 MsgStream log(
msgSvc,
"MilleAlign");
66 log << MSG::INFO <<
"MilleAlign::initialize()" << endreq;
70 m_cgemGeomSvc = cgemGeomSvc;
71 m_cgemFunSvc = cgemFunSvc;
73 for(
int il=0; il<6; il++)
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);
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];
119 for(
int i=0; i<30; i++)
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;
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;
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;
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;
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;
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;
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;
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;
170 m_hresAll_phi =
new TH1F(
"Resi_All_Phi",
"", 200, -2.0, 2.0);
171 m_hlist->Add(m_hresAll_phi);
173 m_hresAll_v =
new TH1F(
"Resi_All_v",
"", 100, -30.0, 30.0);
174 m_hlist->Add(m_hresAll_v);
176 m_hresAll_z =
new TH1F(
"Resi_All_z",
"", 100, -30.0, 30.0);
177 m_hlist->Add(m_hresAll_z);
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]);
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]);
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]);
198 m_hddoca =
new TH1F(
"delt_doca",
"", 200, -1.0, 1.0);
199 m_hlist->Add(m_hddoca);
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]);
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;
233 for(
int i=0; i<
NLAYER; i++)
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);
249 for (
int j=i*m_nlay; j<(i+1)*m_nlay; j++)
256 m_pMilleAlign -> InitMille(&m_dofs[0], &m_sigm[0], m_nlay,m_nglo_on_lay, m_nloc,
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);
269 log << MSG::INFO <<
"MilleAlign::initialize() finished!!" << endreq;
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;
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;
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;
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;
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;
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;
309 std::vector<double> constRX3;
310 std::vector<double> constRY3;
311 std::vector<double> constRX4;
312 std::vector<double> constRY4;
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);
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);
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);
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);
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);
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);
354 constRX3.resize(m_npar);
355 constRY3.resize(m_npar);
356 constRX4.resize(m_npar);
357 constRY4.resize(m_npar);
359 for(i=0; i<m_npar; i++){
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;
449 constDX3DX4[2] = 1.0;
450 constDX3DX4[3] = -1.0;
452 constDZ3DZ4[14] = 1.0;
453 constDZ3DZ4[15] = -1.0;
455 constRX3RX4[20] = 1.0;
456 constRX3RX4[21] = -1.0;
460 constRY3RY4[26] = 1.0;
461 constRY3RY4[27] = -1.0;
465 constRZ3RZ4[32] = 1.0;
466 constRZ3RZ4[33] = -1.0;
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);
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);
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);
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);
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);
500 m_pMilleAlign -> ConstF(&constDX3DX4[0], 0.0);
518 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
519 MsgStream log(
msgSvc,
"MilleAlign");
520 log << MSG::INFO <<
"MilleAlign::fillTree()" << endreq;
544 double hitSigm_Phi=0.0009;
545 double hitSigm_V=0.2239;
559 double phiVUp[3], phiVDown[3];
562 double phiVUp_O[3], phiVDown_O[3];
564 int ntrk =
event -> getNTrk();
566 log << MSG::INFO <<
"Number of tracks : "<<ntrk << endreq;
568 for(itrk=0; itrk<ntrk; itrk++){
569 rectrk =
event->getRecTrk(itrk);
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();
580 nHit = rectrk -> getNcluster();
582 HepVector rechelix = rectrk->
getHelix();
583 HepVector helix = rectrk->
getHelix();
586 if(re_3lay&&nHit<12)
continue;
589 for(ihit=0; ihit<nHit/2; ihit++){
590 rechit = rectrk -> getRecHit(ihit);
592 lay = rechit -> getlayerid();
593 Phi_Meas = rechit -> getrecphi();
594 V_Meas = rechit -> getrecv();
595 Z_Meas = rechit -> getRecZ();
596 Sheet_Meas = rechit -> getsheetid();
598 if(lay==0) Sheet_virtual = (Phi_Meas<0)? 0:1;
599 else Sheet_virtual = Sheet_Meas;
600 lay_virtual = lay*2+Sheet_virtual;
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);
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]);
622 m_middriftplane->
getPointAligned(lay_virtual, trk, posUp, posDown, phiVUp, phiVDown);
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()));
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()));
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]);
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]);
637 resi_z = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? (Z_Meas - posUp.z()) : (Z_Meas - posDown.z());
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());
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]);
647 direction = (fabs(Phi_Meas - phiVUp[0])<fabs(Phi_Meas - phiVDown[0]))? -1 : 1;
653 m_hresAll_phi->Fill(resi_phi);
654 m_hresAll_v->Fill(resi_v);
655 m_hresAll_z->Fill(resi_z);
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);
661 if (fabs(resi_phi) > m_resiCut_phi[lay_virtual])
continue;
662 if (fabs(resi_v) > m_resiCut_v[lay_virtual])
continue;
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]);
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;
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;
681 for(ipar=0; ipar<m_nglo_on_lay; ipar++){
682 iparGB = getAlignParId(lay_virtual, ipar);
683 if( ! getDeriGlo(nTrack_tot, ipar, iparGB, lay_virtual, Sheet_Meas, direction, Phi_Meas, V_Meas, lay_pos, helix, helixErr, deri_Phi, deri_V) )
685 cout <<
"getDeriGlo == false!" << setw(12) << itrk << setw(12) << ipar << endl;
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;
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;
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];
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);
716 bool sc = m_pMilleAlign -> FitLoc(m_pMilleAlign->
GetTrackNumber(), trkparms, 0);
717 if(sc) {m_pMilleAlign -> SetTrackNumber( m_pMilleAlign->
GetTrackNumber()+1 );count[24]++;}