BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
MilleAlign Class Reference

#include <MilleAlign.h>

+ Inheritance diagram for MilleAlign:

Public Member Functions

 MilleAlign ()
 
 ~MilleAlign ()
 
void clear ()
 
void initialize (TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
 
void setParam (MdcAliParams &param)
 
bool fillHist (MdcAliEvent *event)
 
void updateConst (MdcAlignPar *alignPar)
 
- Public Member Functions inherited from MdcAlign
 MdcAlign ()
 
virtual ~MdcAlign ()
 

Public Attributes

std::string fixMomLab
 
- Public Attributes inherited from MdcAlign
std::string fixMomLab
 

Detailed Description

Definition at line 15 of file MilleAlign.h.

Constructor & Destructor Documentation

◆ MilleAlign()

MilleAlign::MilleAlign ( )

Definition at line 28 of file MilleAlign.cxx.

28 {
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}
const int LAYERNMAX
Definition Alignment.h:45

◆ ~MilleAlign()

MilleAlign::~MilleAlign ( )

Definition at line 41 of file MilleAlign.cxx.

41 {
42}

Member Function Documentation

◆ clear()

void MilleAlign::clear ( )
virtual

Implements MdcAlign.

Definition at line 44 of file MilleAlign.cxx.

44 {
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}

◆ fillHist()

bool MilleAlign::fillHist ( MdcAliEvent * event)
virtual

Implements MdcAlign.

Definition at line 191 of file MilleAlign.cxx.

191 {
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}
curve Fill()
IMessageSvc * msgSvc()
virtual double doca(int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const =0
int getLR() const
double getDocaInc() const
int getStat() const
HepSymMatrix getHelixErr() const
HepVector getHelix() const
int GetTrackNumber()
const int NTRKPAR
Definition Alignment.h:49
const int NTRKPARALL
Definition Alignment.h:50

◆ initialize()

void MilleAlign::initialize ( TObjArray * hlist,
IMdcGeomSvc * mdcGeomSvc,
IMdcCalibFunSvc * mdcFunSvc,
IMdcUtilitySvc * mdcUtilitySvc )
virtual

Implements MdcAlign.

Definition at line 56 of file MilleAlign.cxx.

57 {
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,
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}
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)
const double g_res_cut_init
Definition Alignment.h:84
const int NEP
Definition Alignment.h:48
const double g_res_cut
Definition Alignment.h:83
const int NDOFALIGN
Definition Alignment.h:62
const double g_Sigm[3]
Definition Alignment.h:77
const bool g_dofs[3]
Definition Alignment.h:71
const double g_start_chi_cut
Definition Alignment.h:85

◆ setParam()

void MilleAlign::setParam ( MdcAliParams & param)
inlinevirtual

Implements MdcAlign.

Definition at line 81 of file MilleAlign.h.

81 {
82 MdcAlign::setParam(param);
83 m_param = param;
84}
virtual void setParam(MdcAliParams &param)=0
Definition MdcAlign.h:36

◆ updateConst()

void MilleAlign::updateConst ( MdcAlignPar * alignPar)
virtual

Implements MdcAlign.

Definition at line 323 of file MilleAlign.cxx.

323 {
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}
void setErrRz(int iEP, double val)
void setDelDy(int iEP, double val)
void setDelRz(int iEP, double val)
void setErrDx(int iEP, double val)
void setDelDx(int iEP, double val)
void setErrDy(int iEP, double val)

Member Data Documentation

◆ fixMomLab

std::string MilleAlign::fixMomLab

Definition at line 27 of file MilleAlign.h.


The documentation for this class was generated from the following files: