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