CbmRoot
Station3D.cxx
Go to the documentation of this file.
1 /*
2  * To change this license header, choose License Headers in Project Properties.
3  * To change this template file, choose Tools | Templates
4  * and open the template in the editor.
5  */
6 
7 #include <set>
8 
9 #include "Station3D.h"
10 #include "Tracker.h"
11 
12 void CbmBinned3DStation::SearchHits(
13  const CbmTrackParam2& stateVec,
14  Double_t,
15  std::function<void(CbmTBin::HitHolder&)> handleHit) {
16  CbmTrackParam2 minParams = Extrapolate(stateVec, fMinZ);
17  //Double_t minC[21];
18  //minParams.CovMatrix(minC);
19  CbmTrackParam2 maxParams = Extrapolate(stateVec, fMaxZ);
20  //Double_t maxC[21];
21  //maxParams.CovMatrix(maxC);
22  Double_t wXmin =
23  fNofSigmasX * std::sqrt(minParams.GetCovXX() + fDxSq + fScatXSq);
24  Double_t wXmax =
25  fNofSigmasX * std::sqrt(maxParams.GetCovXX() + fDxSq + fScatXSq);
26  Double_t wYmin =
27  fNofSigmasY * std::sqrt(minParams.GetCovYY() + fDySq + fScatYSq);
28  Double_t wYmax =
29  fNofSigmasY * std::sqrt(maxParams.GetCovYY() + fDySq + fScatYSq);
30 
31  Double_t xMin =
32  stateVec.GetTx() > 0 ? minParams.GetX() - wXmin : maxParams.GetX() - wXmax;
33  Double_t xMax =
34  stateVec.GetTx() > 0 ? maxParams.GetX() + wXmax : minParams.GetX() + wXmin;
35 
36  Double_t yMin =
37  stateVec.GetTy() > 0 ? minParams.GetY() - wYmin : maxParams.GetY() - wYmax;
38  Double_t yMax =
39  stateVec.GetTy() > 0 ? maxParams.GetY() + wYmax : minParams.GetY() + wYmin;
40 
41  int lowerXind = GetXInd(xMin);
42  int upperXind = GetXInd(xMax);
43  int lowerYind = GetYInd(yMin);
44  int upperYind = GetYInd(yMax);
45  int lowerTind = GetTInd(0);
46  int upperTind = GetTInd(0);
47 
48  for (int i = lowerYind; i <= upperYind; ++i) {
49  CbmYBin& yBin = fYBins[i];
50 
51  for (int j = lowerXind; j <= upperXind; ++j) {
52  CbmXBin& xBin = yBin[j];
53 
54  for (int k = lowerTind; k <= upperTind; ++k) {
55  CbmTBin& tBin = xBin[k];
56  std::list<CbmTBin::HitHolder>::iterator hitIter = tBin.HitsBegin();
57  std::list<CbmTBin::HitHolder>::iterator hitIterEnd = tBin.HitsEnd();
58 
59  for (; hitIter != hitIterEnd; ++hitIter) {
60  CbmTBin::HitHolder& hitHolder = *hitIter;
61 
62  if (fStage > hitHolder.stage) continue;
63 
64  handleHit(hitHolder);
65  }
66  }
67  }
68  }
69 }
70 
71 void CbmBinned3DStation::SearchHits(Segment& segment,
72  std::function<void(CbmTBin::HitHolder&)>) {
73  const CbmPixelHit* hit1 = segment.begin->hit;
74  Double_t x1 = hit1->GetX();
75  Double_t y1 = hit1->GetY();
76  Double_t dx1Sq = hit1->GetDx() * hit1->GetDx();
77  Double_t dy1Sq = hit1->GetDy() * hit1->GetDy();
78  const CbmPixelHit* hit2 = segment.end->hit;
79  Double_t x2 = hit2->GetX();
80  Double_t y2 = hit2->GetY();
81  Double_t dx2Sq = hit2->GetDx() * hit2->GetDx();
82  Double_t dy2Sq = hit2->GetDy() * hit2->GetDy();
83  Double_t segDeltaZ = hit2->GetZ() - hit1->GetZ();
84  Double_t tx = (x2 - x1) / segDeltaZ;
85  Double_t ty = (y2 - y1) / segDeltaZ;
86  Double_t middleZ = (hit1->GetZ() + hit2->GetZ()) / 2;
87  Double_t deltaZmin = fMinZ - middleZ;
88  Double_t deltaZmax = fMaxZ - middleZ;
89  Double_t coeff1Min = (0.5 - deltaZmin / segDeltaZ);
90  Double_t coeff1MinSq = coeff1Min * coeff1Min;
91  Double_t coeff1Max = (0.5 - deltaZmax / segDeltaZ);
92  Double_t coeff1MaxSq = coeff1Max * coeff1Max;
93  Double_t coeff2Min = (0.5 + deltaZmin / segDeltaZ);
94  Double_t coeff2MinSq = coeff2Min * coeff2Min;
95  Double_t coeff2Max = (0.5 + deltaZmax / segDeltaZ);
96  Double_t coeff2MaxSq = coeff2Max * coeff2Max;
97  Double_t tCoeff = std::sqrt(1 + tx * tx + ty * ty) / cbmBinnedSOL;
98 
99  Double_t searchT;
100  Double_t dtSq;
101  Double_t wT;
102  Double_t tMin;
103  Double_t tMax;
104 
105  if (0 == hit1->GetTimeError()) {
106  searchT = hit2->GetTime();
107  dtSq = hit2->GetTimeError() * hit2->GetTimeError();
108  wT = cbmBinnedSigma * std::sqrt(dtSq + fDtSq);
109  tMin = searchT + tCoeff * (fMinZ - hit2->GetZ()) - wT;
110  tMax = searchT + tCoeff * (fMaxZ - hit2->GetZ()) + wT;
111  } else {
112  searchT = (hit1->GetTime() + hit2->GetTime()) / 2;
113  dtSq = (hit1->GetTimeError() * hit1->GetTimeError()
114  + hit2->GetTimeError() * hit2->GetTimeError())
115  / 2;
116  wT = cbmBinnedSigma * std::sqrt(dtSq + fDtSq);
117  tMin = searchT + tCoeff * deltaZmin - wT;
118  tMax = searchT + tCoeff * deltaZmax + wT;
119  }
120 
121  Double_t wXmin =
122  fNofSigmasX
123  * std::sqrt(coeff1MinSq * dx1Sq + coeff2MinSq * dx2Sq + fDxSq + fScatXSq);
124  Double_t wXmax =
125  fNofSigmasX
126  * std::sqrt(coeff1MaxSq * dx1Sq + coeff2MaxSq * dx2Sq + fDxSq + fScatXSq);
127  Double_t wYmin =
128  fNofSigmasY
129  * std::sqrt(coeff1MinSq * dy1Sq + coeff2MinSq * dy2Sq + fDySq + fScatYSq);
130  Double_t wYmax =
131  fNofSigmasY
132  * std::sqrt(coeff1MaxSq * dy1Sq + coeff2MaxSq * dy2Sq + fDySq + fScatYSq);
133 
134  Double_t xMin = tx > 0 ? coeff1Min * x1 + coeff2Min * x2 - wXmin
135  : coeff1Max * x1 + coeff2Max * x2 - wXmax;
136  Double_t xMax = tx > 0 ? coeff1Max * x1 + coeff2Max * x2 + wXmax
137  : coeff1Min * x1 + coeff2Min * x2 + wXmin;
138 
139  Double_t yMin = ty > 0 ? coeff1Min * y1 + coeff2Min * y2 - wYmin
140  : coeff1Max * y1 + coeff2Max * y2 - wYmax;
141  Double_t yMax = ty > 0 ? coeff1Max * y1 + coeff2Max * y2 + wYmax
142  : coeff1Min * y1 + coeff2Min * y2 + wYmin;
143 
144  int lowerXind = GetXInd(xMin);
145  int upperXind = GetXInd(xMax);
146  int lowerYind = GetYInd(yMin);
147  int upperYind = GetYInd(yMax);
148  int lowerTind = GetTInd(tMin);
149  int upperTind = GetTInd(tMax);
150 
151  for (int i = lowerYind; i <= upperYind; ++i) {
152  CbmYBin& yBin = fYBins[i];
153 
154  for (int j = lowerXind; j <= upperXind; ++j) {
155  CbmXBin& xBin = yBin[j];
156 
157  for (int k = lowerTind; k <= upperTind; ++k) {
158  CbmTBin& tBin = xBin[k];
159  std::list<CbmTBin::HitHolder>::iterator hitIter = tBin.HitsBegin();
160  std::list<CbmTBin::HitHolder>::iterator hitIterEnd = tBin.HitsEnd();
161 
162  for (; hitIter != hitIterEnd; ++hitIter) {
163  const CbmPixelHit* hit = hitIter->hit;
164  Double_t deltaZ = hit->GetZ() - middleZ;
165  Double_t coeff1 = (0.5 - deltaZ / segDeltaZ);
166  Double_t coeff1Sq = coeff1 * coeff1;
167  Double_t coeff2 = (0.5 + deltaZ / segDeltaZ);
168  Double_t coeff2Sq = coeff2 * coeff2;
169  Double_t y = coeff1 * y1 + coeff2 * y2;
170  Double_t deltaY = hit->GetY() - y;
171 
172  if (deltaY * deltaY > fNofSigmasYSq
173  * (coeff1Sq * dy1Sq + coeff2Sq * dy2Sq
174  + hit->GetDy() * hit->GetDy() + fScatYSq))
175  continue;
176 
177  Double_t x = coeff1 * x1 + coeff2 * x2;
178  Double_t deltaX = hit->GetX() - x;
179 
180  if (deltaX * deltaX > fNofSigmasXSq
181  * (coeff1Sq * dx1Sq + coeff2Sq * dx2Sq
182  + hit->GetDx() * hit->GetDx() + fScatXSq))
183  continue;
184 
185  /*Double_t t = 0 == hit1->GetTimeError() ? searchT + tCoeff * (hit->GetZ() - hit2->GetZ()) : searchT + tCoeff * deltaZ;
186  Double_t deltaT = hit->GetTime() - t;
187 
188  if (deltaT * deltaT > cbmBinnedSigmaSq * (dtSq + hit->GetTimeError() * hit->GetTimeError()))
189  continue;*/
190 
191  Segment newSegment(segment.end, &(*hitIter));
192  pair<set<Segment, SegmentComp>::iterator, bool> ir =
193  fSegments.insert(newSegment);
194  segment.children.push_back(const_cast<Segment*>(&(*ir.first)));
195  //handleHit(*hitIter);
196  }
197  }
198  }
199  }
200 }
201 
202 void CbmBinned3DStation::SearchHits(
203  Double_t,
204  Double_t,
205  Double_t minY,
206  Double_t maxY,
207  Double_t minX,
208  Double_t maxX,
209  Double_t minT,
210  Double_t maxT,
211  std::function<void(CbmTBin::HitHolder&)> handleHit) {
212  int lowerXind = GetXInd(minX);
213  int upperXind = GetXInd(maxX);
214  int lowerYind = GetYInd(minY);
215  int upperYind = GetYInd(maxY);
216  int lowerTind = GetTInd(minT);
217  int upperTind = GetTInd(maxT);
218 
219  for (int i = lowerYind; i <= upperYind; ++i) {
220  CbmYBin& yBin = fYBins[i];
221 
222  for (int j = lowerXind; j <= upperXind; ++j) {
223  CbmXBin& xBin = yBin[j];
224 
225  for (int k = lowerTind; k <= upperTind; ++k) {
226  CbmTBin& tBin = xBin[k];
227  std::list<CbmTBin::HitHolder>::iterator hitIter = tBin.HitsBegin();
228  std::list<CbmTBin::HitHolder>::iterator hitIterEnd = tBin.HitsEnd();
229 
230  for (; hitIter != hitIterEnd; ++hitIter) {
231  const CbmPixelHit* hit = hitIter->hit;
232  Double_t y = hit->GetY();
233 
234  if (y < minY || y >= maxY) continue;
235 
236  Double_t x = hit->GetX();
237 
238  if (x < minX || x >= maxX) continue;
239 
240  Double_t t = hit->GetTime();
241 
242  if (t < minT || t >= maxT) continue;
243 
244  handleHit(*hitIter);
245  }
246  }
247  }
248  }
249 }
CbmBinnedStation::fDtSq
Double_t fDtSq
Definition: Station.h:630
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmXBin
Definition: Bins.h:67
CbmBinnedStation::Segment
Definition: Station.h:31
CbmBinnedStation::fStage
char fStage
Definition: Station.h:647
CbmYBin
Definition: Bins.h:100
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmTrackParam2::GetCovYY
Double_t GetCovYY() const
Definition: CbmTrackParam2.h:75
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmBinnedStation::fNofSigmasXSq
Double_t fNofSigmasXSq
Definition: Station.h:644
cbmBinnedSigma
const Double_t cbmBinnedSigma
Definition: Station.h:24
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmBinnedStation::GetXInd
int GetXInd(Double_t v) const
Definition: Station.h:512
CbmBinnedStation::fScatYSq
Double_t fScatYSq
Definition: Station.h:635
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmBinnedStation::fNofSigmasX
Double_t fNofSigmasX
Definition: Station.h:643
CbmBinned3DStation::SearchHits
void SearchHits(const CbmTrackParam2 &stateVec, Double_t stateZ, std::function< void(CbmTBin::HitHolder &)> handleHit)
Definition: Station3D.cxx:12
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmHit::GetTimeError
Double_t GetTimeError() const
Definition: CbmHit.h:76
CbmBinnedStation::Extrapolate
CbmTrackParam2 Extrapolate(const CbmTrackParam2 &parIn, Double_t zOut)
Definition: Station.h:95
CbmTBin::HitHolder::hit
const CbmPixelHit * hit
Definition: Bins.h:29
CbmBinnedStation::fNofSigmasY
Double_t fNofSigmasY
Definition: Station.h:645
CbmBinnedStation::Segment::children
std::list< Segment * > children
Definition: Station.h:35
cbmBinnedSOL
Double_t cbmBinnedSOL
Definition: GeoReader.cxx:107
CbmTBin::HitHolder
Definition: Bins.h:27
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmTBin::HitsEnd
std::list< HitHolder >::iterator HitsEnd()
Definition: Bins.h:51
CbmBinnedStation::fScatXSq
Double_t fScatXSq
Definition: Station.h:633
CbmBinned3DStation::fYBins
CbmYBin * fYBins
Definition: Station3D.h:147
CbmBinnedStation::fDySq
Double_t fDySq
Definition: Station.h:628
CbmBinnedStation::fDxSq
Double_t fDxSq
Definition: Station.h:626
Tracker.h
CbmTrackParam2::GetCovXX
Double_t GetCovXX() const
Definition: CbmTrackParam2.h:73
CbmBinnedStation::fMaxZ
Double_t fMaxZ
Definition: Station.h:612
CbmTrackParam2
Definition: CbmTrackParam2.h:68
CbmBinnedStation::fMinZ
Double_t fMinZ
Definition: Station.h:611
CbmTBin::HitsBegin
std::list< HitHolder >::iterator HitsBegin()
Definition: Bins.h:50
CbmBinnedStation::fNofSigmasYSq
Double_t fNofSigmasYSq
Definition: Station.h:646
CbmBinnedStation::GetTInd
int GetTInd(Double_t v) const
Definition: Station.h:534
CbmBinnedStation::GetYInd
int GetYInd(Double_t v) const
Definition: Station.h:523
CbmBinnedStation::Segment::end
CbmTBin::HitHolder * end
Definition: Station.h:33
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
Station3D.h
CbmPixelHit
Definition: CbmPixelHit.h:21
CbmTBin::HitHolder::stage
char stage
Definition: Bins.h:34
CbmBinnedStation::fSegments
std::set< Segment, SegmentComp > fSegments
Definition: Station.h:638
CbmTBin
Definition: Bins.h:25
CbmBinnedStation::Segment::begin
CbmTBin::HitHolder * begin
Definition: Station.h:32