CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
DotsHelixFitter.cxx
Go to the documentation of this file.
1#include "DotsConnection/DotsHelixFitter.h"
2#include "DotsConnection/TrkFitFun.h"
3//#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/Bootstrap.h"
5#include "RawEvent/RawDataUtil.h"
6#include "Identifier/MdcID.h"
7#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
8#include "MdcGeom/Constants.h"
9
10//using namespace TrkFitFun;
11
12double TrkFitFun::BFIELD;
13
14DotsHelixFitter::DotsHelixFitter():
15 myHelix(NULL),myIsInitialized(false),myMdcGeomSvc(NULL),myCgemGeomSvc(NULL),myMdcCalibFunSvc(NULL),myMdcUtilitySvc(NULL)
16{
17 myHelix_Ea=HepSymMatrix(5,0);
18 myMaxIteration=10;
19 myMinXhitsInCircleFit=3;
20 myMinVhitsInHelixFit=2;
21 myMinHitsInHelixFit=5;
22 myDchi2Converge=5.0;
23 myDchi2Diverge=50.0;
24 myMaxNChi2Increase=2;
25 myChi2Diverge=1000000.0;
26 myFitCircleOnly=false;
27 myFitCircleOriginOnly=false;
28 myUseAxialHitsOnly=false;
29 myUseMdcHitsOnly=false;
30}
31
32
33DotsHelixFitter::~DotsHelixFitter()
34{
35 if(myHelix) delete myHelix;
36}
37
38void DotsHelixFitter::initialize()
39{
40 if(myIsInitialized) return;
41
42 // --- get CgemGeomSvc ---
43 ISvcLocator* svcLocator = Gaudi::svcLocator();
44 ICgemGeomSvc* ISvc;
45 StatusCode Cgem_sc=svcLocator->service("CgemGeomSvc", ISvc);
46 if(Cgem_sc.isSuccess()) myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
47 else cout<<"DotsHelixFitter::initialize(): can not get CgemGeomSvc"<<endl;
48 int nLayerCgem = myCgemGeomSvc->getNumberOfCgemLayer();
49 if(nLayerCgem!=3) {
50 cout<<"DotsHelixFitter::initialize(): nLayerCgem!=3 !!!"<<endl;
51 return;
52 }
53 for(int i=0; i<nLayerCgem; i++)
54 {
55 //myCgemLayer[i]=myCgemGeomSvc->getCgemLayer(i);
56 CgemGeoLayer* CgemLayer = myCgemGeomSvc->getCgemLayer(i);
57 myRmidDGapCgem[i]=(CgemLayer->getMiddleROfGapD())/10.0;//cm
58 myR2midDGapCgem[i]=pow(myRmidDGapCgem[i],2);
59 myRVCgem[i]=(CgemLayer->getInnerROfAnodeCu1())/10.0;//cm
60 myRXCgem[i]=(CgemLayer->getInnerROfAnodeCu2())/10.0;//cm
61 myAngStereoCgem[i]=(CgemLayer->getAngleOfStereo())/180.0*CLHEP::pi;// degree -> radians
62 myNSheets[i]=CgemLayer->getNumberOfSheet();
63 //myRXstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu2();
64 //myRVstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu1();
65 //if(myPrintFlag)
66 cout<<"CGEM layer "<<i<<": "
67 <<" Rmid="<<myRmidDGapCgem[i]<<" cm "
68 <<" RV ="<<myRVCgem[i]<<" cm "
69 <<" RX ="<<myRXCgem[i]<<" cm "
70 <<" , stereo angle = "<<myAngStereoCgem[i]<<" radian "
71 <<myNSheets[i]<<" sheets"
72 //<<", is reversed "<<IsReverse
73 //<<", RX="<<myRXstrips[i]<<", RY="<<myRVstrips[i]
74 <<endl;
75 }
76
77 // --- get MDCGeomSvc ---
78 IMdcGeomSvc* ISvc2;
79 StatusCode mdc_sc=svcLocator->service("MdcGeomSvc", ISvc2);
80 if(mdc_sc.isSuccess()) myMdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc2);
81 else cout<<"DotsHelixFitter::initialize(): can not get MdcGeomSvc"<<endl;
82 int nLayer= myMdcGeomSvc->getLayerSize();
83 if(nLayer!=43) cout<<"DotsHelixFitter::initialize(): MDC nLayer = "<<nLayer<<" !=43 !!!"<<endl;
84 for(int i=0; i<nLayer; i++)
85 {
86 const MdcGeoLayer* aLayer = myMdcGeomSvc->Layer(i);
87 myRWires[i]=0.1*(aLayer->Radius());// cm
88 cout<<"MDC layer "<<i<<" R = "<<myRWires[i]<<endl;
89 double nomShift = myMdcGeomSvc->Layer(i)->nomShift();
90 if(nomShift>0) myWireFlag[i]=1;
91 else if(nomShift<0) myWireFlag[i]=-1;
92 else myWireFlag[i]=0;
93 int nWires=aLayer->NCell();
94 myNWire[i]=nWires;
95 for(int j=0; j<nWires; j++)
96 {
97 const MdcGeoWire* aWire = myMdcGeomSvc->Wire(i,j);
98 HepPoint3D aPointEast = 0.1 * (aWire->Forward());// in cm
99 double* aPointArray = new double[3]{aPointEast.x(),aPointEast.y(),aPointEast.z()};// in cm
100 myWestPosWires[i].push_back(aPointArray);
101 HepPoint3D aPointWest = 0.1*(aWire->Backward());// in cm
102 aPointArray = new double[3]{aPointWest.x(),aPointWest.y(),aPointWest.z()};// in cm
103 myEastPosWires[i].push_back(aPointArray);
104 HepPoint3D aPointMiddle = 0.5*(aPointEast+aPointWest);
105 aPointArray = new double[3]{aPointMiddle.x(),aPointMiddle.y(),aPointMiddle.z()};// in cm
106 myPosWires[i].push_back(aPointArray);
107 //myPhiWires[i].push_back(aWire->Forward().phi());
108 myPhiWires[i].push_back(aPointMiddle.phi());
109 myTensionWires[i].push_back(aWire->Tension());;
110 }
111 }
112
113 // --- MdcCalibFunSvc ---
114 IMdcCalibFunSvc* imdcCalibSvc;
115 StatusCode sc = svcLocator->service("MdcCalibFunSvc", imdcCalibSvc);
116 if ( sc.isSuccess() ){
117 myMdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
118 }
119 else {
120 cout<<"DotsHelixFitter::initialize(): can not get MdcCalibFunSvc"<<endl;
121 }
122
123 // --- get CgemCalibFunSvc ---
124 sc = svcLocator->service("CgemCalibFunSvc", myCgemCalibSvc);
125 if ( !(sc.isSuccess()) )
126 {
127 cout<<"DotsHelixFitter::initialize(): can not get CgemCalibFunSvc"<<endl;
128 }
129
130 // --- set alpha
131 //KalmanFit::Helix aHelix;
132 //myAlpha=aHelix.alpha();
133
134 // --- set alpha and update TrkFitFun::BFIELD
135 updateBField();
136
137 myIsInitialized=true;
138
139 // --- get MdcUtilitySvc
140 sc = svcLocator->service("MdcUtilitySvc", myMdcUtilitySvc);
141 if ( !(sc.isSuccess()) )
142 {
143 cout<<"DotsHelixFitter::initialize(): can not get MdcUtilitySvc"<<endl;
144 }
145
146 return;
147}
148
149void DotsHelixFitter::updateBField()
150{
151 KalmanFit::Helix aHelix;
152 myAlpha=aHelix.alpha();
153 TrkFitFun::BFIELD=fabs(aHelix.bFieldZ())/10.;
154 //cout<<__FILE__<<" "<<__LINE__
155 // <<", myAlpha="<<myAlpha
156 // <<", BFIELD="<<TrkFitFun::BFIELD<<" Tesla"
157 // <<endl;
158}
159
160void DotsHelixFitter::setDChits(vector<const MdcDigi*> aVecMdcDigi, double T0)
161{
162 myVecMdcDigi=aVecMdcDigi;
163 int nDigi = myVecMdcDigi.size();
164 //cout<<"set a vector<const MdcDigi*> with "<<myVecMdcDigi.size()<<" hits"<<endl;
165 myPhiAmbiVecMdcDigi.clear();
166 myChiVecMdcDigi.clear();
167 myAmbiguityMdcDigi.clear();
168 myMdcDigiIsActive.clear();
169 myPhiAmbiVecMdcDigi.resize(nDigi, -9999.);
170 myChiVecMdcDigi.resize(nDigi, 9999);
171 myAmbiguityMdcDigi.resize(nDigi, 0);
172 myMdcDigiIsActive.resize(nDigi, 1);
173 myEventT0 = T0;
174
175 for(int i=0; i<43; i++) myNumMdcDigiPerLayer[i]=0;
176
177 //for(int i=0; i<nDigi; i++)
178 //{
179 // myChiVecMdcDigi[i]=9999;
180 // myAmbiguityMdcDigi[i]=0;
181 // myMdcDigiIsActive[i]=1;
182 //}
183}
184
185/**
186 * @brief set myMdcDigiIsActive from vec
187 *
188 * @param vec
189 * @return true:success;
190 * @return false:fail due to wrong size
191 */
192bool DotsHelixFitter::setDChitsFitFlag(vector<int>& vec)
193{
194 bool flag=false;
195 if(myMdcDigiIsActive.size()==vec.size()) {
196 myMdcDigiIsActive=vec;
197 flag=true;
198 }
199 else cout<<"digi size "<<myMdcDigiIsActive.size()<<" != vec size "<<vec.size()<<endl;
200 return flag;
201}
202
203bool DotsHelixFitter::setDChitsPhiAmbi(const vector<double>& vec)
204{
205 bool flag=false;
206 if(myPhiAmbiVecMdcDigi.size()==vec.size()) {
207 myPhiAmbiVecMdcDigi=vec;
208 flag=true;
209 }
210 else cout<<"DotsHelixFitter::setDChitsPhiAmbi: digi size "<<myPhiAmbiVecMdcDigi.size()<<" != vec size "<<vec.size()<<endl;
211 return flag;
212}
213
214void DotsHelixFitter::setCgemClusters(vector<const RecCgemCluster*> aVecCgemCluster)
215{
216 myVecCgemCluster=aVecCgemCluster;
217 int nCluster1D = myVecCgemCluster.size();
218 myChiVecCgemCluster.resize(nCluster1D,9999);
219 myCgemClusterIsActive.resize(nCluster1D,1);
220 //cout<<"set a vector<RecCgemCluster*> with "<<myVecCgemCluster.size()<<" clusters"<<endl;
221 myVecCgemClusterX.clear();
222 myVecCgemClusterV.clear();
223 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
224 //for(int i=0; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i++)
225 //{
226 // //int layer = (*iter_cluster)->getlayerid();
227 // //int flag = (*iter_cluster)->getflag();
228 // //cout<<"set CGEM cluster "<<i<<", layer "<<layer<<", flag "<<flag<<endl;
229 // myChiVecCgemCluster[i]=9999;
230 // myCgemClusterIsActive[i]=1;
231 //}
232}
233
234void DotsHelixFitter::clearAllHits()
235{
236 myVecCgemCluster.clear();
237 myVecCgemClusterX.clear();
238 myVecCgemClusterV.clear();
239 myChiVecCgemCluster.clear();
240 myCgemClusterIsActive.clear();
241 myVecMdcDigi.clear();
242 myChiVecMdcDigi.clear();
243 myAmbiguityMdcDigi.clear();
244 myMdcDigiIsActive.clear();
245}
246
247void DotsHelixFitter::resetChi2FitFlag()
248{
249 int nCluster1D = myVecCgemCluster.size();
250 myChiVecCgemCluster.resize(nCluster1D,9999);
251 myCgemClusterIsActive.resize(nCluster1D,1);
252 int nDigi = myVecMdcDigi.size();
253 myPhiAmbiVecMdcDigi.resize(nDigi, -9999.);
254 myChiVecMdcDigi.resize(nDigi, 9999);
255 myAmbiguityMdcDigi.resize(nDigi, 0);
256 myMdcDigiIsActive.resize(nDigi, 1);
257}
258
259int DotsHelixFitter::calculateNewHelix()
260{
261 bool debug = false; //debug=true;
262
263 int state = 0;
264
265 double tension = 9999.;
266
267 int nPar=5;
268 if(myFitCircleOnly) nPar=3;
269 if(myFitCircleOriginOnly) nPar=2;
270
271 // --- check hits
272 //if(debug) cout<<"DotsHelixFitter::calculateNewHelix() starts checking hits ... "<<endl;
273 for(int i=0; i<43; i++) myNumMdcDigiPerLayer[i]=0;
274 myNActiveHits=0;
275 myNFittedHits=0;
276 int nXHits[3]={0,0,0}, nVHits[3]={0,0,0};
277 int nDigi=myVecMdcDigi.size();
278 //double* vec_zini = new double[nDigi];
279 vector<double> vec_zini(nDigi);
280 int maxLayer = 0;
281 int i_digi=0;
282 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
283 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
284 {
285 // --- initialize vec_zini ---
286 vec_zini[i_digi]=0;
287 // --- get id, layer, wire ---
288 Identifier id = (*iter_mdcDigi)->identify();
289 int layer = MdcID::layer(id);
290 int wire = MdcID::wire(id);
291 // --- counting active hits
292 if(myMdcDigiIsActive[i_digi]==1)
293 {
294 myNActiveHits++;
295
296 int hitFlag = myWireFlag[layer];
297 if(hitFlag==0)
298 {
299 nXHits[0]++;// total
300 nXHits[1]++;// MDC
301 }
302 else
303 {
304 nVHits[0]++;
305 nVHits[1]++;
306 }
307 }
308 if(maxLayer<layer) maxLayer=layer;
309 int nHits = myNumMdcDigiPerLayer[layer];
310 if(nHits<2) {
311 myIdxMdcDigiNeighbour[layer][nHits]=i_digi;
312 myWireIdMdcDigiNeighbour[layer][nHits]=wire;
313 }
314 myNumMdcDigiPerLayer[layer]++;
315 }
316 //if(debug) cout<<"DotsHelixFitter::calculateNewHelix() ends checking hits ... "<<endl;
317 HepPoint3D farPoint = myHelix->x(M_PI);
318 double maxRTrk = farPoint.perp();
319 for(int i=0; i<43; i++)
320 {
321 if(maxRTrk<myRWires[i]+2) break;// where soft track turning back
322 if(myNumMdcDigiPerLayer[i]==2)
323 {
324 int wire_1 = myWireIdMdcDigiNeighbour[i][0];
325 int wire_2 = myWireIdMdcDigiNeighbour[i][1];
326 int delta_n = abs(wire_1-wire_2);
327 int ambi_1 = 0;
328 int ambi_2 = 0;
329 if(delta_n==1)
330 {
331 ambi_1 = wire_1>wire_2? 1:-1;
332 ambi_2 = wire_1>wire_2? -1:1;
333 }
334 else if(delta_n==myNWire[i]-1)
335 {
336 ambi_1 = wire_1<=1? 1:-1;
337 ambi_2 = wire_2<=1? 1:-1;
338 }
339 //if(debug)
340 // cout<<"layer, wire1, wire2, ambi1, ambi2 = "
341 // <<setw(5)<<i
342 // <<setw(5)<<wire_1
343 // <<setw(5)<<wire_2
344 // <<setw(5)<<ambi_1
345 // <<setw(5)<<ambi_2
346 // <<endl;
347 // --- fix Ambiguity for neighboured mdc hits
348 //myAmbiguityMdcDigi[myIdxMdcDigiNeighbour[i][0]]=ambi_1;
349 //myAmbiguityMdcDigi[myIdxMdcDigiNeighbour[i][1]]=ambi_2;
350 }
351 }
352 i_digi=0;
353 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
354 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
355 {
356 int flag = (*iter_cluster)->getflag();
357 if(flag==0)
358 {
359 nXHits[0]++;// total
360 nXHits[2]++;// CGEM
361 }
362 else if(flag==1)
363 {
364 nVHits[0]++;
365 nVHits[2]++;
366 }
367 myNActiveHits++;
368 }
369 if(myFitCircleOnly)
370 {
371 if(nXHits[0]<myMinXhitsInCircleFit) return 1;
372 }
373 else
374 {
375 if((nXHits[0]+nVHits[0])<myMinHitsInHelixFit) return 2;
376 if(nVHits[0]<myMinVhitsInHelixFit) return 3;
377 }
378
379
380
381
382 // ---> start fit
383 double lastTotChi2 = 999999;
384 double secLastTotChi2 = 999999;// for oscillation check
385 double thirdLastTotChi2 = 999999;// for oscillation check
386 int nIterations = 0;
387 int n_chi2_increase = 0;
388 int n_oscillation=0;
389
390 while(1) // iterations
391 {
392 myNFittedHits=0;
393 myMapFlylenIdx.clear();
394 // --- matrix U P J ---
395 HepSymMatrix U(nPar, 0);
396 //HepMatrix P(5,1,0);
397 //HepMatrix J(5,1,0), JT(1,5,0);
398 //HepMatrix J_dmdx(1,3,0), J_dxda(3,5,0);
399 HepMatrix P(nPar,1,0);
400 HepMatrix J(nPar,1,0), JT(1,nPar,0);
401 HepMatrix J_dmdx(1,3,0), J_dxda(3,nPar,0);
402
403 // --- loop MDC hits ---
404 double totChi2=0;
405 i_digi=0;
406 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
407 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
408 {
409 // --- get id, layer, wire ---
410 Identifier id = (*iter_mdcDigi)->identify();
411 int layer = MdcID::layer(id);
412 int wire = MdcID::wire(id);
413 double charge = (*iter_mdcDigi)->getChargeChannel();
414
415 // --- get doca from track parameters ---
416 //tension=myTensionWires[layer][wire];
417 double wpos[7]={myEastPosWires[layer][wire][0], myEastPosWires[layer][wire][1], myEastPosWires[layer][wire][2]
418 , myWestPosWires[layer][wire][0], myWestPosWires[layer][wire][1], myWestPosWires[layer][wire][2], tension};
419 double doca_trk;
420 double whitPos[3];// approching position on wire (in cm)
421 if(debug)
422 {
423 //cout<<"a = "<<myHelix_aVec<<endl;
424 cout<<"* layer "<<layer<<", wire "<<wire<<", is active "<<myMdcDigiIsActive[i_digi] <<endl;
425 //cout<<"wire positions : ("<<wpos[0]<<", "<<wpos[1]<<", "<<wpos[2]<<") ("<<wpos[3]<<", "<<wpos[4]<<", "<<wpos[5]<<")"<<endl;
426 //cout<<"zini = "<<vec_zini[i_digi]<<endl;
427 }
428 TrkFitFun::getDoca(myHelix_a, wpos, doca_trk, whitPos, vec_zini[i_digi]);
429 //double doca_trk2 = myMdcUtilitySvc->doca(layer, wire, myHelix_aVec, myHelix_Ea, false, false);
430 //if(debug) {
431 // cout<<"doca = "<<doca_trk<<endl;
432 // //cout<<"doca2 = "<<doca_trk2<<endl;
433 //}
434 int leftRight = 2;
435 if(doca_trk!=0) {
436 leftRight=int(doca_trk/fabs(doca_trk));
437 }
438 /*
439 if(myAmbiguityMdcDigi[i_digi]!=0) {
440 leftRight=myAmbiguityMdcDigi[i_digi];
441 //if(debug) cout<<"fix leftRight = "<<leftRight<<endl;
442 }*/
443 int signDoca=1;
444 signDoca=leftRight/fabs(leftRight);
445 // --- conversion of left-right into the BESIII convention
446 // if(leftRight==-1) leftRight=0;
447 if(leftRight==1) leftRight=0;// fixed 2020-11-26
448 else leftRight=1;
449 //if(debug) cout<<"leftRight = "<<leftRight<<endl;
450
451 vec_zini[i_digi] = whitPos[2];// update vec_zini
452
453 // --- get measured doca --- tof in ns, driftTime in ns, T0Walk in ns
454 double rawTime = RawDataUtil::MdcTime((*iter_mdcDigi)->getTimeChannel());
455 double tprop = myMdcCalibFunSvc->getTprop(layer, vec_zini[i_digi]);
456 double T0Walk = myMdcCalibFunSvc->getT0(layer,wire) + myMdcCalibFunSvc->getTimeWalk(layer, charge);
457 //cout<<"line "<<__LINE__<<endl;
458 // --- time of flight (TOF) ---
459 KalmanFit::Helix aHelix = *myHelix;
460 HepPoint3D aNewPivot(whitPos[0],whitPos[1],whitPos[2]);
461 //cout<<"line "<<__LINE__<<endl;
462 aHelix.pivot(aNewPivot);
463 //cout<<"line "<<__LINE__<<endl;
464 double newPhi0 = aHelix.phi0();
465 //cout<<"line "<<__LINE__<<endl;
466 double dphi = newPhi0-myHelix->phi0();
467 //cout<<"line "<<__LINE__<<", dphi="<<dphi<<endl;
468 //if(dphi<-M_PI) dphi+=int(fabs(dphi/(2*M_PI)))*(2*M_PI);
469 //if(dphi> M_PI) dphi-=int(dphi/(2*M_PI))*(2*M_PI);
470 if(dphi<-M_PI||dphi> M_PI) dphi=atan2(sin(dphi),cos(dphi));// into [-pi, pi]
471 double flightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
472 //cout<<"line "<<__LINE__<<endl;
473 myMapFlylenIdx[flightLength]=i_digi;
474 HepLorentzVector p4_pi = myHelix->momentum(dphi, 0.13957);
475 double speed = p4_pi.beta()*TrkFitFun::CC;// cm/second
476 double TOF = flightLength/speed*1.e9;// in ns
477 //cout<<"line "<<__LINE__<<endl;
478 // --- drift time ---
479 double driftT = rawTime - myEventT0 - TOF - T0Walk - tprop;
480 //if(debug) cout<<"myEventT0, driftT = "<<myEventT0<<", "<<driftT<<endl;
481 // --- entrance Angle ---
482 double phiWire = atan2(whitPos[1],whitPos[0]);
483 double phiP = p4_pi.phi();
484 double entranceAngle = phiP-phiWire;
485 //while(entranceAngle<-M_PI) entranceAngle+=2*M_PI;
486 //while(entranceAngle> M_PI) entranceAngle-=2*M_PI;
487 if(entranceAngle<-M_PI||entranceAngle> M_PI) entranceAngle=atan2(sin(entranceAngle),cos(entranceAngle));// into [-pi, pi]
488 // --- measured drift distance ---
489 double doca_hit = 0.1 * myMdcCalibFunSvc->driftTimeToDist(driftT,layer,wire,leftRight, entranceAngle);// in cm
490 // --- get measurement error ---
491 double docaErr_hit = 0.1 * myMdcCalibFunSvc->getSigma(layer, leftRight, doca_hit, entranceAngle, myHelix_a[4], vec_zini[i_digi], charge);// in cm
492 //if(debug) cout<<"error_doca_hit = "<<docaErr_hit<<endl;
493
494 // --- get derivatives, calculate J, update P, U ---
495 doca_hit*=signDoca;
496 //if(debug) cout<<"doca_hit = "<<doca_hit<<endl;
497 double delD = doca_hit-doca_trk;
498 double chi = delD/docaErr_hit;
499 if(debug) cout<<"delta_m, err_m, chi = "<<delD<<", "<<docaErr_hit<<", "<<chi<<endl;
500 myChiVecMdcDigi[i_digi]=chi;
501
502 if(myMdcDigiIsActive[i_digi]!=1) continue;
503 if(myUseAxialHitsOnly && myWireFlag[layer]!=0) continue;
504
505 myNFittedHits++;
506 totChi2+=chi*chi;
507 for(int ipar=0; ipar<nPar; ipar++)
508 {
509 double deri;
510 //if(debug) cout<<"ipar = "<<ipar<<endl;
511 TrkFitFun::getDeriLoc(ipar,myHelix_a, deri, wpos, vec_zini[i_digi]);// in cm
512 J(ipar+1,1) = deri;
513 //if(debug) cout<<"d(doca)/d(a"<<"_"<<ipar<<") = "<<J(ipar+1,1)<<endl;
514 //P(1,ipar+1) += J(1,ipar+1)*(doca_hit-doca_trk)/(docaErr_hit*docaErr_hit);
515 //P(ipar+1,1) += J(ipar+1,1)*chi/docaErr_hit;
516 double scale=1; //if(fabs(chi)>5) scale=fabs(chi*chi);
517 P(ipar+1,1) += J(ipar+1,1)*(doca_hit-doca_trk)/(docaErr_hit*docaErr_hit*scale);// large error for large chi
518 for(int ipar2=0; ipar2<=ipar; ipar2++)
519 {
520 //U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(docaErr_hit*docaErr_hit);
521 U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(docaErr_hit*docaErr_hit*scale);// large error for large chi2
522 }
523 }
524 //if(debug) cout<<"U="<<U<<endl;
525 }// loop MDC hits
526 if(debug) cout<<"end of MDC hits loop"<<endl;
527 // <--- end of MDC hits loop ---
528
529 // --- loop CGEM clusters ---
530 i_digi=0;
531 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
532 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
533 {
534 // --- get layer, cluster flag ---
535 int layer = (*iter_cluster)->getlayerid();
536 int flag = (*iter_cluster)->getflag();
537 if(myUseAxialHitsOnly && flag!=0) continue;
538 if(myUseMdcHitsOnly) continue;
539 int sheet = (*iter_cluster)->getsheetid();
540 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(layer,sheet);
541 if(debug)
542 {
543 cout<<endl<<"layer, sheet, flag = "<<setw(5)<<layer<<setw(5)<<sheet<<setw(5)<<flag<<endl;
544 }
545
546 // --- get position & entrance Angle from track parameters ---
547 double dphi = IntersectCylinder(myRmidDGapCgem[layer]);
548 if(fabs(dphi)<1e-10) {
549 myChiVecCgemCluster[i_digi]=-9999;
550 //cout<<"trk has no intersection with CGEM"<<endl;
551 continue;// track has no intersection with CGEM
552 }
553 double flightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
554 myMapFlylenIdx[flightLength]=-(i_digi+1);
555 HepPoint3D pos = myHelix->x(dphi);
556 if(debug) cout<<"dphi="<<dphi<<", pos="<<pos<<endl;
557 double phi_trk = pos.phi();// in radian
558 double z_trk = pos.z();// in cm
559 J_dxda = dxda_cgem(*myHelix, dphi);
560 Hep3Vector p3_trk = myHelix->momentum(dphi);
561 double phi_p3trk = p3_trk.phi();
562 double incidentPhi = phi_p3trk-phi_trk;
563 //while(incidentPhi>CLHEP::pi) incidentPhi-=CLHEP::twopi;
564 //while(incidentPhi<-CLHEP::pi) incidentPhi+=CLHEP::twopi;
565 if(incidentPhi<-M_PI||incidentPhi> M_PI) incidentPhi=atan2(sin(incidentPhi),cos(incidentPhi));// into [-pi, pi]
566
567 // --- get X-cluster charge, measured X, derivatives, calculate J, update P, U ---
568 double phi_CC(0.0);// in radian
569 double X_CC(0.), V_CC(0.);// in cm
570 double delta_m, err_m;
571 double Q(0.);// in fC
572 double T(100);// not valid at present
573 int mode(2);
574 Q=(*iter_cluster)->getenergydeposit();
575 if(flag==0) {
576 phi_CC = (*iter_cluster)->getrecphi();
577 //phi_CC = (*iter_cluster)->getrecphi_CC();
578 double del_phi = phi_CC-phi_trk;
579 //while(del_phi<-M_PI) del_phi+=2*M_PI;
580 //while(del_phi> M_PI) del_phi-=2*M_PI;
581 if(del_phi<-M_PI||del_phi> M_PI) del_phi=atan2(sin(del_phi),cos(del_phi));// into [-pi, pi]
582 X_CC=phi_CC*myRmidDGapCgem[layer];
583 if(debug) cout<<"phi_trk, phi_rec = "<<phi_trk<<", "<<phi_CC<<endl;
584 //double dX = del_phi*myRmidDGapCgem[layer];
585 //double dX = del_phi;
586 delta_m=del_phi;
587 err_m=myCgemCalibSvc->getSigma(layer,flag,mode,incidentPhi,Q,T)/10.;// in cm
588 err_m/=myRmidDGapCgem[layer];// in radian
589 J_dmdx(1,1)=-1*pos.y()/myR2midDGapCgem[layer];
590 J_dmdx(1,2)=pos.x()/myR2midDGapCgem[layer];
591 J_dmdx(1,3)=0.;
592 }
593 else if(flag==1) {
594 double V_trk = readoutPlane->getVFromPhiZ(phi_trk,z_trk*10,false)/10.0;// in cm
595 V_CC=(*iter_cluster)->getrecv()/10.;// in cm
596 //V_CC=(*iter_cluster)->getrecv_CC()/10.;// in cm
597 if(debug) cout<<"V_trk, V_rec = "<<V_trk<<", "<<V_CC<<endl;
598 delta_m=V_CC-V_trk;// in cm
599 if(fabs(delta_m)>5) {
600 /*int nextSheet = myNSheets[layer]-1-sheet;
601 CgemGeoReadoutPlane* nextReadoutPlane = myCgemGeomSvc->getReadoutPlane(layer,nextSheet);
602 double phiminNext = nextReadoutPlane->getPhimin();
603 double V_trk_next = nextReadoutPlane->getVInNextSheetFromV(V_trk*10, phiminNext)/10.0;// cm*10 -> mm
604 double delta_m_next = V_CC-V_trk_next;
605 if(fabs(delta_m_next)<delta_m) delta_m=delta_m_next;
606 if(debug) cout<<"V_trk_next = "<<V_trk_next<<endl;*/
607 double V_trk_nearPhiMin = readoutPlane->getVFromPhiZ_nearPhiMin(phi_trk,z_trk*10,false)/10.0;// in cm
608 double delta_m_2=V_CC-V_trk_nearPhiMin;// in cm
609 if(fabs(delta_m_2)<fabs(delta_m)) delta_m=delta_m_2;
610 if(debug) cout<<"V_trk_nearPhiMin= "<<V_trk_nearPhiMin<<endl;
611 }
612 err_m=myCgemCalibSvc->getSigma(layer,flag,mode,incidentPhi,Q,T)/10.;// in cm
613 J_dmdx(1,1)=-myRVCgem[layer]*cos(myAngStereoCgem[layer])*pos.y()/myR2midDGapCgem[layer];
614 J_dmdx(1,2)= myRVCgem[layer]*cos(myAngStereoCgem[layer])*pos.x()/myR2midDGapCgem[layer];
615 J_dmdx(1,3)= -sin(myAngStereoCgem[layer]);
616 }
617 else {
618 cout<<"flag ="<<flag<<", DotsHelixFitter::calculateNewHelix() is not ready for it!"<<endl;
619 continue;// next cluster
620 }
621 JT=J_dmdx*J_dxda;
622 J=JT.T();
623 //if(debug) {
624 // cout<<"J_dmdx="<<J_dmdx<<endl;
625 // cout<<"J_dxda="<<J_dxda<<endl;
626 // cout<<"J="<<J<<endl;
627 //}
628 double chi = delta_m/err_m;
629 myChiVecCgemCluster[i_digi]=chi;
630 if(debug) cout<<"delta_m, err_m, chi = "<<delta_m<<", "<<err_m<<", "<<chi<<endl;
631 totChi2+=chi*chi;
632 myNFittedHits++;
633 for(int ipar=0; ipar<nPar; ipar++)
634 {
635 //if(debug)
636 // cout<<"d(doca)/d(a"<<"_"<<ipar<<") = "<<J(ipar+1,1)<<endl;
637 P(ipar+1,1) += J(ipar+1,1)*chi/err_m;
638 for(int ipar2=0; ipar2<=ipar; ipar2++)
639 {
640 U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(err_m*err_m);
641 }
642 }
643 if(debug) cout<<"U="<<U<<endl;
644 }// loop myVecCgemCluster
645 if(debug)
646 cout<<"end of CGEM cluster loop"<<endl;
647 // <--- end of CGEM clusters loop ---
648
649 // --- delta a ---
650 //if(debug) cout<<"U="<<U<<endl;
651 int ifail=99;
652 U.invert(ifail);// also the error matrix
653 //if(debug)
654 //{
655 // cout<<"ifail = "<<ifail<<endl;
656 // cout<<"U^-1="<<U<<endl;
657 // cout<<"P="<<P<<endl;
658 //}
659 HepVector da = U*P;// in cm
660 //if(debug) cout<<"da = "<<da<<endl;
661
662 // --- update helix ---
663 for(int ipar=0; ipar<nPar; ipar++)
664 {
665 myHelix_a[ipar]+=da[ipar];
666 if(ipar==1&&(myHelix_a[ipar]>2*M_PI||myHelix_a[ipar]<0)) // put phi0 in (0, 2pi), added by wangll, 2022-05-12
667 {
668 double aphi=myHelix_a[ipar];
669 aphi=atan2(sin(aphi),cos(aphi));
670 if(aphi<0) aphi+=2*M_PI;
671 myHelix_a[ipar]=aphi;
672 }
673 myHelix_aVec[ipar]=myHelix_a[ipar];
674 }
675 myHelix->a(myHelix_aVec);
676 if(debug)
677 {
678 cout<<"aNew = "<<myHelix_aVec
679 <<" chi2 = "<<totChi2
680 <<endl;
681 }
682 myChi2=totChi2;
683 //if(nIterations==0) cout<<"before fit chi2 = "<<totChi2<<endl;
684
685 // --- converge or too many iterations
686 nIterations++;
687 if(debug) cout<<" iteration "<<nIterations<<", chi2="<<totChi2<<endl;
688 if(fabs(lastTotChi2-totChi2)<myDchi2Converge) {
689 if(debug)
690 cout<<"DotsHelixFitter::calculateNewHelix(): converge after "<<nIterations<<" iterations"
691 <<" with lastTotChi2="<<lastTotChi2<<", totChi2="<<totChi2
692 <<endl;
693 if(nPar==5) {
694 myHelix_Ea = U;
695 myHelix->Ea(U);
696 }
697 return 0; // good
698 //cout<<"myHelix_Ea updated "<<endl;
699 //break;
700 }
701 else if(totChi2-lastTotChi2>myDchi2Diverge) {
702 n_chi2_increase++;
703 if(n_chi2_increase>myMaxNChi2Increase||totChi2>myChi2Diverge)
704 return 5;// divergence
705 }
706 if( nIterations>3 && fabs(secLastTotChi2-totChi2)<(0.1*myDchi2Converge) && fabs(thirdLastTotChi2-lastTotChi2)<(0.1*myDchi2Converge) ) {
707 n_oscillation++;
708 //cout<<" n_oscillation="<<n_oscillation<<", thirdLastTotChi2, secLastTotChi2, lastTotChi2 = "<<thirdLastTotChi2<<", "<<secLastTotChi2<<", "<<lastTotChi2<<endl;
709 if(n_oscillation>0&&totChi2<lastTotChi2) {
710 if(nPar==5) {
711 myHelix_Ea = U;
712 myHelix->Ea(U);
713 }
714 if(debug) cout<<"DotsHelixFitter::calculateNewHelix(): oscillation happened, thirdLastTotChi2, secLastTotChi2, lastTotChi2 = "<<thirdLastTotChi2<<", "<<secLastTotChi2<<", "<<lastTotChi2<<endl;
715 return 0;// break from oscillation
716 }
717 }
718 if(nIterations>myMaxIteration) {
719 if(debug) cout<<"DotsHelixFitter::calculateNewHelix(): "<<nIterations
720 <<" iterations, break with lastTotChi2, totChi2 = "<<lastTotChi2<<", "<<totChi2<<endl;
721 return 4;// too many iterations
722 //break;
723 }
724 thirdLastTotChi2=secLastTotChi2;
725 secLastTotChi2=lastTotChi2;//for oscillation check
726 lastTotChi2=totChi2;
727 }// <--- iterations
728
729 // --- delete vec_zini --- // delete []vec_zini;
730
731 //if(nIterations>myMaxIteration) return 4;
732 //else return 0;
733}
734
735
736int DotsHelixFitter::deactiveHits(double chi_cut, int nMax, int layerMin, int layerMax)
737{
738 map<double, int> map_abschi_idx;
739 //double maxChi2=0;
740 //vector<const MdcDigi*>::iterator iter_mdcDigi_maxChi2=NULL;
741 //int i_digi_maxChi2=-1;
742
743 int i_digi=0;
744 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
745 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
746 {
747 if(myMdcDigiIsActive[i_digi]==1)
748 {
749 Identifier id = (*iter_mdcDigi)->identify();
750 int layer = MdcID::layer(id);
751 if(layer<layerMin||layer>layerMax) continue;
752 double abs_chi=fabs(myChiVecMdcDigi[i_digi]);
753
754 if( abs_chi>chi_cut ) // && maxChi2<(abs_chi*abs_chi)
755 {
756 //maxChi2=abs_chi*abs_chi;
757 //iter_mdcDigi_maxChi2=iter_mdcDigi;
758 //i_digi_maxChi2=i_digi;
759 map_abschi_idx[abs_chi]=i_digi;
760 }
761 }
762 }
763
764 int nDeactive=0;
765 if(nMax>0&&map_abschi_idx.size()>0)
766 {
767 map<double, int>::reverse_iterator it_map;
768 for(it_map=map_abschi_idx.rbegin(); it_map!=map_abschi_idx.rend(); it_map++)
769 {
770 myMdcDigiIsActive[it_map->second]=-1;
771 //cout<<"a hit with |chi|="<<it_map->first<<" deactived"<<endl;
772 nDeactive++;
773 if(nDeactive>=nMax) break;
774 }
775 }
776 //if(i_digi_maxChi2!=-1)
777 //{
778 // myMdcDigiIsActive[i_digi_maxChi2] = 0;
779 // //cout<<"the worst hit "<<i_digi_maxChi2<<" (chi2="<<maxChi2<<") is deactive"<<endl;
780 // return 1;
781 //}
782 //else {
783 // //cout<<"no bad hits found!"<<endl;
784 // return 0;
785 //}
786
787 if(nDeactive<nMax) {
788 i_digi=0;
789 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
790 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
791 {
792 if(myCgemClusterIsActive[i_digi]!=1) continue;
793 int layer = (*iter_cluster)->getlayerid();
794 if(layer<layerMin||layer>layerMax) continue;
795 double abs_chi=fabs(myChiVecCgemCluster[i_digi]);
796 if( abs_chi>chi_cut ) map_abschi_idx[abs_chi]=-(i_digi+1);
797 }
798 if(nMax>0&&map_abschi_idx.size()>0)
799 {
800 map<double, int>::reverse_iterator it_map;
801 for(it_map=map_abschi_idx.rbegin(); it_map!=map_abschi_idx.rend(); it_map++)
802 {
803 int digiIdx=it_map->second;
804 if(digiIdx>=0) continue;
805 digiIdx=-1*digiIdx-1;
806 myCgemClusterIsActive[digiIdx]=-1;
807 //cout<<"a CGEM cluster with |chi|="<<it_map->first<<" deactived"<<endl;
808 nDeactive++;
809 if(nDeactive>=nMax) break;
810 }
811 }
812 }
813
814
815 return nDeactive;
816}
817
818/**
819 * @brief activate myMdcDigiIsActive's hit based on myChiVecMdcDigi
820 *
821 * @param chi_cut require |myChiVecMdcDigi[i]|<chi_cut
822 * @param sel 0:inactive and uncertain; 1:only uncertain
823 * @return int
824 */
825int DotsHelixFitter::activeHits(double chi_cut, int sel)
826{
827 int nActived=0;
828 int i_digi=0;
829 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
830 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
831 {
832 if((sel==0&&myMdcDigiIsActive[i_digi]<=0)||(sel==1&&myMdcDigiIsActive[i_digi]==0))
833 {
834 double abs_chi=fabs(myChiVecMdcDigi[i_digi]);
835
836 if( abs_chi<chi_cut )
837 {
838 myMdcDigiIsActive[i_digi] = 1;
839 nActived++;
840 }
841 }
842 }
843 return nActived;
844}
845
846
847int DotsHelixFitter::deactiveHits(int layer_max, int nHit_max)
848{
849 int nActive=0;
850 int nDeactive=0;
851
852 vector<const MdcDigi*>::iterator iter_mdcDigi_begin = myVecMdcDigi.begin();
853 map<double, int>::iterator it_map;
854 //cout<<"digi size, flyLen size = "<<myVecMdcDigi.size()<<", "<<myMapFlylenIdx.size()<<endl;
855 for(it_map=myMapFlylenIdx.begin(); it_map!=myMapFlylenIdx.end(); it_map++)
856 {
857 //cout<<"nActive = "<<nActive<<endl;
858 int i_hit = it_map->second;
859 if(i_hit>=0) // for MDC digi
860 {
861 vector<const MdcDigi*>::iterator iter_mdcDigi = iter_mdcDigi_begin+(i_hit);
862 if(myMdcDigiIsActive[i_hit]==1)
863 {
864 if(nActive>=nHit_max) myMdcDigiIsActive[i_hit]=-1;
865 Identifier id = (*iter_mdcDigi)->identify();
866 int layer = MdcID::layer(id);
867 if(layer>=layer_max) myMdcDigiIsActive[i_hit]=-1;
868 if(myMdcDigiIsActive[i_hit]==1) nActive++;
869 else {
870 nDeactive++;
871 //cout<<"a hit on layer "<<layer<<" is deactived"<<endl;
872 }
873 }
874 }
875 }
876 //cout<<nDeactive<<" hits deactived "<<endl;
877
878 return nDeactive;
879}
880
881void DotsHelixFitter::setInitialHelix(KalmanFit::Helix aHelix)
882{
883 HepPoint3D pivot(0,0,0);
884 myHelix_aVec = aHelix.a();
885 if(myHelix!=NULL) delete myHelix;
886 myHelix = new KalmanFit::Helix(pivot, myHelix_aVec);
887 //cout<<"alpha="<<myHelix->alpha()<<endl;
888 for(int i=0; i<5; i++) myHelix_a[i]=myHelix_aVec[i];
889
890}
891
892void DotsHelixFitter::loadOneDcDigi(const MdcDigi* aDcDigi)
893{
894 // --- get id, layer, wire ---
895 Identifier id = aDcDigi->identify();
896 myLayer = MdcID::layer(id);
897 myWire = MdcID::wire(id);
898 myCharge= aDcDigi->getChargeChannel();
899
900 // --- get myWirePos ---
901 double tension = 9999.;
902 //tension=myTensionWires[myLayer][myWire];
903 myWirePos[0]=myEastPosWires[myLayer][myWire][0];
904 myWirePos[1]=myEastPosWires[myLayer][myWire][1];
905 myWirePos[2]=myEastPosWires[myLayer][myWire][2];
906 myWirePos[3]=myWestPosWires[myLayer][myWire][0];
907 myWirePos[4]=myWestPosWires[myLayer][myWire][1];
908 myWirePos[5]=myWestPosWires[myLayer][myWire][2];
909 myWirePos[6]=tension;
910
911}
912
913void DotsHelixFitter::calculateRoughDD(const MdcDigi* aMdcDigi)
914{
915 loadOneDcDigi(aMdcDigi);
916
917 double rawTime = RawDataUtil::MdcTime(aMdcDigi->getTimeChannel());
918 double tprop = myMdcCalibFunSvc->getTprop(myLayer, 0.);
919 double T0Walk = myMdcCalibFunSvc->getT0(myLayer,myWire) + myMdcCalibFunSvc->getTimeWalk(myLayer, myCharge);
920 double TOF = 10*myRWires[myLayer]/Constants::c*1.0e9;// in ns
921 myDriftTime = rawTime - myEventT0 - TOF - T0Walk - tprop;
922
923 int lrCalib=2;
924 double entranceAngle = 0;
925 myDocaFromDigi = 0.1 * myMdcCalibFunSvc->driftTimeToDist(myDriftTime,myLayer,myWire,lrCalib,entranceAngle);// in cm
926}
927
928void DotsHelixFitter::calculateDocaFromTrk(const MdcDigi* aDcDigi)
929{
930 bool debug = false; //debug=true;
931
932 loadOneDcDigi(aDcDigi);
933
934 // --- get doca from track parameters ---
935 TrkFitFun::getDoca(myHelix_a, myWirePos, myDocaFromTrk, myPosOnWire, myPosOnWire[2]);
936 if(debug)
937 {
938 cout<<"DotsHelixFitter::UpdateDcDigiInfo(): "<<endl;
939 cout<<"doca = "<<myDocaFromTrk<<endl;
940 cout<<"point on wire: "<<myPosOnWire[0]<<", "<<myPosOnWire[1]<<", "<<myPosOnWire[2]<<endl;
941 }
942
943 // --- flight length ---
944 KalmanFit::Helix aHelix = *myHelix;
945 HepPoint3D aNewPivot(myPosOnWire[0],myPosOnWire[1],myPosOnWire[2]);
946 aHelix.pivot(aNewPivot);
947 double newPhi0 = aHelix.phi0();
948 double dphi = newPhi0-myHelix->phi0();
949 //while(dphi<-M_PI) dphi+=2*M_PI;
950 //while(dphi> M_PI) dphi-=2*M_PI;
951 if(dphi<-M_PI||dphi> M_PI) dphi=atan2(sin(dphi),cos(dphi));// into [-pi, pi]
952 myFlightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
953}
954
955/**
956 * @brief based on input aDcDigi, modifies myFlightLength, myLeftRight, myPosOnWire, myDocaFromDigi, myDriftDistErr, myDcChi, ...
957 *
958 * @param aDcDigi
959 */
960void DotsHelixFitter::updateDcDigiInfo(const MdcDigi* aDcDigi)
961{
962 bool debug = false;
963
964 loadOneDcDigi(aDcDigi);
965
966 // --- get doca from track parameters ---
967 myPosOnWire[2]=0;
968 TrkFitFun::getDoca(myHelix_a, myWirePos, myDocaFromTrk, myPosOnWire, myPosOnWire[2]);
969 double phiWire = atan2(myPosOnWire[1],myPosOnWire[0]);
970 if(debug)
971 {
972 cout<<"DotsHelixFitter::UpdateDcDigiInfo(): "<<endl;
973 cout<<"doca = "<<myDocaFromTrk<<endl;
974 cout<<"point on wire: "<<myPosOnWire[0]<<", "<<myPosOnWire[1]<<", "<<myPosOnWire[2]<<endl;
975 }
976
977 // --- flight length ---
978 KalmanFit::Helix aHelix = *myHelix;
979 HepPoint3D aNewPivot(myPosOnWire[0],myPosOnWire[1],myPosOnWire[2]);
980 aHelix.pivot(aNewPivot);
981 double newPhi0 = aHelix.phi0();
982 double dphi = newPhi0-myHelix->phi0();
983 //while(dphi<-M_PI) dphi+=2*M_PI;
984 //while(dphi> M_PI) dphi-=2*M_PI;
985 if(dphi<-M_PI||dphi> M_PI) dphi=atan2(sin(dphi),cos(dphi));// into [-pi, pi]
986 myFlightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
987
988 // --- approaching point on Helix
989 HepPoint3D posOnTrk = aHelix.x();
990 myPosOnTrk[0]=posOnTrk.x();
991 myPosOnTrk[1]=posOnTrk.y();
992 myPosOnTrk[2]=posOnTrk.z();
993 double phiPosOnTrk=atan2(myPosOnTrk[1],myPosOnTrk[0]);
994
995 // --- left or right
996 myLeftRight = 2;
997 int signDoca=1;
998 if(myDocaFromTrk!=0)
999 {
1000 myLeftRight=int(myDocaFromTrk/fabs(myDocaFromTrk));
1001 signDoca=myLeftRight;
1002
1003 // ---> conversion of left-right into the BESIII convention
1004 //if(myLeftRight==-1) myLeftRight=0;
1005 if(myLeftRight==1) myLeftRight=0;// fixed 2020-11-26
1006 else myLeftRight=1;
1007 double dphiLR=phiWire-phiPosOnTrk;
1008 //while(dphiLR<-M_PI) dphiLR+=2*M_PI;
1009 //while(dphiLR> M_PI) dphiLR-=2*M_PI;
1010 if(dphiLR<-M_PI||dphiLR> M_PI) dphiLR=atan2(sin(dphiLR),cos(dphiLR));// into [-pi, pi]
1011 //if(dphiLR>0) myLeftRight=0;
1012 //else if(dphiLR<0) myLeftRight=1;
1013 //cout<<"myLeftRight, dphiLR = "<<myLeftRight<<", "<<dphiLR<<endl;
1014 }
1015 if(debug) cout<<"myLeftRight = "<<myLeftRight<<endl;
1016
1017 // --- time of flight (TOF) ---
1018 HepLorentzVector p4_pi = myHelix->momentum(dphi, 0.13957);
1019 double speed = p4_pi.beta()*TrkFitFun::CC;// cm/second
1020 double TOF = myFlightLength/speed*1.e9;// in ns
1021
1022 // --- get measured doca --- tof in ns, driftTime in ns, T0Walk in ns
1023 double rawTime = RawDataUtil::MdcTime(aDcDigi->getTimeChannel());
1024 double tprop = myMdcCalibFunSvc->getTprop(myLayer, myPosOnWire[2]);
1025 double T0Walk = myMdcCalibFunSvc->getT0(myLayer,myWire) + myMdcCalibFunSvc->getTimeWalk(myLayer, myCharge);
1026
1027 // --- drift time ---
1028 myDriftTime = rawTime - myEventT0 - TOF - T0Walk - tprop;
1029 if(debug) cout<<"driftT = "<<myDriftTime<<endl;
1030 // --- entrance Angle ---
1031 double phiP = p4_pi.phi();
1032 double entranceAngle = phiP-phiWire;
1033 //while(entranceAngle<-M_PI) entranceAngle+=2*M_PI;
1034 //while(entranceAngle> M_PI) entranceAngle-=2*M_PI;
1035 if(entranceAngle<-M_PI||entranceAngle> M_PI) entranceAngle=atan2(sin(entranceAngle),cos(entranceAngle));// into [-pi, pi]
1036 myEntranceAngle = entranceAngle;
1037 // --- measured drift distance ---
1038 myDocaFromDigi = 0.1 * myMdcCalibFunSvc->driftTimeToDist(myDriftTime,myLayer,myWire,myLeftRight,entranceAngle);// in cm
1039 myDriftDist[myLeftRight]=myDocaFromDigi;
1040 myDriftDist[1-myLeftRight]=0.1 * myMdcCalibFunSvc->driftTimeToDist(myDriftTime,myLayer,myWire,1-myLeftRight,entranceAngle);// in cm
1041 // --- get measurement error ---
1042 double docaErr_hit = 0.1 * myMdcCalibFunSvc->getSigma(myLayer, myLeftRight, myDocaFromDigi, myEntranceAngle, myHelix_a[4], myPosOnWire[2], myCharge);// in cm
1043 //if(debug) cout<<"error_doca_hit = "<<docaErr_hit<<endl;
1044 myDriftDistErr[myLeftRight]=docaErr_hit;
1045 myDriftDistErr[1-myLeftRight] = 0.1 * myMdcCalibFunSvc->getSigma(myLayer, 1-myLeftRight, myDocaFromDigi, myEntranceAngle, myHelix_a[4], myPosOnWire[2], myCharge);// in cm
1046
1047 // --- get chi ---
1048 myDocaFromDigi*=signDoca;
1049 if(debug) cout<<"doca_hit = "<<myDocaFromDigi<<endl;
1050 double delD = myDocaFromDigi-myDocaFromTrk;
1051 myDcChi = delD/docaErr_hit;
1052}
1053
1054//vector<pair<double, double> > DotsHelixFitter::getSZ(const MdcDigi* aDcDigi)
1055//{
1056// // --- center and radius of the circle
1057// double xc = myHelix->center().x();
1058// double yc = myHelix->center().y();
1059// double rTrack = myHelix->radius();// signed FIXME
1060//
1061// // --- drift distance
1062// double drift = 0;
1063//
1064// // --- wire line segment on xy plane
1065// //loadOneDcDigi(aDcDigi);
1066// updateDcDigiInfo(aDcDigi);
1067// drift=0.5*(myDriftDist[0]+myDriftDist[1]);
1068// double west2east = sqrt((myWirePos[0]-myWirePos[3])*(myWirePos[0]-myWirePos[3])+(myWirePos[1]-myWirePos[4])*(myWirePos[1]-myWirePos[4]));
1069// double slope = (myWirePos[1]-myWirePos[4])/(myWirePos[0]-myWirePos[3]);
1070// double intercept = (myWirePos[4]-slope*myWirePos[3]+myWirePos[1]-slope*myWirePos[0])/2;
1071//
1072// // --- calculate sz
1073// vector<pair<double, double> > vec_sz;
1074// double s(0),z(0);
1075// double a = 1+slope*slope;
1076// double b = -2*(xc+slope*yc-slope*intercept);
1077// cout<<"layer "<<myLayer<<": ";
1078// for(int LR=-1; LR<2; LR+=2)
1079// {
1080// double c=xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+LR*drift)*(rTrack+LR*drift);
1081// double delta=b*b-4*a*c;
1082// if(delta>=0)
1083// {
1084// for(int sign=-1; sign<2; sign+=2)
1085// {
1086// double x = (-b+sign*sqrt(delta))/(2.0*a);
1087// double y = slope*x+intercept;
1088// if((x-myWirePos[3])*(x-myWirePos[0])<0)
1089// {
1090// HepPoint3D point(x,y,0);
1091// s = myHelix->flightArc(point);
1092// double l = sqrt((x-myWirePos[3])*(x-myWirePos[3])+(y-myWirePos[4])*(y-myWirePos[4]));
1093// z = myWirePos[5] + l/west2east*fabs((myWirePos[2]-myWirePos[5]));
1094// cout<<"("<<s<<","<<z<<") ";
1095// pair<double,double> sz(s,z);
1096// vec_sz.push_back(sz);
1097// }// if the position is on wire
1098// }// loop sign
1099// }// if delta>=0
1100// }// loop LR
1101//
1102// return vec_sz;
1103//}
1104
1105/**
1106 * @brief return square of x, aka pow(x,2) or x*x
1107 */
1108inline static double squared(double x){
1109 return x*x;
1110}
1111
1112/**
1113 * @brief return S, the track length in XY; and Z of a hit on a stereo wire
1114 *
1115 * Current basics:
1116 * line function
1117 * X_vec2 = lambda * (Xb_vec2-Xa_vec2) + Xa_vec2;
1118 * and circle function
1119 * (X-Xc)**2 + (Y-Yc)**2 = (R +-r_d)**2
1120 * solve system to yield lambda here,
1121 * so like X_vec2, we get z from X_vec3
1122 * place X_vec3 into myHelix->flightArc() for s
1123 *
1124 * @param aDcDigi
1125 * @return vector<pair<double, double> >
1126 */
1127vector<pair<double, double> > DotsHelixFitter::getSZ(const MdcDigi *aDcDigi) {
1128 bool debug=false; //debug=true;
1129
1130 vector<pair<double, double> > vec_sz;
1131 // --- center and radius of the circle
1132 double xc = myHelix->center().x();
1133 double yc = myHelix->center().y();
1134 double rTrack = myHelix->radius(); // signed FIXME
1135
1136 // --- drift distance
1137 updateDcDigiInfo(aDcDigi);
1138 double r_drift = 0.5 * (myDriftDist[0] + myDriftDist[1]);
1139
1140 // --- wire line segment
1141 HepPoint3D Xa(myWirePos[0], myWirePos[1], myWirePos[2]);
1142 HepPoint3D Xb(myWirePos[3], myWirePos[4], myWirePos[5]);
1143 HepPoint3D Xdelta = Xb - Xa;
1144
1145 // function
1146 // A*lambda **2 + B*lambda + C ==0
1147 double A = squared(Xdelta.x() ) + squared(Xdelta.y());
1148 double B = 2 * Xdelta.x() * (Xa.x() - xc) + 2 * Xdelta.y() * (Xa.y() - yc);
1149 double C_part = squared(Xa.x() - xc) + squared(Xa.y() - yc);
1150 double C1 = C_part - squared(rTrack + r_drift);
1151 double C2 = C_part - squared(rTrack - r_drift);
1152
1153 // for each of the two C
1154
1155 if(debug) cout<<"layer "<<myLayer<<": ";
1156 for (int _ = 0; _ < 2; _++) {
1157 double C = _ ? C2 : C1;
1158 double Delta = squared(B) - 4 * A * C;
1159
1160 if (Delta < 0) continue;// require Delta >=0
1161
1162 for (int sign = -1; sign < 2; sign += 2) {// if Delta >0, we have 2 roots
1163 double lambda = (-B + sign * sqrt(Delta)) / 2 / A;
1164 if (lambda > 1.2 or lambda < -0.2) continue; // require local; +/-20% for safety
1165
1166 HepPoint3D Xhit = Xa + lambda * Xdelta;
1167 double s = myHelix->flightArc(Xhit);
1168 if(debug) cout<<"("<<s<<","<<Xhit.z()<<") ";
1169 pair<double, double> sz(s, Xhit.z());
1170 vec_sz.push_back(sz);
1171 if (Delta == 0) break; // if Delta==0, run only once.
1172 }
1173
1174 }
1175 if(debug) cout<<"; ";
1176 return vec_sz;
1177}
1178
1179
1180RecMdcHit DotsHelixFitter::makeRecMdcHit(const MdcDigi* aDcDigi)
1181{
1182 updateDcDigiInfo(aDcDigi);
1183 RecMdcHit aRecMdcHit; // = new RecMdcHit;
1184 aRecMdcHit.setDriftDistLeft(-1*myDriftDist[0]);
1185 aRecMdcHit.setDriftDistRight(myDriftDist[1]);
1186 aRecMdcHit.setErrDriftDistLeft(myDriftDistErr[0]);
1187 aRecMdcHit.setErrDriftDistRight(myDriftDistErr[1]);
1188 aRecMdcHit.setChisqAdd(myDcChi*myDcChi);
1189 aRecMdcHit.setFlagLR(myLeftRight);
1190 //aRecMdcHit.setStat(stat);
1191 Identifier mdcId = aDcDigi->identify();
1192 aRecMdcHit.setMdcId(mdcId);
1193 double tdc = aDcDigi->getTimeChannel();
1194 double adc = aDcDigi->getChargeChannel();
1195 aRecMdcHit.setTdc(tdc);
1196 aRecMdcHit.setAdc(adc);
1197 aRecMdcHit.setDriftT(myDriftTime);
1198 // --- doca
1199 double doca=fabs(myDocaFromTrk);
1200 if(myLeftRight==0) doca*=-1;
1201 aRecMdcHit.setDoca(doca);
1202 aRecMdcHit.setEntra(myEntranceAngle);
1203 aRecMdcHit.setZhit(myPosOnWire[2]);
1204 aRecMdcHit.setFltLen(myFlightLength);
1205 return aRecMdcHit;
1206}
1207
1208RecMdcHit DotsHelixFitter::makeRecMdcHit(const MdcDigi* aDcDigi, KalmanFit::Helix aHelix)
1209{
1210 setInitialHelix(aHelix);
1211 return makeRecMdcHit(aDcDigi);
1212}
1213
1214
1215vector<RecMdcHit> DotsHelixFitter::makeRecMdcHitVec(int sel)
1216{
1217 vector<RecMdcHit> aRecMdcHitVec;
1218 int i_digi=0;
1219 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
1220 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
1221 {
1222 if(sel==1&&myMdcDigiIsActive[i_digi]!=1) continue;// skip inactive hits
1223 aRecMdcHitVec.push_back(makeRecMdcHit(*iter_mdcDigi));
1224 }
1225 return aRecMdcHitVec;
1226}
1227
1228
1229double DotsHelixFitter::IntersectCylinder(double r)
1230{
1231 double m_rad = myHelix->radius();
1232 double l = myHelix->center().perp();
1233
1234 double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
1235
1236 if(cosPhi < -1 || cosPhi > 1) return 0;
1237
1238 double dPhi = myHelix->center().phi() - acos(cosPhi) - myHelix->phi0();
1239 //cout<<"DotsHelixFitter::IntersectCylinder: center phi="<<myHelix->center().phi()<<", Phi="<<acos(cosPhi)<<", phi0="<<myHelix->phi0()<<", dphi="<<dPhi<<endl;
1240
1241 if(dPhi < -M_PI) dPhi += 2 * M_PI;
1242
1243 return dPhi;
1244}
1245
1246HepMatrix DotsHelixFitter::dxda_cgem(KalmanFit::Helix a, double phi)
1247{
1248 HepMatrix dXDA(3,5,0);
1249
1250 double alpha = a.alpha();
1251 double dr = a.dr();
1252 double phi0 = a.phi0();
1253 double kappa = a.kappa();
1254 double dz = a.dz();
1255 double tanl = a.tanl();
1256
1257 HepPoint3D pos = a.x(phi);
1258 double x = pos.x();
1259 double y = pos.y();
1260 double z = pos.z();
1261 double r2 = x*x+y*y;
1262
1263 double cosPhi=cos(phi);
1264 double sinPhi=sin(phi);
1265
1266 double dPhiDdr = -(kappa/alpha-cosPhi/(dr+alpha/kappa))/sinPhi;
1267 //double dPhiDdr2= -kappa*kappa*(2.0*alpha*dr+kappa*(dr*dr+r2))/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
1268 //cout<<"dPhiDdr = "<<dPhiDdr<<endl;
1269 //cout<<"dPhiDdr2= "<<dPhiDdr2<<endl;
1270 double dPhiDkappa = -kappa*(dr*dr-r2)*(2*alpha+dr*kappa)/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
1271
1272 double dxDdr = cos(phi0)+alpha/kappa*sin(phi0+phi)*dPhiDdr;
1273 double dyDdr = sin(phi0)-alpha/kappa*cos(phi0+phi)*dPhiDdr;
1274 double dxDphi0 = -dr*sin(phi0)+alpha/kappa*(-sin(phi0)+sin(phi0+phi));
1275 double dyDphi0 = dr*cos(phi0)+alpha/kappa*( cos(phi0)-cos(phi0+phi));
1276 double dxDkappa = -alpha/(kappa*kappa)*(cos(phi0)-cos(phi0+phi))+alpha/kappa*sin(phi0+phi)*dPhiDkappa;
1277 double dyDkappa = -alpha/(kappa*kappa)*(sin(phi0)-sin(phi0+phi))-alpha/kappa*cos(phi0+phi)*dPhiDkappa;
1278 double dzDdr = -alpha/kappa*tanl*dPhiDdr;
1279 double dzDkappa = alpha/(kappa*kappa)*tanl*phi-alpha/kappa*tanl*dPhiDkappa;
1280 double dzDtanl = -alpha/kappa*phi;
1281
1282 dXDA(1,1) = dxDdr;
1283 dXDA(1,2) = dxDphi0;
1284 dXDA(1,3) = dxDkappa;
1285 dXDA(2,1) = dyDdr;
1286 dXDA(2,2) = dyDphi0;
1287 dXDA(2,3) = dyDkappa;
1288 dXDA(3,1) = dzDdr;
1289 dXDA(3,3) = dzDkappa;
1290 dXDA(3,4) = 1.0;
1291 dXDA(3,5) = dzDtanl;
1292
1293 return dXDA;
1294}
1295
1296vector<double> DotsHelixFitter::CircleFitByTaubin(const vector< pair<double, double> >& pos_hits)
1297{
1298
1299 int nHits = pos_hits.size();
1300 // --- for Taubin fit
1301 int i, iter, IterMAX=99;
1302 double Xi,Yi,Zi;
1303 double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
1304 double A0,A1,A2,A22,A3,A33;
1305 double Dy,xnew,x,ynew,y;
1306 double DET,Xcenter,Ycenter;
1307
1308
1309 // --- Compute x- and y- sample means (via a function in the class "data")
1310 double meanX(0.), meanY(0.);
1311 for (i=0; i<nHits; i++)
1312 {
1313 meanX+=pos_hits[i].first;
1314 meanY+=pos_hits[i].second;
1315 }
1316 meanX/=nHits;
1317 meanY/=nHits;
1318
1319
1320 // --- computing moments
1321 Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
1322 for (i=0; i<nHits; i++)
1323 {
1324 Xi = pos_hits[i].first - meanX; // centered x-coordinates
1325 Yi = pos_hits[i].second - meanY; // centered y-coordinates
1326 Zi = Xi*Xi + Yi*Yi;
1327
1328 Mxy += Xi*Yi;
1329 Mxx += Xi*Xi;
1330 Myy += Yi*Yi;
1331 Mxz += Xi*Zi;
1332 Myz += Yi*Zi;
1333 Mzz += Zi*Zi;
1334 }
1335 Mxx /= nHits;
1336 Myy /= nHits;
1337 Mxy /= nHits;
1338 Mxz /= nHits;
1339 Myz /= nHits;
1340 Mzz /= nHits;
1341
1342 // --- computing coefficients of the characteristic polynomial
1343 Mz = Mxx + Myy;
1344 Cov_xy = Mxx*Myy - Mxy*Mxy;
1345 Var_z = Mzz - Mz*Mz;
1346 A3 = 4*Mz;
1347 A2 = -3*Mz*Mz - Mzz;
1348 A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
1349 A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
1350 A22 = A2 + A2;
1351 A33 = A3 + A3 + A3;
1352
1353 // finding the root of the characteristic polynomial
1354 // using Newton's method starting at x=0
1355 // (it is guaranteed to converge to the right root)
1356 for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
1357 {
1358 Dy = A1 + x*(A22 + A33*x);
1359 xnew = x - y/Dy;
1360 if ((xnew == x)||(!isfinite(xnew))) break;
1361 ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
1362 if (abs(ynew)>=abs(y)) break;
1363 x = xnew; y = ynew;
1364 }
1365 //if(iter==IterMAX) cout<<"IterMAX reached!"<<endl;
1366
1367 // --- computing paramters of the fitting circle
1368 DET = x*x - x*Mz + Cov_xy;
1369 Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
1370 Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
1371
1372 // --- assembling the output
1373 double xc_fit = Xcenter + meanX;
1374 double yc_fit = Ycenter + meanY;
1375 double r_fit = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
1376 //cout<<"fitted trk par: xc, yc, r = "<<xc_fit<<" ,"<<yc_fit<<", "<<r_fit<<endl;
1377 //cout<<"after "<<iter<<" iterations! "<<endl;
1378 //cout<<"CircleFitByTaubin "<<iter<<" iterations! "<<endl;
1379
1380 vector<double> output;
1381 output.push_back(xc_fit);
1382 output.push_back(yc_fit);
1383 output.push_back(r_fit);
1384 output.push_back(iter);
1385 output.push_back(meanX);
1386 output.push_back(meanY);
1387
1388 return output;
1389}
1390
1391vector<double> DotsHelixFitter::circleParConversion(const vector<double>& vec_par_Taubin)
1392{
1393 const vector<double>& par_new = vec_par_Taubin;
1394 double xc_fit = par_new[0];
1395 double yc_fit = par_new[1];
1396 double r_fit = par_new[2];
1397 double dr_fit = sqrt(xc_fit*xc_fit+yc_fit*yc_fit)-r_fit;
1398 double phi0_fit = atan2(yc_fit,xc_fit);
1399 double kappa_fit= myAlpha/r_fit;
1400 // --- charge
1401 double xmean = par_new[4];
1402 double ymean = par_new[5];
1403 Hep3Vector pos_mean(xmean, ymean);
1404 Hep3Vector pos_center(xc_fit, yc_fit);
1405 Hep3Vector vec_z=pos_mean.cross(pos_center);
1406 double charge=vec_z.z()/fabs(vec_z.z());
1407 if(charge>0)
1408 {
1409 dr_fit = -dr_fit;
1410 phi0_fit+=M_PI;
1411 kappa_fit = -kappa_fit;
1412 }
1413 while(phi0_fit< 0 ) phi0_fit+=2*M_PI;
1414 while(phi0_fit> 2*M_PI) phi0_fit-=2*M_PI;
1415
1416 vector<double> par;
1417 par.push_back(dr_fit);
1418 par.push_back(phi0_fit);
1419 par.push_back(kappa_fit);
1420
1421 return par;
1422}
1423
1424
1425vector<double> DotsHelixFitter::calculateCirclePar_Taubin(bool useIniHelix, int layerMin, int layerMax, bool useExtraPoints, int debug)
1426{
1427 //bool debug=false; //debug=true;
1428 //cout<<"calculateCirclePar_Taubin: "<<endl;
1429 int n_iter=0;
1430
1431 // --- fitted parameters
1432 double dr_fit = myHelix_a[0];
1433 double phi0_fit = myHelix_a[1];
1434 double kappa_fit= myHelix_a[2];
1435 //double dr_fit_last = 9999.;
1436 //double phi0_fit_last = 9999.;
1437 //double kappa_fit_last= 9999.;
1438 double totChi2=0;
1439 myNActiveHits=0;
1440 myNInnerXHits=0;
1441 myNOuterXHits=0;
1442
1443
1444 // --- prepare (x,y) pairs
1445 vector< pair<double, double> > pos_hits;
1446 vector< pair<double, double> > CGEM_pos_hits;
1447
1448 // --- CGEM
1449 int i_digi=0;
1450 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
1451 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
1452 {
1453 int layer = (*iter_cluster)->getlayerid();
1454 int flag = (*iter_cluster)->getflag();
1455 //cout<<"CGEM cluster "<<i_digi<<", layer "<<layer<<", flag "<<flag<<endl;
1456 if(flag!=0) continue;
1457 if(myCgemClusterIsActive[i_digi]!=1) continue;
1458 double phi_cluster = (*iter_cluster)->getrecphi();
1459 pair<double, double> aHitPos(myRmidDGapCgem[layer]*cos(phi_cluster), myRmidDGapCgem[layer]*sin(phi_cluster));
1460 pos_hits.push_back(aHitPos);
1461 //cout<<"add one position pair"<<endl;
1462 }
1463 //cout<<"filled pos CGEM"<<endl;
1464 CGEM_pos_hits=pos_hits;
1465 //cout<<"set CGEM_pos_hits"<<endl;
1466 //cout<<"CGEM_pos_hits.size = "<<CGEM_pos_hits.size()<<endl;
1467
1468 // --- layer, wire vector for MDC digi
1469 vector<int> vecLayer, vecWire;
1470 int nDigi = myVecMdcDigi.size();
1471 vecLayer.resize(nDigi);
1472 vecWire.resize(nDigi);
1473 bool layerFilled=false;
1474 double converged=-1.0;
1475 double lastChi2=-99.0;
1476 double last2Chi2=-99.0;
1477 while(1) {
1478 myNActiveOuterXHits=0;
1479 // --- reset pos_hits
1480 pos_hits.clear();
1481 pos_hits=CGEM_pos_hits;
1482
1483 for(int i=0; i<43; i++) myNumMdcDigiPerLayer[i]=0;
1484 totChi2=0;
1485 myMapFlylenIdx.clear();
1486 //myMapFlylenChi.clear();
1487
1488 // --- update helix parameters
1489 HepPoint3D pivot(0,0,0);
1490 HepVector a(5,0);
1491 a(1)=dr_fit;
1492 a(2)=phi0_fit;
1493 a(3)=kappa_fit;
1494 //a(4)=myHelix_a[3];
1495 //a(5)=myHelix_a[4];
1496 KalmanFit::Helix aHelix(pivot,a);
1497 setInitialHelix(aHelix);
1498
1499 // --- update chi of CGEM clusters
1500 HepPoint3D pos;
1501 double phi_trk, phi_clu;
1502 i_digi=0;
1503 iter_cluster = myVecCgemCluster.begin();
1504 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
1505 {
1506 int layer = (*iter_cluster)->getlayerid();
1507 int flag = (*iter_cluster)->getflag();
1508 //cout<<"CGEM cluster "<<i_digi<<", layer "<<layer<<", flag "<<flag<<endl;
1509 if(flag!=0) continue;
1510 if(myCgemClusterIsActive[i_digi]!=1) continue;
1511 double dphi = IntersectCylinder(myRmidDGapCgem[layer]);
1512 if(fabs(dphi)<1e-10) {
1513 myChiVecCgemCluster[i_digi]=-9999;
1514 if(debug) cout<<"trk has no intersection with CGEM layer "<<layer<<endl;
1515 continue;// track has no intersection with CGEM
1516 }
1517 pos = myHelix->x(dphi);
1518 myFlightLength = fabs(dphi*myHelix->radius());// in cm
1519 myMapFlylenIdx[myFlightLength]=-(i_digi+1);
1520 //if(debug) cout<<"dphi="<<dphi<<", pos="<<pos<<endl;
1521 phi_trk = pos.phi();// in radian
1522 phi_clu = (*iter_cluster)->getrecphi();
1523 double del_phi = phi_clu-phi_trk;
1524 if(del_phi<-M_PI||del_phi> M_PI) del_phi=atan2(sin(del_phi),cos(del_phi));
1525 double del_X=del_phi*getRmidGapCgem(layer);// in cm
1526 double Q(0.);// in fC
1527 double T(100);// not valid at present
1528 int mode(2);
1529 Q=(*iter_cluster)->getenergydeposit();
1530 Hep3Vector p3_trk = myHelix->momentum(dphi);
1531 double phi_p3trk = p3_trk.phi();
1532 double incidentPhi = phi_p3trk-phi_trk;
1533 //while(incidentPhi>CLHEP::pi) incidentPhi-=CLHEP::twopi;
1534 //while(incidentPhi<-CLHEP::pi) incidentPhi+=CLHEP::twopi;
1535 if(incidentPhi<-M_PI||incidentPhi> M_PI) incidentPhi=atan2(sin(incidentPhi),cos(incidentPhi));// into [-pi, pi]
1536 double err_m=myCgemCalibSvc->getSigma(layer,flag,mode,incidentPhi,Q,T)/10.;// in cm
1537 double chi_clu = del_X/err_m;
1538 if(debug) cout<<"CGEM cluster "<<i_digi<<", layer "<<layer
1539 <<", phi_clu, phi_trk = "<<phi_clu<<", "<<phi_trk
1540 <<", del_phi="<<del_phi<<", chi="<<chi_clu<<endl;
1541 myChiVecCgemCluster[i_digi]=chi_clu;
1542 totChi2+=chi_clu*chi_clu;
1543 }
1544
1545 // --- MDC
1546 i_digi=0;
1547 std::vector<Hep2Vector> wirePos;
1548 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
1549 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
1550 {
1551 //cout<<__FILE__<<": "<<__FUNCTION__<<": i_digi="<<i_digi<<endl;
1552 // --- get id, layer, wire ---
1553 Identifier id = (*iter_mdcDigi)->identify();
1554 int layer = MdcID::layer(id);
1555 int wire = MdcID::wire(id);
1556 if(!layerFilled) {
1557 vecLayer[i_digi]=layer;
1558 vecWire[i_digi]=wire;
1559 }
1560 //double charge = (*iter_mdcDigi)->getChargeChannel();
1561
1562 myNumMdcDigiPerLayer[layer]++;
1563
1564 // --- position from MDC digi
1565 double xpos(0), ypos(0);
1566 if(!useIniHelix&&n_iter==0) // --- use wire position
1567 {
1568 xpos=myEastPosWires[layer][wire][0];
1569 ypos=myEastPosWires[layer][wire][1];
1570 //pair<double, double> aHitPos(myEastPosWires[layer][wire][0], myEastPosWires[layer][wire][1]);
1571 if(myUsePhiAmbi) {
1572 double phiLR = myPhiAmbiVecMdcDigi.at(i_digi);
1573 if(phiLR>-1000) {
1574 calculateRoughDD(*iter_mdcDigi);
1575 xpos += myDocaFromDigi*cos(phiLR);
1576 ypos += myDocaFromDigi*sin(phiLR);
1577 }
1578 }
1579 }
1580 else {// --- hit position calculated iteratively
1581 updateDcDigiInfo(*iter_mdcDigi);
1582 myMapFlylenIdx[myFlightLength]=i_digi;
1583 myChiVecMdcDigi[i_digi]=myDcChi;
1584 HepPoint3D newPivot(myEastPosWires[layer][wire][0],myEastPosWires[layer][wire][1],0);
1585 aHelix.pivot(newPivot);
1586 double phi0_new = aHelix.phi0();
1587 double dr_new = aHelix.dr();
1588 double d_measured = myDriftDist[myLeftRight];
1589 if(myUsePhiAmbi) {
1590 double phiLR = myPhiAmbiVecMdcDigi.at(i_digi);
1591 double dphi = dPhi(phi0_new,phiLR);// compare phi_LR
1592 if(fabs(dphi)>0.5*M_PI) phi0_new+=M_PI;// update phi0_new
1593 }
1594 else {
1595 d_measured=dr_new/fabs(dr_new)*d_measured;// Left-right ambiguityfrom helix and wire position
1596 }
1597 xpos = myEastPosWires[layer][wire][0]+d_measured*cos(phi0_new);
1598 ypos = myEastPosWires[layer][wire][1]+d_measured*sin(phi0_new);
1599 }
1600 if(myWireFlag[layer]!=0) continue;
1601 if(layer<layerMin||layer>layerMax) continue;
1602 if(n_iter==0) {
1603 if(layer>=36) myNOuterXHits++;
1604 else myNInnerXHits++;
1605 }
1606 if(myMdcDigiIsActive[i_digi]!=1) continue;
1607
1608 if(debug){
1609 wirePos.push_back(Hep2Vector(myEastPosWires[layer][wire][0],myEastPosWires[layer][wire][1]));
1610 // usedPos.push_back(Hep2Vector(xpos,ypos));
1611 }
1612
1613 pair<double, double> aHitPos(xpos,ypos);
1614 pos_hits.push_back(aHitPos);
1615 totChi2+=myDcChi*myDcChi;
1616 if(layer>=36) myNActiveOuterXHits++;
1617 }
1618 layerFilled=true;
1619
1620 //cout<<"pos_hits.size = "<<pos_hits.size()<<endl;
1621 if(pos_hits.size()<3) {
1622 if(debug) cout<<"DotsHelixFitter::"<<__FUNCTION__<<": N_xhits<3"<<endl;
1623 converged=-2.;
1624 break;
1625 }
1626 if(useExtraPoints&&(n_iter==0||myWireCrossPointPersistent)) {
1627 pos_hits.insert(pos_hits.end(), myExtraPoints.begin(), myExtraPoints.end() );
1628 }
1629 if(n_iter>20) {
1630 //if(debug)
1631 {
1632 cout<<" Taubin not converge after "<< n_iter<<" iterations "<<endl;
1633 cout<<" Circle par from last Taubin fit: "<<a[0]
1634 <<", "<<a[1]
1635 <<", "<<a[2]
1636 <<endl;
1637 }
1638 converged=-3.;
1639 break;
1640 }
1641 if(n_iter>1) {
1642 double dPhi0=phi0_fit-a(2);
1643 dPhi0=fabs(atan2(sin(dPhi0),cos(dPhi0)));
1644 double dKappa=fabs((kappa_fit-a(3))/kappa_fit);
1645 //cout<<"diff = "<<fabs(dr_fit-a(1))<<", "<<dPhi0<<", "<<dKappa<<endl;
1646 if(n_iter>90) if(debug) cout<<"diff: "<<dr_fit-a(1)<<", "<<dPhi0<<", "<<dKappa<<", "<<lastChi2-totChi2<<endl;
1647 if(fabs(dr_fit-a(1))<0.0001&&dPhi0<0.0001&&dKappa<0.0001&&(fabs(lastChi2-totChi2)<0.2||fabs(last2Chi2-totChi2)<0.2))
1648 {
1649 if(debug) cout<<"Taubin converges at iteration "<<n_iter<<endl;
1650 converged=1.;
1651 break;
1652 }
1653 }
1654 last2Chi2=lastChi2;
1655 lastChi2=totChi2;
1656
1657 // --- update circle par from Taubin fit
1658 vector<double> par_new = CircleFitByTaubin(pos_hits);
1659 double xc_fit = par_new[0];
1660 double yc_fit = par_new[1];
1661 double r_fit = par_new[2];
1662 dr_fit = sqrt(xc_fit*xc_fit+yc_fit*yc_fit)-r_fit;
1663 phi0_fit = atan2(yc_fit,xc_fit);
1664 kappa_fit= myAlpha/r_fit;
1665 // --- charge
1666 double xmean = par_new[4];
1667 double ymean = par_new[5];
1668 Hep3Vector pos_mean(xmean, ymean);
1669 Hep3Vector pos_center(xc_fit, yc_fit);
1670 Hep3Vector vec_z=pos_mean.cross(pos_center);
1671 double charge=vec_z.z()/fabs(vec_z.z());
1672 if(charge>0)
1673 {
1674 dr_fit = -dr_fit;
1675 phi0_fit+=M_PI;
1676 kappa_fit = -kappa_fit;
1677 }
1678 while(phi0_fit< 0 ) phi0_fit+=2*M_PI;
1679 while(phi0_fit> 2*M_PI) phi0_fit-=2*M_PI;
1680 if(debug) cout<<"iter "<<n_iter<<", circlr par = "<<dr_fit<<", "<<phi0_fit<<", "<<kappa_fit<<", chi2 = "<<totChi2<<" , circle aka center=("<<xc_fit<<","<<yc_fit<<"), radius = "<<r_fit<<endl;
1681 extern bool global_show_this_event_detail;
1682 if(debug&&global_show_this_event_detail){
1683 cout<<"\twirePos: ";
1684 for (size_t i=0; i<wirePos.size(); i++){
1685 cout<<wirePos[i]<<" ";
1686 }
1687 cout<<"\n";
1688 cout<<"\tusedPos: ";
1689 for (size_t i=0; i<pos_hits.size();i++){
1690 cout<<"("<<pos_hits[i].first<<","<<pos_hits[i].second<<") ";
1691 }
1692 cout<<"\n";
1693 }
1694
1695 n_iter++;
1696 }// iterations
1697 myNActiveHits=pos_hits.size();
1698 myChi2=totChi2;
1699
1700 if(debug) cout<<setw(20)<<"hit chi (fitted): ";
1701 // --- print chi of CGEM clusters
1702 i_digi=0;
1703 iter_cluster = myVecCgemCluster.begin();
1704 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
1705 {
1706 int layer = (*iter_cluster)->getlayerid();
1707 int flag = (*iter_cluster)->getflag();
1708 //cout<<"CGEM cluster "<<i_digi<<", layer "<<layer<<", flag "<<flag<<endl;
1709 if(flag!=0) continue;
1710 if(debug) cout<<" "<<myChiVecCgemCluster[i_digi]<<"("<<myCgemClusterIsActive[i_digi]<<")Cgem"<<layer;
1711 }
1712
1713 // --- print chi of MDC digis
1714 int i_idx = 0;
1715 int n_idx=myMapFlylenIdx.size();
1716 map<double, int>::iterator iter_idx = myMapFlylenIdx.begin();
1717 int gapSize=0;
1718 int gapMax=0;
1719 double Smax=9999.0;
1720 double Stmp=0;
1721 for(; iter_idx!=myMapFlylenIdx.end(); iter_idx++, i_idx++)
1722 {
1723 if(debug) if(i_idx>0&&i_idx%4==0) cout<<endl<<setw(20)<<" ";
1724 int i_hit = iter_idx->second;
1725 if(i_hit>=0) {
1726 if(debug) cout<<" (L"<<vecLayer[i_hit]<<",W"<<vecWire[i_hit]<<")Chi"<<myChiVecMdcDigi[i_hit]<<"("<<myMdcDigiIsActive[i_hit]<<")S"<<iter_idx->first;
1727 if(myMdcDigiIsActive[i_hit]<=0) gapSize++;
1728 else {
1729 if(gapMax<gapSize) gapMax=gapSize;
1730 if(gapSize>=2&&Stmp<Smax) Smax=Stmp;
1731 Stmp=iter_idx->first;
1732 gapSize=0;
1733 }
1734 }
1735 }
1736 if(gapMax<gapSize) gapMax=gapSize;
1737 if(gapSize>=2&&Stmp<Smax) Smax=Stmp;
1738 //if(i_idx%10!=1)
1739 if(debug){
1740 cout<<endl;
1741 cout<<"Smax="<<Smax<<endl;
1742 }
1743 //iter_idx = myMapFlylenIdx.find(Smax);
1744 //if(iter_idx!=myMapFlylenIdx.end()) iter_idx++;
1745 //for(; iter_idx!=myMapFlylenIdx.end(); iter_idx++)
1746 //{
1747 // int i_hit = iter_idx->second;
1748 // if(i_hit>=0) {
1749 // myMdcDigiIsActive[i_hit]=-1;
1750 // cout<<i_hit<<"th X-hit deactived, ";
1751 // }
1752 //}
1753 //cout<<endl;
1754
1755 // --- fill parameters for helix
1756 vector<double> par_fit;
1757 par_fit.push_back(dr_fit);
1758 par_fit.push_back(phi0_fit);
1759 par_fit.push_back(kappa_fit);
1760 par_fit.push_back(converged);
1761 return par_fit;
1762}
1763
1764bool DotsHelixFitter::getPosAtCgem(int layer, HepPoint3D& pos)
1765{
1766 bool sc=true;
1767 if(layer>=0&&layer<3) {
1768 double dphi = IntersectCylinder(myRmidDGapCgem[layer]);
1769 if(fabs(dphi)<1e-10) sc=false;
1770 else {
1771 pos = myHelix->x(dphi);
1772 }
1773 }
1774 else sc=false;
1775
1776 return sc;
1777}
1778
1779
1780void DotsHelixFitter::printActiveHits()
1781{
1782
1783}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:87
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
global_show_this_event_detail
bool global_show_this_event_detail
Definition DotsConnection.cxx:25
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
alpha
const double alpha
Definition FastVertexFit.cxx:4
output
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in output
Definition FoamA.h:89
s
XmlRpcServer s
Definition HelloServer.cpp:11
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
M_PI
#define M_PI
Definition TConstant.h:4
_
#define _(A, B)
Definition cfortran.h:44
CgemGeoLayer::getMiddleROfGapD
double getMiddleROfGapD() const
Definition CgemGeoLayer.h:105
CgemGeoLayer::getInnerROfAnodeCu1
double getInnerROfAnodeCu1() const
Definition CgemGeoLayer.h:163
CgemGeoLayer::getNumberOfSheet
int getNumberOfSheet() const
Definition CgemGeoLayer.h:215
CgemGeoLayer::getAngleOfStereo
double getAngleOfStereo() const
Definition CgemGeoLayer.h:218
CgemGeoLayer::getInnerROfAnodeCu2
double getInnerROfAnodeCu2() const
Definition CgemGeoLayer.h:167
CgemGeoReadoutPlane::getVFromPhiZ
double getVFromPhiZ(double phi, double z, bool checkRange=true) const
Definition CgemGeoReadoutPlane.h:178
CgemGeoReadoutPlane::getVFromPhiZ_nearPhiMin
double getVFromPhiZ_nearPhiMin(double phi, double z, bool checkRange=true) const
Definition CgemGeoReadoutPlane.h:189
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
CgemGeomSvc::getNumberOfCgemLayer
int getNumberOfCgemLayer() const
Definition CgemGeomSvc.h:45
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition CgemGeomSvc.h:140
Constants::c
static const double c
Definition Constants.h:43
DotsHelixFitter::setInitialHelix
void setInitialHelix(KalmanFit::Helix aHelix)
Definition DotsHelixFitter.cxx:881
DotsHelixFitter::activeHits
int activeHits(double chi_cut=10, int sel=0)
activate myMdcDigiIsActive's hit based on myChiVecMdcDigi
Definition DotsHelixFitter.cxx:825
DotsHelixFitter::setDChits
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
Definition DotsHelixFitter.cxx:160
DotsHelixFitter::dPhi
double dPhi(double phi1, double phi2)
Definition DotsHelixFitter.h:51
DotsHelixFitter::setCgemClusters
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
Definition DotsHelixFitter.cxx:214
DotsHelixFitter::getPosAtCgem
bool getPosAtCgem(int layer, HepPoint3D &pos)
Definition DotsHelixFitter.cxx:1764
DotsHelixFitter::deactiveHits
int deactiveHits(double chi_cut=10, int nMax=1, int layerMin=0, int layerMax=50)
Definition DotsHelixFitter.cxx:736
DotsHelixFitter::dxda_cgem
HepMatrix dxda_cgem(KalmanFit::Helix a, double phi)
Definition DotsHelixFitter.cxx:1246
DotsHelixFitter::setDChitsPhiAmbi
bool setDChitsPhiAmbi(const vector< double > &vec)
Definition DotsHelixFitter.cxx:203
DotsHelixFitter::calculateDocaFromTrk
void calculateDocaFromTrk(const MdcDigi *aDcDigi)
Definition DotsHelixFitter.cxx:928
DotsHelixFitter::circleParConversion
vector< double > circleParConversion(const vector< double > &vec_par_Taubin)
Definition DotsHelixFitter.cxx:1391
DotsHelixFitter::getRmidGapCgem
double getRmidGapCgem(int i)
Definition DotsHelixFitter.h:112
DotsHelixFitter::getSZ
vector< pair< double, double > > getSZ(const MdcDigi *aDcDigi)
return S, the track length in XY; and Z of a hit on a stereo wire
Definition DotsHelixFitter.cxx:1127
DotsHelixFitter::setDChitsFitFlag
bool setDChitsFitFlag(vector< int > &vec)
set myMdcDigiIsActive from vec
Definition DotsHelixFitter.cxx:192
DotsHelixFitter::CircleFitByTaubin
vector< double > CircleFitByTaubin(const vector< pair< double, double > > &pos_hits)
Definition DotsHelixFitter.cxx:1296
DotsHelixFitter::loadOneDcDigi
void loadOneDcDigi(const MdcDigi *aDcDigi)
Definition DotsHelixFitter.cxx:892
DotsHelixFitter::calculateCirclePar_Taubin
vector< double > calculateCirclePar_Taubin(bool useIniHelix=false, int layerMin=-1, int layerMax=50, bool useExtraPoints=false, int debug=0)
Definition DotsHelixFitter.cxx:1425
DotsHelixFitter::makeRecMdcHitVec
vector< RecMdcHit > makeRecMdcHitVec(int sel=1)
Definition DotsHelixFitter.cxx:1215
DotsHelixFitter::calculateRoughDD
void calculateRoughDD(const MdcDigi *aMdcDigi)
Definition DotsHelixFitter.cxx:913
DotsHelixFitter::makeRecMdcHit
RecMdcHit makeRecMdcHit(const MdcDigi *aDcDigi)
Definition DotsHelixFitter.cxx:1180
DotsHelixFitter::IntersectCylinder
double IntersectCylinder(double r)
Definition DotsHelixFitter.cxx:1229
DotsHelixFitter::updateDcDigiInfo
void updateDcDigiInfo(const MdcDigi *aDcDigi)
based on input aDcDigi, modifies myFlightLength, myLeftRight, myPosOnWire, myDocaFromDigi,...
Definition DotsHelixFitter.cxx:960
ICgemCalibFunSvc::getSigma
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
KalmanFit::Helix::bFieldZ
double bFieldZ(double)
sets/returns z componet of the magnetic field.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:291
KalmanFit::Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:205
KalmanFit::Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/src/helix/Helix.cxx:206
KalmanFit::Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:271
KalmanFit::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/src/helix/Helix.cxx:151
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:229
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:265
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:211
MdcCalibFunSvc::getSigma
double getSigma(int layid, int lr, double dist, double entrance=0.0, double tanlam=0.0, double z=0.0, double Q=1000.0) const
Definition MdcCalibFunSvc.cxx:278
MdcCalibFunSvc::getT0
double getT0(int layid, int cellid) const
Definition MdcCalibFunSvc.cxx:507
MdcCalibFunSvc::driftTimeToDist
double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const
Definition MdcCalibFunSvc.cxx:149
MdcCalibFunSvc::getTprop
double getTprop(int lay, double z) const
Definition MdcCalibFunSvc.cxx:143
MdcCalibFunSvc::getTimeWalk
double getTimeWalk(int layid, double Q) const
Definition MdcCalibFunSvc.cxx:514
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:770
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:786
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:54
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition RawDataUtil.h:8
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition RawData.cxx:45
RawData::getTimeChannel
unsigned int getTimeChannel() const
Definition RawData.cxx:40
RecMdcHit::setMdcId
void setMdcId(Identifier mdcid)
Definition RecMdcHit.h:67
RecMdcHit::setErrDriftDistRight
void setErrDriftDistRight(double erddr)
Definition RecMdcHit.h:63
RecMdcHit::setFltLen
void setFltLen(double fltLen)
Definition RecMdcHit.h:74
RecMdcHit::setErrDriftDistLeft
void setErrDriftDistLeft(double erddl)
Definition RecMdcHit.h:62
RecMdcHit::setDriftDistLeft
void setDriftDistLeft(double ddl)
Definition RecMdcHit.h:60
RecMdcHit::setDoca
void setDoca(double doca)
Definition RecMdcHit.h:71
RecMdcHit::setTdc
void setTdc(double tdc)
Definition RecMdcHit.h:68
RecMdcHit::setAdc
void setAdc(double adc)
Definition RecMdcHit.h:69
RecMdcHit::setFlagLR
void setFlagLR(int lr)
Definition RecMdcHit.h:65
RecMdcHit::setChisqAdd
void setChisqAdd(double pChisq)
Definition RecMdcHit.h:64
RecMdcHit::setZhit
void setZhit(double zhit)
Definition RecMdcHit.h:73
RecMdcHit::setDriftT
void setDriftT(double driftT)
Definition RecMdcHit.h:70
RecMdcHit::setDriftDistRight
void setDriftDistRight(double ddr)
Definition RecMdcHit.h:61
RecMdcHit::setEntra
void setEntra(double entra)
Definition RecMdcHit.h:72
TrkFitFun::getDeriLoc
bool getDeriLoc(int ipar, double helix[], double &deri, double wpos[], double zini)
Definition TrkFitFun.cxx:309
TrkFitFun::getDoca
bool getDoca(double trkpar[], double wpos[], double &doca, double whitPos[], double zini)
doca: distance of closest approach.
Definition TrkFitFun.cxx:180
TrkFitFun::CC
const double CC
Definition TrkFitFun.h:26
vec
dble_vec_t vec[12]
Definition ranlxd.c:372