CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
DotsHelixFitter Class Reference

#include <DotsHelixFitter.h>

Public Member Functions

 DotsHelixFitter ()
 
 ~DotsHelixFitter ()
 
 DotsHelixFitter (KalmanFit::Helix iniHelix, vector< const MdcDigi * > vecMdcDigi, vector< const RecCgemCluster * > vecCgemCluster)
 
void initialize ()
 
void setMaxIterations (int n=10)
 
void setMinXhitsInCircleFit (int n=3)
 
void setMinVhitsInHelixFit (int n=2)
 
void setMinHitsInHelixFit (int n=5)
 
void setDchi2Converge (double cut=5)
 
void setDchi2Diverge (double cut=50)
 
void setMaxNChi2Increase (int n=2)
 
void setChi2Diverge (double cut=1000000)
 
void setInitialHelix (KalmanFit::Helix aHelix)
 
void setT0 (double T0)
 
void setCgemClusters (vector< const RecCgemCluster * > aVecCgemCluster)
 
void setDChits (vector< const MdcDigi * > aVecMdcDigi, double T0)
 
void fitCircleOnly (bool x=true)
 
void fitModeCircle ()
 
void fitModeHelix ()
 
void useAxialHitsOnly (bool x=true)
 
KalmanFit::Helix getClassHelix ()
 
HepVector getHelix ()
 
const HepSymMatrix & getEa ()
 
double getChi2 ()
 
void fit ()
 
void calculateInitialHelix ()
 
void calcuMeasDeriv ()
 
int calculateNewHelix ()
 
int deactiveHits (double chi_cut=10, int nMax=1)
 
int activeHits (double chi_cut=10)
 
int deactiveHits (int layer_max=43, int nHit_max=1000)
 
int getNActiveHits ()
 
void updateChi2 ()
 
void loadOneDcDigi (const MdcDigi *aDcDigi)
 
void calculateDocaFromTrk (const MdcDigi *aDcDigi)
 
void updateDcDigiInfo (const MdcDigi *aDcDigi)
 
double getFlightLength ()
 
double getFlightLength (const MdcDigi *aDcDigi)
 
int getLayer (const MdcDigi *aDcDigi)
 
double getDocaFromTrk ()
 
double getDocaFromTrk (const MdcDigi *aDcDigi)
 
RecMdcHit makeRecMdcHit (const MdcDigi *aDcDigi)
 
RecMdcHit makeRecMdcHit (const MdcDigi *aDcDigi, KalmanFit::Helix aHelix)
 
vector< RecMdcHit > makeRecMdcHitVec (int sel=1)
 
vector< const MdcDigi * > getVecMdcDigi ()
 
vector< double > getVecChiMdcDigi ()
 
vector< int > getVecMdcDigiIsAct ()
 
vector< const RecCgemCluster * > getVecCgemCluster ()
 
vector< int > getVecIsActiveCgemCluster ()
 
vector< double > getVecChiCgemCluster ()
 
double IntersectCylinder (double r)
 
HepMatrix dxda_cgem (KalmanFit::Helix a, double phi)
 
double getRmidGapCgem (int i)
 

Detailed Description

Definition at line 16 of file DotsHelixFitter.h.

Constructor & Destructor Documentation

◆ DotsHelixFitter() [1/2]

DotsHelixFitter::DotsHelixFitter ( )

Definition at line 12 of file DotsHelixFitter.cxx.

12 :
13 myHelix(NULL),myIsInitialized(false),myMdcGeomSvc(NULL),myCgemGeomSvc(NULL),myMdcCalibFunSvc(NULL),myMdcUtilitySvc(NULL)
14{
15 myHelix_Ea=HepSymMatrix(5,0);
16 myMaxIteration=10;
17 myMinXhitsInCircleFit=3;
18 myMinVhitsInHelixFit=2;
19 myMinHitsInHelixFit=5;
20 myDchi2Converge=5.0;
21 myDchi2Diverge=50.0;
22 myMaxNChi2Increase=2;
23 myChi2Diverge=1000000.0;
24 myFitCircleOnly=false;
25 myUseAxialHitsOnly=false;
26}

◆ ~DotsHelixFitter()

DotsHelixFitter::~DotsHelixFitter ( )

Definition at line 29 of file DotsHelixFitter.cxx.

30{
31 if(myHelix) delete myHelix;
32}

◆ DotsHelixFitter() [2/2]

DotsHelixFitter::DotsHelixFitter ( KalmanFit::Helix  iniHelix,
vector< const MdcDigi * >  vecMdcDigi,
vector< const RecCgemCluster * >  vecCgemCluster 
)

Member Function Documentation

◆ activeHits()

int DotsHelixFitter::activeHits ( double  chi_cut = 10)

Definition at line 650 of file DotsHelixFitter.cxx.

651{
652 int nActived=0;
653 int i_digi=0;
654 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
655 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
656 {
657 if(!myMdcDigiIsActive[i_digi])
658 {
659 double abs_chi=fabs(myChiVecMdcDigi[i_digi]);
660
661 if( abs_chi<chi_cut )
662 {
663 myMdcDigiIsActive[i_digi] = 1;
664 nActived++;
665 }
666 }
667 }
668 return nActived;
669}

◆ calculateDocaFromTrk()

void DotsHelixFitter::calculateDocaFromTrk ( const MdcDigi *  aDcDigi)

Definition at line 734 of file DotsHelixFitter.cxx.

735{
736 bool debug = false;
737
738 loadOneDcDigi(aDcDigi);
739
740 // --- get doca from track parameters ---
741 getDoca(myHelix_a, myWirePos, myDocaFromTrk, myPosOnWire, myPosOnWire[2]);
742 if(debug)
743 {
744 cout<<"DotsHelixFitter::UpdateDcDigiInfo(): "<<endl;
745 cout<<"doca = "<<myDocaFromTrk<<endl;
746 cout<<"point on wire: "<<myPosOnWire[0]<<", "<<myPosOnWire[1]<<", "<<myPosOnWire[2]<<endl;
747 }
748
749 // --- flight length ---
750 KalmanFit::Helix aHelix = *myHelix;
751 HepPoint3D aNewPivot(myPosOnWire[0],myPosOnWire[1],myPosOnWire[2]);
752 aHelix.pivot(aNewPivot);
753 double newPhi0 = aHelix.phi0();
754 double dphi = newPhi0-myHelix->phi0();
755 while(dphi<-M_PI) dphi+=2*M_PI;
756 while(dphi> M_PI) dphi-=2*M_PI;
757 myFlightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
758}
M_PI
#define M_PI
Definition: TConstant.h:4
DotsHelixFitter::loadOneDcDigi
void loadOneDcDigi(const MdcDigi *aDcDigi)
Definition: DotsHelixFitter.cxx:713
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:211
TrkFitFun::getDoca
bool getDoca(double trkpar[], double wpos[], double &doca, double whitPos[], double zini)
Definition: TrkFitFun.cxx:169

Referenced by getDocaFromTrk(), and getFlightLength().

◆ calculateInitialHelix()

void DotsHelixFitter::calculateInitialHelix ( )

◆ calculateNewHelix()

int DotsHelixFitter::calculateNewHelix ( )

Definition at line 177 of file DotsHelixFitter.cxx.

178{
179 bool debug = false; //debug=true;
180
181 int state = 0;
182
183 double tension = 9999.;
184
185 int nPar=5;
186 if(myFitCircleOnly) nPar=3;
187
188 // --- check hits
189 //if(debug) cout<<"DotsHelixFitter::calculateNewHelix() starts checking hits ... "<<endl;
190 for(int i=0; i<43; i++) myNumMdcDigiPerLayer[i]=0;
191 myNActiveHits=0;
192 int nXHits[3]={0,0,0}, nVHits[3]={0,0,0};
193 int nDigi=myVecMdcDigi.size();
194 double* vec_zini = new double[nDigi];
195 int maxLayer = 0;
196 int i_digi=0;
197 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
198 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
199 {
200 // --- initialize vec_zini ---
201 vec_zini[i_digi]=0;
202 // --- get id, layer, wire ---
203 Identifier id = (*iter_mdcDigi)->identify();
204 int layer = MdcID::layer(id);
205 int wire = MdcID::wire(id);
206 // --- counting active hits
207 if(myMdcDigiIsActive[i_digi])
208 {
209 myNActiveHits++;
210
211 int hitFlag = myWireFlag[layer];
212 if(hitFlag==0)
213 {
214 nXHits[0]++;// total
215 nXHits[1]++;// MDC
216 }
217 else
218 {
219 nVHits[0]++;
220 nVHits[1]++;
221 }
222 }
223 if(maxLayer<layer) maxLayer=layer;
224 int nHits = myNumMdcDigiPerLayer[layer];
225 if(nHits<2) {
226 myIdxMdcDigiNeighbour[layer][nHits]=i_digi;
227 myWireIdMdcDigiNeighbour[layer][nHits]=wire;
228 }
229 myNumMdcDigiPerLayer[layer]++;
230 }
231 //if(debug) cout<<"DotsHelixFitter::calculateNewHelix() ends checking hits ... "<<endl;
232 HepPoint3D farPoint = myHelix->x(M_PI);
233 double maxRTrk = farPoint.perp();
234 for(int i=0; i<43; i++)
235 {
236 if(maxRTrk<myRWires[i]+2) break;// where soft track turning back
237 if(myNumMdcDigiPerLayer[i]==2)
238 {
239 int wire_1 = myWireIdMdcDigiNeighbour[i][0];
240 int wire_2 = myWireIdMdcDigiNeighbour[i][1];
241 int delta_n = abs(wire_1-wire_2);
242 int ambi_1 = 0;
243 int ambi_2 = 0;
244 if(delta_n==1)
245 {
246 ambi_1 = wire_1>wire_2? 1:-1;
247 ambi_2 = wire_1>wire_2? -1:1;
248 }
249 else if(delta_n==myNWire[i]-1)
250 {
251 ambi_1 = wire_1<=1? 1:-1;
252 ambi_2 = wire_2<=1? 1:-1;
253 }
254 if(debug)
255 cout<<"layer, wire1, wire2, ambi1, ambi2 = "
256 <<setw(5)<<i
257 <<setw(5)<<wire_1
258 <<setw(5)<<wire_2
259 <<setw(5)<<ambi_1
260 <<setw(5)<<ambi_2
261 <<endl;
262 // --- fix Ambiguity for neighboured mdc hits
263 //myAmbiguityMdcDigi[myIdxMdcDigiNeighbour[i][0]]=ambi_1;
264 //myAmbiguityMdcDigi[myIdxMdcDigiNeighbour[i][1]]=ambi_2;
265 }
266 }
267 i_digi=0;
268 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
269 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
270 {
271 int flag = (*iter_cluster)->getflag();
272 if(flag==0)
273 {
274 nXHits[0]++;// total
275 nXHits[2]++;// CGEM
276 }
277 else if(flag==1)
278 {
279 nVHits[0]++;
280 nVHits[2]++;
281 }
282 myNActiveHits++;
283 }
284 if(myFitCircleOnly)
285 {
286 if(nXHits[0]<myMinXhitsInCircleFit) return 1;
287 }
288 else
289 {
290 if((nXHits[0]+nVHits[0])<myMinHitsInHelixFit) return 2;
291 if(nVHits[0]<myMinVhitsInHelixFit) return 3;
292 }
293
294
295
296
297 // ---> start fit
298 double lastTotChi2 = 999999;
299 int nIterations = 0;
300 int n_chi2_increase = 0;
301
302 while(1) // iterations
303 {
304 myMapFlylenIdx.clear();
305 // --- matrix U P J ---
306 HepSymMatrix U(nPar, 0);
307 //HepMatrix P(5,1,0);
308 //HepMatrix J(5,1,0), JT(1,5,0);
309 //HepMatrix J_dmdx(1,3,0), J_dxda(3,5,0);
310 HepMatrix P(nPar,1,0);
311 HepMatrix J(nPar,1,0), JT(1,nPar,0);
312 HepMatrix J_dmdx(1,3,0), J_dxda(3,nPar,0);
313
314 // --- loop MDC hits ---
315 double totChi2=0;
316 i_digi=0;
317 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
318 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
319 {
320 // --- get id, layer, wire ---
321 Identifier id = (*iter_mdcDigi)->identify();
322 int layer = MdcID::layer(id);
323 int wire = MdcID::wire(id);
324 double charge = (*iter_mdcDigi)->getChargeChannel();
325
326 // --- get doca from track parameters ---
327 //tension=myTensionWires[layer][wire];
328 double wpos[7]={myEastPosWires[layer][wire][0], myEastPosWires[layer][wire][1], myEastPosWires[layer][wire][2]
329 , myWestPosWires[layer][wire][0], myWestPosWires[layer][wire][1], myWestPosWires[layer][wire][2], tension};
330 double doca_trk;
331 double whitPos[3];// approching position on wire (in cm)
332 if(debug)
333 {
334 //cout<<"a = "<<myHelix_aVec<<endl;
335 cout<<"* layer "<<layer<<", wire "<<wire<<", is active "<<myMdcDigiIsActive[i_digi] <<endl;
336 //cout<<"wire positions : ("<<wpos[0]<<", "<<wpos[1]<<", "<<wpos[2]<<") ("<<wpos[3]<<", "<<wpos[4]<<", "<<wpos[5]<<")"<<endl;
337 //cout<<"zini = "<<vec_zini[i_digi]<<endl;
338 }
339 getDoca(myHelix_a, wpos, doca_trk, whitPos, vec_zini[i_digi]);
340 //double doca_trk2 = myMdcUtilitySvc->doca(layer, wire, myHelix_aVec, myHelix_Ea, false, false);
341 if(debug) {
342 cout<<"doca = "<<doca_trk<<endl;
343 //cout<<"doca2 = "<<doca_trk2<<endl;
344 }
345 int leftRight = 2;
346 if(doca_trk!=0) {
347 leftRight=int(doca_trk/fabs(doca_trk));
348 }
349 /*
350 if(myAmbiguityMdcDigi[i_digi]!=0) {
351 leftRight=myAmbiguityMdcDigi[i_digi];
352 //if(debug) cout<<"fix leftRight = "<<leftRight<<endl;
353 }*/
354 int signDoca=1;
355 signDoca=leftRight/fabs(leftRight);
356 // --- conversion of left-right into the BESIII convention
357 // if(leftRight==-1) leftRight=0;
358 if(leftRight==1) leftRight=0;// fixed 2020-11-26
359 else leftRight=1;
360 //if(debug) cout<<"leftRight = "<<leftRight<<endl;
361
362 vec_zini[i_digi] = whitPos[2];// update vec_zini
363
364 // --- get measured doca --- tof in ns, driftTime in ns, T0Walk in ns
365 double rawTime = RawDataUtil::MdcTime((*iter_mdcDigi)->getTimeChannel());
366 double tprop = myMdcCalibFunSvc->getTprop(layer, vec_zini[i_digi]);
367 double T0Walk = myMdcCalibFunSvc->getT0(layer,wire) + myMdcCalibFunSvc->getTimeWalk(layer, charge);
368 // --- time of flight (TOF) ---
369 KalmanFit::Helix aHelix = *myHelix;
370 HepPoint3D aNewPivot(whitPos[0],whitPos[1],whitPos[2]);
371 aHelix.pivot(aNewPivot);
372 double newPhi0 = aHelix.phi0();
373 double dphi = newPhi0-myHelix->phi0();
374 while(dphi<-M_PI) dphi+=2*M_PI;
375 while(dphi> M_PI) dphi-=2*M_PI;
376 double flightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
377 myMapFlylenIdx[flightLength]=i_digi;
378 HepLorentzVector p4_pi = myHelix->momentum(dphi, 0.13957);
379 double speed = p4_pi.beta()*CC;// cm/second
380 double TOF = flightLength/speed*1.e9;// in ns
381 // --- drift time ---
382 double driftT = rawTime - myEventT0 - TOF - T0Walk - tprop;
383 //if(debug) cout<<"myEventT0, driftT = "<<myEventT0<<", "<<driftT<<endl;
384 // --- entrance Angle ---
385 double phiWire = atan2(whitPos[1],whitPos[0]);
386 double phiP = p4_pi.phi();
387 double entranceAngle = phiP-phiWire;
388 while(entranceAngle<-M_PI) entranceAngle+=2*M_PI;
389 while(entranceAngle> M_PI) entranceAngle-=2*M_PI;
390 // --- measured drift distance ---
391 double doca_hit = 0.1 * myMdcCalibFunSvc->driftTimeToDist(driftT,layer,wire,leftRight, entranceAngle);// in cm
392 // --- get measurement error ---
393 double docaErr_hit = 0.1 * myMdcCalibFunSvc->getSigma(layer, leftRight, doca_hit, entranceAngle, myHelix_a[4], vec_zini[i_digi], charge);// in cm
394 //if(debug) cout<<"error_doca_hit = "<<docaErr_hit<<endl;
395
396 // --- get derivatives, calculate J, update P, U ---
397 doca_hit*=signDoca;
398 //if(debug) cout<<"doca_hit = "<<doca_hit<<endl;
399 double delD = doca_hit-doca_trk;
400 double chi = delD/docaErr_hit;
401 if(debug)
402 cout<<"delta_m, err_m, chi = "<<delD<<", "<<docaErr_hit<<", "<<chi<<endl;
403 myChiVecMdcDigi[i_digi]=chi;
404
405 if(!myMdcDigiIsActive[i_digi]) continue;
406 if(myUseAxialHitsOnly && myWireFlag[layer]!=0) continue;
407
408 totChi2+=chi*chi;
409 for(int ipar=0; ipar<nPar; ipar++)
410 {
411 double deri;
412 //if(debug) cout<<"ipar = "<<ipar<<endl;
413 getDeriLoc(ipar,myHelix_a, deri, wpos, vec_zini[i_digi]);// in cm
414 J(ipar+1,1) = deri;
415 //if(debug) cout<<"d(doca)/d(a"<<"_"<<ipar<<") = "<<J(ipar+1,1)<<endl;
416 //P(1,ipar+1) += J(1,ipar+1)*(doca_hit-doca_trk)/(docaErr_hit*docaErr_hit);
417 //P(ipar+1,1) += J(ipar+1,1)*chi/docaErr_hit;
418 double scale=1; //if(fabs(chi)>5) scale=fabs(chi*chi);
419 P(ipar+1,1) += J(ipar+1,1)*(doca_hit-doca_trk)/(docaErr_hit*docaErr_hit*scale);// large error for large chi
420 for(int ipar2=0; ipar2<=ipar; ipar2++)
421 {
422 //U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(docaErr_hit*docaErr_hit);
423 U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(docaErr_hit*docaErr_hit*scale);// large error for large chi2
424 }
425 }
426 }// loop MDC hits
427 if(debug) cout<<"end of MDC hits loop"<<endl;
428 // <--- end of MDC hits loop ---
429
430 // --- loop CGEM clusters ---
431 i_digi=0;
432 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
433 for(; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i_digi++)
434 {
435 // --- get layer, cluster flag ---
436 int layer = (*iter_cluster)->getlayerid();
437 int flag = (*iter_cluster)->getflag();
438 if(myUseAxialHitsOnly && flag!=0) continue;
439 int sheet = (*iter_cluster)->getsheetid();
440 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(layer,sheet);
441 if(debug)
442 {
443 cout<<endl<<"layer, sheet, flag = "<<setw(5)<<layer<<setw(5)<<sheet<<setw(5)<<flag<<endl;
444 }
445
446 // --- get position & entrance Angle from track parameters ---
447 double dphi = IntersectCylinder(myRmidDGapCgem[layer]);
448 HepPoint3D pos = myHelix->x(dphi);
449 //cout<<"dphi="<<dphi<<", pos="<<pos<<endl;
450 double phi_trk = pos.phi();// in radian
451 double z_trk = pos.z();// in cm
452 J_dxda = dxda_cgem(*myHelix, dphi);
453 Hep3Vector p3_trk = myHelix->momentum(dphi);
454 double phi_p3trk = p3_trk.phi();
455 double incidentPhi = phi_p3trk-phi_trk;
456 while(incidentPhi>CLHEP::pi) incidentPhi-=CLHEP::twopi;
457 while(incidentPhi<-CLHEP::pi) incidentPhi+=CLHEP::twopi;
458
459 // --- get X-cluster charge, measured X, derivatives, calculate J, update P, U ---
460 double phi_CC(0.0);// in radian
461 double X_CC(0.), V_CC(0.);// in cm
462 double delta_m, err_m;
463 double Q(0.);// in fC
464 double T(100);// not valid at present
465 int mode(2);
466 Q=(*iter_cluster)->getenergydeposit();
467 if(flag==0) {
468 phi_CC = (*iter_cluster)->getrecphi();
469 //phi_CC = (*iter_cluster)->getrecphi_CC();
470 double del_phi = phi_CC-phi_trk;
471 while(del_phi<-M_PI) del_phi+=2*M_PI;
472 while(del_phi> M_PI) del_phi-=2*M_PI;
473 X_CC=phi_CC*myRmidDGapCgem[layer];
474 if(debug) cout<<"phi_trk, phi_rec = "<<phi_trk<<", "<<phi_CC<<endl;
475 //double dX = del_phi*myRmidDGapCgem[layer];
476 //double dX = del_phi;
477 delta_m=del_phi;
478 err_m=myCgemCalibSvc->getSigma(layer,flag,mode,incidentPhi,Q,T)/10.;// in cm
479 err_m/=myRmidDGapCgem[layer];// in radian
480 J_dmdx(1,1)=-1*pos.y()/myR2midDGapCgem[layer];
481 J_dmdx(1,2)=pos.x()/myR2midDGapCgem[layer];
482 J_dmdx(1,3)=0.;
483 }
484 else if(flag==1) {
485 double V_trk = readoutPlane->getVFromPhiZ(phi_trk,z_trk*10,false)/10.0;// in cm
486 V_CC=(*iter_cluster)->getrecv()/10.;// in cm
487 //V_CC=(*iter_cluster)->getrecv_CC()/10.;// in cm
488 if(debug) cout<<"V_trk, V_rec = "<<V_trk<<", "<<V_CC<<endl;
489 delta_m=V_CC-V_trk;// in cm
490 if(fabs(delta_m)>5) {
491 /*int nextSheet = myNSheets[layer]-1-sheet;
492 CgemGeoReadoutPlane* nextReadoutPlane = myCgemGeomSvc->getReadoutPlane(layer,nextSheet);
493 double phiminNext = nextReadoutPlane->getPhimin();
494 double V_trk_next = nextReadoutPlane->getVInNextSheetFromV(V_trk*10, phiminNext)/10.0;// cm*10 -> mm
495 double delta_m_next = V_CC-V_trk_next;
496 if(fabs(delta_m_next)<delta_m) delta_m=delta_m_next;
497 if(debug) cout<<"V_trk_next = "<<V_trk_next<<endl;*/
498 double V_trk_nearPhiMin = readoutPlane->getVFromPhiZ_nearPhiMin(phi_trk,z_trk*10,false)/10.0;// in cm
499 double delta_m_2=V_CC-V_trk_nearPhiMin;// in cm
500 if(fabs(delta_m_2)<fabs(delta_m)) delta_m=delta_m_2;
501 if(debug) cout<<"V_trk_nearPhiMin= "<<V_trk_nearPhiMin<<endl;
502 }
503 err_m=myCgemCalibSvc->getSigma(layer,flag,mode,incidentPhi,Q,T)/10.;// in cm
504 J_dmdx(1,1)=-myRVCgem[layer]*cos(myAngStereoCgem[layer])*pos.y()/myR2midDGapCgem[layer];
505 J_dmdx(1,2)= myRVCgem[layer]*cos(myAngStereoCgem[layer])*pos.x()/myR2midDGapCgem[layer];
506 J_dmdx(1,3)= -sin(myAngStereoCgem[layer]);
507 }
508 else {
509 cout<<"flag ="<<flag<<", DotsHelixFitter::calculateNewHelix() is not ready for it!"<<endl;
510 continue;// next cluster
511 }
512 JT=J_dmdx*J_dxda;
513 J=JT.T();
514 double chi = delta_m/err_m;
515 myChiVecCgemCluster[i_digi]=chi;
516 if(debug)
517 cout<<"delta_m, err_m, chi = "<<delta_m<<", "<<err_m<<", "<<chi<<endl;
518 totChi2+=chi*chi;
519 for(int ipar=0; ipar<nPar; ipar++)
520 {
521 //if(debug)
522 // cout<<"d(doca)/d(a"<<"_"<<ipar<<") = "<<J(ipar+1,1)<<endl;
523 P(ipar+1,1) += J(ipar+1,1)*chi/err_m;
524 for(int ipar2=0; ipar2<=ipar; ipar2++)
525 {
526 U(ipar+1,ipar2+1)+=J(ipar+1,1)*J(ipar2+1,1)/(err_m*err_m);
527 }
528 }
529 }// loop myVecCgemCluster
530 if(debug)
531 cout<<"end of CGEM cluster loop"<<endl;
532 // <--- end of CGEM clusters loop ---
533
534 // --- delta a ---
535 //if(debug)
536 // cout<<"U="<<U<<endl;
537 int ifail=99;
538 U.invert(ifail);// also the error matrix
539 //if(debug)
540 //{
541 // cout<<"ifail = "<<ifail<<endl;
542 // cout<<"U^-1="<<U<<endl;
543 // cout<<"P="<<P<<endl;
544 //}
545 HepVector da = U*P;// in cm
546 //if(debug)
547 // cout<<"da = "<<da<<endl;
548
549 // --- update helix ---
550 for(int ipar=0; ipar<nPar; ipar++)
551 {
552 myHelix_a[ipar]+=da[ipar];
553 myHelix_aVec[ipar]=myHelix_a[ipar];
554 }
555 myHelix->a(myHelix_aVec);
556 if(debug)
557 {
558 cout<<"aNew = "<<myHelix_aVec
559 <<" chi2 = "<<totChi2
560 <<endl;
561 }
562 myChi2=totChi2;
563 //if(nIterations==0) cout<<"before fit chi2 = "<<totChi2<<endl;
564
565 // --- converge or too many iterations
566 nIterations++;
567 if(fabs(lastTotChi2-totChi2)<myDchi2Converge) {
568 if(debug)
569 cout<<"DotsHelixFitter::calculateNewHelix(): converge after "<<nIterations<<" iterations"
570 <<" with lastTotChi2="<<lastTotChi2<<", totChi2="<<totChi2
571 <<endl;
572 if(nPar==5) {
573 myHelix_Ea = U;
574 myHelix->Ea(U);
575 }
576 //cout<<"myHelix_Ea updated "<<endl;
577 break;
578 }
579 else if(totChi2-lastTotChi2>myDchi2Diverge) {
580 n_chi2_increase++;
581 if(n_chi2_increase>myMaxNChi2Increase||totChi2>myChi2Diverge)
582 return 5;// divergence
583 }
584 if(nIterations>myMaxIteration) {
585 if(debug) cout<<"DotsHelixFitter::calculateNewHelix(): "<<nIterations
586 <<" iterations, break with lastTotChi2, totChi2 = "<<lastTotChi2<<", "<<totChi2<<endl;
587 break;
588 }
589 lastTotChi2=totChi2;
590 }// <--- iterations
591
592 // --- delete vec_zini ---
593 delete []vec_zini;
594
595 if(nIterations>myMaxIteration) return 4;
596 else return 0;
597}
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
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
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::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition: CgemGeomSvc.h:140
DotsHelixFitter::dxda_cgem
HepMatrix dxda_cgem(KalmanFit::Helix a, double phi)
Definition: DotsHelixFitter.cxx:915
DotsHelixFitter::IntersectCylinder
double IntersectCylinder(double r)
Definition: DotsHelixFitter.cxx:899
ICgemCalibFunSvc::getSigma
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
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-14/src/helix/Helix.cxx:206
KalmanFit::Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/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-14/src/helix/Helix.cxx:151
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:265
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
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
TrkFitFun::getDeriLoc
bool getDeriLoc(int ipar, double helix[], double &deri, double wpos[], double zini)
Definition: TrkFitFun.cxx:298

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ calcuMeasDeriv()

void DotsHelixFitter::calcuMeasDeriv ( )

◆ deactiveHits() [1/2]

int DotsHelixFitter::deactiveHits ( double  chi_cut = 10,
int  nMax = 1 
)

Definition at line 600 of file DotsHelixFitter.cxx.

601{
602 map<double, int> map_abschi_idx;
603 //double maxChi2=0;
604 //vector<const MdcDigi*>::iterator iter_mdcDigi_maxChi2=NULL;
605 //int i_digi_maxChi2=-1;
606
607 int i_digi=0;
608 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
609 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
610 {
611 if(myMdcDigiIsActive[i_digi])
612 {
613 double abs_chi=fabs(myChiVecMdcDigi[i_digi]);
614
615 if( abs_chi>chi_cut ) // && maxChi2<(abs_chi*abs_chi)
616 {
617 //maxChi2=abs_chi*abs_chi;
618 //iter_mdcDigi_maxChi2=iter_mdcDigi;
619 //i_digi_maxChi2=i_digi;
620 map_abschi_idx[abs_chi]=i_digi;
621 }
622 }
623 }
624 int nDeactive=0;
625 if(nMax>0&&map_abschi_idx.size()>0)
626 {
627 map<double, int>::reverse_iterator it_map;
628 for(it_map=map_abschi_idx.rbegin(); it_map!=map_abschi_idx.rend(); it_map++)
629 {
630 myMdcDigiIsActive[it_map->second]=0;
631 //cout<<"a hit with |chi|="<<it_map->first<<" deactived"<<endl;
632 nDeactive++;
633 if(nDeactive>=nMax) break;
634 }
635 }
636 //if(i_digi_maxChi2!=-1)
637 //{
638 // myMdcDigiIsActive[i_digi_maxChi2] = 0;
639 // //cout<<"the worst hit "<<i_digi_maxChi2<<" (chi2="<<maxChi2<<") is deactive"<<endl;
640 // return 1;
641 //}
642 //else {
643 // //cout<<"no bad hits found!"<<endl;
644 // return 0;
645 //}
646 return nDeactive;
647}

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ deactiveHits() [2/2]

int DotsHelixFitter::deactiveHits ( int  layer_max = 43,
int  nHit_max = 1000 
)

Definition at line 672 of file DotsHelixFitter.cxx.

673{
674 int nActive=0;
675 int nDeactive=0;
676
677 vector<const MdcDigi*>::iterator iter_mdcDigi_begin = myVecMdcDigi.begin();
678 map<double, int>::iterator it_map;
679 //cout<<"digi size, flyLen size = "<<myVecMdcDigi.size()<<", "<<myMapFlylenIdx.size()<<endl;
680 for(it_map=myMapFlylenIdx.begin(); it_map!=myMapFlylenIdx.end(); it_map++)
681 {
682 //cout<<"nActive = "<<nActive<<endl;
683 vector<const MdcDigi*>::iterator iter_mdcDigi = iter_mdcDigi_begin+(it_map->second);
684 if(myMdcDigiIsActive[it_map->second])
685 {
686 if(nActive>=nHit_max) myMdcDigiIsActive[it_map->second]=0;
687 Identifier id = (*iter_mdcDigi)->identify();
688 int layer = MdcID::layer(id);
689 if(layer>=layer_max) myMdcDigiIsActive[it_map->second]=0;
690 if(myMdcDigiIsActive[it_map->second]) nActive++;
691 else {
692 nDeactive++;
693 //cout<<"a hit on layer "<<layer<<" is deactived"<<endl;
694 }
695 }
696 }
697 //cout<<nDeactive<<" hits deactived "<<endl;
698
699 return nDeactive;
700}

◆ dxda_cgem()

HepMatrix DotsHelixFitter::dxda_cgem ( KalmanFit::Helix  a,
double  phi 
)

Definition at line 915 of file DotsHelixFitter.cxx.

916{
917 HepMatrix dXDA(3,5,0);
918
919 double alpha = a.alpha();
920 double dr = a.dr();
921 double phi0 = a.phi0();
922 double kappa = a.kappa();
923 double dz = a.dz();
924 double tanl = a.tanl();
925
926 HepPoint3D pos = a.x(phi);
927 double x = pos.x();
928 double y = pos.y();
929 double z = pos.z();
930 double r2 = x*x+y*y;
931
932 double cosPhi=cos(phi);
933 double sinPhi=sin(phi);
934
935 double dPhiDdr = -(kappa/alpha-cosPhi/(dr+alpha/kappa))/sinPhi;
936 //double dPhiDdr2= -kappa*kappa*(2.0*alpha*dr+kappa*(dr*dr+r2))/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
937 //cout<<"dPhiDdr = "<<dPhiDdr<<endl;
938 //cout<<"dPhiDdr2= "<<dPhiDdr2<<endl;
939 double dPhiDkappa = -kappa*(dr*dr-r2)*(2*alpha+dr*kappa)/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
940
941 double dxDdr = cos(phi0)+alpha/kappa*sin(phi0+phi)*dPhiDdr;
942 double dyDdr = sin(phi0)-alpha/kappa*cos(phi0+phi)*dPhiDdr;
943 double dxDphi0 = -dr*sin(phi0)+alpha/kappa*(-sin(phi0)+sin(phi0+phi));
944 double dyDphi0 = dr*cos(phi0)+alpha/kappa*( cos(phi0)-cos(phi0+phi));
945 double dxDkappa = -alpha/(kappa*kappa)*(cos(phi0)-cos(phi0+phi))+alpha/kappa*sin(phi0+phi)*dPhiDkappa;
946 double dyDkappa = -alpha/(kappa*kappa)*(sin(phi0)-sin(phi0+phi))-alpha/kappa*cos(phi0+phi)*dPhiDkappa;
947 double dzDdr = -alpha/kappa*tanl*dPhiDdr;
948 double dzDkappa = alpha/(kappa*kappa)*tanl*phi-alpha/kappa*tanl*dPhiDkappa;
949 double dzDtanl = -alpha/kappa*phi;
950
951 dXDA(1,1) = dxDdr;
952 dXDA(1,2) = dxDphi0;
953 dXDA(1,3) = dxDkappa;
954 dXDA(2,1) = dyDdr;
955 dXDA(2,2) = dyDphi0;
956 dXDA(2,3) = dyDkappa;
957 dXDA(3,1) = dzDdr;
958 dXDA(3,3) = dzDkappa;
959 dXDA(3,4) = 1.0;
960 dXDA(3,5) = dzDtanl;
961
962 return dXDA;
963}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
alpha
const double alpha
Definition: FastVertexFit.cxx:4
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:229

Referenced by calculateNewHelix().

◆ fit()

void DotsHelixFitter::fit ( )

◆ fitCircleOnly()

void DotsHelixFitter::fitCircleOnly ( bool  x = true)
inline

Definition at line 37 of file DotsHelixFitter.h.

37{myFitCircleOnly=x;};

◆ fitModeCircle()

void DotsHelixFitter::fitModeCircle ( )
inline

Definition at line 38 of file DotsHelixFitter.h.

38{myFitCircleOnly=true;}

Referenced by HoughTrack::fitCircle().

◆ fitModeHelix()

void DotsHelixFitter::fitModeHelix ( )
inline

Definition at line 39 of file DotsHelixFitter.h.

39{myFitCircleOnly=false; myUseAxialHitsOnly=false;}

Referenced by HoughTrack::fitHelix().

◆ getChi2()

double DotsHelixFitter::getChi2 ( )
inline

Definition at line 45 of file DotsHelixFitter.h.

45{return myChi2;};

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ getClassHelix()

KalmanFit::Helix DotsHelixFitter::getClassHelix ( )
inline

Definition at line 42 of file DotsHelixFitter.h.

42{return *myHelix;} // in cm, radian, GeV/c

◆ getDocaFromTrk() [1/2]

double DotsHelixFitter::getDocaFromTrk ( )
inline

Definition at line 64 of file DotsHelixFitter.h.

64{return myDocaFromTrk;};

◆ getDocaFromTrk() [2/2]

double DotsHelixFitter::getDocaFromTrk ( const MdcDigi *  aDcDigi)
inline

Definition at line 65 of file DotsHelixFitter.h.

65{ calculateDocaFromTrk(aDcDigi); return myDocaFromTrk;};
DotsHelixFitter::calculateDocaFromTrk
void calculateDocaFromTrk(const MdcDigi *aDcDigi)
Definition: DotsHelixFitter.cxx:734

◆ getEa()

const HepSymMatrix & DotsHelixFitter::getEa ( )
inline

Definition at line 44 of file DotsHelixFitter.h.

44{return myHelix_Ea;};// error matrix

Referenced by HoughTrack::fitHelix().

◆ getFlightLength() [1/2]

double DotsHelixFitter::getFlightLength ( )
inline

Definition at line 61 of file DotsHelixFitter.h.

61{ return myFlightLength; };

◆ getFlightLength() [2/2]

double DotsHelixFitter::getFlightLength ( const MdcDigi *  aDcDigi)
inline

Definition at line 62 of file DotsHelixFitter.h.

62{ calculateDocaFromTrk(aDcDigi); return myFlightLength; };

◆ getHelix()

HepVector DotsHelixFitter::getHelix ( )
inline

Definition at line 43 of file DotsHelixFitter.h.

43{return myHelix_aVec;};

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ getLayer()

int DotsHelixFitter::getLayer ( const MdcDigi *  aDcDigi)
inline

Definition at line 63 of file DotsHelixFitter.h.

63{ loadOneDcDigi(aDcDigi); return myLayer; };

◆ getNActiveHits()

int DotsHelixFitter::getNActiveHits ( )
inline

Definition at line 55 of file DotsHelixFitter.h.

55{return myNActiveHits;};

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ getRmidGapCgem()

double DotsHelixFitter::getRmidGapCgem ( int  i)
inline

Definition at line 82 of file DotsHelixFitter.h.

83 {
84 if(i>2) i=2;
85 if(i<0) i=0;
86 return myRmidDGapCgem[i];
87 }

◆ getVecCgemCluster()

vector< const RecCgemCluster * > DotsHelixFitter::getVecCgemCluster ( )
inline

Definition at line 75 of file DotsHelixFitter.h.

75{ return myVecCgemCluster;}

◆ getVecChiCgemCluster()

vector< double > DotsHelixFitter::getVecChiCgemCluster ( )
inline

Definition at line 77 of file DotsHelixFitter.h.

77{ return myChiVecCgemCluster;}

◆ getVecChiMdcDigi()

vector< double > DotsHelixFitter::getVecChiMdcDigi ( )
inline

Definition at line 72 of file DotsHelixFitter.h.

72{return myChiVecMdcDigi;}

◆ getVecIsActiveCgemCluster()

vector< int > DotsHelixFitter::getVecIsActiveCgemCluster ( )
inline

Definition at line 76 of file DotsHelixFitter.h.

76{ return myCgemClusterIsActive;}

◆ getVecMdcDigi()

vector< const MdcDigi * > DotsHelixFitter::getVecMdcDigi ( )
inline

Definition at line 71 of file DotsHelixFitter.h.

71{return myVecMdcDigi;}

◆ getVecMdcDigiIsAct()

vector< int > DotsHelixFitter::getVecMdcDigiIsAct ( )
inline

Definition at line 73 of file DotsHelixFitter.h.

73{return myMdcDigiIsActive;}

◆ initialize()

void DotsHelixFitter::initialize ( )

Definition at line 34 of file DotsHelixFitter.cxx.

35{
36 if(myIsInitialized) return;
37
38 // --- get CgemGeomSvc ---
39 ISvcLocator* svcLocator = Gaudi::svcLocator();
40 ICgemGeomSvc* ISvc;
41 StatusCode Cgem_sc=svcLocator->service("CgemGeomSvc", ISvc);
42 if(Cgem_sc.isSuccess()) myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
43 else cout<<"DotsHelixFitter::initialize(): can not get CgemGeomSvc"<<endl;
44 int nLayerCgem = myCgemGeomSvc->getNumberOfCgemLayer();
45 if(nLayerCgem!=3) {
46 cout<<"DotsHelixFitter::initialize(): nLayerCgem!=3 !!!"<<endl;
47 return;
48 }
49 for(int i=0; i<nLayerCgem; i++)
50 {
51 //myCgemLayer[i]=myCgemGeomSvc->getCgemLayer(i);
52 CgemGeoLayer* CgemLayer = myCgemGeomSvc->getCgemLayer(i);
53 myRmidDGapCgem[i]=(CgemLayer->getMiddleROfGapD())/10.0;//cm
54 myR2midDGapCgem[i]=pow(myRmidDGapCgem[i],2);
55 myRVCgem[i]=(CgemLayer->getInnerROfAnodeCu1())/10.0;//cm
56 myRXCgem[i]=(CgemLayer->getInnerROfAnodeCu2())/10.0;//cm
57 myAngStereoCgem[i]=(CgemLayer->getAngleOfStereo())/180.0*CLHEP::pi;// degree -> radians
58 myNSheets[i]=CgemLayer->getNumberOfSheet();
59 //myRXstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu2();
60 //myRVstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu1();
61 //if(myPrintFlag)
62 cout<<"CGEM layer "<<i<<": "
63 <<" Rmid="<<myRmidDGapCgem[i]<<" cm "
64 <<" RV ="<<myRVCgem[i]<<" cm "
65 <<" RX ="<<myRXCgem[i]<<" cm "
66 <<" , stereo angle = "<<myAngStereoCgem[i]<<" radian "
67 <<myNSheets[i]<<" sheets"
68 //<<", is reversed "<<IsReverse
69 //<<", RX="<<myRXstrips[i]<<", RY="<<myRVstrips[i]
70 <<endl;
71 }
72
73 // --- get MDCGeomSvc ---
74 IMdcGeomSvc* ISvc2;
75 StatusCode mdc_sc=svcLocator->service("MdcGeomSvc", ISvc2);
76 if(mdc_sc.isSuccess()) myMdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc2);
77 else cout<<"DotsHelixFitter::initialize(): can not get MdcGeomSvc"<<endl;
78 int nLayer= myMdcGeomSvc->getLayerSize();
79 if(nLayer!=43) cout<<"DotsHelixFitter::initialize(): MDC nLayer = "<<nLayer<<" !=43 !!!"<<endl;
80 for(int i=0; i<nLayer; i++)
81 {
82 const MdcGeoLayer* aLayer = myMdcGeomSvc->Layer(i);
83 myRWires[i]=0.1*(aLayer->Radius());// cm
84 cout<<"MDC layer "<<i<<" R = "<<myRWires[i]<<endl;
85 double nomShift = myMdcGeomSvc->Layer(i)->nomShift();
86 if(nomShift>0) myWireFlag[i]=1;
87 else if(nomShift<0) myWireFlag[i]=-1;
88 else myWireFlag[i]=0;
89 int nWires=aLayer->NCell();
90 myNWire[i]=nWires;
91 for(int j=0; j<nWires; j++)
92 {
93 const MdcGeoWire* aWire = myMdcGeomSvc->Wire(i,j);
94 HepPoint3D aPointEast = 0.1 * (aWire->Forward());// in cm
95 double* aPointArray = new double[3]{aPointEast.x(),aPointEast.y(),aPointEast.z()};// in cm
96 myWestPosWires[i].push_back(aPointArray);
97 HepPoint3D aPointWest = 0.1*(aWire->Backward());// in cm
98 aPointArray = new double[3]{aPointWest.x(),aPointWest.y(),aPointWest.z()};// in cm
99 myEastPosWires[i].push_back(aPointArray);
100 HepPoint3D aPointMiddle = 0.5*(aPointEast+aPointWest);
101 aPointArray = new double[3]{aPointMiddle.x(),aPointMiddle.y(),aPointMiddle.z()};// in cm
102 myPosWires[i].push_back(aPointArray);
103 //myPhiWires[i].push_back(aWire->Forward().phi());
104 myPhiWires[i].push_back(aPointMiddle.phi());
105 myTensionWires[i].push_back(aWire->Tension());;
106 }
107 }
108
109 // --- MdcCalibFunSvc ---
110 IMdcCalibFunSvc* imdcCalibSvc;
111 StatusCode sc = svcLocator->service("MdcCalibFunSvc", imdcCalibSvc);
112 if ( sc.isSuccess() ){
113 myMdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
114 }
115 else {
116 cout<<"DotsHelixFitter::initialize(): can not get MdcCalibFunSvc"<<endl;
117 }
118
119 // --- get CgemCalibFunSvc ---
120 sc = svcLocator->service("CgemCalibFunSvc", myCgemCalibSvc);
121 if ( !(sc.isSuccess()) )
122 {
123 cout<<"DotsHelixFitter::initialize(): can not get CgemCalibFunSvc"<<endl;
124 }
125
126 myIsInitialized=true;
127
128 // --- get MdcUtilitySvc
129 sc = svcLocator->service("MdcUtilitySvc", myMdcUtilitySvc);
130 if ( !(sc.isSuccess()) )
131 {
132 cout<<"DotsHelixFitter::initialize(): can not get MdcUtilitySvc"<<endl;
133 }
134
135 return;
136}
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
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition: CgemGeomSvc.h:48
CgemGeomSvc::getNumberOfCgemLayer
int getNumberOfCgemLayer() const
Definition: CgemGeomSvc.h:45
MdcGeoLayer::Radius
double Radius(void) const
Definition: MdcGeoLayer.h:160
MdcGeoLayer::nomShift
double nomShift(void) const
Definition: MdcGeoLayer.h:169
MdcGeoLayer::NCell
int NCell(void) const
Definition: MdcGeoLayer.h:165
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
MdcGeoWire::Tension
double Tension(void) const
Definition: MdcGeoWire.h:136
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
MdcGeomSvc::getLayerSize
const int getLayerSize()
Definition: MdcGeomSvc.cxx:676

Referenced by DotsConnection::initialize(), and HoughFinder::initialize().

◆ IntersectCylinder()

double DotsHelixFitter::IntersectCylinder ( double  r)

Definition at line 899 of file DotsHelixFitter.cxx.

900{
901 double m_rad = myHelix->radius();
902 double l = myHelix->center().perp();
903
904 double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
905
906 if(cosPhi < -1 || cosPhi > 1) return 0;
907
908 double dPhi = myHelix->center().phi() - acos(cosPhi) - myHelix->phi0();
909
910 if(dPhi < -M_PI) dPhi += 2 * M_PI;
911
912 return dPhi;
913}
KalmanFit::Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:205

Referenced by calculateNewHelix().

◆ loadOneDcDigi()

void DotsHelixFitter::loadOneDcDigi ( const MdcDigi *  aDcDigi)

Definition at line 713 of file DotsHelixFitter.cxx.

714{
715 // --- get id, layer, wire ---
716 Identifier id = aDcDigi->identify();
717 myLayer = MdcID::layer(id);
718 myWire = MdcID::wire(id);
719 myCharge= aDcDigi->getChargeChannel();
720
721 // --- get myWirePos ---
722 double tension = 9999.;
723 //tension=myTensionWires[myLayer][myWire];
724 myWirePos[0]=myEastPosWires[myLayer][myWire][0];
725 myWirePos[1]=myEastPosWires[myLayer][myWire][1];
726 myWirePos[2]=myEastPosWires[myLayer][myWire][2];
727 myWirePos[3]=myWestPosWires[myLayer][myWire][0];
728 myWirePos[4]=myWestPosWires[myLayer][myWire][1];
729 myWirePos[5]=myWestPosWires[myLayer][myWire][2];
730 myWirePos[6]=tension;
731
732}
RawData::identify
virtual Identifier identify() const
Definition: RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition: RawData.cxx:45

Referenced by calculateDocaFromTrk(), getLayer(), and updateDcDigiInfo().

◆ makeRecMdcHit() [1/2]

RecMdcHit DotsHelixFitter::makeRecMdcHit ( const MdcDigi *  aDcDigi)

Definition at line 850 of file DotsHelixFitter.cxx.

851{
852 updateDcDigiInfo(aDcDigi);
853 RecMdcHit aRecMdcHit; // = new RecMdcHit;
854 aRecMdcHit.setDriftDistLeft(-1*myDriftDist[0]);
855 aRecMdcHit.setDriftDistRight(myDriftDist[1]);
856 aRecMdcHit.setErrDriftDistLeft(myDriftDistErr[0]);
857 aRecMdcHit.setErrDriftDistRight(myDriftDistErr[1]);
858 aRecMdcHit.setChisqAdd(myDcChi*myDcChi);
859 aRecMdcHit.setFlagLR(myLeftRight);
860 //aRecMdcHit.setStat(stat);
861 Identifier mdcId = aDcDigi->identify();
862 aRecMdcHit.setMdcId(mdcId);
863 double tdc = aDcDigi->getTimeChannel();
864 double adc = aDcDigi->getChargeChannel();
865 aRecMdcHit.setTdc(tdc);
866 aRecMdcHit.setAdc(adc);
867 aRecMdcHit.setDriftT(myDriftTime);
868 // --- doca
869 double doca=fabs(myDocaFromTrk);
870 if(myLeftRight==0) doca*=-1;
871 aRecMdcHit.setDoca(doca);
872 aRecMdcHit.setEntra(myEntranceAngle);
873 aRecMdcHit.setZhit(myPosOnWire[2]);
874 aRecMdcHit.setFltLen(myFlightLength);
875 return aRecMdcHit;
876}
DotsHelixFitter::updateDcDigiInfo
void updateDcDigiInfo(const MdcDigi *aDcDigi)
Definition: DotsHelixFitter.cxx:760
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

Referenced by makeRecMdcHit(), and makeRecMdcHitVec().

◆ makeRecMdcHit() [2/2]

RecMdcHit DotsHelixFitter::makeRecMdcHit ( const MdcDigi *  aDcDigi,
KalmanFit::Helix  aHelix 
)

Definition at line 878 of file DotsHelixFitter.cxx.

879{
880 setInitialHelix(aHelix);
881 return makeRecMdcHit(aDcDigi);
882}
DotsHelixFitter::setInitialHelix
void setInitialHelix(KalmanFit::Helix aHelix)
Definition: DotsHelixFitter.cxx:702
DotsHelixFitter::makeRecMdcHit
RecMdcHit makeRecMdcHit(const MdcDigi *aDcDigi)
Definition: DotsHelixFitter.cxx:850

◆ makeRecMdcHitVec()

vector< RecMdcHit > DotsHelixFitter::makeRecMdcHitVec ( int  sel = 1)

Definition at line 885 of file DotsHelixFitter.cxx.

886{
887 vector<RecMdcHit> aRecMdcHitVec;
888 int i_digi=0;
889 vector<const MdcDigi*>::iterator iter_mdcDigi = myVecMdcDigi.begin();
890 for(; iter_mdcDigi!=myVecMdcDigi.end(); iter_mdcDigi++, i_digi++)
891 {
892 if(sel==1&&!myMdcDigiIsActive[i_digi]) continue;// skip inactive hits
893 aRecMdcHitVec.push_back(makeRecMdcHit(*iter_mdcDigi));
894 }
895 return aRecMdcHitVec;
896}

◆ setCgemClusters()

void DotsHelixFitter::setCgemClusters ( vector< const RecCgemCluster * >  aVecCgemCluster)

Definition at line 159 of file DotsHelixFitter.cxx.

160{
161 myVecCgemCluster=aVecCgemCluster;
162 int nCluster1D = myVecCgemCluster.size();
163 myChiVecCgemCluster.resize(nCluster1D);
164 myCgemClusterIsActive.resize(nCluster1D);
165 //cout<<"set a vector<RecCgemCluster*> with "<<myVecCgemCluster.size()<<" clusters"<<endl;
166 myVecCgemClusterX.clear();
167 myVecCgemClusterV.clear();
168 vector<const RecCgemCluster*>::iterator iter_cluster = myVecCgemCluster.begin();
169 for(int i=0; iter_cluster!=myVecCgemCluster.end(); iter_cluster++, i++)
170 {
171 int flag = (*iter_cluster)->getflag();
172 myChiVecCgemCluster[i]=9999;
173 myCgemClusterIsActive[i]=1;
174 }
175}

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ setChi2Diverge()

void DotsHelixFitter::setChi2Diverge ( double  cut = 1000000)
inline

Definition at line 30 of file DotsHelixFitter.h.

30{myChi2Diverge=cut;};
cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition: KarFin.h:27

Referenced by DotsConnection::initialize().

◆ setDchi2Converge()

void DotsHelixFitter::setDchi2Converge ( double  cut = 5)
inline

Definition at line 27 of file DotsHelixFitter.h.

27{myDchi2Converge=cut;};

◆ setDchi2Diverge()

void DotsHelixFitter::setDchi2Diverge ( double  cut = 50)
inline

Definition at line 28 of file DotsHelixFitter.h.

28{myDchi2Diverge=cut;};

◆ setDChits()

void DotsHelixFitter::setDChits ( vector< const MdcDigi * >  aVecMdcDigi,
double  T0 
)

Definition at line 139 of file DotsHelixFitter.cxx.

140{
141 myVecMdcDigi=aVecMdcDigi;
142 //cout<<"set a vector<const MdcDigi*> with "<<myVecMdcDigi.size()<<" hits"<<endl;
143 int nDigi = myVecMdcDigi.size();
144 myChiVecMdcDigi.resize(nDigi);
145 myAmbiguityMdcDigi.resize(nDigi);
146 myMdcDigiIsActive.resize(nDigi);
147 myEventT0 = T0;
148
149 for(int i=0; i<43; i++) myNumMdcDigiPerLayer[i]=0;
150
151 for(int i=0; i<nDigi; i++)
152 {
153 myChiVecMdcDigi[i]=9999;
154 myAmbiguityMdcDigi[i]=0;
155 myMdcDigiIsActive[i]=1;
156 }
157}

Referenced by HoughTrack::fitCircle(), and HoughTrack::fitHelix().

◆ setInitialHelix()

void DotsHelixFitter::setInitialHelix ( KalmanFit::Helix  aHelix)

Definition at line 702 of file DotsHelixFitter.cxx.

703{
704 HepPoint3D pivot(0,0,0);
705 myHelix_aVec = aHelix.a();
706 if(myHelix!=NULL) delete myHelix;
707 myHelix = new KalmanFit::Helix(pivot, myHelix_aVec);
708 //cout<<"alpha="<<myHelix->alpha()<<endl;
709 for(int i=0; i<5; i++) myHelix_a[i]=myHelix_aVec[i];
710
711}

Referenced by HoughTrack::fitCircle(), HoughTrack::fitHelix(), and makeRecMdcHit().

◆ setMaxIterations()

void DotsHelixFitter::setMaxIterations ( int  n = 10)
inline

Definition at line 23 of file DotsHelixFitter.h.

23{myMaxIteration=n;};
n
const Int_t n
Definition: DataBase/tau_mode.c:65

◆ setMaxNChi2Increase()

void DotsHelixFitter::setMaxNChi2Increase ( int  n = 2)
inline

Definition at line 29 of file DotsHelixFitter.h.

29{myMaxNChi2Increase=n;};

◆ setMinHitsInHelixFit()

void DotsHelixFitter::setMinHitsInHelixFit ( int  n = 5)
inline

Definition at line 26 of file DotsHelixFitter.h.

26{myMinHitsInHelixFit=n;};

◆ setMinVhitsInHelixFit()

void DotsHelixFitter::setMinVhitsInHelixFit ( int  n = 2)
inline

Definition at line 25 of file DotsHelixFitter.h.

25{myMinVhitsInHelixFit=n;};

◆ setMinXhitsInCircleFit()

void DotsHelixFitter::setMinXhitsInCircleFit ( int  n = 3)
inline

Definition at line 24 of file DotsHelixFitter.h.

24{myMinXhitsInCircleFit=n;};

◆ setT0()

void DotsHelixFitter::setT0 ( double  T0)
inline

Definition at line 33 of file DotsHelixFitter.h.

33{myEventT0=T0;};

◆ updateChi2()

void DotsHelixFitter::updateChi2 ( )

◆ updateDcDigiInfo()

void DotsHelixFitter::updateDcDigiInfo ( const MdcDigi *  aDcDigi)

Definition at line 760 of file DotsHelixFitter.cxx.

761{
762 bool debug = false;
763
764 loadOneDcDigi(aDcDigi);
765
766 // --- get doca from track parameters ---
767 getDoca(myHelix_a, myWirePos, myDocaFromTrk, myPosOnWire, myPosOnWire[2]);
768 double phiWire = atan2(myPosOnWire[1],myPosOnWire[0]);
769 if(debug)
770 {
771 cout<<"DotsHelixFitter::UpdateDcDigiInfo(): "<<endl;
772 cout<<"doca = "<<myDocaFromTrk<<endl;
773 cout<<"point on wire: "<<myPosOnWire[0]<<", "<<myPosOnWire[1]<<", "<<myPosOnWire[2]<<endl;
774 }
775
776 // --- flight length ---
777 KalmanFit::Helix aHelix = *myHelix;
778 HepPoint3D aNewPivot(myPosOnWire[0],myPosOnWire[1],myPosOnWire[2]);
779 aHelix.pivot(aNewPivot);
780 double newPhi0 = aHelix.phi0();
781 double dphi = newPhi0-myHelix->phi0();
782 while(dphi<-M_PI) dphi+=2*M_PI;
783 while(dphi> M_PI) dphi-=2*M_PI;
784 myFlightLength = fabs(dphi*myHelix->radius())*sqrt(1+myHelix->tanl()*myHelix->tanl());// in cm
785
786 // --- approaching point on Helix
787 HepPoint3D posOnTrk = aHelix.x();
788 myPosOnTrk[0]=posOnTrk.x();
789 myPosOnTrk[1]=posOnTrk.y();
790 myPosOnTrk[2]=posOnTrk.z();
791 double phiPosOnTrk=atan2(myPosOnTrk[1],myPosOnTrk[0]);
792
793 // --- left or right
794 myLeftRight = 2;
795 int signDoca=1;
796 if(myDocaFromTrk!=0)
797 {
798 myLeftRight=int(myDocaFromTrk/fabs(myDocaFromTrk));
799 signDoca=myLeftRight;
800
801 // ---> conversion of left-right into the BESIII convention
802 //if(myLeftRight==-1) myLeftRight=0;
803 if(myLeftRight==1) myLeftRight=0;// fixed 2020-11-26
804 else myLeftRight=1;
805 double dphiLR=phiWire-phiPosOnTrk;
806 while(dphiLR<-M_PI) dphiLR+=2*M_PI;
807 while(dphiLR> M_PI) dphiLR-=2*M_PI;
808 //if(dphiLR>0) myLeftRight=0;
809 //else if(dphiLR<0) myLeftRight=1;
810 //cout<<"myLeftRight, dphiLR = "<<myLeftRight<<", "<<dphiLR<<endl;
811 }
812 if(debug) cout<<"myLeftRight = "<<myLeftRight<<endl;
813
814 // --- time of flight (TOF) ---
815 HepLorentzVector p4_pi = myHelix->momentum(dphi, 0.13957);
816 double speed = p4_pi.beta()*CC;// cm/second
817 double TOF = myFlightLength/speed*1.e9;// in ns
818
819 // --- get measured doca --- tof in ns, driftTime in ns, T0Walk in ns
820 double rawTime = RawDataUtil::MdcTime(aDcDigi->getTimeChannel());
821 double tprop = myMdcCalibFunSvc->getTprop(myLayer, myPosOnWire[2]);
822 double T0Walk = myMdcCalibFunSvc->getT0(myLayer,myWire) + myMdcCalibFunSvc->getTimeWalk(myLayer, myCharge);
823
824 // --- drift time ---
825 myDriftTime = rawTime - myEventT0 - TOF - T0Walk - tprop;
826 if(debug) cout<<"driftT = "<<myDriftTime<<endl;
827 // --- entrance Angle ---
828 double phiP = p4_pi.phi();
829 double entranceAngle = phiP-phiWire;
830 while(entranceAngle<-M_PI) entranceAngle+=2*M_PI;
831 while(entranceAngle> M_PI) entranceAngle-=2*M_PI;
832 myEntranceAngle = entranceAngle;
833 // --- measured drift distance ---
834 myDocaFromDigi = 0.1 * myMdcCalibFunSvc->driftTimeToDist(myDriftTime,myLayer,myWire,myLeftRight,entranceAngle);// in cm
835 myDriftDist[myLeftRight]=myDocaFromDigi;
836 myDriftDist[1-myLeftRight]=0.1 * myMdcCalibFunSvc->driftTimeToDist(myDriftTime,myLayer,myWire,1-myLeftRight,entranceAngle);// in cm
837 // --- get measurement error ---
838 double docaErr_hit = 0.1 * myMdcCalibFunSvc->getSigma(myLayer, myLeftRight, myDocaFromDigi, myEntranceAngle, myHelix_a[4], myPosOnWire[2], myCharge);// in cm
839 //if(debug) cout<<"error_doca_hit = "<<docaErr_hit<<endl;
840 myDriftDistErr[myLeftRight]=docaErr_hit;
841 myDriftDistErr[1-myLeftRight] = 0.1 * myMdcCalibFunSvc->getSigma(myLayer, 1-myLeftRight, myDocaFromDigi, myEntranceAngle, myHelix_a[4], myPosOnWire[2], myCharge);// in cm
842
843 // --- get chi ---
844 myDocaFromDigi*=signDoca;
845 if(debug) cout<<"doca_hit = "<<myDocaFromDigi<<endl;
846 double delD = myDocaFromDigi-myDocaFromTrk;
847 myDcChi = delD/docaErr_hit;
848}

Referenced by makeRecMdcHit().

◆ useAxialHitsOnly()

void DotsHelixFitter::useAxialHitsOnly ( bool  x = true)
inline

Definition at line 40 of file DotsHelixFitter.h.

40{myUseAxialHitsOnly=x;}
The documentation for this class was generated from the following files: