13#define HEP_SHORT_NAMES
16#include "CLHEP/Matrix/Vector.h"
17#include "CLHEP/Matrix/SymMatrix.h"
18#include "CLHEP/Matrix/DiagMatrix.h"
19#include "CLHEP/Matrix/Matrix.h"
22#include "TrkReco/TCosmicFitter.h"
23#include "TrkReco/TTrack.h"
25#include "TrkReco/TrkReco.h"
27#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
28#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
30#include "GaudiKernel/MsgStream.h"
31#include "GaudiKernel/AlgFactory.h"
32#include "GaudiKernel/ISvcLocator.h"
33#include "GaudiKernel/IMessageSvc.h"
34#include "GaudiKernel/IDataProviderSvc.h"
35#include "GaudiKernel/IDataManagerSvc.h"
36#include "GaudiKernel/SmartDataPtr.h"
37#include "GaudiKernel/IDataProviderSvc.h"
38#include "GaudiKernel/PropertyMgr.h"
39#include "GaudiKernel/IJobOptionsSvc.h"
40#include "GaudiKernel/IMessageSvc.h"
41#include "GaudiKernel/Bootstrap.h"
42#include "GaudiKernel/StatusCode.h"
44#include "GaudiKernel/ContainedObject.h"
45#include "GaudiKernel/SmartRef.h"
46#include "GaudiKernel/SmartRefVector.h"
47#include "GaudiKernel/ObjectVector.h"
48#include "EventModel/EventModel.h"
50#include "EvTimeEvent/RecEsTime.h"
52#include "CLHEP/Matrix/Matrix.h"
53#include "GaudiKernel/StatusCode.h"
54#include "GaudiKernel/IInterface.h"
55#include "GaudiKernel/Kernel.h"
56#include "GaudiKernel/Service.h"
57#include "GaudiKernel/ISvcLocator.h"
58#include "GaudiKernel/SvcFactory.h"
59#include "GaudiKernel/IDataProviderSvc.h"
60#include "GaudiKernel/Bootstrap.h"
61#include "GaudiKernel/MsgStream.h"
62#include "GaudiKernel/SmartDataPtr.h"
63#include "GaudiKernel/IMessageSvc.h"
66using CLHEP::HepVector;
67using CLHEP::HepSymMatrix;
68using CLHEP::HepDiagMatrix;
69using CLHEP::HepMatrix;
73#define CAL_TOF_HELIX 0
76#define Convergence 1.0e-5
86 StatusCode scmgn = Gaudi::svcLocator()->service(
"MagneticFieldSvc",m_pmgnIMF);
87 if(scmgn!=StatusCode::SUCCESS) {
88 std::cout<<
"Unable to open Magnetic field service"<<std::endl;
106 int err =
fit(b, 0.);
115#ifdef TRKRECO_DEBUG_DETAIL
116 std::cout <<
" TCosmicFitter::fit ..." << std::endl;
120 Gaudi::svcLocator() -> service(
"MdcCalibFunSvc", l_mdcCalFunSvc);
125 for (
unsigned i = 0; i < b.
links().length(); i++) {
126 unsigned state = b.
links()[i]->hit()->state();
138 double _t0_bes = -1.;
140 Gaudi::svcLocator()->service(
"MessageSvc",
msgSvc);
141 MsgStream log(
msgSvc,
"TCosmicFitter");
143 IDataProviderSvc* eventSvc = NULL;
144 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
146 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc,
"/Event/Recon/RecEsTimeCol");
148 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
149 _t0_bes = (*iter_evt)->getTest();
152 log << MSG::WARNING <<
"Could not find RecEsTimeCol" << endreq;
170 bool allAxial =
true;
177 for (
unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX);
185 for (
unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
190 const float Rad = 81;
191 float dRho =
t.helix().a()[0];
192 float Lproj = sqrt(Rad*Rad - dRho*dRho);
193 float tlmd =
t.helix().a()[4];
194 float fct = sqrt(1. + tlmd * tlmd);
199 while (
TMLink * l =
t.links()[i++]) {
211 int doSagCorrection = 0;
213 t.approach(* l, doSagCorrection);
214 double dPhi = l->dPhi();
215 const HepPoint3D & onTrack = l->positionOnTrack();
216 const HepPoint3D & onWire = l->positionOnWire();
218#ifdef TRKRECO_DEBUG_DETAIL
225 if (onWire.cross(onTrack).z() < 0.) leftRight =
WireHitLeft;
226 double distance = h.
drift(leftRight);
227 double eDistance = h.
dDrift(leftRight);
229 if(nTrial && !allAxial){
230 int side = leftRight;
232 double Tfly = _t0_bes/220.*(110.-onWire.y());
234 float Rad_wire = sqrt(
x[0]*
x[0]+
x[1]*
x[1]);
235 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho);
236 float flyLength = Lproj - L_wire;
237 if (
x[1]<0) flyLength = Lproj + L_wire;
238 Tfly = flyLength/29.98*sqrt(1.0+(0.105/0.5)*(0.105/0.5))*fct;
240 float time = l->hit()->reccdc()->tdc - Tfly;
244 float dist = distance;
245 float edist = eDistance;
246 double T0 = l_mdcCalFunSvc->
getT0(layerId,wire);
248 rawadc =l->hit()->reccdc()->adc;
250 double timeWalk = l_mdcCalFunSvc->
getTimeWalk(layerId, rawadc);
251 double drifttime =
time -T0-timeWalk;
252 dist = l_mdcCalFunSvc->
driftTimeToDist(drifttime,layerId, wire, side, 0.0);
253 edist = l_mdcCalFunSvc->
getSigma(layerId, side, dist, 0.0);
260 distance = (double) dist;
261 eDistance = (double) edist;
263 l->drift(distance,0);
264 l->drift(distance,1);
265 l->dDrift(eDistance,0);
266 l->dDrift(eDistance,1);
269 double eDistance2 = eDistance * eDistance;
273 double vmag =
v.mag();
274 double dDistance = vmag - distance;
277 this->dxda(*l,
t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection);
283 ? ((
v.x() * (1. - vw.x() * vw.x()) -
284 v.y() * vw.x() * vw.y() -
v.z() * vw.x() * vw.z())
286 (
v.y() * (1. - vw.y() * vw.y()) -
287 v.z() * vw.y() * vw.z() -
v.x() * vw.y() * vw.x())
289 (
v.z() * (1. - vw.z() * vw.z()) -
290 v.x() * vw.z() * vw.x() -
v.y() * vw.z() * vw.y())
294 std::cout <<
" in fit " << onTrack <<
", " << onWire;
297 dchi2da += (dDistance / eDistance2) * dDda;
298 d2chi2d2a += vT_times_v(dDda) / eDistance2;
299 double pChi2 = dDistance * dDistance / eDistance2;
303 l->update(onTrack, onWire, leftRight, pChi2);
307 double change = chi2Old - chi2;
308 if (fabs(change) < convergence)
break;
312#ifdef TRKRECO_DEBUG_DETAIL
313 std::cout <<
"chi2Old, chi2=" << chi2Old <<
" "<< chi2 << std::endl;
321 for (
unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
326 while (
TMLink * l =
t.links()[i++]) {
338 int doSagCorrection = 0;
340 t.approach(* l, doSagCorrection);
341 double dPhi = l->dPhi();
342 const HepPoint3D & onTrack = l->positionOnTrack();
343 const HepPoint3D & onWire = l->positionOnWire();
345#ifdef TRKRECO_DEBUG_DETAIL
352 if (onWire.cross(onTrack).z() < 0.) leftRight =
WireHitLeft;
353 double distance = h.
drift(leftRight);
354 double eDistance = h.
dDrift(leftRight);
355 if(nTrial && !allAxial){
356 int side = leftRight;
358 double Tfly = _t0_bes/220.*(110.-onWire.y());
360 float Rad_wire = sqrt(
x[0]*
x[0]+
x[1]*
x[1]);
361 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho);
362 float flyLength = Lproj - L_wire;
363 if (
x[1]<0) flyLength = Lproj + L_wire;
364 Tfly = flyLength/29.98*sqrt(1.0+(0.105/10)*(0.105/10))*fct;
367 float time = l->hit()->reccdc()->tdc - Tfly;
371 float dist= distance;
372 float edist = eDistance;
375 double T0 = l_mdcCalFunSvc->
getT0(layerId,wire);
377 rawadc =l->hit()->reccdc()->adc;
378 double timeWalk = l_mdcCalFunSvc->
getTimeWalk(layerId, rawadc);
379 double drifttime =
time -T0-timeWalk;
380 dist = l_mdcCalFunSvc->
driftTimeToDist(drifttime,layerId, wire, side, 0.0);
381 edist = l_mdcCalFunSvc->
getSigma(layerId, side, dist, 0.0);
388 distance = (double) dist;
389 eDistance = (double) edist;
391 l->drift(distance,0);
392 l->drift(distance,1);
393 l->dDrift(eDistance,0);
394 l->dDrift(eDistance,1);
397 double eDistance2 = eDistance * eDistance;
401 double vmag =
v.mag();
402 double dDistance = vmag - distance;
405 this->dxda(*l,
t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection);
411 ? ((
v.x() * (1. - vw.x() * vw.x()) -
412 v.y() * vw.x() * vw.y() -
v.z() * vw.x() * vw.z())
414 (
v.y() * (1. - vw.y() * vw.y()) -
415 v.z() * vw.y() * vw.z() -
v.x() * vw.y() * vw.x())
417 (
v.z() * (1. - vw.z() * vw.z()) -
418 v.x() * vw.z() * vw.x() -
v.y() * vw.z() * vw.y())
422 std::cout <<
" in fit " << onTrack <<
", " << onWire;
425 dchi2da += (dDistance / eDistance2) * dDda;
426 d2chi2d2a += vT_times_v(dDda) / eDistance2;
427 double pChi2 = dDistance * dDistance / eDistance2;
431 l->update(onTrack, onWire, leftRight, pChi2);
435#ifdef TRKRECO_DEBUG_DETAIL
436 std::cout <<
"factor = " << factor << std::endl;
437 std::cout <<
"chi2 = " << chi2 << std::endl;
439 if(factor < 0.01)
break;
446 f = dchi2da.sub(1, 3);
447 e = d2chi2d2a.sub(1, 3);
456 da = solve(d2chi2d2a, dchi2da);
459#ifdef TRKRECO_DEBUG_DETAIL
477 Eb = d2chi2d2a.sub(1, 3);
478 Ec = Eb.inverse(err);
488 Ea = d2chi2d2a.inverse(err);
500 t._ndf = nValid - dim;
953TCosmicFitter::dxda(
const TMLink & link,
959 int doSagCorrection)
const {
969 double rho = 333.564095/(-1000*(m_pmgnIMF->
getReferField())) / kappa;
970cout<<
"TCosmicFitter::rho: "<<333.564095/(-1000*(m_pmgnIMF->
getReferField()))<<endl;
971 double tanLambda = a[4];
973 double sinPhi0 =
sin(phi0);
974 double cosPhi0 =
cos(phi0);
975 double sinPhi0dPhi =
sin(phi0 + dPhi);
976 double cosPhi0dPhi =
cos(phi0 + dPhi);
980 HepPoint3D wireBackwardPosition =
w.backwardPosition();
983 if(doSagCorrection ){
984 cout<<
"doSagCorrection in TCosmicFitter!"<<endl;
986 float wirePosition[3];
991 if(wireZ >
w.backwardPosition().z() &&
992 wireZ <
w.forwardPosition().z() ){
997 xw.setX( (
double) wirePosition[0]);
998 xw.setY( (
double) wirePosition[1]);
999 xw.setZ( (
double) wirePosition[2]);
1001 wireBackwardPosition.setY((
double) ybSag);
1004 HepVector3D v_aux(
w.forwardPosition().x() -
w.backwardPosition().x(),
1006 w.forwardPosition().z() -
w.backwardPosition().z());
1015 double dmag2 = d.mag2();
1017 dphida[0] = (sinPhi0 * d.x() - cosPhi0 * d.y()) / dmag2;
1018 dphida[1] = (dRho + rho) * (cosPhi0 * d.x() + sinPhi0 * d.y())
1020 dphida[2] = (- rho / kappa) * (sinPhi0 * d.x() - cosPhi0 * d.y())
1031 c =
HepPoint3D(wireBackwardPosition - (
v * wireBackwardPosition) *
v);
1037 dxdphi[0] = rho * sinPhi0dPhi;
1038 dxdphi[1] = - rho * cosPhi0dPhi;
1039 dxdphi[2] = - rho * tanLambda;
1042 d2xdphi2[0] = rho * cosPhi0dPhi;
1043 d2xdphi2[1] = rho * sinPhi0dPhi;
1046 double dxdphi_dot_v = dxdphi[0]*
v.x() + dxdphi[1]*
v.y() + dxdphi[2]*
v.z();
1047 double x_dot_v =
x[0]*
v.x() +
x[1]*
v.y() +
x[2]*
v.z();
1049 double dfdphi = - (dxdphi[0] - dxdphi_dot_v*
v.x()) * dxdphi[0]
1050 - (dxdphi[1] - dxdphi_dot_v*
v.y()) * dxdphi[1]
1051 - (dxdphi[2] - dxdphi_dot_v*
v.z()) * dxdphi[2]
1052 - (x[0] - c[0] - x_dot_v*
v.x()) * d2xdphi2[0]
1053 - (x[1] - c[1] - x_dot_v*
v.y()) * d2xdphi2[1];
1059 dxda_phi[0] = cosPhi0;
1060 dxda_phi[1] = -(dRho + rho)*sinPhi0 + rho*sinPhi0dPhi;
1061 dxda_phi[2] = -(rho/kappa)*( cosPhi0 - cosPhi0dPhi );
1066 dyda_phi[0] = sinPhi0;
1067 dyda_phi[1] = (dRho + rho)*cosPhi0 - rho*cosPhi0dPhi;
1068 dyda_phi[2] = -(rho/kappa)*( sinPhi0 - sinPhi0dPhi );
1075 dzda_phi[2] = (rho/kappa)*tanLambda*dPhi;
1077 dzda_phi[4] = -rho*dPhi;
1081 d2xdphida[1] = rho*cosPhi0dPhi;
1082 d2xdphida[2] = -(rho/kappa)*sinPhi0dPhi;
1088 d2ydphida[1] = rho*sinPhi0dPhi;
1089 d2ydphida[2] = (rho/kappa)*cosPhi0dPhi;
1096 d2zdphida[2] = (rho/kappa)*tanLambda;
1098 d2zdphida[4] = -rho;
1101 for(
int i = 0; i < 5; i++ ){
1102 double d_dot_v =
v.x()*dxda_phi[i]
1104 +
v.z()*dzda_phi[i];
1105 dfda[i] = - (dxda_phi[i] - d_dot_v*
v.x()) * dxdphi[0]
1106 - (dyda_phi[i] - d_dot_v*
v.y()) * dxdphi[1]
1107 - (dzda_phi[i] - d_dot_v*
v.z()) * dxdphi[2]
1108 - (x[0] - c[0] - x_dot_v*
v.x()) * d2xdphida[i]
1109 - (x[1] - c[1] - x_dot_v*
v.y()) * d2ydphida[i]
1110 - (x[2] - c[2] - x_dot_v*
v.z()) * d2zdphida[i];
1111 dphida[i] = - dfda[i] /dfdphi;
1115 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
1116 dxda[1] = - (dRho + rho) * sinPhi0 + rho * sinPhi0dPhi * (1. + dphida[1]);
1117 dxda[2] = - rho / kappa * (cosPhi0 - cosPhi0dPhi)
1118 + rho * sinPhi0dPhi * dphida[2];
1119 dxda[3] = rho * sinPhi0dPhi * dphida[3];
1120 dxda[4] = rho * sinPhi0dPhi * dphida[4];
1122 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
1123 dyda[1] = (dRho + rho) * cosPhi0 - rho * cosPhi0dPhi * (1. + dphida[1]);
1124 dyda[2] = - rho / kappa * (sinPhi0 - sinPhi0dPhi)
1125 - rho * cosPhi0dPhi * dphida[2];
1126 dyda[3] = - rho * cosPhi0dPhi * dphida[3];
1127 dyda[4] = - rho * cosPhi0dPhi * dphida[4];
1129 dzda[0] = - rho * tanLambda * dphida[0];
1130 dzda[1] = - rho * tanLambda * dphida[1];
1131 dzda[2] = rho / kappa * tanLambda * dPhi - rho * tanLambda * dphida[2];
1132 dzda[3] = 1. - rho * tanLambda * dphida[3];
1133 dzda[4] = - rho * dPhi - rho * tanLambda * dphida[4];
HepGeom::Point3D< double > HepPoint3D
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.
A class to fit a TTrackBase object.
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
returns a list of masked TMLinks assigned to this track. 'mask' will be applied if mask is not 0.
virtual unsigned objectType(void) const
returns object type.
A class to represent a track in tracking.
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")