CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
MilleAlign.cxx
Go to the documentation of this file.
1#include "MdcAlignAlg/MilleAlign.h"
2#include "MdcAlignAlg/MdcAlignAlg.h"
3#include "MdcAlignAlg/Alignment.h"
4
5#include "GaudiKernel/IMessageSvc.h"
6#include "GaudiKernel/StatusCode.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/Bootstrap.h"
9
10#include "Identifier/MdcID.h"
11
12#include <iostream>
13#include <fstream>
14#include <iomanip>
15#include <string>
16#include <cstring>
17#include <cmath>
18
19#ifndef ENABLE_BACKWARDS_COMPATIBILITY
20// backwards compatblty wll be enabled ONLY n CLHEP 1.9
21typedef HepGeom::Point3D<double> HepPoint3D;
22#endif
23
24#include "TSpline.h"
25
26using namespace std;
27
28MilleAlign::MilleAlign(){
29 for(int lay=0; lay<LAYERNMAX; lay++){
30 m_resiCut[lay] = 1.5;
31 if(lay<8){
32 m_docaCut[lay][0] = 0.5;
33 m_docaCut[lay][1] = 5.5;
34 } else{
35 m_docaCut[lay][0] = 0.5;
36 m_docaCut[lay][1] = 7.5;
37 }
38 }
39}
40
41MilleAlign::~MilleAlign(){
42}
43
44void MilleAlign::clear(){
45 delete m_hresAll;
46 delete m_hresInn;
47 delete m_hresStp;
48 delete m_hresOut;
49 for(int lay=0; lay<LAYERNMAX; lay++) delete m_hresLay[lay];
50 delete m_hresAllRec;
51 for(int lay=0; lay<LAYERNMAX; lay++) delete m_hresLayRec[lay];
52 delete m_hddoca;
53 delete m_pMilleAlign;
54}
55
56void MilleAlign::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
57 IMdcCalibFunSvc* mdcFunSvc){
58 IMessageSvc* msgSvc;
59 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
60 MsgStream log(msgSvc, "MilleAlign");
61 log << MSG::INFO << "MilleAlign::initialize()" << endreq;
62
63 // Initialze MdcUtilitySvc
64 IMdcUtilitySvc* imdcUtilitySvc;
65 StatusCode sc = Gaudi::svcLocator() -> service ("MdcUtilitySvc",imdcUtilitySvc);
66 m_mdcUtilitySvc= dynamic_cast<MdcUtilitySvc*> (imdcUtilitySvc);
67 if ( sc.isFailure() ){
68 log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endreq;
69 }
70
71 m_hlist = hlist;
72 m_mdcGeomSvc = mdcGeomSvc;
73 m_mdcFunSvc = mdcFunSvc;
74
75 // initialize hitograms
76 m_hresAll = new TH1F("HResAllInc", "", 200, -1.0, 1.0);
77 m_hlist->Add(m_hresAll);
78
79 m_hresInn = new TH1F("HResInnInc", "", 200, -1.0, 1.0);
80 m_hlist->Add(m_hresInn);
81
82 m_hresStp = new TH1F("HResStpInc", "", 200, -1.0, 1.0);
83 m_hlist->Add(m_hresStp);
84
85 m_hresOut = new TH1F("HResOutInc", "", 200, -1.0, 1.0);
86 m_hlist->Add(m_hresOut);
87
88 char hname[200];
89 for(int lay=0; lay<LAYERNMAX; lay++){
90 sprintf(hname, "Res_Layer%02d", lay);
91 m_hresLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
92 m_hlist->Add(m_hresLay[lay]);
93 }
94
95 m_hresAllRec = new TH1F("HResAllRecInc", "", 200, -1.0, 1.0);
96 m_hlist->Add(m_hresAllRec);
97 for(int lay=0; lay<LAYERNMAX; lay++){
98 sprintf(hname, "Res_LayerRec%02d", lay);
99 m_hresLayRec[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
100 m_hlist->Add(m_hresLayRec[lay]);
101 }
102
103 // for debug
104 m_hddoca = new TH1F("delt_doca", "", 200, -1.0, 1.0);
105 m_hlist->Add(m_hddoca);
106
107 for(int lay=0; lay<LAYERNMAX; lay++){
108 sprintf(hname, "delt_docaLay%02d", lay);
109 m_hddocaLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
110 m_hlist->Add(m_hddocaLay[lay]);
111 }
112
113 // initialize millepede
114 m_nglo = NEP;
115 m_nloc = NTRKPAR;
116 m_nGloHit = 2 * NDOFALIGN;
117 m_npar = NDOFALIGN * m_nglo;
118
119 int i;
120 for(i=0; i<NDOFALIGN; i++){
121 m_dofs[i] = g_dofs[i];
122 m_sigm[i] = g_Sigm[i];
123 }
124
125 m_pMilleAlign = new Millepede();
126 m_pMilleAlign -> InitMille(&m_dofs[0], &m_sigm[0], m_nglo, m_nloc,
127 g_start_chi_cut, 3, g_res_cut, g_res_cut_init);
128
129 m_derGB.resize(m_npar);
130 m_derNonLin.resize(m_npar);
131 m_par.resize(m_npar);
132 m_error.resize(m_npar);
133 m_pull.resize(m_npar);
134
135 m_derLC.resize(m_nloc);
136
137 // contraints
138 std::vector<double> constTX;
139 std::vector<double> constTY;
140 std::vector<double> constRZ;
141
142 std::vector<double> constTXE;
143 std::vector<double> constTXW;
144 std::vector<double> constTYE;
145 std::vector<double> constTYW;
146 std::vector<double> constRZE;
147 std::vector<double> constRZW;
148
149 constTX.resize(m_npar);
150 constTY.resize(m_npar);
151 constRZ.resize(m_npar);
152
153 constTXE.resize(m_npar);
154 constTXW.resize(m_npar);
155 constTYE.resize(m_npar);
156 constTYW.resize(m_npar);
157 constRZE.resize(m_npar);
158 constRZW.resize(m_npar);
159
160 for(i=0; i<m_npar; i++){
161 constTX[i] = 0.0;
162 constTY[i] = 0.0;
163 constRZ[i] = 0.0;
164
165 constTXE[i] = 0.0;
166 constTXW[i] = 0.0;
167 constTYE[i] = 0.0;
168 constTYW[i] = 0.0;
169 constRZE[i] = 0.0;
170 constRZW[i] = 0.0;
171 }
172 constTX[7] = 1.0;
173 constTX[15] = 1.0;
174 constTY[23] = 1.0;
175 constTY[31] = 1.0;
176 constRZ[39] = 1.0;
177 constRZ[47] = 1.0;
178
179 constTXE[7] = 1.0;
180 constTXW[15] = 1.0;
181 constTYE[23] = 1.0;
182 constTYW[31] = 1.0;
183 constRZE[39] = 1.0;
184 constRZW[47] = 1.0;
185
186 //m_pMilleAlign -> ConstF(&constTX[0], 0.0);
187 //m_pMilleAlign -> ConstF(&constTY[0], 0.0);
188// m_pMilleAlign -> ConstF(&constRZ[0], 0.0);
189
190 m_pMilleAlign -> ConstF(&constTXE[0], 0.0);
191 m_pMilleAlign -> ConstF(&constTXW[0], 0.0);
192 m_pMilleAlign -> ConstF(&constTYE[0], 0.0);
193 m_pMilleAlign -> ConstF(&constTYW[0], 0.0);
194 m_pMilleAlign -> ConstF(&constRZE[0], 0.0);
195 m_pMilleAlign -> ConstF(&constRZW[0], 0.0);
196}
197
198bool MilleAlign::fillHist(MdcAliEvent* event){
199 IMessageSvc* msgSvc;
200 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
201 MsgStream log(msgSvc, "MilleAlign");
202 log << MSG::DEBUG << "MilleAlign::fillTree()" << endreq;
203
204 int recFlag;
205 int itrk;
206 int ihit;
207 int fgGetDoca;
208 int lr;
209 int lay;
210 int cel;
211 double doca;
212 double dmeas;
213 double zhit;
214 double resi;
215 double resiRec;
216 double deri;
217 double hitSigm;
218 double hitpos[3];
219 double wpos[7]; // wpos[6] is wire tension
220 const MdcGeoWire* pWire;
221
222 double trkpar[NTRKPAR];
223 double trkparms[NTRKPARALL]; // track parameters and errors
224
225 // numerical derivative
226 int ipar;
227 int iparGB;
228
229 MdcAliRecTrk* rectrk;
230 MdcAliRecHit* rechit;
231
232 int ntrk = event -> getNTrk();
233 if( (ntrk<m_param.nTrkCut[0]) || (ntrk>m_param.nTrkCut[1])) return false;
234
235 for(itrk=0; itrk<ntrk; itrk++){
236 rectrk = event->getRecTrk(itrk);
237 recFlag = rectrk->getStat();
238
239 trkpar[0] = rectrk -> getDr();
240 trkpar[1] = rectrk -> getPhi0();
241 trkpar[2] = rectrk -> getKappa();
242 trkpar[3] = rectrk -> getDz();
243 trkpar[4] = rectrk -> getTanLamda();
244
245 int nHit = rectrk -> getNHits();
246 if(nHit < m_param.nHitCut) continue;
247 if(fabs(trkpar[0]) > m_param.drCut) continue;
248 if(fabs(trkpar[3]) > m_param.dzCut) continue;
249
250 HepVector rechelix = rectrk->getHelix();
251 HepVector helix = rectrk->getHelix();
252 HepSymMatrix helixErr = rectrk->getHelixErr();
253
254 int nHitUse = 0;
255 for(ihit=0; ihit<nHit; ihit++){
256 rechit = rectrk -> getRecHit(ihit);
257 lr = rechit->getLR();
258 lay = rechit -> getLayid();
259 cel = rechit -> getCellid();
260 pWire = m_mdcGeomSvc -> Wire(lay, cel);
261 dmeas = rechit -> getDmeas();
262 zhit = rechit -> getZhit();
263 hitSigm = rechit -> getErrDmeas();
264
265 wpos[0] = pWire -> Backward().x(); // east end
266 wpos[1] = pWire -> Backward().y();
267 wpos[2] = pWire -> Backward().z();
268 wpos[3] = pWire -> Forward().x(); // west end
269 wpos[4] = pWire -> Forward().y();
270 wpos[5] = pWire -> Forward().z();
271 wpos[6] = pWire -> Tension();
272
273 double docaRec = rechit->getDocaInc();
274 double doca = (m_mdcUtilitySvc->doca(lay, cel, helix, helixErr))*10.0;
275
276 resi = -1.0*dmeas - doca;
277 if ((fabs(doca) < m_docaCut[lay][0]) || (fabs(doca) > m_docaCut[lay][1]) ||
278 (fabs(resi) > m_resiCut[lay])) continue;
279 nHitUse++;
280
281 resiRec = rechit -> getResiIncLR();
282
283 double dd = fabs(doca) - fabs(rechit->getDocaInc());
284 m_hddoca -> Fill(dd);
285 m_hddocaLay[lay] -> Fill(dd);
286
287 // fill histograms
288 m_hresAll->Fill(resi);
289 if(lay < 8) m_hresInn->Fill(resi);
290 else if(lay < 20) m_hresStp->Fill(resi);
291 else m_hresOut->Fill(resi);
292 m_hresLay[lay]->Fill(resi);
293
294 m_hresAllRec->Fill(resiRec);
295 m_hresLayRec[lay]->Fill(resiRec);
296
297 // reset the derivatives arrays
298 m_pMilleAlign -> ZerLoc(&m_derGB[0], &m_derLC[0], &m_derNonLin[0]);
299
300 // derivatives of local parameters
301 for(ipar=0; ipar<m_nloc; ipar++){
302 if( ! getDeriLoc(ipar, lay, cel ,rechelix, helixErr, deri) ){
303 cout << "getDeriLoc == false!" << setw(12) << itrk << setw(12) << ipar << endl;
304 return false;
305 }
306 m_derLC[ipar] = deri;
307 }
308
309 // derivatives of global parameters
310 // ipar 0 - 5: dx_east, dx_west, dy_east, dy_west, rz_east, rz_west
311 for(ipar=0; ipar<m_nGloHit; ipar++){
312 iparGB = getAlignParId(lay, ipar);
313 if( ! getDeriGlo(ipar, iparGB, lay, cel, helix, helixErr, wpos, deri) )
314 {
315 cout << "getDeriGlo == false!" << setw(12) << itrk << setw(12) << ipar << endl;
316 return false;
317 }
318 m_derGB[iparGB] = deri;
319 }
320 m_pMilleAlign -> EquLoc(&m_derGB[0], &m_derLC[0], &m_derNonLin[0], resi, hitSigm);
321 } // loop of nhit
322
323 // local fit in Millepede
324 bool sc = m_pMilleAlign -> FitLoc(m_pMilleAlign->GetTrackNumber(), trkparms, 0);
325 if(sc) m_pMilleAlign -> SetTrackNumber( m_pMilleAlign->GetTrackNumber()+1 );
326 } // track loop
327 return true;
328}
329
330
331void MilleAlign::updateConst(MdcAlignPar* alignPar){
332 IMessageSvc* msgSvc;
333 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
334 MsgStream log(msgSvc, "MilleAlign");
335 log << MSG::INFO << "MilleAlign::updateConst()" << endreq;
336
337 m_pMilleAlign -> MakeGlobalFit(&m_par[0], &m_error[0], &m_pull[0]);
338
339 int iEP;
340 int ipar;
341 double val;
342 double err;
343 for(int i=0; i<NDOFALIGN; i++){
344 for(iEP=0; iEP<NEP; iEP++){
345 ipar = i * NEP + iEP;
346 if(m_dofs[i]){
347 val = m_par[ipar];
348 err = m_error[ipar];
349 } else{
350 val = 0.0;
351 err = 0.0;
352 }
353
354 if(0 == i){
355 alignPar->setDelDx(iEP, val);
356 alignPar->setErrDx(iEP, err);
357 } else if(1 == i){
358 alignPar->setDelDy(iEP, val);
359 alignPar->setErrDy(iEP, err);
360 } else if(2 == i){
361 alignPar->setDelRz(iEP, val/1000.0); // mrad -> rad
362 alignPar->setErrRz(iEP, err/1000.0); // mrad -> rad
363 }
364 }
365 }
366}
367
368int MilleAlign::getAlignParId(int lay, int iparHit){
369 int ip;
370 if(lay < 8) ip = 0;
371 else if(lay < 10) ip = 1;
372 else if(lay < 12) ip = 2;
373 else if(lay < 14) ip = 3;
374 else if(lay < 16) ip = 4;
375 else if(lay < 18) ip = 5;
376 else if(lay < 20) ip = 6;
377 else ip = 7;
378
379 // iparHit 0 - 5: dx_east, dx_west, dy_east, dy_west, rz_east, rz_west
380 int ipar = iparHit * 8 + ip;
381 return ipar;
382}
383
384bool MilleAlign::getDeriLoc(int ipar, int lay, int cel, HepVector rechelix, HepSymMatrix &helixErr, double &deri){
385 int i;
386 double doca;
387 HepVector sampar(NTRKPAR, 0);
388 double xxLC[gNsamLC];
389 double yyLC[gNsamLC];
390
391 for(i=0; i<m_nloc; i++) sampar[i] = rechelix[i];
392 double startpar = rechelix[ipar] - 0.5*gStepLC[ipar]*(double)gNsamLC;
393
394 for(i=0; i<gNsamLC; i++){
395 sampar[ipar] = startpar + (double)i * gStepLC[ipar];
396 xxLC[i] = sampar[ipar];
397 if(0==ipar || 3==ipar) xxLC[i] *= 10.; // cm -> mm
398
399 HepVector helix = sampar;
400 bool passCellRequired = false;
401 doca = (m_mdcUtilitySvc->doca(lay, cel, helix, helixErr,passCellRequired))*10.0;
402
403 if(NULL == doca){
404// cout << "in getDeriLoc, doca = " << doca << endl;
405 return false;
406 }
407 yyLC[i] = doca;
408 }
409
410 if (0==ipar || 3==ipar) rechelix[ipar] *= 10.; // cm -> mm
411 TSpline3* pSpline3 = new TSpline3("deri", xxLC, yyLC, gNsamLC);
412 deri = pSpline3->Derivative(rechelix[ipar]);
413 delete pSpline3;
414 return true;
415}
416
417bool MilleAlign::getDeriGlo(int iparHit, int iparGB, int lay, int cel, HepVector helix,
418 HepSymMatrix &helixErr, double wpos[], double &deri){
419 int i;
420 double doca;
421 double xxGB[gNsamGB];
422 double yyGB[gNsamGB];
423 double dAlignPar;
424 double dAlignParini = 0.0;
425 double startpar = dAlignParini - 0.5*gStepGB[iparGB]*(double)gNsamGB;
426 double wposSam[7];
427 for(i=0; i<7; i++) wposSam[i] = wpos[i]; // 0-2:east; 3-5:west
428
429 for(i=0; i<gNsamGB; i++){
430 dAlignPar = startpar + (double)i * gStepGB[iparGB];
431 xxGB[i] = dAlignPar;
432 if(0 == iparHit){ // dx_east
433 wposSam[0] = wpos[0] + dAlignPar;
434 } else if(1 == iparHit){ // dx_west
435 wposSam[3] = wpos[3] + dAlignPar;
436 } else if(2 == iparHit){ // dy_east
437 wposSam[1] = wpos[1] + dAlignPar;
438 } else if(3 == iparHit){ // dy_west
439 wposSam[4] = wpos[4] + dAlignPar;
440 } else if(4 == iparHit){ // rz_east
441 wposSam[0] = wpos[0] - (wpos[1] * dAlignPar * 0.001);
442 wposSam[1] = wpos[1] + (wpos[0] * dAlignPar * 0.001);
443 } else if(5 == iparHit){ // rz_west
444 wposSam[3] = wpos[3] - (wpos[4] * dAlignPar * 0.001);
445 wposSam[4] = wpos[4] + (wpos[3] * dAlignPar * 0.001);
446 }
447
448 HepPoint3D eastP(wposSam[0]/10., wposSam[1]/10., wposSam[2]/10.);
449 HepPoint3D westP(wposSam[3]/10., wposSam[4]/10., wposSam[5]/10.);
450 doca = (m_mdcUtilitySvc->doca(lay, cel, eastP, westP, helix, helixErr))*10.0; // cm->mm
451
452 if(NULL == doca) return false;
453
454 yyGB[i] = doca;
455 }
456
457 TSpline3* pSpline3 = new TSpline3("deri", xxGB, yyGB, gNsamGB);
458 deri = pSpline3 -> Derivative(dAlignParini);
459 delete pSpline3;
460 return true;
461}
462
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: MilleAlign.cxx:21
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
MdcAliParams::dzCut
double dzCut
Definition: MdcAliParams.h:29
MdcAliParams::drCut
double drCut
Definition: MdcAliParams.h:28
MdcAliParams::nTrkCut
int nTrkCut[2]
Definition: MdcAliParams.h:21
MdcAliRecHit::getLR
int getLR() const
Definition: MdcAliRecHit.h:22
MdcAliRecHit::getDocaInc
double getDocaInc() const
Definition: MdcAliRecHit.h:24
MdcAliRecTrk::getStat
int getStat() const
Definition: MdcAliRecTrk.h:26
MdcAliRecTrk::getHelixErr
HepSymMatrix getHelixErr() const
Definition: MdcAliRecTrk.h:33
MdcAliRecTrk::getHelix
HepVector getHelix() const
Definition: MdcAliRecTrk.h:32
MdcAlignPar::setErrRz
void setErrRz(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:40
MdcAlignPar::setDelDy
void setDelDy(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:29
MdcAlignPar::setDelRz
void setDelRz(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:33
MdcAlignPar::setErrDx
void setErrDx(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:35
MdcAlignPar::setDelDx
void setDelDx(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:28
MdcAlignPar::setErrDy
void setErrDy(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:36
MdcUtilitySvc::doca
double doca(int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const
Definition: MdcUtilitySvc.cxx:216
MilleAlign::clear
void clear()
Definition: MilleAlign.cxx:44
MilleAlign::initialize
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc)
Definition: MilleAlign.cxx:56
MilleAlign::fillHist
bool fillHist(MdcAliEvent *event)
Definition: MilleAlign.cxx:198
MilleAlign::updateConst
void updateConst(MdcAlignPar *alignPar)
Definition: MilleAlign.cxx:331
Alignment::g_res_cut_init
const double g_res_cut_init
Definition: Alignment.h:84
Alignment::gNsamGB
const int gNsamGB
Definition: Alignment.h:88
Alignment::LAYERNMAX
const int LAYERNMAX
Definition: Alignment.h:45
Alignment::gStepGB
const double gStepGB[48]
Definition: Alignment.h:91
Alignment::NEP
const int NEP
Definition: Alignment.h:48
Alignment::NTRKPAR
const int NTRKPAR
Definition: Alignment.h:49
Alignment::gStepLC
const double gStepLC[5]
Definition: Alignment.h:90
Alignment::NTRKPARALL
const int NTRKPARALL
Definition: Alignment.h:50
Alignment::g_res_cut
const double g_res_cut
Definition: Alignment.h:83
Alignment::gNsamLC
const int gNsamLC
Definition: Alignment.h:87
Alignment::NDOFALIGN
const int NDOFALIGN
Definition: Alignment.h:62
Alignment::g_Sigm
const double g_Sigm[3]
Definition: Alignment.h:77
Alignment::g_dofs
const bool g_dofs[3]
Definition: Alignment.h:71
Alignment::g_start_chi_cut
const double g_start_chi_cut
Definition: Alignment.h:85