CbmRoot
CbmTofTrackletTools.cxx
Go to the documentation of this file.
1 
7 #include "CbmTofTrackletTools.h"
8 #include "CbmTofHit.h"
9 #include "CbmTofTracklet.h"
10 #include "LKFMinuit.h"
11 
12 #include "FairLogger.h"
13 
14 #include "Rtypes.h" // for Double_t, Double32_t, Int_t, etc
15 #include "TDecompSVD.h"
16 #include "TMatrixD.h"
17 #include "TMatrixFSymfwd.h" // for TMatrixFSym
18 #include "TVectorD.h"
19 
20 using std::vector;
21 LKFMinuit fMinuit;
22 
23 CbmTofTrackletTools::CbmTofTrackletTools() {
24  fMinuit = *LKFMinuit::Instance();
25  //if( &fMinuit == NULL ) fMinuit.Initialize();
26 }
27 
28 CbmTofTrackletTools::~CbmTofTrackletTools() {}
29 
30 Double_t* CbmTofTrackletTools::Line3DFit(CbmTofTracklet* pTrk, Int_t iDetId) {
31  TGraph2DErrors* gr = new TGraph2DErrors();
32  Int_t NHit = 0;
33  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
34  if (iDetId == pTrk->GetTofDetIndex(iHit))
35  continue; // skip this hit of tracklet
36  gr->SetPoint(NHit,
37  pTrk->GetTofHitPointer(iHit)->GetX(),
38  pTrk->GetTofHitPointer(iHit)->GetY(),
39  pTrk->GetTofHitPointer(iHit)->GetZ());
40  gr->SetPointError(NHit,
41  pTrk->GetTofHitPointer(iHit)->GetDx(),
42  pTrk->GetTofHitPointer(iHit)->GetDy(),
43  pTrk->GetTofHitPointer(iHit)->GetDz());
44  NHit++;
45  }
46  // fit the graph now
47  Double_t pStart[4] = {0., 0., 0., 0.};
48  pStart[0] = pTrk->GetFitX(0.);
49  pStart[1] = (pTrk->GetTrackParameter())->GetTx();
50  pStart[2] = pTrk->GetFitY(0.);
51  pStart[3] = (pTrk->GetTrackParameter())->GetTy();
52  LOG(debug) << "Line3DFit init: X0 " << pStart[0] << ", TX " << pStart[1]
53  << ", Y0 " << pStart[2] << ", TY " << pStart[3];
54  fMinuit.DoFit(gr, pStart);
55  gr->Delete();
56  Double_t* dRes;
57  // LOG(info) << "Get fit results ";
58  dRes = fMinuit.GetParFit();
59  LOG(debug) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0]
60  << ", TX " << dRes[1] << ", Y0 " << dRes[2] << ", TY " << dRes[3]
61  << ", Chi2DoF: " << fMinuit.GetChi2DoF();
62  return dRes;
63 }
64 
65 Double_t CbmTofTrackletTools::FitTt(CbmTofTracklet* pTrk, Int_t iDetId) {
66  Int_t nValidHits = 0;
67  Int_t iHit1 = 0;
68  Double_t dDist1 = 100.;
69  //LOG(info) << "FitTt: " << pTrk->GetNofHits() << " hits, exclude " << Form("0x%08x",iDetId);
70  Double_t aR[pTrk->GetNofHits() - 1];
71  Double_t at[pTrk->GetNofHits() - 1];
72  Double_t ae[pTrk->GetNofHits() - 1];
73  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
74  if (iDetId == pTrk->GetTofDetIndex(iHit)) continue;
75  if (nValidHits == 0) {
76  iHit1 = iHit;
77  //LOG(info) << "Init Dist1 with " << iHit1;
78  Double_t* dPar = Line3DFit(pTrk, iDetId); // spatial fit without iDetId
79  dDist1 = pTrk->GetTofHitPointer(iHit1)->GetZ()
80  * TMath::Sqrt(1. + dPar[1] * dPar[1] + dPar[3] * dPar[3]);
81  //LOG(info) << "Dist1 = " << dDist1;
82  }
83  Double_t dSign = 1.;
84  if (pTrk->GetTofHitPointer(iHit)->GetZ()
85  < pTrk->GetTofHitPointer(iHit1)->GetZ())
86  dSign = -1.;
87  aR[nValidHits] = dDist1
88  + dSign
89  * pTrk->Dist3D(pTrk->GetTofHitPointer(iHit),
90  pTrk->GetTofHitPointer(iHit1));
91  at[nValidHits] = pTrk->GetTofHitPointer(iHit)->GetTime();
92  ae[nValidHits] = 0.1; // const timing error, FIXME (?)
93  //LOG(info) << "Init vector index " << nValidHits;
94  nValidHits++;
95  }
96  if (nValidHits < 2) return 0.;
97  if (nValidHits == 2) return (at[0] - at[1]) / (aR[0] - aR[1]);
98 
99  Double_t RRsum = 0; // Values will follow this procedure:
100  Double_t Rsum = 0; // $Rsum=\sum_{i}^{nValidHits}\frac{R_i}{e_i^2}$
101  Double_t tsum =
102  0; // where e_i will always be the error on the t measurement
103  Double_t esum =
104  0; // RR=R^2 in numerator, e denotes 1 in numerator , Rt= R*t in numerator
105  Double_t Rtsum = 0; //
106  Double_t sig_weight = 0; // ae^2
107  Double_t yoffset =
108  at[0] - 10; // T0 time offset to scale times to ns regime and not 10^10ns
109  for (Int_t i = 0; i < nValidHits; i++) {
110  at[i] -= yoffset; // substract offset
111  sig_weight = (ae[i] * ae[i]);
112  Rsum += (aR[i] / sig_weight);
113  tsum += (at[i] / sig_weight);
114  RRsum += (aR[i] * aR[i] / sig_weight);
115  Rtsum += (aR[i] * at[i] / sig_weight);
116  esum += (1 / sig_weight);
117  }
118  Double_t det_cov_mat =
119  esum * RRsum
120  - Rsum
121  * Rsum; // Determinant of inverted Covariance Matrix -> 1/det is common Faktor of Cavariance Matrix
122  //fT0=(RRsum*tsum-Rsum*Rtsum)/det_cov_mat; // Best Guess for time at origin
123  return (-Rsum * tsum + esum * Rtsum)
124  / det_cov_mat; // Best guess for inverted velocity
125  //fT0Err=TMath::Sqrt(RRsum/det_cov_mat); // sqrt of (0,0) in Covariance matrix -> error on fT0
126  //fTtErr=TMath::Sqrt(esum/det_cov_mat); // sqrt of (1,1) in Covariance Matrix -> error on fTt
127  //fT0TtCov=-Rsum/det_cov_mat; // (0,1)=(1,0) in Covariance Matrix -> cov(fT0,fTt)
128 
129  //fT0+=yoffset; // Adding yoffset again
130 }
131 
132 Double_t CbmTofTrackletTools::GetXdif(CbmTofTracklet* pTrk,
133  Int_t iDetId,
134  CbmTofHit* pHit) {
135  Double_t dXref = 0.;
136  Int_t iNref = 0;
137  Double_t dTx = 0;
138 
139  if (1) {
140  for (Int_t iHL = 0; iHL < pTrk->GetNofHits() - 1; iHL++) {
141  if (iDetId == pTrk->GetTofDetIndex(iHL) || 0 == pTrk->GetTofDetIndex(iHL))
142  continue; // exclude faked hits
143  for (Int_t iHH = iHL + 1; iHH < pTrk->GetNofHits(); iHH++) {
144  if (iDetId == pTrk->GetTofDetIndex(iHH)
145  || 0 == pTrk->GetTofDetIndex(iHH))
146  continue; // exclude faked hits
147  //dTt+=(pTrk->GetTofHitPointer(iHH)->GetTime()-pTrk->GetTofHitPointer(iHL)->GetTime())/(pTrk->GetTofHitPointer(iHH)->GetR()-pTrk->GetTofHitPointer(iHL)->GetR()); // for projective geometries only !!!
148  dTx += (pTrk->GetTofHitPointer(iHH)->GetX()
149  - pTrk->GetTofHitPointer(iHL)->GetX())
150  / (pTrk->GetTofHitPointer(iHH)->GetZ()
151  - pTrk->GetTofHitPointer(iHL)->GetZ());
152  iNref++;
153  }
154  }
155  dTx /= iNref;
156  } else {
157  dTx = pTrk->GetTrackParameter()->GetTx();
158  }
159 
160  iNref = 0;
161  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
162  if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit))
163  continue;
164 
165  Double_t dDZ = pHit->GetZ() - pTrk->GetTofHitPointer(iHit)->GetZ();
166  dXref += pTrk->GetTofHitPointer(iHit)->GetX() + dTx * dDZ;
167  iNref++;
168  }
169 
170  if (iNref == 0) {
171  LOG(error) << "DetId " << iDetId << ", Nref " << iNref << " sizes "
172  << pTrk->GetNofHits() << ", " << pTrk->GetNofHits();
173  return 1.E20;
174  }
175 
176  dXref /= iNref;
177 
178  return pHit->GetX() - dXref;
179 }
180 
181 Double_t CbmTofTrackletTools::GetYdif(CbmTofTracklet* pTrk,
182  Int_t iDetId,
183  CbmTofHit* pHit) {
184  Double_t dYref = 0.;
185  Int_t iNref = 0;
186  Double_t dTy = 0;
187 
188  if (1) {
189  for (Int_t iHL = 0; iHL < pTrk->GetNofHits() - 1; iHL++) {
190  if (iDetId == pTrk->GetTofDetIndex(iHL) || 0 == pTrk->GetTofDetIndex(iHL))
191  continue; // exclude faked hits
192  for (Int_t iHH = iHL + 1; iHH < pTrk->GetNofHits(); iHH++) {
193  if (iDetId == pTrk->GetTofDetIndex(iHH)
194  || 0 == pTrk->GetTofDetIndex(iHH))
195  continue; // exclude faked hits
196  dTy += (pTrk->GetTofHitPointer(iHH)->GetY()
197  - pTrk->GetTofHitPointer(iHL)->GetY())
198  / (pTrk->GetTofHitPointer(iHH)->GetZ()
199  - pTrk->GetTofHitPointer(iHL)->GetZ());
200  iNref++;
201  }
202  }
203  dTy /= iNref;
204  } else {
205  dTy = pTrk->GetTrackParameter()->GetTy();
206  }
207 
208  iNref = 0;
209  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
210  if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit))
211  continue;
212 
213  Double_t dDZ = pHit->GetZ() - pTrk->GetTofHitPointer(iHit)->GetZ();
214  dYref += pTrk->GetTofHitPointer(iHit)->GetY() + dTy * dDZ;
215  iNref++;
216  }
217 
218  if (iNref == 0) {
219  LOG(error) << "DetId " << iDetId << ", Nref " << iNref << " sizes "
220  << pTrk->GetNofHits() << ", " << pTrk->GetNofHits();
221  return 1.E20;
222  }
223 
224  dYref /= iNref;
225 
226  return pHit->GetY() - dYref;
227 }
228 
229 Double_t CbmTofTrackletTools::GetTdif(CbmTofTracklet* pTrk,
230  Int_t iDetId,
231  CbmTofHit* pHit) {
232  Double_t dTref = 0.;
233  Double_t Nref = 0;
234  Double_t dTt = 0.;
235  Int_t iNt = 0;
236 
237  if (0) {
238  for (Int_t iHL = 0; iHL < pTrk->GetNofHits() - 1; iHL++) {
239  if (iDetId == pTrk->GetTofDetIndex(iHL) || 0 == pTrk->GetTofDetIndex(iHL))
240  continue; // exclude faked hits
241  for (Int_t iHH = iHL + 1; iHH < pTrk->GetNofHits(); iHH++) {
242  if (iDetId == pTrk->GetTofDetIndex(iHH)
243  || 0 == pTrk->GetTofDetIndex(iHH))
244  continue; // exclude faked hits
245  //dTt+=(pTrk->GetTofHitPointer(iHH)->GetTime()-pTrk->GetTofHitPointer(iHL)->GetTime())/(pTrk->GetTofHitPointer(iHH)->GetR()-pTrk->GetTofHitPointer(iHL)->GetR()); // for projective geometries only !!!
246  Double_t dSign = 1.;
247  if (pTrk->GetTofHitPointer(iHH)->GetZ()
248  < pTrk->GetTofHitPointer(iHL)->GetZ())
249  dSign = -1.;
250  dTt += (pTrk->GetTofHitPointer(iHH)->GetTime()
251  - pTrk->GetTofHitPointer(iHL)->GetTime())
252  / pTrk->Dist3D(pTrk->GetTofHitPointer(iHH),
253  pTrk->GetTofHitPointer(iHL))
254  * dSign;
255  iNt++;
256  }
257  }
258 
259  if (iNt == 0) {
260  LOG(error) << "No valid hit pair ";
261  return 1.E20;
262  }
263  dTt /= (Double_t) iNt;
264  } else {
265  dTt = FitTt(pTrk, iDetId);
266  }
267 
268  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
269  if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit))
270  continue;
271  //dTref += pTrk->GetTofHitPointer(iHit)->GetTime() - dTt*(pTrk->GetTofHitPointer(iHit)->GetR()-pHit->GetR());
272  Double_t dSign = 1.;
273  if (pTrk->GetTofHitPointer(iHit)->GetZ() < pHit->GetZ()) dSign = -1;
274  dTref += pTrk->GetTofHitPointer(iHit)->GetTime()
275  - dTt * dSign * pTrk->Dist3D(pTrk->GetTofHitPointer(iHit), pHit);
276  Nref++;
277  }
278  if (Nref == 0) {
279  LOG(error) << "DetId " << iDetId << ", Nref " << Nref << " sizes "
280  << pTrk->GetNofHits();
281  return 1.E20;
282  }
283  dTref /= (Double_t) Nref;
284  Double_t dTdif = pHit->GetTime() - dTref;
285  // LOG(debug) << "iSt "<< iSt<<" DetId "<<iDetId<<", Nref "<<Nref<<" Tdif
286  // "<<dTdif;
287  return dTdif;
288 }
289 
290 Double_t CbmTofTrackletTools::GetTexpected(CbmTofTracklet* pTrk,
291  Int_t iDetId,
292  CbmTofHit* pHit) {
293  Double_t dTref = 0.;
294  Double_t Nref = 0;
295  Double_t dTt = 0.;
296 
297  dTt = FitTt(pTrk, iDetId);
298 
299  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
300  if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit))
301  continue;
302  Double_t dSign = 1.;
303  if (pTrk->GetTofHitPointer(iHit)->GetZ() < pHit->GetZ()) dSign = -1;
304  dTref += pTrk->GetTofHitPointer(iHit)->GetTime()
305  - dTt * dSign * pTrk->Dist3D(pTrk->GetTofHitPointer(iHit), pHit);
306  Nref++;
307  }
308  if (Nref == 0) {
309  LOG(error) << "DetId " << iDetId << ", Nref " << Nref << " sizes "
310  << pTrk->GetNofHits();
311  return 1.E20;
312  }
313  dTref /= (Double_t) Nref;
314  return dTref;
315 }
316 
317 ClassImp(CbmTofTrackletTools)
CbmTofTrackletTools::CbmTofTrackletTools
CbmTofTrackletTools()
Definition: CbmTofTrackletTools.cxx:23
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmTofTrackletTools
contains fits and resolution functions
Definition: CbmTofTrackletTools.h:21
LKFMinuit::GetParFit
double * GetParFit()
Definition: LKFMinuit.h:37
CbmTofTracklet::GetTrackParameter
CbmTofTrackletParam * GetTrackParameter()
Definition: CbmTofTracklet.h:96
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmTofTracklet::GetNofHits
Int_t GetNofHits() const
Definition: CbmTofTracklet.h:48
CbmTofTrackletTools::Line3DFit
Double_t * Line3DFit(CbmTofTracklet *pTrk, Int_t iDetId)
Definition: CbmTofTrackletTools.cxx:30
CbmTofTrackletParam::GetTy
Double_t GetTy() const
Definition: CbmTofTrackletParam.h:53
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmTofTrackletParam::GetTx
Double_t GetTx() const
Definition: CbmTofTrackletParam.h:52
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
LKFMinuit.h
CbmTofTrackletTools::GetTexpected
virtual Double_t GetTexpected(CbmTofTracklet *pTrk, Int_t iDetId, CbmTofHit *pHit)
Definition: CbmTofTrackletTools.cxx:290
CbmTofTracklet
Provides information on attaching a TofHit to a TofTrack.
Definition: CbmTofTracklet.h:25
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
gr
Char_t * gr
Definition: CbmEvDisTracks.cxx:289
LKFMinuit
Definition: LKFMinuit.h:19
fMinuit
LKFMinuit fMinuit
Definition: CbmTofTrackletTools.cxx:21
CbmTofTracklet::Dist3D
virtual Double_t Dist3D(CbmTofHit *pHit0, CbmTofHit *pHit1)
Definition: CbmTofTracklet.cxx:557
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmTofTracklet.h
CbmTofTrackletTools::GetYdif
virtual Double_t GetYdif(CbmTofTracklet *pTrk, Int_t iDetId, CbmTofHit *pHit)
Definition: CbmTofTrackletTools.cxx:181
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmTofTrackletTools::FitTt
Double_t FitTt(CbmTofTracklet *pTrk, Int_t iDetId)
Definition: CbmTofTrackletTools.cxx:65
LKFMinuit::DoFit
int DoFit(TGraph2DErrors *gr, double pStart[])
Definition: LKFMinuit.cxx:28
CbmTofTrackletTools::GetXdif
virtual Double_t GetXdif(CbmTofTracklet *pTrk, Int_t iDetId, CbmTofHit *pHit)
Definition: CbmTofTrackletTools.cxx:132
LKFMinuit::GetChi2DoF
double GetChi2DoF()
Definition: LKFMinuit.h:39
CbmTofTracklet::GetFitX
Double_t GetFitX(Double_t Z)
Definition: CbmTofTracklet.cxx:516
CbmTofTracklet::GetFitY
Double_t GetFitY(Double_t Z)
Definition: CbmTofTracklet.cxx:520
CbmTofTrackletTools.h
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
CbmTofTrackletTools::~CbmTofTrackletTools
virtual ~CbmTofTrackletTools()
Definition: CbmTofTrackletTools.cxx:28
CbmTofTracklet::GetTofDetIndex
Int_t GetTofDetIndex(Int_t ind) const
Definition: CbmTofTracklet.h:73
LKFMinuit::Instance
static LKFMinuit * Instance()
Definition: LKFMinuit.h:21
CbmTofTracklet::GetTofHitPointer
CbmTofHit * GetTofHitPointer(Int_t ind)
Definition: CbmTofTracklet.h:72
CbmHit::GetDz
Double_t GetDz() const
Definition: CbmHit.h:71
CbmTofTrackletTools::GetTdif
virtual Double_t GetTdif(CbmTofTracklet *pTrk, Int_t iDetId, CbmTofHit *pHit)
Definition: CbmTofTrackletTools.cxx:229