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