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