BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
L1TrackExtender.cxx
Go to the documentation of this file.
1#include "L1AddMaterial.h"
2#include "L1Algo.h"
3#include "L1Branch.h"
4#include "L1Extrapolation.h"
5#include "L1Filtration.h"
6#include "L1HitArea.h"
7#include "L1HitPoint.h"
8#include "L1Track.h"
9#include "L1TrackPar.h"
10
11#include <iostream>
12
13// using namespace std;
14using std::cout;
15using std::endl;
16using std::vector;
17
18/// Fit track
19/// t - track with hits
20/// T - track params
21/// dir - 0 - forward, 1 - backward
22/// qp0 - momentum for extrapolation
23/// initialize - should be params ititialized. 1 - yes.
24void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0, const bool initParams)
25{
26 // cout << " L1Algo::BranchFitterFast " << endl; exit(1); //AZ-220223
27 L1_assert(t.StsHits.size() >= 3);
28
29 // get hits of current track
30 const std::vector<THitI>& hits = t.StsHits; // array of indeses of hits of current track
31 const int nHits = t.StsHits.size();
32
33 const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
34 const int iFirstHit = (dir) ? nHits - 1 : 0;
35 const int iLastHit = (dir) ? 0 : nHits - 1;
36
37 L1StsHit& hit0 = vStsHits[hits[iFirstHit]];
38 L1StsHit& hit1 = vStsHits[hits[iFirstHit + step]];
39 L1StsHit& hit2 = vStsHits[hits[iFirstHit + 2 * step]];
40
41 int ista0 = GetFStation(vSFlag[hit0.f]);
42 int ista1 = GetFStation(vSFlag[hit1.f]);
43 int ista2 = GetFStation(vSFlag[hit2.f]);
44
45 L1Station& sta0 = vStations[ista0];
46 L1Station& sta1 = vStations[ista1];
47 L1Station& sta2 = vStations[ista2];
48
49 fvec x0, y0;
50 fvec z0 = vStsZPos[hit0.iz];
51#ifdef AZ
52 // AZ-230223
53 int ihp = vStsHits2Unus[t.StsHits[iFirstHit]];
54 x0 = (*vStsHitPointsUnused)[ihp].X();
55 y0 = (*vStsHitPointsUnused)[ihp].Y();
56#else
57 fvec u0 = static_cast<fscal>(vStsStrips[hit0.f]);
58 fvec v0 = static_cast<fscal>(vStsStripsB[hit0.b]);
59 StripsToCoor(u0, v0, x0, y0, sta0);
60#endif
61
62 fvec x1, y1;
63 fvec z1 = vStsZPos[hit1.iz];
64#ifdef AZ
65 // AZ-230223
66 ihp = vStsHits2Unus[t.StsHits[iFirstHit + step]];
67 x1 = (*vStsHitPointsUnused)[ihp].X();
68 y1 = (*vStsHitPointsUnused)[ihp].Y();
69#else
70 fvec u1 = static_cast<fscal>(vStsStrips[hit1.f]);
71 fvec v1 = static_cast<fscal>(vStsStripsB[hit1.b]);
72 StripsToCoor(u1, v1, x1, y1, sta1);
73#endif
74
75 fvec x2, y2;
76#ifdef AZ
77 // AZ-230223
78 ihp = vStsHits2Unus[t.StsHits[iFirstHit + 2 * step]];
79 x2 = (*vStsHitPointsUnused)[ihp].X();
80 y2 = (*vStsHitPointsUnused)[ihp].Y();
81#else
82 fvec u2 = static_cast<fscal>(vStsStrips[hit2.f]);
83 fvec v2 = static_cast<fscal>(vStsStripsB[hit2.b]);
84 StripsToCoor(u2, v2, x2, y2, sta2);
85#endif
86 // fvec z2 = vStsZPos[hit2.iz];
87
88 fvec dzi = 1. / (z1 - z0);
89
90 const fvec vINF = .1;
91 T.x = x0;
92 T.y = y0;
93 if (initParams) {
94 T.tx = (x1 - x0) * dzi;
95 T.ty = (y1 - y0) * dzi;
96 T.qp = qp0;
97 }
98
99 T.z = z0;
100 T.chi2 = 0.;
101 T.NDF = 2.;
102 T.C00 = sta0.XYInfo.C00;
103 T.C10 = sta0.XYInfo.C10;
104 T.C11 = sta0.XYInfo.C11;
105
106 T.C20 = T.C21 = 0;
107 T.C30 = T.C31 = T.C32 = 0;
108 T.C40 = T.C41 = T.C42 = T.C43 = 0;
109 T.C22 = T.C33 = vINF;
110 T.C44 = 1.;
111
112 L1FieldValue fB0, fB1, fB2 _fvecalignment;
114 fvec fz0 = sta1.z; // suppose field is smoth
115 fvec fz1 = sta2.z;
116 fvec fz2 = sta0.z;
117
118 sta1.fieldSlice.GetFieldValue(x1, y1, fB0);
119 sta2.fieldSlice.GetFieldValue(x2, y2, fB1);
120 sta0.fieldSlice.GetFieldValue(x0, y0, fB2);
121
122 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
123
124 int ista_prev = ista1;
125 int ista = ista2;
126
127 for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
128 L1StsHit& hit = vStsHits[hits[i]];
129 ista_prev = ista;
130 ista = GetFStation(vSFlag[hit.f]);
131
132 L1Station& sta = vStations[ista];
133
134 L1Extrapolate(T, vStsZPos[hit.iz], qp0, fld);
135 L1AddMaterial(T, sta.materialInfo, qp0);
136 if ((step * ista <= step * (NMvdStations + (step + 1) / 2 - 1))
137 && (step * ista_prev >= step * (NMvdStations + (step + 1) / 2 - 1 - step)))
138 L1AddPipeMaterial(T, qp0);
139
140 fvec x, y;
141#ifdef AZ
142 // AZ-260223
143 ihp = vStsHits2Unus[t.StsHits[i]];
144 L1HitPoint& point = (*vStsHitPointsUnused)[ihp];
145 x = point.X();
146 y = point.Y();
147 fvec u = static_cast<fscal>(point.U());
148 fvec v = static_cast<fscal>(point.V());
149 fvec du = static_cast<fscal>(point.DU());
150 fvec dv = static_cast<fscal>(point.DV());
151#else
152 fvec u = static_cast<fscal>(vStsStrips[hit.f]);
153 fvec v = static_cast<fscal>(vStsStripsB[hit.b]);
154 fvec du = static_cast<fscal>(vStsStrips[hit.f].errf()); // AZ
155 fvec dv = static_cast<fscal>(vStsStripsB[hit.b].errf()); // AZ
156#endif
157 // L1Filter( T, sta.frontInfo, u );
158 // L1Filter( T, sta.backInfo, v );
159 fvec w = 1.0; // AZ
160 du *= du; // AZ
161 dv *= dv; // AZ
162 L1Filter(T, sta.frontInfo, u, w, &du); // AZ
163 L1Filter(T, sta.backInfo, v, w, &dv); // AZ
164 fB0 = fB1;
165 fB1 = fB2;
166 fz0 = fz1;
167 fz1 = fz2;
168 // AZ-260223 fvec x,y;
169#ifndef AZ
170 StripsToCoor(u, v, x, y, sta);
171#endif
172 sta.fieldSlice.GetFieldValue(x, y, fB2);
173
174 fz2 = sta.z;
175 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
176 } // i
177
178} // void L1Algo::BranchFitterFast
179
180/// like BranchFitterFast but more precise
181void L1Algo::BranchFitter(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0, const bool initParams)
182{
183 BranchFitterFast(t, T, dir, qp0, initParams);
184 for (int i = 0; i < 1; i++) {
185 BranchFitterFast(t, T, !dir, T.qp, false);
186 BranchFitterFast(t, T, dir, T.qp, false);
187 }
188} // void L1Algo::BranchFitter
189
190//----------------------------------------------------------------------------------------------
191
192/// Find additional hits for existing track
193/// t - track with hits
194/// T - track params
195/// dir - 0 - forward, 1 - backward
196/// qp0 - momentum for extrapolation
197/// initialize - should be params ititialized. 1 - yes.
198void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0) // TODO take into account pipe
199{
200 // cout << " L1Algo::FindMoreHits " << endl; exit(1); //AZ-220223
201 std::vector<THitI> newHits;
202 newHits.clear();
203
204 const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
205 const int iFirstHit = (dir) ? 2 : t.StsHits.size() - 3;
206 // int ista = GetFStation(vSFlag[vStsHits[t.StsHits[iFirstHit]].f]) + 2*step; // current station. set to the end of
207 // track
208
209 const L1StsHit& hit0 = vStsHits[t.StsHits[iFirstHit]]; // optimize
210 const L1StsHit& hit1 = vStsHits[t.StsHits[iFirstHit + step]];
211 const L1StsHit& hit2 = vStsHits[t.StsHits[iFirstHit + 2 * step]];
212
213 const int ista0 = GetFStation(vSFlag[hit0.f]);
214 const int ista1 = GetFStation(vSFlag[hit1.f]);
215 const int ista2 = GetFStation(vSFlag[hit2.f]);
216
217 const L1Station& sta0 = vStations[ista0];
218 const L1Station& sta1 = vStations[ista1];
219 const L1Station& sta2 = vStations[ista2];
220
221 fvec x0, y0;
222#ifdef AZ
223 // AZ-230223
224 int ihp = vStsHits2Unus[t.StsHits[iFirstHit]];
225 x0 = (*vStsHitPointsUnused)[ihp].X();
226 y0 = (*vStsHitPointsUnused)[ihp].Y();
227#else
228 fvec u0 = static_cast<fscal>(vStsStrips[hit0.f]);
229 fvec v0 = static_cast<fscal>(vStsStripsB[hit0.b]);
230 StripsToCoor(u0, v0, x0, y0, sta0);
231#endif
232 // fvec z0 = vStsZPos[hit0.iz];
233
234 fvec x1, y1;
235#ifdef AZ
236 // AZ-230223
237 ihp = vStsHits2Unus[t.StsHits[iFirstHit + step]];
238 x1 = (*vStsHitPointsUnused)[ihp].X();
239 y1 = (*vStsHitPointsUnused)[ihp].Y();
240#else
241 fvec u1 = static_cast<fscal>(vStsStrips[hit1.f]);
242 fvec v1 = static_cast<fscal>(vStsStripsB[hit1.b]);
243 StripsToCoor(u1, v1, x1, y1, sta1);
244#endif
245 // fvec z1 = vStsZPos[hit1.iz];
246
247 fvec x2, y2;
248#ifdef AZ
249 // AZ-230223
250 ihp = vStsHits2Unus[t.StsHits[iFirstHit + 2 * step]];
251 x2 = (*vStsHitPointsUnused)[ihp].X();
252 y2 = (*vStsHitPointsUnused)[ihp].Y();
253#else
254 fvec u2 = static_cast<fscal>(vStsStrips[hit2.f]);
255 fvec v2 = static_cast<fscal>(vStsStripsB[hit2.b]);
256 StripsToCoor(u2, v2, x2, y2, sta2);
257#endif
258 // fvec z2 = vStsZPos[hit2.iz];
259
260 // fvec dzi = 1./(z1-z0);
261
262 L1FieldValue fB0, fB1, fB2 _fvecalignment;
264 fvec fz0 = sta1.z;
265 fvec fz1 = sta2.z;
266 fvec fz2 = sta0.z;
267
268 sta1.fieldSlice.GetFieldValue(x1, y1, fB0);
269 sta2.fieldSlice.GetFieldValue(x2, y2, fB1);
270 sta0.fieldSlice.GetFieldValue(x0, y0, fB2);
271
272 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
273
274 int ista = ista2 + 2 * step; // skip one station. if there would be hit it has to be found on previous stap
275 if (ista2 == FIRSTCASTATION)
276 ista = ista2 + step;
277
278 const fvec Pick_gather2 = Pick_gather * Pick_gather;
279 for (; (ista < NStations) && (ista >= 0); ista += step) { // CHECKME why ista2?
280
281 L1Station& sta = vStations[ista];
282
283 L1Extrapolate(T, sta.z, qp0, fld);
284
285 fscal r2_best = 1e8; // best distance to hit
286 int iHit_best = -1; // index of the best hit
287
288 const fscal iz = 1 / T.z[0];
289 L1HitArea area(vGrid[ista], T.x[0] * iz, T.y[0] * iz,
290 (sqrt(Pick_gather * (T.C00 + sta.XYInfo.C00)) + MaxDZ * fabs(T.tx))[0] * iz,
291 (sqrt(Pick_gather * (T.C11 + sta.XYInfo.C11)) + MaxDZ * fabs(T.ty))[0] * iz);
292 THitI ih = 0;
293 while (area.GetNext(ih)) {
294 ih += StsHitsUnusedStartIndex[ista];
295 L1StsHit& hit = (*vStsHitsUnused)[ih];
296
297 if (GetFUsed(vSFlag[hit.f] | vSFlagB[hit.b]))
298 continue; // if used
299
300 fscal xh, yh, zh;
301#ifdef AZ
302 // AZ-230223
303 L1HitPoint& point = (*vStsHitPointsUnused)[ih];
304 xh = point.X(), yh = point.Y(), zh = point.Z();
305#else
306 GetHitCoor(hit, xh, yh, zh, sta); // faster
307 // L1HitPoint &point = (*vStsHitPointsUnused)[ih];
308 // fscal xh = point.x, yh = point.y, zh = point.z;
309#endif
310 fvec y, C11;
311 L1ExtrapolateYC11Line(T, zh, y, C11);
312
313 fscal dym_est = (Pick_gather * sqrt(fabs(C11 + sta.XYInfo.C11)))[0];
314 fscal y_minus_new = y[0] - dym_est;
315 if (yh < y_minus_new)
316 continue; // CHECKME take into account overlaping?
317
318 fvec x, C00;
319 L1ExtrapolateXC00Line(T, zh, x, C00);
320
321 fscal dx = xh - x[0];
322 fscal dy = yh - y[0];
323 fscal d2 = dx * dx + dy * dy;
324 if (d2 > r2_best)
325 continue;
326
327 fscal dxm_est2 = (Pick_gather2 * (fabs(C00 + sta.XYInfo.C00)))[0];
328 if (dx * dx > dxm_est2)
329 continue;
330
331 r2_best = d2;
332 iHit_best = ih;
333 }
334 if (iHit_best < 0)
335 break;
336
337 newHits.push_back(RealIHit[iHit_best]);
338
339 fvec x, y, z;
340 L1StsHit& hit = (*vStsHitsUnused)[iHit_best];
341#ifdef AZ
342 // AZ-260223
343 L1HitPoint& point = (*vStsHitPointsUnused)[iHit_best];
344 x = point.X(), y = point.Y();
345 fvec u = static_cast<fscal>(point.U());
346 fvec v = static_cast<fscal>(point.V());
347 fvec du = static_cast<fscal>(point.DU());
348 fvec dv = static_cast<fscal>(point.DV());
349#else
350 fvec u = static_cast<fvec>(vStsStrips[hit.f]);
351 fvec v = static_cast<fvec>(vStsStripsB[hit.b]);
352 fvec du = static_cast<fvec>(vStsStrips[hit.f].errf()); // AZ
353 fvec dv = static_cast<fvec>(vStsStripsB[hit.b].errf()); // AZ
354 StripsToCoor(u, v, x, y, sta);
355#endif
356 z = vStsZPos[hit.iz];
357
358 L1ExtrapolateLine(T, z);
359 L1AddMaterial(T, sta.materialInfo, qp0);
360 // L1Filter( T, sta.frontInfo, u );
361 // L1Filter( T, sta.backInfo, v );
362 fvec w = 1.0; // AZ
363 du *= du; // AZ
364 dv *= dv; // AZ
365 L1Filter(T, sta.frontInfo, u, w, &du); // AZ
366 L1Filter(T, sta.backInfo, v, w, &dv); // AZ
367
368 fB0 = fB1;
369 fB1 = fB2;
370 fz0 = fz1;
371 fz1 = fz2;
372 sta.fieldSlice.GetFieldValue(x, y, fB2);
373 fz2 = sta.z;
374 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
375 }
376
377 // save hits
378 if (dir) { // backward
379 const unsigned int NOldHits = t.StsHits.size();
380 const unsigned int NNewHits = newHits.size();
381 t.StsHits.resize(NNewHits + NOldHits);
382 for (int i = NOldHits - 1; i >= 0; i--) {
383 t.StsHits[NNewHits + i] = t.StsHits[i];
384 }
385 for (unsigned int i = 0, ii = NNewHits - 1; i < NNewHits; i++, ii--) {
386 t.StsHits[i] = newHits[ii];
387 }
388 } else { // forward
389 const unsigned int NOldHits = t.StsHits.size();
390 t.StsHits.resize(newHits.size() + NOldHits);
391 for (unsigned int i = 0; i < newHits.size(); i++) {
392 t.StsHits[NOldHits + i] = (newHits[i]);
393 }
394 }
395
396} // void L1Algo::FindMoreHits
397
398//----------------------------------------------------------------------------------------------
399
400/// Try to extrapolate and find additional hits on other stations
402{
403 // const unsigned int minNHits = 3;
404
405 L1TrackPar T;
406
407 // forward
408 bool dir = 0;
409
410 BranchFitter(t, T, dir);
411 // BranchFitterFast (t, T, dir, 0, true);
412
413 // if (t.StsHits.size() < minNHits) return T.chi2[0];
414 FindMoreHits(t, T, dir, T.qp);
415
416 // backward
417 dir = 1;
418 BranchFitterFast(t, T, dir, T.qp, false); // 577
419
420 FindMoreHits(t, T, dir, T.qp);
421
422 return T.chi2[0];
423}
#define L1_assert(v)
Definition CbmL1Def.h:49
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
void L1ExtrapolateLine(L1TrackPar &T, fvec z_out)
void L1ExtrapolateYC11Line(const L1TrackPar &T, fvec z_out, fvec &y, fvec &C11)
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, L1FieldRegion &F, fvec *w=0)
void L1ExtrapolateXC00Line(const L1TrackPar &T, fvec z_out, fvec &x, fvec &C00)
void L1Filter(L1TrackPar &T, L1UMeasurementInfo &info, fvec u, fvec w=1., fvec *du2=0)
unsigned int THitI
Definition L1StsHit.h:6
float fscal
Definition P4_F32vec4.h:232
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
__m128 v
Definition P4_F32vec4.h:1
int i
Definition P4_F32vec4.h:22
F32vec4 fvec
Definition P4_F32vec4.h:231
void StripsToCoor(const fscal &u, const fscal &v, fscal &x, fscal &y, const L1Station &sta) const
convert strip positions to coordinates
Definition L1Algo.cxx:185
vector< L1StsHit > vStsHits
Definition L1Algo.h:162
THitI * RealIHit
Definition L1Algo.h:174
L1Station vStations[MaxNStations] _fvecalignment
Definition L1Algo.h:155
int NStations
Definition L1Algo.h:153
vector< fscal > vStsZPos
Definition L1Algo.h:161
vector< L1Strip > vStsStripsB
Definition L1Algo.h:160
unsigned char GetFStation(unsigned char flag)
– Flags routines –
Definition L1Algo.h:271
vector< unsigned char > vSFlagB
Definition L1Algo.h:165
THitI StsHitsUnusedStartIndex[MaxNStations+1]
Definition L1Algo.h:175
fvec MaxDZ
Definition L1Algo.h:526
bool GetFUsed(unsigned char flag)
Definition L1Algo.h:272
vector< unsigned char > vSFlag
Definition L1Algo.h:164
fscal BranchExtender(L1Branch &t)
Find additional hits for existing track return chi2.
int NMvdStations
Definition L1Algo.h:154
void BranchFitterFast(const L1Branch &t, L1TrackPar &T, const bool dir, const fvec qp0=0., const bool initParams=true)
Fit track t - track with hits T - track params dir - 0 - forward, 1 - backward qp0 - momentum for ext...
fvec Pick_gather
parameters which are different for different iterations. Set in the begin of CAL1TrackFinder
Definition L1Algo.h:530
vector< L1Strip > vStsStrips
Definition L1Algo.h:159
vector< int > vStsHits2Unus
Definition L1Algo.h:177
int FIRSTCASTATION
================================= DATA PART =================================
Definition L1Algo.h:476
void BranchFitter(const L1Branch &t, L1TrackPar &T, const bool dir, const fvec qp0=0., const bool initParams=true)
Fit track. more precise than FitterFast.
void FindMoreHits(L1Branch &t, L1TrackPar &T, const bool dir, const fvec qp0=0.0)
Find additional hits for existing track t - track with hits T - track params dir - 0 - forward,...
L1Grid vGrid[MaxNStations]
Definition L1Algo.h:163
void GetHitCoor(const L1StsHit &_h, fscal &_x, fscal &_y, fscal &_z, const L1Station &sta)
--— Hit-point-strips conversion routines ---—
Definition L1Algo.cxx:173
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition L1Field.h:47
bool GetNext(THitI &i)
look up the next hit in the requested area.
Definition L1HitArea.h:86
L1XYMeasurementInfo XYInfo
Definition L1Station.h:31
L1UMeasurementInfo backInfo
Definition L1Station.h:30
L1FieldSlice fieldSlice
Definition L1Station.h:29
L1MaterialInfo materialInfo
Definition L1Station.h:28
L1UMeasurementInfo frontInfo
Definition L1Station.h:30
TZPosI iz
Definition L1StsHit.h:13
TStripI f
Definition L1StsHit.h:12
TStripI b
Definition L1StsHit.h:12
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
std::vector< THitI > StsHits
Definition L1Branch.h:39
contain strips positions and coordinates of hit
Definition L1HitPoint.h:9
fscal V() const
Definition L1HitPoint.h:43
fscal DU() const
Definition L1HitPoint.h:44
fscal U() const
Definition L1HitPoint.h:42
fscal Z() const
Definition L1HitPoint.h:41
fscal Y() const
Definition L1HitPoint.h:40
fscal X() const
Definition L1HitPoint.h:39
fscal DV() const
Definition L1HitPoint.h:45