1 #include <CLHEP/Matrix/Matrix.h>
2 #include <CLHEP/Matrix/Vector.h>
4 #define SPEED_OF_LIGHT 29.9792548
15 const double c = 29.9792548;
20 vector<const DCDCTrackHit*> *trackHits,
24 n_fdc(pseudopoints->
size()), n_cdc(trackHits->
size()), ppPtr(pseudopoints),
25 trkhitPtr(trackHits), trajPtr(trajectory), delta(trajPtr->getDelta()),
26 debug_level(level), lorentz_def(lorentz_def_in), storeDetails(false),
27 innerResidFrac(1.0), rCDC(trackHits, trajectory, level),
28 rFDC(pseudopoints, trajectory, lorentz_def_in, level)
38 cout <<
"combinedResidFunc::resid: resid called\n";
39 cout <<
" params: " << *
x;
45 double thisChiSquared = 0.0;
51 vector<double> residsF;
54 vector<HepVector> pointF;
55 vector<double> docasF, errorsF;
56 vector<HepLorentzVector> pocasF;
58 for (
unsigned int ir = 0; ir <
n_fdc; ir++) {
59 thisResid = residsF[ir];
60 (*f)(ir + 1) = thisResid;
62 thisChiSquared += thisResid*thisResid;
64 FDCHitDetailsPtr->
doca = docasF[ir];
65 FDCHitDetailsPtr->
poca = pocasF[ir];
66 FDCHitDetailsPtr->
rCorr = pointF[ir];
75 vector<double> residsC;
78 vector<double> docasC, distsC, errorsC;
79 vector<HepLorentzVector> pocasC;
80 vector<HepVector> posWiresC;
82 for (
unsigned int ir = 0; ir <
n_cdc; ir++) {
83 thisResid = residsC[ir];
84 (*f)(n_fdc + ir + 1) = thisResid;
86 thisChiSquared += thisResid*thisResid;
88 CDCHitDetailsPtr->
doca = docasC[ir];
89 CDCHitDetailsPtr->
poca = pocasC[ir];
90 CDCHitDetailsPtr->
dist = distsC[ir];
91 CDCHitDetailsPtr->
posWire = posWiresC[ir];
96 if (
debug_level > 2) cout <<
"combinedResidFunc::resid: resids:" << *
f;
113 double delta_x = 0.0, delta_y = 0.0;
118 point(1) = ppoint.
xy.X() + delta_x;
119 point(2) = ppoint.
xy.Y() + delta_y;
121 point(1) = ppoint.
xy.X() - delta_x;
122 point(2) = ppoint.
xy.Y() - delta_y;
126 cout <<
"combinedResidFunc::pseudo2HepVector, x = " << x <<
" y = " << y <<
" z = " << z <<
" ct = " << ct <<
" delta_x = " << delta_x <<
" delta_y = " << delta_y <<
" ispos = " << ispos <<
" point = " << point << endl;
132 double xWire = ppoint.
wire->
udir(0);
133 double yWire = ppoint.
wire->
udir(1);
134 double wireCrossTraj = xWire*(y - ppoint.
xy.Y()) - yWire*(x - ppoint.
xy.X());
135 double wireCrossDelta = xWire*delta_y - yWire*delta_x;
136 bool isposTraj = wireCrossTraj > 0?
true:
false;
137 bool isposDelta = wireCrossDelta > 0?
true:
false;
138 bool ispos = !(isposTraj ^ isposDelta);
140 <<
"combinedResidFunc::getCorrectionSign,"
143 <<
" ppx = " << ppoint.
xy.X()
144 <<
" ppy = " << ppoint.
xy.Y()
145 <<
" dx = " << x - ppoint.
xy.X()
146 <<
" dy = " << y - ppoint.
xy.Y()
147 <<
" delta_x = " << delta_x
148 <<
" delta_y = " << delta_y
149 <<
" xWire = " << xWire
150 <<
" yWire = " << yWire
151 <<
" wireCrossTraj = " << wireCrossTraj
152 <<
" wireCrossDelta = " << wireCrossDelta
153 <<
" isposTraj = " << isposTraj
154 <<
" isposDelta = " << isposDelta
155 <<
" ispos = " << ispos
163 double cosangle = ppoint.
wire->
udir(1);
164 double sinangle= ppoint.
wire->
udir(0);
166 if (
debug_level > 3) cout <<
"x, y, z, ct = " << x <<
' ' << y <<
' ' << z <<
' ' << ct <<
" time, ct/c = " << ppoint.
time <<
" " << ct/
c <<
" driftdist = " << driftDist << endl;
168 double ds = -lorentzShift;
169 double dw = driftDist;
170 if (
debug_level > 3) cout <<
"ds, dw" << ds <<
" " << dw << endl;
171 delta_x = dw*cosangle + ds*sinangle;
172 delta_y = -dw*sinangle + ds*cosangle;
180 double z = wire->
origin.Z();
184 double theta = acos(wire->
udir.z());
185 double phi = atan2(wire->
udir.y(), wire->
udir.x());
190 DLine line(x, y, z, theta, phi);
197 details.
rCorr = point;
202 HepLorentzVector poca;
230 vector<double> residsC;
234 vector<double> residsF;
239 for (
unsigned int i = 0; i <
n_fdc; i++) {
240 residsBoth[i] = residsF[i];
243 for (
unsigned int i = 0; i <
n_cdc; i++) {
244 residsBoth[n_fdc + i] = residsC[i];
252 HepVector paramsCentral = *params, paramsThis;
255 int iHep = 0, jHep = 0;
257 for (
unsigned int j = 0; j < nParams; j++) {
259 cout <<
"central params: jHep = " << jHep <<
", values:" << paramsCentral(jHep) << endl;
265 vector<double> residsCentral, residsThis;
268 cout <<
"combinedResidFunc::deriv2: central resids, ";
269 for (
int k = 0; k < nResids; k++) {
270 cout << k <<
"=" << residsCentral[k] <<
" ";
276 for (
unsigned int j = 0; j < nParams; j++) {
279 paramsThis = paramsCentral;
280 paramsThis(jHep) = paramsCentral(jHep) +
delta[j];
281 if (
debug_level > 2) cout <<
"perturbed params: jHep = " << jHep <<
", values:" << paramsThis << endl;
287 cout <<
"combinedResidFunc::deriv2: resids, ";
288 for (
int k = 0; k < nResids; k++) {
289 cout << k <<
"=" << residsThis[k] <<
" ";
294 for (
int i = 0; i < nResids; i++) {
296 (*Jacobian)(iHep, jHep) = (residsThis[i] - residsCentral[i])/
delta[j];
int getXYT(double z, double &x, double &y, double &ct) const
DVector2 xy
rough x,y coordinates in lab coordinate system
vector< FDCHitDetails * > FDCDetails
const DFDCWire * wire
DFDCWire for this wire.
void getDetails(vector< HepVector > &points, vector< double > &docasRef, vector< double > &errorsRef, vector< HepLorentzVector > &pocasRef)
double doca(C &spaceObject, HepLorentzVector &poca) const
class DFDCPseudo: definition for a reconstructed point in the FDC
void setStoreDetails(bool storeDetailsValue)
bool getCorrectionSign(const DFDCPseudo &pseudopoint, double x, double y, double deltaX, double deltaY)
void getCorrectionValue(const DFDCPseudo &pseudopoint, double x, double y, double z, double t, double &delta_x, double &delta_y)
void getResids(vector< double > &residsRef)
HepVector poca(HepVector point)
void getResidsBoth(vector< double > &residsBoth)
void deriv(const HepVector *x, void *data, HepMatrix *J)
DLine trackhit2line(const DCDCTrackHit &trackhit)
vector< CDCHitDetails * > CDCDetails
void residAndDeriv(const HepVector *x, void *data, HepVector *f, HepMatrix *J)
virtual unsigned int getNumberOfParams()
double time
time corresponding to this pseudopoint.
void resid(const HepVector *x, void *data, HepVector *f)
double GetLorentzCorrection(double x, double y, double z, double alpha, double dx) const
void setInnerResidFrac(double innerResidFracIn)
virtual void swim(HepVector startingPoint, double theta, double phi)
void setInnerResidFrac(double innerResidFracIn)
void getDetails(vector< double > &docasRef, vector< double > &distsRef, vector< double > &errorsRef, vector< HepLorentzVector > &pocasRef, vector< HepVector > &posWiresRef)
HepVector pseudo2HepVector(const DFDCPseudo &pseudopoint)
FDCHitDetails getDetails(const DFDCPseudo *ppoint, HepVector point)
void getResids(vector< double > &residsRef)
const DLorentzDeflections * lorentz_def
combinedResidFunc(vector< const DFDCPseudo * > *pseudopoints, vector< const DCDCTrackHit * > *trackhits, MyTrajectory *trajectory, const DLorentzDeflections *lorentz_def, int level=1)