40 #define COVEXC "cov_is_zero" 60 for (
unsigned int i = 0; i < nreps; ++i) {
76 std::vector<std::vector<int>*> planes;
78 int nPlanes = planes.size();
81 for (
unsigned int iBeta = 0; iBeta <
fBeta.size(); iBeta++) {
84 for (
int ipl = 0; ipl < nPlanes; ++ipl) {
86 std::vector<GFAbsRecoHit*>
hits;
87 unsigned int nhitsInPlane = planes.at(ipl)->size();
88 for (
unsigned int ihitInPl = 0; ihitInPl < nhitsInPlane; ++ihitInPl) {
89 hits.push_back(trk->
getHit(planes.at(ipl)->at(ihitInPl)));
93 for (
unsigned int irep = 0; irep < nreps; ++irep) {
98 if (bk->
hitFailed(planes.at(ipl)->at(0)) > 0)
continue;
101 int repDim = rep->
getDim();
102 TMatrixT<Double_t> state(repDim, 1);
103 TMatrixT<Double_t> cov(repDim, repDim);
106 if (ipl > 0 || (ipl == 0 && iBeta == 0)) {
111 pl = hits.at(0)->getDetPlane(rep);
115 if (cov[0][0] < 1.
E-50) {
121 std::cerr << exc.
what();
123 for (
unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
138 TMatrixT<Double_t> H = hits.at(0)->getHMatrix(rep);
139 bk->
setMatrix(
"fPredSt", planes.at(ipl)->at(0), state);
140 bk->
setMatrix(
"fPredCov", planes.at(ipl)->at(0), cov);
142 TMatrixT<Double_t> covInv;
145 TMatrixT<Double_t> Htransp(H);
148 bk->
setMatrix(
"H", planes.at(ipl)->at(0), H);
151 TMatrixT<Double_t> stMod(state.GetNrows(), 1);
154 for (
unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
156 bk->
getNumber(
"p", planes.at(ipl)->at(ihit), pki);
159 TMatrixT<Double_t> V = hits.at(0)->getHitCov(pl);
160 TMatrixT<Double_t> Vinv;
162 bk->
setMatrix(
"V", planes.at(ipl)->at(0), V);
163 for (
unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
166 TMatrixT<Double_t> m = hits.at(ihit)->getHitCoord(pl);
167 bk->
setMatrix(
"m", planes.at(ipl)->at(ihit), m);
170 bk->
getNumber(
"p", planes.at(ipl)->at(ihit), pki);
171 TMatrixT<Double_t> Gain =
calcGain(cov, V, H, sumPk);
173 stMod += pki * Gain * (m - H * state);
183 for (
int ipl = nPlanes - 1; ipl >= 0; --ipl) {
185 std::vector<GFAbsRecoHit*>
hits;
186 unsigned int nhitsInPlane = planes.at(ipl)->size();
187 for (
unsigned int ihitInPl = 0; ihitInPl < nhitsInPlane; ++ihitInPl) {
188 hits.push_back(trk->
getHit(planes.at(ipl)->at(ihitInPl)));
192 for (
unsigned int irep = 0; irep < nreps; ++irep) {
197 if (bk->
hitFailed(planes.at(ipl)->at(0)) > 0)
continue;
199 int repDim = rep->
getDim();
200 TMatrixT<Double_t> state(repDim, 1);
201 TMatrixT<Double_t> cov(repDim, repDim);
204 TMatrixT<Double_t> H;
205 bk->
getMatrix(
"H", planes.at(ipl)->at(0), H);
206 TMatrixT<Double_t> Htransp(H);
211 if (ipl < (nPlanes - 1)) {
215 if (cov[0][0] < 1.
E-50) {
222 std::cerr << exc.
what();
225 for (
unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
235 TMatrixT<Double_t> covInv;
238 bk->
setMatrix(
"bPredSt", planes.at(ipl)->at(0), state);
239 bk->
setMatrix(
"bPredCov", planes.at(ipl)->at(0), cov);
241 TMatrixT<Double_t> stMod(state.GetNrows(), 1);
244 for (
unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
246 bk->
getNumber(
"p", planes.at(ipl)->at(ihit), pki);
250 TMatrixT<Double_t> V;
251 bk->
getMatrix(
"V", planes.at(ipl)->at(0), V);
252 TMatrixT<Double_t> Vinv;
255 for (
unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
257 TMatrixT<Double_t> m;
258 bk->
getMatrix(
"m", planes.at(ipl)->at(ihit), m);
261 bk->
getNumber(
"p", planes.at(ipl)->at(ihit), pki);
262 TMatrixT<Double_t> Gain =
calcGain(cov, V, H, sumPk);
265 stMod += pki * Gain * (m - H * state);
274 TMatrixT<Double_t> fSt, fCov, bSt, bCov, smooSt, smooCov, fCovInv, bCovInv, m, H, smooResid;
275 for (
int ipl = 0; ipl < nPlanes; ++ipl) {
276 for (
unsigned int irep = 0; irep < nreps; ++irep) {
279 if (bk->
hitFailed(planes.at(ipl)->at(0)) > 0)
continue;
280 bk->
getMatrix(
"fPredSt", planes.at(ipl)->at(0), fSt);
281 bk->
getMatrix(
"bPredSt", planes.at(ipl)->at(0), bSt);
282 bk->
getMatrix(
"fPredCov", planes.at(ipl)->at(0), fCov);
283 bk->
getMatrix(
"bPredCov", planes.at(ipl)->at(0), bCov);
284 fCovInv.ResizeTo(fCov);
285 bCovInv.ResizeTo(bCov);
288 smooCov.ResizeTo(fCov);
290 smooSt.ResizeTo(fSt.GetNrows(), 1);
291 smooSt = smooCov * (fCovInv * fSt + bCovInv * bSt);
293 bk->
getMatrix(
"H", planes.at(ipl)->at(0), H);
294 for (
unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
296 bk->
getMatrix(
"m", planes.at(ipl)->at(ihit), m);
297 smooResid.ResizeTo(m.GetNrows(), 1);
298 smooResid = m - H * smooSt;
301 bk->
setMatrix(
"smooResid", planes.at(ipl)->at(ihit), smooResid);
307 if (iBeta !=
fBeta.size() - 1) {
308 for (
unsigned int irep = 0; irep < nreps; ++irep) {
310 for (
int ipl = 0; ipl < nPlanes; ++ipl) {
312 if (bk->
hitFailed(planes.at(ipl)->at(0)) > 0)
continue;
314 double sumPhiCut = 0.;
315 std::vector<double> phi;
316 TMatrixT<Double_t> V;
317 bk->
getMatrix(
"V", planes.at(ipl)->at(0), V);
318 V *=
fBeta.at(iBeta);
319 TMatrixT<Double_t> Vinv;
321 for (
unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
323 bk->
getMatrix(
"smooResid", planes.at(ipl)->at(ihit), smooResid);
324 TMatrixT<Double_t> smooResidTransp(smooResid);
326 int dimV = V.GetNrows();
327 double detV = V.Determinant();
328 TMatrixT<Double_t> expArg = smooResidTransp * Vinv * smooResid;
329 if (expArg.GetNrows() != 1)
330 throw std::runtime_error(
"genf::GFDaf::processTrack(): wrong matrix dimensions!");
332 1. / (pow(2. * TMath::Pi(), dimV / 2.) * sqrt(detV)) * exp(-0.5 * expArg[0][0]);
333 phi.push_back(thisPhi);
338 throw std::runtime_error(
"genf::GFDaf::processTrack(): chi2 cut too low");
339 sumPhiCut += 1. / (2 * TMath::Pi() * sqrt(detV)) * exp(-0.5 * cutVal /
fBeta.at(iBeta));
341 for (
unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
342 bk->
setNumber(
"p", planes.at(ipl)->at(ihit), phi.at(ihit) / (sumPhi + sumPhiCut));
356 for (
int irep = 0; irep < nreps; ++irep) {
360 TMatrixT<Double_t> cov = arep->
getCov();
361 for (
int i = 0; i < cov.GetNrows(); ++i) {
362 for (
int j = 0; j < cov.GetNcols(); ++j) {
382 if (TMath::IsNaN(det)) {
383 GFException e(
"Daf Gain: det of matrix is nan", __LINE__, __FILE__);
388 GFException e(
"cannot invert matrix in Daf - det=0", __LINE__, __FILE__);
395 const TMatrixT<Double_t>& V,
396 const TMatrixT<Double_t>& H,
401 TMatrixT<Double_t> Cinv;
403 TMatrixT<Double_t> Vinv;
406 TMatrixT<Double_t> Htransp(H);
409 TMatrixT<Double_t> covsum = Cinv + (p * Htransp * Vinv * H);
410 TMatrixT<Double_t> covsumInv;
413 return (covsumInv * Htransp * Vinv);
419 if (fabs(val - 0.01) < 1.
E-10) {
426 else if (fabs(val - 0.005) < 1.
E-10) {
433 else if (fabs(val - 0.001) < 1.
E-10) {
442 "genf::GFTrackCand::GFDafsetProbCut(): value not supported", __LINE__, __FILE__)
489 if (b1 <= 0)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b1 <= 0");
492 if (b2 <= 0 || b2 > b1)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b2 in wrong range");
496 if (b3 > b2)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b3 > b2");
500 if (b4 > b3)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b4 > b3");
504 if (b5 > b4)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b5 > b4");
508 if (b6 > b5)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b6 > b5");
512 if (b7 > b6)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b7 > b6");
516 if (b8 > b7)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b8 > b7");
520 if (b9 > b8)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b9 > b8");
523 if (b10 < 0.)
return;
524 if (b10 > b9)
throw std::runtime_error(
"genf::GFDaf::setBetas(): b10 > b9");
525 fBeta.push_back(b10);
void setDetPlane(std::string key, unsigned int index, const GFDetPlane &pl)
GFException & setNumbers(std::string, const std::vector< double > &)
set list of numbers with description
GFDaf()
Standard CTOR. Sets default values for annealing scheme and probablity cut.
const GFDetPlane & getReferencePlane() const
void invertMatrix(const TMatrixT< Double_t > &, TMatrixT< Double_t > &)
invert a matrix. First argument is matrix to be inverted, second is return by ref.
void setBetas(double b1, double b2, double b3=-1., double b4=-1., double b5=-1., double b6=-1., double b7=-1., double b8=-1., double b9=-1., double b10=-1.)
Configure the annealing scheme. In the current implementation you need to provide at least two temper...
GFAbsTrackRep * getTrackRep(int id) const
Accessor for track representations.
void setNumber(std::string key, unsigned int index, const double &num)
const TMatrixT< Double_t > & getState() const
void processTrack(GFTrack *)
Performs DAF fit on all track representations in a GFTrack.
Base Class for genfit track representations. Defines interface for track parameterizations.
void setCov(const TMatrixT< Double_t > &aCov)
void setMatrix(std::string key, unsigned int index, const TMatrixT< Double_t > &mat)
void bookGFDetPlanes(std::string key)
unsigned int getNumReps() const
Get number of track represenatations.
std::vector< double > fBeta
virtual void setData(const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
void addFailedHit(unsigned int irep, unsigned int id)
bool isFatal()
get fatal flag.
GFAbsRecoHit * getHit(int id) const
TMatrixT< Double_t > calcGain(const TMatrixT< Double_t > &cov, const TMatrixT< Double_t > &HitCov, const TMatrixT< Double_t > &H, const double &p)
Calculate Kalman Gain.
virtual double extrapolate(const GFDetPlane &plane, TMatrixT< Double_t > &statePred)
returns the tracklength spanned in this extrapolation
unsigned int hitFailed(unsigned int)
const TMatrixT< Double_t > & getCov() const
GFBookkeeping * getBK(int index=-1)
get GFBookKeeping object for particular track rep (default is cardinal rep)
std::map< int, double > chi2Cuts
void setProbCut(double val)
Set the probabilty cut for the weight calculation for the hits. Currently supported are the values 0...
void info()
print information in the exception object
const char * what() const noexcept override
standard error message handling for exceptions. use like "std::cerr << e.what();" ...
bool getNumber(std::string key, unsigned int index, double &num) const
unsigned int getNumHits() const
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
void setStatusFlag(int _val)
void bookMatrices(std::string key)
bool getMatrix(std::string key, unsigned int index, TMatrixT< Double_t > &mat) const
void bookNumbers(std::string key, double val=0.)
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
void getHitsByPlane(std::vector< std::vector< int > * > &retVal)
unsigned int getDim() const
returns dimension of state vector
void blowUpCovs(GFTrack *trk)
This is needed to blow up the covariance matrix before a fitting pass. The method drops off-diagonal ...
bool getDetPlane(std::string key, unsigned int index, GFDetPlane &pl) const