CbmRoot
CbmStsAlgoFindHits.cxx
Go to the documentation of this file.
1 
5 #include "CbmStsAlgoFindHits.h"
6 
7 #include <iostream>
8 
9 #include "CbmStsCluster.h"
10 #include "CbmStsHit.h"
11 #include "CbmStsParSensor.h"
12 #include "CbmStsParSensorCond.h"
13 #include <TGeoMatrix.h>
14 #include <TMath.h>
15 
16 #include "CbmDigiManager.h"
17 #include "CbmStsDigi.h"
18 
19 using std::pair;
20 using std::vector;
21 
22 
23 // ----- Constructor ---------------------------------------------------
25 // -------------------------------------------------------------------------
26 
27 
28 // ----- Create a new hit ----------------------------------------------
29 void CbmStsAlgoFindHits::CreateHit(Double_t xLocal,
30  Double_t yLocal,
31  Double_t varX,
32  Double_t varY,
33  Double_t varXY,
34  const CbmStsCluster& clusterF,
35  const CbmStsCluster& clusterB,
36  UInt_t indexF,
37  UInt_t indexB,
38  Double_t du,
39  Double_t dv) {
40 
41  // --- Check output array
42  assert(fHits);
43 
44  // --- Transform into global coordinate system
45  Double_t local[3] = {xLocal, yLocal, 0.};
46  Double_t global[3];
47  if (fMatrix)
48  fMatrix->LocalToMaster(local, global);
49  else {
50  global[0] = local[0];
51  global[1] = local[1];
52  global[2] = local[2];
53  } //? no transformation matrix available
54 
55  // We assume here that the local-to-global transformations is only translation
56  // plus maybe rotation upside down or front-side back. In that case, the
57  // global covariance matrix is the same as the local one.
58  // TODO: Proper transformation of covariance matrix
59  Double_t error[3] = {TMath::Sqrt(varX), TMath::Sqrt(varY), 0.};
60 
61  // --- Calculate hit time (average of cluster times)
62  Double_t hitTime = 0.5 * (clusterF.GetTime() + clusterB.GetTime());
63  Double_t etF = clusterF.GetTimeError();
64  Double_t etB = clusterB.GetTimeError();
65  Double_t hitTimeError = 0.5 * TMath::Sqrt(etF * etF + etB * etB);
66 
67  // --- Create hit
68  fHits->emplace_back(fAddress,
69  global,
70  error,
71  varXY,
72  indexF,
73  indexB,
74  hitTime,
75  hitTimeError,
76  du,
77  dv);
78 
79  return;
80 }
81 // -------------------------------------------------------------------------
82 
83 
84 // ----- Algorithm execution -------------------------------------------
85 Long64_t CbmStsAlgoFindHits::Exec(const vector<CbmStsCluster>& clustersF,
86  const vector<CbmStsCluster>& clustersB,
87  vector<CbmStsHit>& hits,
88  UInt_t address,
89  Double_t timeCutSig,
90  Double_t timeCutAbs,
91  Double_t dY,
92  UInt_t nStrips,
93  Double_t pitch,
94  Double_t stereoF,
95  Double_t stereoB,
96  Double_t lorentzF,
97  Double_t lorentzB,
98  TGeoHMatrix* matrix) {
99 
100  // Assert validity of parameters
101  assert(nStrips > 0);
102  assert(pitch > 0.);
103  assert(TMath::Abs(stereoF - stereoB) > 0.1);
104 
105  // Set internal parameters
106  fAddress = address;
107  fNofStrips = nStrips;
108  fDy = dY;
109  fPitch = pitch;
110  fStereoF = stereoF;
111  fStereoB = stereoB;
112  fLorentzF = lorentzF;
113  fLorentzB = lorentzB;
114  fMatrix = matrix;
115  fHits = &hits;
116  fHits->clear();
117  fTanStereoF = TMath::Tan(fStereoF * TMath::DegToRad());
118  fTanStereoB = TMath::Tan(fStereoB * TMath::DegToRad());
120  fDx = Double_t(fNofStrips) * fPitch;
121 
122  // Determine the maximum cluster time errors
123  Double_t maxTerrF = 0.;
124  for (auto& cluster : clustersF) {
125  if (cluster.GetTimeError() > maxTerrF) maxTerrF = cluster.GetTimeError();
126  }
127  Double_t maxTerrB = 0.;
128  for (auto& cluster : clustersB) {
129  if (cluster.GetTimeError() > maxTerrB) maxTerrB = cluster.GetTimeError();
130  }
131  const Double_t maxTerrF2 = maxTerrF * maxTerrF;
132  const Double_t maxTerrB2 = maxTerrB * maxTerrB;
133  const Double_t max_sigma_both = 4. * TMath::Sqrt(maxTerrF2 + maxTerrB2);
134 
135  // --- Loop over front and back side clusters
136  Long64_t nHits = 0;
137  Int_t startB = 0;
138  for (UInt_t iClusterF = 0; iClusterF < clustersF.size(); iClusterF++) {
139  const CbmStsCluster& clusterF = clustersF[iClusterF];
140 
141  Double_t tF = clusterF.GetTime();
142  Double_t tFerr2 = clusterF.GetTimeError() * clusterF.GetTimeError();
143  Double_t max_sigma = 4. * TMath::Sqrt(tFerr2 + maxTerrB2);
144 
145  for (UInt_t iClusterB = startB; iClusterB < clustersB.size(); iClusterB++) {
146  const CbmStsCluster& clusterB = clustersB[iClusterB];
147 
148  Double_t timeDiff = tF - clusterB.GetTime();
149  if ((timeDiff > 0) && (timeDiff > max_sigma_both)) {
150  startB++;
151  continue;
152  } else if ((timeDiff > 0) && (timeDiff > max_sigma)) {
153  continue;
154  } else if ((timeDiff < 0) && (fabs(timeDiff) > max_sigma))
155  break;
156 
157  // Cut on time difference of the two clusters
158  Double_t timeCut = -1.;
159  if (timeCutAbs > 0.)
160  timeCut = timeCutAbs; // absolute cut value
161  else {
162  if (timeCutSig > 0.) {
163  Double_t eF = clusterF.GetTimeError();
164  Double_t eB = clusterB.GetTimeError();
165  timeCut = timeCutSig * TMath::Sqrt(eF * eF + eB * eB);
166  } //? cut calculated from cluster errors
167  }
168  if (fabs(clusterF.GetTime() - clusterB.GetTime()) > timeCut) continue;
169 
170  // --- Calculate intersection points
171  Int_t nOfHits =
172  IntersectClusters(clusterF, clusterB, iClusterF, iClusterB);
173  nHits += nOfHits;
174 
175  } //# clusters back side
176 
177  } //# clusters front side
178 
179  return nHits;
180 }
181 // -------------------------------------------------------------------------
182 
183 
184 // ----- Get cluster position at read-out edge -------------------------
186  Double_t& xCluster,
187  Int_t& side) {
188 
189  // Take integer channel
190  Int_t iChannel = Int_t(centre);
191  Double_t xDif = centre - Double_t(iChannel);
192 
193  // Calculate corresponding strip on sensor
194  Int_t iStrip = -1;
195  pair<Int_t, Int_t> stripSide = GetStrip(iChannel);
196  iStrip = stripSide.first;
197  side = stripSide.second;
198 
199  // Re-add difference to integer channel. Convert channel to
200  // coordinate
201  xCluster = (Double_t(iStrip) + xDif + 0.5) * fPitch;
202 
203  // Correct for Lorentz-Shift
204  // Simplification: The correction uses only the y component of the
205  // magnetic field. The shift is calculated using the mid-plane of the
206  // sensor, which is not correct for tracks not traversing the entire
207  // sensor thickness (i.e., are created or stopped somewhere in the sensor).
208  // However, this is the best one can do in reconstruction.
209  if (side == 0)
210  xCluster -= fLorentzF;
211  else
212  xCluster -= fLorentzB;
213 
214  return;
215 }
216 // -------------------------------------------------------------------------
217 
218 
219 // ----- Get strip and side from channel number ------------------------
220 pair<Int_t, Int_t> CbmStsAlgoFindHits::GetStrip(UInt_t channel) const {
221 
222  Int_t stripNr = -1;
223  Int_t side = -1;
224 
225  // --- Determine front or back side
226  if (channel < fNofStrips) { // front side
227  side = 0;
228  stripNr = channel;
229  } else { // back side
230  side = 1;
231  stripNr = channel - fNofStrips;
232  }
233 
234  // --- Horizontal cross-connection
235  while (stripNr < 0)
236  stripNr += fNofStrips;
237  while (stripNr >= Int_t(fNofStrips))
238  stripNr -= fNofStrips;
239 
240  return (pair<Int_t, Int_t>(stripNr, side));
241 }
242 // -------------------------------------------------------------------------
243 
244 
245 // ----- Intersection of two lines along the strips --------------------
246 Bool_t CbmStsAlgoFindHits::Intersect(Double_t xF,
247  Double_t exF,
248  Double_t xB,
249  Double_t exB,
250  Double_t& x,
251  Double_t& y,
252  Double_t& varX,
253  Double_t& varY,
254  Double_t& varXY) {
255 
256  // In the coordinate system with origin at the bottom left corner,
257  // a line along the strips with coordinate x0 at the top edge is
258  // given by the function y(x) = Dy - ( x - x0 ) / tan(phi), if
259  // the stereo angle phi does not vanish. Two lines yF(x), yB(x) with top
260  // edge coordinates xF, xB intersect at
261  // x = ( tan(phiB)*xF - tan(phiF)*xB ) / (tan(phiB) - tan(phiF)
262  // y = Dy + ( xB - xF ) / ( tan(phiB) - tan(phiF) )
263  // For the case that one of the stereo angles vanish (vertical strips),
264  // the calculation of the intersection is straightforward.
265 
266  // --- First check whether stereo angles are different. Else there is
267  // --- no intersection.
268  if (TMath::Abs(fStereoF - fStereoB) < 0.5) {
269  x = -1000.;
270  y = -1000.;
271  return kFALSE;
272  }
273 
274  // --- Now treat vertical front strips
275  if (TMath::Abs(fStereoF) < 0.001) {
276  x = xF;
277  y = fDy - (xF - xB) / fTanStereoB;
278  varX = exF * exF;
279  varY = (exF * exF + exB * exB) / fTanStereoB / fTanStereoB;
280  varXY = -1. * exF * exF / fTanStereoB;
281  return IsInside(x - fDx / 2., y - fDy / 2.);
282  }
283 
284  // --- Maybe the back side has vertical strips?
285  if (TMath::Abs(fStereoB) < 0.001) {
286  x = xB;
287  y = fDy - (xB - xF) / fTanStereoF;
288  varX = exB * exB;
289  varY = (exF * exF + exB * exB) / fTanStereoF / fTanStereoF;
290  varXY = -1. * exB * exB / fTanStereoF;
291  return IsInside(x - fDx / 2., y - fDy / 2.);
292  }
293 
294  // --- OK, both sides have stereo angle
295  x = (fTanStereoB * xF - fTanStereoF * xB) / (fTanStereoB - fTanStereoF);
296  y = fDy + (xB - xF) / (fTanStereoB - fTanStereoF);
297  varX = fErrorFac
298  * (exF * exF * fTanStereoB * fTanStereoB
299  + exB * exB * fTanStereoF * fTanStereoF);
300  varY = fErrorFac * (exF * exF + exB * exB);
301  varXY = -1. * fErrorFac * (exF * exF * fTanStereoB + exB * exB * fTanStereoF);
302 
303  // --- Check for being in active area.
304  return IsInside(x - fDx / 2., y - fDy / 2.);
305 }
306 // -------------------------------------------------------------------------
307 
308 
309 // ----- Intersect two clusters ----------------------------------------
311  const CbmStsCluster& clusterB,
312  UInt_t indexF,
313  UInt_t indexB) {
314 
315  // --- Calculate cluster centre position on readout edge
316  Int_t side = -1;
317  Double_t xF = -1.;
318  Double_t xB = -1.;
319  GetClusterPosition(clusterF.GetPosition(), xF, side);
320  //std::cout << "Cluster position " << clusterF.GetPosition() << ": x "
321  // << xF << " side " << side;
322  if (side != 0) {
323  std::cout << "Cluster position " << clusterF.GetPosition() << ": x " << xF
324  << " side " << side << std::endl;
325  }
326  assert(side == 0);
327  Double_t exF = clusterF.GetPositionError() * fPitch;
328  Double_t du = exF * TMath::Cos(TMath::DegToRad() * fStereoF);
329  GetClusterPosition(clusterB.GetPosition(), xB, side);
330  assert(side == 1);
331  Double_t exB = clusterB.GetPositionError() * fPitch;
332  Double_t dv = exB * TMath::Cos(TMath::DegToRad() * fStereoB);
333 
334  // --- Should be inside active area
335  if (!(xF >= 0. || xF <= fDx)) return 0;
336  if (!(xB >= 0. || xB <= fDx)) return 0;
337 
338  // --- Hit counter
339  Int_t nHits = 0;
340 
341  // --- Calculate number of line segments due to horizontal
342  // --- cross-connection. If x(y=0) does not fall on the bottom edge,
343  // --- the strip is connected to the one corresponding to the line
344  // --- with top edge coordinate xF' = xF +/- Dx. For odd combinations
345  // --- of stereo angle and sensor dimensions, this could even happen
346  // --- multiple times. For each of these lines, the intersection with
347  // --- the line on the other side is calculated. If inside the active area,
348  // --- a hit is created.
349  Int_t nF = Int_t((xF + fDy * fTanStereoF) / fDx);
350  Int_t nB = Int_t((xB + fDy * fTanStereoB) / fDx);
351 
352  // --- If n is positive, all lines from 0 to n must be considered,
353  // --- if it is negative (phi negative), all lines from n to 0.
354  Int_t nF1 = TMath::Min(0, nF);
355  Int_t nF2 = TMath::Max(0, nF);
356  Int_t nB1 = TMath::Min(0, nB);
357  Int_t nB2 = TMath::Max(0, nB);
358 
359  // --- Double loop over possible lines
360  Double_t xC = -1.; // x coordinate of intersection point
361  Double_t yC = -1.; // y coordinate of intersection point
362  Double_t varX = 0.; // variance of xC
363  Double_t varY = 0.; // variance of yC
364  Double_t varXY = 0.; // covariance xC-yC
365  for (Int_t iF = nF1; iF <= nF2; iF++) {
366  Double_t xFi = xF - Double_t(iF) * fDx;
367  for (Int_t iB = nB1; iB <= nB2; iB++) {
368  Double_t xBi = xB - Double_t(iB) * fDx;
369 
370  // --- Intersect the two lines
371  Bool_t found = Intersect(xFi, exF, xBi, exB, xC, yC, varX, varY, varXY);
372  if (found) {
373 
374  // --- Transform into sensor system with origin at sensor centre
375  xC -= 0.5 * fDx;
376  yC -= 0.5 * fDy;
377  // --- Create the hit
378  CreateHit(xC,
379  yC,
380  varX,
381  varY,
382  varXY,
383  clusterF,
384  clusterB,
385  indexF,
386  indexB,
387  du,
388  dv);
389  nHits++;
390 
391  } //? Intersection of lines
392  } // lines on back side
393  } // lines on front side
394 
395  return nHits;
396 
397  return 0;
398 }
399 // -------------------------------------------------------------------------
400 
401 
402 // ----- Check whether a point is inside the active area ---------------
403 Bool_t CbmStsAlgoFindHits::IsInside(Double_t x, Double_t y) {
404  if (x < -fDx / 2.) return kFALSE;
405  if (x > fDx / 2.) return kFALSE;
406  if (y < -fDy / 2.) return kFALSE;
407  if (y > fDy / 2.) return kFALSE;
408  return kTRUE;
409 }
410 // -------------------------------------------------------------------------
411 
412 
CbmStsAlgoFindHits::fTanStereoF
Double_t fTanStereoF
///< Transformation matrix to global C.S.
Definition: CbmStsAlgoFindHits.h:215
CbmStsAlgoFindHits::fDx
Double_t fDx
Active size in x [cm].
Definition: CbmStsAlgoFindHits.h:207
CbmStsAlgoFindHits::IsInside
Bool_t IsInside(Double_t x, Double_t y)
Check whether a point (x,y) is inside the active area.
Definition: CbmStsAlgoFindHits.cxx:403
CbmStsCluster::GetPosition
Double_t GetPosition() const
Cluster position @value Cluster position in channel number units.
Definition: CbmStsCluster.h:67
CbmStsCluster
Data class for STS clusters.
Definition: CbmStsCluster.h:31
CbmStsAlgoFindHits::CreateHit
void CreateHit(Double_t xLocal, Double_t yLocal, Double_t varX, Double_t varY, Double_t varXY, const CbmStsCluster &clusterF, const CbmStsCluster &clusterB, UInt_t indexF, UInt_t indexB, Double_t du=0., Double_t dv=0.)
Create a new hit in the output array.
Definition: CbmStsAlgoFindHits.cxx:29
CbmStsAlgoFindHits::fDy
Double_t fDy
Active size in y [cm].
Definition: CbmStsAlgoFindHits.h:208
CbmStsAlgoFindHits::fPitch
Double_t fPitch
Strip pitch [cm].
Definition: CbmStsAlgoFindHits.h:209
CbmStsAlgoFindHits::fStereoF
Double_t fStereoF
Stereo angle front side [deg].
Definition: CbmStsAlgoFindHits.h:210
CbmStsAlgoFindHits::fLorentzB
Double_t fLorentzB
Lorentz shift correction back side [cm].
Definition: CbmStsAlgoFindHits.h:213
CbmStsParSensor.h
CbmStsAlgoFindHits::fTanStereoB
Double_t fTanStereoB
Tangent of stereo angle back side.
Definition: CbmStsAlgoFindHits.h:216
CbmStsParSensorCond.h
CbmStsAlgoFindHits::IntersectClusters
Int_t IntersectClusters(const CbmStsCluster &clusterF, const CbmStsCluster &clusterB, UInt_t indexF, UInt_t indexB)
Find the intersection points of two clusters.
Definition: CbmStsAlgoFindHits.cxx:310
CbmStsAlgoFindHits::GetStrip
std::pair< Int_t, Int_t > GetStrip(UInt_t channel) const
Get strip and side from module channel.
Definition: CbmStsAlgoFindHits.cxx:220
CbmStsCluster::GetTime
Double_t GetTime() const
Get cluster time.
Definition: CbmStsCluster.h:90
CbmStsDigi.h
CbmStsAlgoFindHits
Algorithm for hit finding in the sensors of the CBM-STS.
Definition: CbmStsAlgoFindHits.h:34
CbmStsAlgoFindHits::fNofStrips
UInt_t fNofStrips
Number of strips.
Definition: CbmStsAlgoFindHits.h:206
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmStsAlgoFindHits::fMatrix
TGeoHMatrix * fMatrix
Definition: CbmStsAlgoFindHits.h:214
CbmStsAlgoFindHits::fErrorFac
Double_t fErrorFac
For calculation of the hit error.
Definition: CbmStsAlgoFindHits.h:217
CbmStsCluster::GetPositionError
Double_t GetPositionError() const
Cluster position error @value Error (r.m.s.) of cluster position in channel number units.
Definition: CbmStsCluster.h:73
CbmDigiManager.h
CbmStsAlgoFindHits::fLorentzF
Double_t fLorentzF
Lorentz shift correction front side [cm].
Definition: CbmStsAlgoFindHits.h:212
CbmStsAlgoFindHits::fStereoB
Double_t fStereoB
Stereo angle back side [deg].
Definition: CbmStsAlgoFindHits.h:211
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
CbmStsAlgoFindHits::Intersect
Bool_t Intersect(Double_t xF, Double_t exF, Double_t xB, Double_t exB, Double_t &x, Double_t &y, Double_t &varX, Double_t &varY, Double_t &varXY)
Intersection point of two strips / cluster centres.
Definition: CbmStsAlgoFindHits.cxx:246
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
hits
static vector< vector< QAHit > > hits
Definition: CbmTofHitFinderTBQA.cxx:114
CbmStsAlgoFindHits::fHits
std::vector< CbmStsHit > * fHits
Output vector of hits.
Definition: CbmStsAlgoFindHits.h:220
CbmStsAlgoFindHits::GetClusterPosition
void GetClusterPosition(Double_t ClusterCentre, Double_t &xCluster, Int_t &side)
Definition: CbmStsAlgoFindHits.cxx:185
CbmStsAlgoFindHits.h
CbmStsCluster::GetTimeError
Double_t GetTimeError() const
Get error of cluster time.
Definition: CbmStsCluster.h:96
CbmStsCluster.h
Data class for STS clusters.
CbmStsAlgoFindHits::fAddress
UInt_t fAddress
Unique address for hits (sensor)
Definition: CbmStsAlgoFindHits.h:203
CbmStsAlgoFindHits::CbmStsAlgoFindHits
CbmStsAlgoFindHits()
Constructor.
Definition: CbmStsAlgoFindHits.cxx:24
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsAlgoFindHits::Exec
Long64_t Exec(const std::vector< CbmStsCluster > &clustersF, const std::vector< CbmStsCluster > &clustersB, std::vector< CbmStsHit > &hits, UInt_t address, Double_t timeCutSig, Double_t timeCutAbs, Double_t dY, UInt_t nStrips, Double_t pitch, Double_t stereoF, Double_t stereoB, Double_t lorentzF, Double_t lorentzB, TGeoHMatrix *matrix)
Execute algorithm.
Definition: CbmStsAlgoFindHits.cxx:85