FairRoot/PandaRoot
GFTools.cxx
Go to the documentation of this file.
1 
2 #include <GFTools.h>
3 #include <typeinfo>
4 TMatrixT<double> GFTools::getSmoothedPos(GFTrack* trk, unsigned int irep, unsigned int ihit) {
5 
6  TMatrixT<double> smoothed_state;
7  TMatrixT<double> smoothed_cov;
8  TMatrixT<double> pos;
9 
10  if(GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov)) {
11 
12  TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
13  TMatrixT<double> pos_tmp(H * smoothed_state);
14  pos.ResizeTo(pos_tmp);
15  pos = pos_tmp;
16 
17  }
18  return pos;
19 
20 }
21 
22 TMatrixT<double> GFTools::getBiasedSmoothedPos(GFTrack* trk, unsigned int irep, unsigned int ihit) {
23 
24 TMatrixT<double> smoothed_state;
25  TMatrixT<double> smoothed_cov;
26  TMatrixT<double> pos;
27 
28  if(GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov)) {
29 
30  TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
31  H.Print();smoothed_state.Print();
32  TMatrixT<double> pos_tmp(H * smoothed_state);
33  pos.ResizeTo(pos_tmp);
34  pos = pos_tmp;
35 
36  }
37  return pos;
38 
39 }
40 
41 TMatrixT<double> GFTools::getSmoothedCov(GFTrack* trk, unsigned int irep, unsigned int ihit) {
42 
43  TMatrixT<double> smoothed_state;
44  TMatrixT<double> smoothed_cov;
45  GFDetPlane pl;
46 
47  GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl);
48  TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
49 
50  TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
51  TMatrixT<double> cov(H,TMatrixT<double>::kMult,cov_tmp);
52 
53  return cov;
54 
55 }
56 
57 TMatrixT<double> GFTools::getBiasedSmoothedCov(GFTrack* trk, unsigned int irep, unsigned int ihit) {
58 
59  TMatrixT<double> smoothed_state;
60  TMatrixT<double> smoothed_cov;
61  GFDetPlane pl;
62 
63  GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl);
64  TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
65 
66  TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
67  TMatrixT<double> cov(H,TMatrixT<double>::kMult,cov_tmp);
68 
69  return cov;
70 
71 }
72 
73 bool GFTools::getSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov) {
74 
75  GFDetPlane pl;
76  TMatrixT<double> auxInfo;
77  return GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl, auxInfo);
78 
79 }
80 
81 bool GFTools::getBiasedSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov) {
82 
83  GFDetPlane pl;
84  TMatrixT<double> auxInfo;
85  return GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl, auxInfo);
86 
87 }
88 
89 bool GFTools::getSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov, GFDetPlane& smoothing_plane) {
90 
91  TMatrixT<double> auxInfo;
92  return GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, smoothing_plane, auxInfo);
93 
94 }
95 
96 bool GFTools::getBiasedSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov, GFDetPlane& smoothing_plane) {
97 
98  TMatrixT<double> auxInfo;
99  return GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, smoothing_plane, auxInfo);
100 
101 }
102 
103 bool GFTools::getSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov, GFDetPlane& smoothing_plane, TMatrixT<double>& auxInfo) {
104 
105  if(!trk->getSmoothing()) {
106  std::cout << "Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
107  TMatrixT<double> failed(1,1);
108  return false;
109  }
110 
111  if(ihit >= trk->getNumHits()) {
112  std::cout << "Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
113  TMatrixT<double> failed(1,1);
114  failed(0,0) = 0;
115  return false;
116  }
117 
118  TMatrixT<double> fUpSt;
119  TMatrixT<double> fUpCov;
120  TMatrixT<double> fAuxInfo;
121  TMatrixT<double>* fAuxInfoP;
122  TMatrixT<double> bUpSt;
123  TMatrixT<double> bUpCov;
124  TMatrixT<double> bAuxInfo;
125  TMatrixT<double>* bAuxInfoP;
126 
127  GFDetPlane fPl;
128  GFDetPlane bPl;
129 
130  GFAbsTrackRep* rep = trk->getTrackRep(irep)->clone();
131 
132  if(trk->getTrackRep(irep)->hasAuxInfo()) {
133  trk->getBK(irep)->getMatrix("fAuxInfo",ihit,auxInfo);
134  }
135 
136  if(ihit == 0) {
137  trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
138  trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
139  trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
140  trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);
141  if(trk->getTrackRep(irep)->hasAuxInfo()) {
142  trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
143  bAuxInfoP = &bAuxInfo;
144  } else {
145  bAuxInfoP = NULL;
146  }
147  if(bUpSt.GetNrows() == 0) return false;
148  rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
149  rep->extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
150  return true;
151  }
152 
153  if(ihit == trk->getNumHits()-1) {
154  trk->getBK(irep)->getMatrix("fUpSt",ihit-1,fUpSt);
155  trk->getBK(irep)->getMatrix("fUpCov",ihit-1,fUpCov);
156  trk->getBK(irep)->getDetPlane("fPl",ihit-1,fPl);
157  trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
158  if(trk->getTrackRep(irep)->hasAuxInfo()) {
159  trk->getBK(irep)->getMatrix("fAuxInfo",ihit-1,fAuxInfo);
160  fAuxInfoP = &fAuxInfo;
161  } else {
162  fAuxInfoP = NULL;
163  }
164  if(fUpSt.GetNrows() == 0) return false;
165  rep->setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
166  rep->extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
167  return true;
168  }
169 
170  trk->getBK(irep)->getMatrix("fUpSt",ihit-1,fUpSt);
171  trk->getBK(irep)->getMatrix("fUpCov",ihit-1,fUpCov);
172  trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
173  trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
174  if(trk->getTrackRep(irep)->hasAuxInfo()) {
175  trk->getBK(irep)->getMatrix("fAuxInfo",ihit-1,fAuxInfo);
176  trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
177  fAuxInfoP = &fAuxInfo;
178  bAuxInfoP = &bAuxInfo;
179  } else {
180  fAuxInfoP = bAuxInfoP = NULL;
181  }
182  trk->getBK(irep)->getDetPlane("fPl",ihit-1,fPl);
183  trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
184  trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);
185 
186 
187  TMatrixT<double> fSt;
188  TMatrixT<double> fCov;
189  TMatrixT<double> bSt;
190  TMatrixT<double> bCov;
191 
192  if(fUpSt.GetNrows() == 0 || bUpSt.GetNrows() == 0) return false;
193 
194  rep->setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
195  rep->extrapolate(smoothing_plane,fSt,fCov);
196  rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
197  rep->extrapolate(smoothing_plane,bSt,bCov);
198 
199  delete rep;
200 
201  TMatrixT<double> fCovInvert;
202  TMatrixT<double> bCovInvert;
203 
204  GFTools::invertMatrix(fCov, fCovInvert);
205  GFTools::invertMatrix(bCov, bCovInvert);
206 
207  GFTools::invertMatrix(fCovInvert + bCovInvert, smoothed_cov);
208 
209  smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
210  smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);
211 
212  return true;
213 
214 }
215 
216 bool GFTools::getBiasedSmoothedData(GFTrack* trk, unsigned int irep, unsigned int ihit, TMatrixT<double>& smoothed_state, TMatrixT<double>& smoothed_cov, GFDetPlane& smoothing_plane, TMatrixT<double>& auxInfo) {
217 
218  if(!trk->getSmoothing()) {
219  std::cout << "Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
220  TMatrixT<double> failed(1,1);
221  return false;
222  }
223 
224  if(ihit >= trk->getNumHits()) {
225  std::cout << "Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
226  TMatrixT<double> failed(1,1);
227  failed(0,0) = 0;
228  return false;
229  }
230 
231  TMatrixT<double> bUpSt;
232  TMatrixT<double> bUpCov;
233  TMatrixT<double> bAuxInfo;
234  TMatrixT<double>* bAuxInfoP;
235  GFDetPlane bPl;
236 
237  GFAbsTrackRep* rep = trk->getTrackRep(irep)->clone();
238 
239  if(trk->getTrackRep(irep)->hasAuxInfo()) {
240  trk->getBK(irep)->getMatrix("fAuxInfo",ihit,auxInfo);
241  }
242 
243  if(ihit == 0) {
244  trk->getBK(irep)->getMatrix("bUpSt",ihit,smoothed_state);
245  trk->getBK(irep)->getMatrix("bUpCov",ihit,smoothed_cov);
246  trk->getBK(irep)->getDetPlane("bPl",ihit,smoothing_plane);
247  return true;
248  }
249 
250  if(ihit == trk->getNumHits()-1) {
251  trk->getBK(irep)->getMatrix("fUpSt",ihit,smoothed_state);
252  trk->getBK(irep)->getMatrix("fUpCov",ihit,smoothed_cov);
253  trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
254  return true;
255  }
256 
257  TMatrixT<double> fSt;
258  TMatrixT<double> fCov;
259 
260  trk->getBK(irep)->getMatrix("fUpSt",ihit,fSt);
261  trk->getBK(irep)->getMatrix("fUpCov",ihit,fCov);
262  trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
263  trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
264  if(trk->getTrackRep(irep)->hasAuxInfo()) {
265  trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
266  bAuxInfoP = &bAuxInfo;
267  } else {
268  bAuxInfoP = NULL;
269  }
270  trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
271  trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);
272 
273 
274  TMatrixT<double> bSt;
275  TMatrixT<double> bCov;
276 
277  if(bUpSt.GetNrows() == 0) return false;
278 
279  rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
280  rep->extrapolate(smoothing_plane,bSt,bCov);
281 
282  delete rep;
283 
284  TMatrixT<double> fCovInvert;
285  TMatrixT<double> bCovInvert;
286 
287  GFTools::invertMatrix(fCov, fCovInvert);
288  GFTools::invertMatrix(bCov, bCovInvert);
289 
290  GFTools::invertMatrix(fCovInvert + bCovInvert, smoothed_cov);
291 
292  smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
293  smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);
294 
295  return true;
296 
297 }
298 
299 GFDetPlane GFTools::getSmoothingPlane(GFTrack* trk, unsigned int irep, unsigned int ihit) {
300 
301  GFDetPlane pl;
302  trk->getBK(irep)->getDetPlane("fPl",ihit,pl);
303  return pl;
304 
305 }
306 
307 
308 /*
309 void GFTools::invertMatrix(const TMatrixT<double>& mat, TMatrixT<double>& inv){
310  inv.ResizeTo(mat);
311  // since ROOT has issues with absolute matrix entries that are below 10**-16 we force the 11 element to be one
312  //double factor = 1./mat[0][0];
313  //inv = (factor*mat);
314  inv = mat;
315  inv.SetTol(1E-150);
316  double det=0.;
317  inv.Invert(&det);
318  if(TMath::IsNaN(det)) {
319  GFException e("GFTools::invertMatrix(): det of matrix is nan",__LINE__,__FILE__);
320  e.setFatal();
321  throw e;
322  }
323  if(det==0){
324  GFException e("cannot invert matrix GFTools::invertMatrix() - det=0",
325  __LINE__,__FILE__);
326  e.setFatal();
327  throw e;
328  }
329  //recover multiplication with factor at the beginning
330  //inv *= factor;
331 }
332 */
333 
334 void GFTools::invertMatrix(const TMatrixT<double>& mat, TMatrixT<double>& inv){
335  inv.ResizeTo(mat);
336  bool status = 0;
337  TDecompSVD invertAlgo(mat);
338 
339  // check if numerical limits are reached (i.e at least one entry < 1E-100 and/or at least one entry > 1E100)
340  if (!(mat<1.E100) || !(mat>-1.E100)){
341  GFException e("cannot invert matrix GFTools::invertMatrix(), entries too big (>1E100)",
342  __LINE__,__FILE__);
343  e.setFatal();
344  throw e;
345  }
346 
347  invertAlgo.SetTol(1E-50); //this is a hack because a tolerance of 1E-22 does not make any sense for doubles only for floats
348 
349  inv = invertAlgo.Invert(status);
350  if(status == 0){
351  GFException e("cannot invert matrix GFTools::invertMatrix(), status = 0",
352  __LINE__,__FILE__);
353  e.setFatal();
354  throw e;
355  }
356 }
357 
358 double GFTools::getSmoothedChiSqu(GFTrack* const trk, unsigned int irep, unsigned int ihit){
359  TMatrixT<double> smoothed_state;
360  TMatrixT<double> smoothed_cov;
361  GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov);
362  TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
363  TMatrixT<double> HT(TMatrixT<double>::kTransposed, H);
364  TMatrixT<double> pos = H * smoothed_state;
365  TMatrixT<double> m = trk->getHit(ihit)->getRawHitCoord(); //measurement of hit
366  TMatrixT<double> V = trk->getHit(ihit)->getRawHitCov(); //covariance matrix of hit
367  TMatrixT<double> res = m - pos;
368  TMatrixT<double> resT(TMatrixT<double>::kTransposed, res);
369  TMatrixT<double> R = V - H*smoothed_cov*HT;
370  TMatrixT<double> invR;
371  invertMatrix(R,invR);
372  TMatrixT<double> smoothedChiSqu = resT*invR*res;
373  return smoothedChiSqu[0][0];
374 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
GFTools::getBiasedSmoothedData
bool getBiasedSmoothedData(GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov)
Get biased smoothed state vector and state covariance.
Definition: GFTools.cxx:81
GFAbsTrackRep
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:80
GFTrack::getBK
GFBookkeeping * getBK(int index=-1)
get GFBookKeeping object for particular track rep (default is cardinal rep)
Definition: GFTrack.h:326
GFTools::getBiasedSmoothedPos
TMatrixT< double > getBiasedSmoothedPos(GFTrack *trk, unsigned int irep, unsigned int ihit)
Get biased smoothed track position in plane coordinates.
Definition: GFTools.cxx:22
GFTrack::getNumHits
unsigned int getNumHits() const
Definition: GFTrack.h:148
res
Int_t res
Definition: anadigi.C:166
m
__m128 m
Definition: P4_F32vec4.h:28
GFTools::getSmoothingPlane
GFDetPlane getSmoothingPlane(GFTrack *trk, unsigned int irep, unsigned int ihit)
Get smoothing plane.
Definition: GFTools.cxx:299
GFDetPlane
Detector plane genfit geometry class.
Definition: GFDetPlane.h:59
GFTrack
Track object for genfit. genfit algorithms work on these objects.
Definition: GFTrack.h:60
GFTrack::getSmoothing
bool getSmoothing()
Read back if smoothing is/was turned on or off for this track.
Definition: GFTrack.h:405
GFTools::getSmoothedData
bool getSmoothedData(GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov)
Get smoothed state vector and state covariance.
Definition: GFTools.cxx:73
GFTools::getSmoothedPos
TMatrixT< double > getSmoothedPos(GFTrack *trk, unsigned int irep, unsigned int ihit)
Get smoothed track position in plane coordinates.
Definition: GFTools.cxx:4
GFBookkeeping::getDetPlane
bool getDetPlane(std::string key, unsigned int index, GFDetPlane &pl)
Definition: GFBookkeeping.cxx:323
GFTools::invertMatrix
void invertMatrix(const TMatrixT< double > &mat, TMatrixT< double > &inv)
Invert a matrix, throwing GFException when inversion fails.
Definition: GFTools.cxx:334
GFAbsTrackRep::hasAuxInfo
virtual bool hasAuxInfo()
See if the track representation has auxillary information stored.
Definition: GFAbsTrackRep.h:329
GFAbsRecoHit::getHMatrix
virtual TMatrixT< double > getHMatrix(const GFAbsTrackRep *stateVector)=0
Get transformation matrix. Transformation between hit coordinates and track representation coordinate...
GFAbsRecoHit::getRawHitCoord
TMatrixT< double > getRawHitCoord() const
Get raw hit coordinates.
Definition: GFAbsRecoHit.h:158
GFAbsTrackRep::extrapolate
virtual double extrapolate(const GFDetPlane &plane, TMatrixT< double > &statePred)
returns the tracklength spanned in this extrapolation
Definition: GFAbsTrackRep.cxx:42
GFTools::getSmoothedChiSqu
double getSmoothedChiSqu(GFTrack *const trk, unsigned int irep, unsigned int ihit)
Get smoothed chi2 for a specific hit (ihit).
Definition: GFTools.cxx:358
GFTrack::getTrackRep
GFAbsTrackRep * getTrackRep(int id) const
Accessor for track representations.
Definition: GFTrack.h:186
trk
GFTrack * trk
Definition: checkgenfit.C:13
GFBookkeeping::getMatrix
bool getMatrix(std::string key, unsigned int index, TMatrixT< double > &mat)
Definition: GFBookkeeping.cxx:303
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
GFTrack::getHit
GFAbsRecoHit * getHit(int id) const
Definition: GFTrack.h:144
GFException::setFatal
void setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:77
GFAbsTrackRep::setData
virtual void setData(const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
Puts the track representation in a given state.
Definition: GFAbsTrackRep.h:267
GFAbsRecoHit::getRawHitCov
TMatrixT< double > getRawHitCov() const
Get raw hit covariances.
Definition: GFAbsRecoHit.h:153
R
Double_t R
Definition: checkhelixhit.C:61
GFTools::getBiasedSmoothedCov
TMatrixT< double > getBiasedSmoothedCov(GFTrack *trk, unsigned int irep, unsigned int ihit)
Get biased smoothed track covariance in plane coordinates.
Definition: GFTools.cxx:57
GFAbsTrackRep::clone
virtual GFAbsTrackRep * clone() const =0
status
int status[10]
Definition: f_Init.h:28
GFTools::getSmoothedCov
TMatrixT< double > getSmoothedCov(GFTrack *trk, unsigned int irep, unsigned int ihit)
Get smoothed track covariance in plane coordinates.
Definition: GFTools.cxx:41
Generated on Wed Apr 3 2019 10:02:46 for FairRoot/PandaRoot by   doxygen 1.8.5