16#include "CLHEP/Matrix/DiagMatrix.h"
17#include "CLHEP/Matrix/Matrix.h"
18#include "CLHEP/Matrix/SymMatrix.h"
19#include "CLHEP/Matrix/Vector.h"
22#include "TrkReco/TCosmicFitter.h"
23#include "TrkReco/TTrack.h"
27#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
30#include "GaudiKernel/Bootstrap.h"
31#include "GaudiKernel/IDataManagerSvc.h"
32#include "GaudiKernel/IDataProviderSvc.h"
33#include "GaudiKernel/IMessageSvc.h"
34#include "GaudiKernel/ISvcLocator.h"
35#include "GaudiKernel/MsgStream.h"
36#include "GaudiKernel/PropertyMgr.h"
37#include "GaudiKernel/SmartDataPtr.h"
38#include "GaudiKernel/StatusCode.h"
40#include "EventModel/EventModel.h"
41#include "GaudiKernel/ContainedObject.h"
42#include "GaudiKernel/ObjectVector.h"
43#include "GaudiKernel/SmartRef.h"
44#include "GaudiKernel/SmartRefVector.h"
46#include "EvTimeEvent/RecEsTime.h"
48#include "CLHEP/Matrix/Matrix.h"
49#include "GaudiKernel/Bootstrap.h"
50#include "GaudiKernel/IDataProviderSvc.h"
51#include "GaudiKernel/IInterface.h"
52#include "GaudiKernel/IMessageSvc.h"
53#include "GaudiKernel/ISvcLocator.h"
54#include "GaudiKernel/Kernel.h"
55#include "GaudiKernel/MsgStream.h"
56#include "GaudiKernel/Service.h"
57#include "GaudiKernel/SmartDataPtr.h"
58#include "GaudiKernel/StatusCode.h"
60using CLHEP::HepDiagMatrix;
61using CLHEP::HepMatrix;
62using CLHEP::HepSymMatrix;
63using CLHEP::HepVector;
65typedef CLHEP::HepMatrix
Matrix;
67#define CAL_TOF_HELIX 0
70#define Convergence 1.0e-5
92 int err =
fit( b, 0. );
100#ifdef TRKRECO_DEBUG_DETAIL
101 std::cout <<
" TCosmicFitter::fit ..." << std::endl;
105 Gaudi::svcLocator()->service(
"MdcCalibFunSvc", l_mdcCalFunSvc );
110 for (
unsigned i = 0; i < b.
links().length(); i++ )
112 unsigned state = b.
links()[i]->hit()->state();
124 double _t0_bes = -1.;
126 Gaudi::svcLocator()->service(
"MessageSvc",
msgSvc );
127 MsgStream log(
msgSvc,
"TCosmicFitter" );
129 IDataProviderSvc* eventSvc = NULL;
130 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc );
132 SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc,
"/Event/Recon/RecEsTimeCol" );
135 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
136 _t0_bes = ( *iter_evt )->getTest();
139 else { log << MSG::WARNING <<
"Could not find RecEsTimeCol" << endmsg; }
156 bool allAxial =
true;
163 for (
unsigned i = 0; i < 5; i++ ) maxDouble[i] = ( FLT_MAX );
172 for (
unsigned j = 0; j < 5; j++ ) dchi2da[j] = 0.;
177 const float Rad = 81;
178 float dRho =
t.helix().a()[0];
179 float Lproj = sqrt( Rad * Rad - dRho * dRho );
180 float tlmd =
t.helix().a()[4];
181 float fct = sqrt( 1. + tlmd * tlmd );
186 while (
TMLink* l =
t.links()[i++] )
199 int doSagCorrection = 0;
201 t.approach( *l, doSagCorrection );
202 double dPhi = l->dPhi();
203 const HepPoint3D& onTrack = l->positionOnTrack();
204 const HepPoint3D& onWire = l->positionOnWire();
206#ifdef TRKRECO_DEBUG_DETAIL
213 if ( onWire.cross( onTrack ).z() < 0. ) leftRight =
WireHitLeft;
214 double distance = h.
drift( leftRight );
215 double eDistance = h.
dDrift( leftRight );
217 if ( nTrial && !allAxial )
219 int side = leftRight;
221 double Tfly = _t0_bes / 220. * ( 110. - onWire.y() );
223 float Rad_wire = sqrt( x[0] * x[0] + x[1] * x[1] );
224 float L_wire = sqrt( Rad_wire * Rad_wire - dRho * dRho );
225 float flyLength = Lproj - L_wire;
226 if ( x[1] < 0 ) flyLength = Lproj + L_wire;
227 Tfly = flyLength / 29.98 * sqrt( 1.0 + ( 0.105 / 0.5 ) * ( 0.105 / 0.5 ) ) *
230 float time = l->hit()->reccdc()->tdc - Tfly;
234 float dist = distance;
235 float edist = eDistance;
236 double T0 = l_mdcCalFunSvc->
getT0( layerId, wire );
238 rawadc = l->hit()->reccdc()->adc;
240 double timeWalk = l_mdcCalFunSvc->
getTimeWalk( layerId, rawadc );
241 double drifttime =
time - T0 - timeWalk;
242 dist = l_mdcCalFunSvc->
driftTimeToDist( drifttime, layerId, wire, side, 0.0 );
243 edist = l_mdcCalFunSvc->
getSigma( layerId, side, dist, 0.0 );
245 edist = edist / 10.0;
250 distance = (double)dist;
251 eDistance = (double)edist;
253 l->drift( distance, 0 );
254 l->drift( distance, 1 );
255 l->dDrift( eDistance, 0 );
256 l->dDrift( eDistance, 1 );
259 double eDistance2 = eDistance * eDistance;
263 double vmag =
v.mag();
264 double dDistance = vmag - distance;
267 this->dxda( *l,
t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection );
272 dDda = ( vmag > 0. ) ? ( (
v.x() * ( 1. - vw.x() * vw.x() ) -
v.y() * vw.x() * vw.y() -
273 v.z() * vw.x() * vw.z() ) *
275 (
v.y() * ( 1. - vw.y() * vw.y() ) -
v.z() * vw.y() * vw.z() -
276 v.x() * vw.y() * vw.x() ) *
278 (
v.z() * ( 1. - vw.z() * vw.z() ) -
v.x() * vw.z() * vw.x() -
279 v.y() * vw.z() * vw.y() ) *
285 std::cout <<
" in fit " << onTrack <<
", " << onWire;
288 dchi2da += ( dDistance / eDistance2 ) * dDda;
289 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
290 double pChi2 = dDistance * dDistance / eDistance2;
294 l->update( onTrack, onWire, leftRight, pChi2 );
298 double change = chi2Old - chi2;
299 if ( fabs( change ) < convergence )
break;
304#ifdef TRKRECO_DEBUG_DETAIL
305 std::cout <<
"chi2Old, chi2=" << chi2Old <<
" " << chi2 << std::endl;
313 for (
unsigned j = 0; j < 5; j++ ) dchi2da[j] = 0.;
318 while (
TMLink* l =
t.links()[i++] )
331 int doSagCorrection = 0;
333 t.approach( *l, doSagCorrection );
334 double dPhi = l->dPhi();
335 const HepPoint3D& onTrack = l->positionOnTrack();
336 const HepPoint3D& onWire = l->positionOnWire();
338#ifdef TRKRECO_DEBUG_DETAIL
345 if ( onWire.cross( onTrack ).z() < 0. ) leftRight =
WireHitLeft;
346 double distance = h.
drift( leftRight );
347 double eDistance = h.
dDrift( leftRight );
348 if ( nTrial && !allAxial )
350 int side = leftRight;
352 double Tfly = _t0_bes / 220. * ( 110. - onWire.y() );
354 float Rad_wire = sqrt( x[0] * x[0] + x[1] * x[1] );
355 float L_wire = sqrt( Rad_wire * Rad_wire - dRho * dRho );
356 float flyLength = Lproj - L_wire;
357 if ( x[1] < 0 ) flyLength = Lproj + L_wire;
358 Tfly = flyLength / 29.98 * sqrt( 1.0 + ( 0.105 / 10 ) * ( 0.105 / 10 ) ) *
362 float time = l->hit()->reccdc()->tdc - Tfly;
366 float dist = distance;
367 float edist = eDistance;
370 double T0 = l_mdcCalFunSvc->
getT0( layerId, wire );
372 rawadc = l->hit()->reccdc()->adc;
373 double timeWalk = l_mdcCalFunSvc->
getTimeWalk( layerId, rawadc );
374 double drifttime =
time - T0 - timeWalk;
375 dist = l_mdcCalFunSvc->
driftTimeToDist( drifttime, layerId, wire, side, 0.0 );
376 edist = l_mdcCalFunSvc->
getSigma( layerId, side, dist, 0.0 );
378 edist = edist / 10.0;
383 distance = (double)dist;
384 eDistance = (double)edist;
386 l->drift( distance, 0 );
387 l->drift( distance, 1 );
388 l->dDrift( eDistance, 0 );
389 l->dDrift( eDistance, 1 );
392 double eDistance2 = eDistance * eDistance;
396 double vmag =
v.mag();
397 double dDistance = vmag - distance;
400 this->dxda( *l,
t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection );
405 dDda = ( vmag > 0. ) ? ( (
v.x() * ( 1. - vw.x() * vw.x() ) -
v.y() * vw.x() * vw.y() -
406 v.z() * vw.x() * vw.z() ) *
408 (
v.y() * ( 1. - vw.y() * vw.y() ) -
v.z() * vw.y() * vw.z() -
409 v.x() * vw.y() * vw.x() ) *
411 (
v.z() * ( 1. - vw.z() * vw.z() ) -
v.x() * vw.z() * vw.x() -
412 v.y() * vw.z() * vw.y() ) *
418 std::cout <<
" in fit " << onTrack <<
", " << onWire;
421 dchi2da += ( dDistance / eDistance2 ) * dDda;
422 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
423 double pChi2 = dDistance * dDistance / eDistance2;
427 l->update( onTrack, onWire, leftRight, pChi2 );
431#ifdef TRKRECO_DEBUG_DETAIL
432 std::cout <<
"factor = " << factor << std::endl;
433 std::cout <<
"chi2 = " << chi2 << std::endl;
435 if ( factor < 0.01 )
break;
443 f = dchi2da.sub( 1, 3 );
444 e = d2chi2d2a.sub( 1, 3 );
452 else { da = solve( d2chi2d2a, dchi2da ); }
454#ifdef TRKRECO_DEBUG_DETAIL
473 Eb = d2chi2d2a.sub( 1, 3 );
474 Ec = Eb.inverse( err );
485 Ea = d2chi2d2a.inverse( err );
496 t._ndf = nValid - dim;
947int TCosmicFitter::dxda(
const TMLink& link,
const Helix& h,
double dPhi,
Vector& dxda,
948 Vector& dyda,
Vector& dzda,
int doSagCorrection )
const {
958 double rho = 333.564095 / ( -1000 * ( getPmgnIMF()->
getReferField() ) ) / kappa;
959 cout <<
"TCosmicFitter::rho: " << 333.564095 / ( -1000 * ( getPmgnIMF()->
getReferField() ) )
961 double tanLambda = a[4];
963 double sinPhi0 =
sin( phi0 );
964 double cosPhi0 =
cos( phi0 );
965 double sinPhi0dPhi =
sin( phi0 + dPhi );
966 double cosPhi0dPhi =
cos( phi0 + dPhi );
970 HepPoint3D wireBackwardPosition =
w.backwardPosition();
973 if ( doSagCorrection )
975 cout <<
"doSagCorrection in TCosmicFitter!" << endl;
977 float wirePosition[3];
978 wirePosition[0] = wirePosition[1] = wirePosition[2] = 0.;
983 if ( wireZ >
w.backwardPosition().z() && wireZ <
w.forwardPosition().z() )
989 xw.setX( (
double)wirePosition[0] );
990 xw.setY( (
double)wirePosition[1] );
991 xw.setZ( (
double)wirePosition[2] );
993 wireBackwardPosition.setY( (
double)ybSag );
996 HepVector3D v_aux(
w.forwardPosition().x() -
w.backwardPosition().x(), yfSag - ybSag,
997 w.forwardPosition().z() -
w.backwardPosition().z() );
1007 double dmag2 = d.mag2();
1009 dphida[0] = ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
1010 dphida[1] = ( dRho + rho ) * ( cosPhi0 * d.x() + sinPhi0 * d.y() ) / dmag2 - 1.;
1011 dphida[2] = ( -rho / kappa ) * ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
1022 c =
HepPoint3D( wireBackwardPosition - (
v * wireBackwardPosition ) *
v );
1028 dxdphi[0] = rho * sinPhi0dPhi;
1029 dxdphi[1] = -rho * cosPhi0dPhi;
1030 dxdphi[2] = -rho * tanLambda;
1033 d2xdphi2[0] = rho * cosPhi0dPhi;
1034 d2xdphi2[1] = rho * sinPhi0dPhi;
1037 double dxdphi_dot_v = dxdphi[0] *
v.x() + dxdphi[1] *
v.y() + dxdphi[2] *
v.z();
1038 double x_dot_v =
x[0] *
v.x() +
x[1] *
v.y() +
x[2] *
v.z();
1040 double dfdphi = -( dxdphi[0] - dxdphi_dot_v *
v.x() ) * dxdphi[0] -
1041 ( dxdphi[1] - dxdphi_dot_v *
v.y() ) * dxdphi[1] -
1042 ( dxdphi[2] - dxdphi_dot_v *
v.z() ) * dxdphi[2] -
1043 (
x[0] - c[0] - x_dot_v *
v.x() ) * d2xdphi2[0] -
1044 (
x[1] - c[1] - x_dot_v *
v.y() ) * d2xdphi2[1];
1049 dxda_phi[0] = cosPhi0;
1050 dxda_phi[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi;
1051 dxda_phi[2] = -( rho / kappa ) * ( cosPhi0 - cosPhi0dPhi );
1056 dyda_phi[0] = sinPhi0;
1057 dyda_phi[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi;
1058 dyda_phi[2] = -( rho / kappa ) * ( sinPhi0 - sinPhi0dPhi );
1065 dzda_phi[2] = ( rho / kappa ) * tanLambda * dPhi;
1067 dzda_phi[4] = -rho * dPhi;
1071 d2xdphida[1] = rho * cosPhi0dPhi;
1072 d2xdphida[2] = -( rho / kappa ) * sinPhi0dPhi;
1078 d2ydphida[1] = rho * sinPhi0dPhi;
1079 d2ydphida[2] = ( rho / kappa ) * cosPhi0dPhi;
1086 d2zdphida[2] = ( rho / kappa ) * tanLambda;
1088 d2zdphida[4] = -rho;
1091 for (
int i = 0; i < 5; i++ )
1093 double d_dot_v =
v.x() * dxda_phi[i] +
v.y() * dyda_phi[i] +
v.z() * dzda_phi[i];
1094 dfda[i] = -( dxda_phi[i] - d_dot_v *
v.x() ) * dxdphi[0] -
1095 ( dyda_phi[i] - d_dot_v *
v.y() ) * dxdphi[1] -
1096 ( dzda_phi[i] - d_dot_v *
v.z() ) * dxdphi[2] -
1097 (
x[0] - c[0] - x_dot_v *
v.x() ) * d2xdphida[i] -
1098 (
x[1] - c[1] - x_dot_v *
v.y() ) * d2ydphida[i] -
1099 (
x[2] - c[2] - x_dot_v *
v.z() ) * d2zdphida[i];
1100 dphida[i] = -dfda[i] / dfdphi;
1104 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
1105 dxda[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi * ( 1. + dphida[1] );
1106 dxda[2] = -rho / kappa * ( cosPhi0 - cosPhi0dPhi ) + rho * sinPhi0dPhi * dphida[2];
1107 dxda[3] = rho * sinPhi0dPhi * dphida[3];
1108 dxda[4] = rho * sinPhi0dPhi * dphida[4];
1110 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
1111 dyda[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi * ( 1. + dphida[1] );
1112 dyda[2] = -rho / kappa * ( sinPhi0 - sinPhi0dPhi ) - rho * cosPhi0dPhi * dphida[2];
1113 dyda[3] = -rho * cosPhi0dPhi * dphida[3];
1114 dyda[4] = -rho * cosPhi0dPhi * dphida[4];
1116 dzda[0] = -rho * tanLambda * dphida[0];
1117 dzda[1] = -rho * tanLambda * dphida[1];
1118 dzda[2] = rho / kappa * tanLambda * dPhi - rho * tanLambda * dphida[2];
1119 dzda[3] = 1. - rho * tanLambda * dphida[3];
1120 dzda[4] = -rho * dPhi - rho * tanLambda * dphida[4];
1126 if (
nullptr == m_pmgnIMF )
1128 StatusCode sc = Gaudi::svcLocator()->service(
"MagneticFieldSvc", m_pmgnIMF );
1129 if ( sc.isFailure() )
1130 { std::cout <<
"Unable to open Magnetic field service" << std::endl; }
HepGeom::Vector3D< double > HepVector3D
HepGeom::Point3D< double > HepPoint3D
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
double sin(const BesAngle a)
double cos(const BesAngle a)
CLHEP::HepSymMatrix SymMatrix
#define WireHitFittingValid
#define WireHitInvalidForFit
#define TFitAlreadyFitted
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
const HepVector & a(void) const
returns helix parameters.
virtual double getReferField()=0
virtual double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const =0
virtual 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 =0
virtual double getT0(int layid, int cellid) const =0
virtual double getTimeWalk(int layid, double Q) const =0
virtual ~TCosmicFitter()
Destructor.
TCosmicFitter(const std::string &name)
Constructor.
int fit(TTrackBase &) const
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
unsigned state(void) const
returns state.
float dDrift(unsigned) const
returns drift distance error.
float drift(unsigned) const
returns drift distance.
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
A class to represent a wire in MDC.
unsigned localId(void) const
returns local id in a wire layer.
unsigned layerId(void) const
returns layer id.
const HepVector3D & direction(void) const
returns direction vector of the wire.
bool stereo(void) const
returns true if this wire is in a stereo layer.
const std::string & name(void) const
returns name.
TMFitter(const std::string &name)
Constructor.
A class to relate TMDCWireHit and TTrack objects.
const HepPoint3D & positionOnTrack(void) const
returns the closest point on track to wire.
const TMDCWire *const wire(void) const
returns a pointer to a wire.
A virtual class for a track class in tracking.
bool fitted(void) const
returns true if fitted.
const AList< TMLink > & links(unsigned mask=0) const
virtual unsigned objectType(void) const
returns object type.
A class to represent a track in tracking.