CbmRoot
CbmStsAlgoFindHitsOrtho.cxx
Go to the documentation of this file.
1 
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 ----------------------------------------------
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 CbmStsAlgoFindHitsOrtho::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  UInt_t nStripsF,
92  UInt_t nStripsB,
93  Double_t pitchF,
94  Double_t pitchB,
95  Double_t lorentzF,
96  Double_t lorentzB,
97  TGeoHMatrix* matrix) {
98 
99  // Assert validity of parameters
100  assert(nStripsF > 0);
101  assert(nStripsB > 0);
102  assert(pitchF > 0.);
103  assert(pitchB > 0.);
104 
105  // Set internal parameters
106  fAddress = address;
107  fNofStripsF = nStripsF;
108  fNofStripsB = nStripsB;
109  fPitchF = pitchF;
110  fPitchB = pitchB;
111  fLorentzF = lorentzF;
112  fLorentzB = lorentzB;
113  fMatrix = matrix;
114  fHits = &hits;
115  fHits->clear();
116  fDx = Double_t(fNofStripsF) * fPitchF;
117  fDy = Double_t(fNofStripsB) * fPitchB;
118 
119  // Determine the maximum cluster time errors
120  Double_t maxTerrF = 0.;
121  for (auto& cluster : clustersF) {
122  if (cluster.GetTimeError() > maxTerrF) maxTerrF = cluster.GetTimeError();
123  }
124  Double_t maxTerrB = 0.;
125  for (auto& cluster : clustersB) {
126  if (cluster.GetTimeError() > maxTerrB) maxTerrB = cluster.GetTimeError();
127  }
128  const Double_t maxTerrF2 = maxTerrF * maxTerrF;
129  const Double_t maxTerrB2 = maxTerrB * maxTerrB;
130  const Double_t max_sigma_both = 4. * TMath::Sqrt(maxTerrF2 + maxTerrB2);
131 
132  // --- Loop over front and back side clusters
133  Long64_t nHits = 0;
134  Int_t startB = 0;
135  for (UInt_t iClusterF = 0; iClusterF < clustersF.size(); iClusterF++) {
136  const CbmStsCluster& clusterF = clustersF[iClusterF];
137 
138  Double_t tF = clusterF.GetTime();
139  Double_t tFerr2 = clusterF.GetTimeError() * clusterF.GetTimeError();
140  Double_t max_sigma = 4. * TMath::Sqrt(tFerr2 + maxTerrB2);
141 
142  for (UInt_t iClusterB = startB; iClusterB < clustersB.size(); iClusterB++) {
143  const CbmStsCluster& clusterB = clustersB[iClusterB];
144 
145  Double_t timeDiff = tF - clusterB.GetTime();
146  if ((timeDiff > 0) && (timeDiff > max_sigma_both)) {
147  startB++;
148  continue;
149  } else if ((timeDiff > 0) && (timeDiff > max_sigma)) {
150  continue;
151  } else if ((timeDiff < 0) && (fabs(timeDiff) > max_sigma))
152  break;
153 
154  // Cut on time difference of the two clusters
155  Double_t timeCut = -1.;
156  if (timeCutAbs > 0.)
157  timeCut = timeCutAbs; // absolute cut value
158  else {
159  if (timeCutSig > 0.) {
160  Double_t eF = clusterF.GetTimeError();
161  Double_t eB = clusterB.GetTimeError();
162  timeCut = timeCutSig * TMath::Sqrt(eF * eF + eB * eB);
163  } //? cut calculated from cluster errors
164  }
165  if (fabs(clusterF.GetTime() - clusterB.GetTime()) > timeCut) continue;
166 
167  // --- Calculate intersection points
168  Int_t nOfHits =
169  IntersectClusters(clusterF, clusterB, iClusterF, iClusterB);
170  nHits += nOfHits;
171 
172  } //# clusters back side
173 
174  } //# clusters front side
175 
176  return nHits;
177 }
178 // -------------------------------------------------------------------------
179 
180 
181 // ----- Get cluster position at read-out edge -------------------------
183  Double_t& xCluster,
184  Int_t& side) {
185 
186  // Take integer channel
187  Int_t iChannel = Int_t(centre);
188  Double_t xDif = centre - Double_t(iChannel);
189 
190  // Calculate corresponding strip on sensor
191  Int_t iStrip = -1;
192  pair<Int_t, Int_t> stripSide = GetStrip(iChannel);
193  iStrip = stripSide.first;
194  side = stripSide.second;
195  assert(side == 0 || side == 1);
196 
197  // Re-add difference to integer channel. Convert channel to
198  // coordinate
199  if (side == 0)
200  xCluster = (Double_t(iStrip) + xDif + 0.5) * fPitchF;
201  else
202  xCluster = (Double_t(iStrip) + xDif + 0.5) * fPitchB;
203 
204  // Correct for Lorentz-Shift
205  // Simplification: The correction uses only the y component of the
206  // magnetic field. The shift is calculated using the mid-plane of the
207  // sensor, which is not correct for tracks not traversing the entire
208  // sensor thickness (i.e., are created or stopped somewhere in the sensor).
209  // However, this is the best one can do in reconstruction.
210  if (side == 0)
211  xCluster -= fLorentzF;
212  else
213  xCluster -= fLorentzB;
214 
215  return;
216 }
217 // -------------------------------------------------------------------------
218 
219 
220 // ----- Get strip and side from channel number ------------------------
221 pair<Int_t, Int_t> CbmStsAlgoFindHitsOrtho::GetStrip(UInt_t channel) const {
222 
223  Int_t stripNr = -1;
224  Int_t side = -1;
225 
226  // --- Determine front or back side
227  if (channel < fNofStripsF) { // front side
228  side = 0;
229  stripNr = channel;
230  } else { // back side
231  side = 1;
232  stripNr = channel - fNofStripsF;
233  }
234 
235  return (pair<Int_t, Int_t>(stripNr, side));
236 }
237 // -------------------------------------------------------------------------
238 
239 
240 // ----- Intersect two clusters ----------------------------------------
242  const CbmStsCluster& clusterB,
243  UInt_t indexF,
244  UInt_t indexB) {
245 
246  // --- Calculate cluster centre position on readout edge
247  Int_t side = -1;
248  Double_t xF = -1.;
249  Double_t xB = -1.;
250  GetClusterPosition(clusterF.GetPosition(), xF, side);
251  assert(side == 0);
252  Double_t exF = clusterF.GetPositionError() * fPitchF;
253  Double_t du = exF;
254  GetClusterPosition(clusterB.GetPosition(), xB, side);
255  assert(side == 1);
256  Double_t exB = clusterB.GetPositionError() * fPitchB;
257  Double_t dv = exB;
258 
259  // --- Should be inside active area
260  if (!(xF >= 0. || xF <= fDx)) return 0;
261  if (!(xB >= 0. || xB <= fDy)) return 0;
262 
263  // --- Hit counter
264  Int_t nHits = 0;
265 
266  // In orthogonal sensor, all pairs of (front, back) cluster have
267  // a single intersection
268  // => exactly one hit!
269 
270  // --- Prepare hit coordinates and errors
271  // In the coordinate system with origin at the bottom left corner,
272  // the coordinates in the orthogonal sensor are straightforward.
273  Double_t xC = xF; // x coordinate of intersection point
274  Double_t yC = xB; // y coordinate of intersection point
275  Double_t varX = exF * exF; // variance of xC
276  Double_t varY = exB * exB; // variance of yC
277  Double_t varXY = 0.; // covariance xC-yC => independent variables!
278 
279  // --- Transform into sensor system with origin at sensor centre
280  xC -= 0.5 * fDx;
281  yC -= 0.5 * fDy;
282 
283  // --- Create the hit
284  CreateHit(
285  xC, yC, varX, varY, varXY, clusterF, clusterB, indexF, indexB, du, dv);
286  nHits++;
287 
288  return nHits;
289 }
290 // -------------------------------------------------------------------------
291 
292 
293 // ----- Check whether a point is inside the active area ---------------
294 Bool_t CbmStsAlgoFindHitsOrtho::IsInside(Double_t x, Double_t y) {
295  if (x < -fDx / 2.) return kFALSE;
296  if (x > fDx / 2.) return kFALSE;
297  if (y < -fDy / 2.) return kFALSE;
298  if (y > fDy / 2.) return kFALSE;
299  return kTRUE;
300 }
301 // -------------------------------------------------------------------------
302 
303 
CbmStsAlgoFindHitsOrtho::GetClusterPosition
void GetClusterPosition(Double_t ClusterCentre, Double_t &xCluster, Int_t &side)
Definition: CbmStsAlgoFindHitsOrtho.cxx:182
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
CbmStsAlgoFindHitsOrtho::fDx
Double_t fDx
Active size in x [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:179
CbmStsParSensor.h
CbmStsAlgoFindHitsOrtho::fLorentzB
Double_t fLorentzB
Lorentz shift correction back side [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:184
CbmStsParSensorCond.h
CbmStsAlgoFindHitsOrtho::fHits
std::vector< CbmStsHit > * fHits
///< Transformation matrix to global C.S.
Definition: CbmStsAlgoFindHitsOrtho.h:188
CbmStsCluster::GetTime
Double_t GetTime() const
Get cluster time.
Definition: CbmStsCluster.h:90
CbmStsAlgoFindHitsOrtho::fAddress
UInt_t fAddress
Unique address for hits (sensor)
Definition: CbmStsAlgoFindHitsOrtho.h:174
CbmStsAlgoFindHitsOrtho::fNofStripsB
UInt_t fNofStripsB
Number of strips backs side.
Definition: CbmStsAlgoFindHitsOrtho.h:178
CbmStsDigi.h
CbmStsAlgoFindHitsOrtho::fMatrix
TGeoHMatrix * fMatrix
Definition: CbmStsAlgoFindHitsOrtho.h:185
CbmStsAlgoFindHitsOrtho::GetStrip
std::pair< Int_t, Int_t > GetStrip(UInt_t channel) const
Get strip and side from module channel.
Definition: CbmStsAlgoFindHitsOrtho.cxx:221
CbmStsAlgoFindHitsOrtho::fNofStripsF
UInt_t fNofStripsF
Number of strips front side.
Definition: CbmStsAlgoFindHitsOrtho.h:177
CbmStsAlgoFindHitsOrtho.h
CbmStsAlgoFindHitsOrtho::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: CbmStsAlgoFindHitsOrtho.cxx:29
CbmStsAlgoFindHitsOrtho::IntersectClusters
Int_t IntersectClusters(const CbmStsCluster &clusterF, const CbmStsCluster &clusterB, UInt_t indexF, UInt_t indexB)
Find the intersection points of two clusters.
Definition: CbmStsAlgoFindHitsOrtho.cxx:241
CbmStsAlgoFindHitsOrtho
Algorithm for hit finding in sensors with orthogonal strips.
Definition: CbmStsAlgoFindHitsOrtho.h:35
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmStsAlgoFindHitsOrtho::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, UInt_t nStripsF, UInt_t nStripsB, Double_t pitchF, Double_t pitchB, Double_t lorentzF, Double_t lorentzB, TGeoHMatrix *matrix)
Execute algorithm.
Definition: CbmStsAlgoFindHitsOrtho.cxx:85
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
CbmStsAlgoFindHitsOrtho::fDy
Double_t fDy
Active size in y [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:180
CbmStsAlgoFindHitsOrtho::fPitchB
Double_t fPitchB
Strip pitch back side [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:182
CbmDigiManager.h
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmStsAlgoFindHitsOrtho::fLorentzF
Double_t fLorentzF
Lorentz shift correction front side [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:183
hits
static vector< vector< QAHit > > hits
Definition: CbmTofHitFinderTBQA.cxx:114
CbmStsAlgoFindHitsOrtho::CbmStsAlgoFindHitsOrtho
CbmStsAlgoFindHitsOrtho()
Constructor.
Definition: CbmStsAlgoFindHitsOrtho.cxx:24
CbmStsCluster::GetTimeError
Double_t GetTimeError() const
Get error of cluster time.
Definition: CbmStsCluster.h:96
CbmStsCluster.h
Data class for STS clusters.
CbmStsAlgoFindHitsOrtho::fPitchF
Double_t fPitchF
Strip pitch front side [cm].
Definition: CbmStsAlgoFindHitsOrtho.h:181
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsAlgoFindHitsOrtho::IsInside
Bool_t IsInside(Double_t x, Double_t y)
Check whether a point (x,y) is inside the active area.
Definition: CbmStsAlgoFindHitsOrtho.cxx:294