Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFPTrackVector.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFPTrackVector.cxx
1 /*
2  * This file is part of KFParticle package
3  * Copyright (C) 2007-2019 FIAS Frankfurt Institute for Advanced Studies
4  * 2007-2019 Goethe University of Frankfurt
5  * 2007-2019 Ivan Kisel <I.Kisel@compeng.uni-frankfurt.de>
6  * 2007-2019 Maksym Zyzak
7  *
8  * KFParticle is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * KFParticle is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program. If not, see <https://www.gnu.org/licenses/>.
20  */
21 
22 #include "KFPTrackVector.h"
23 #include <iostream>
24 
25 void KFPTrackVector::SetParameter(const float_v& value, int iP, int iTr)
26 {
33 // gather caused errors at XeonPhi, temporarly replaced with the simple copying
34 // if( (iTr+float_vLen) < Size())
35 // reinterpret_cast<float_v&>(fP[iP][iTr]) = value;
36 // else
37 // {
38 // const uint_v index(uint_v::IndexesFromZero());
39 // (reinterpret_cast<float_v&>(fP[iP][iTr])).gather(reinterpret_cast<const float*>(&value), index, float_m(index<(Size() - iTr)));
40 // }
41 
42  if( (iTr+float_vLen) < Size())
43  reinterpret_cast<float_v&>(fP[iP][iTr]) = value;
44  else
45  for(int i=0; i<float_v::Size; i++)
46  {
47  if(iTr + i >= Size()) continue;
48  fP[iP][iTr+i] = value[i];
49  }
50 }
51 void KFPTrackVector::SetCovariance(const float_v& value, int iC, int iTr)
52 {
59 // gather caused errors at XeonPhi, temporarly replaced with the simple copying
60 // if( (iTr+float_vLen) < Size())
61 // reinterpret_cast<float_v&>(fC[iC][iTr]) = value;
62 // else
63 // {
64 // const uint_v index(uint_v::IndexesFromZero());
65 // (reinterpret_cast<float_v&>(fC[iC][iTr])).gather(reinterpret_cast<const float*>(&value), index, float_m(index<(Size() - iTr)));
66 // }
67 
68  if( (iTr+float_vLen) < Size())
69  reinterpret_cast<float_v&>(fC[iC][iTr]) = value;
70  else
71  for(int i=0; i<float_v::Size; i++)
72  {
73  if(iTr + i >= Size()) continue;
74  fC[iC][iTr+i] = value[i];
75  }
76 }
77 
78 
79 void KFPTrackVector::Resize(const int n)
80 {
84  for(int i=0; i<6; i++)
85  fP[i].resize(n);
86  for(int i=0; i<21; i++)
87  fC[i].resize(n);
88 #ifdef NonhomogeneousField
89  for(int i=0; i<10; i++)
90  fField[i].resize(n);
91 #endif
92 // fChi2.resize(n);
93 // fNDF.resize(n);
94  fId.resize(n);
95  fPDG.resize(n);
96  fQ.resize(n);
97  fPVIndex.resize(n);
98  fNPixelHits.resize(n);
99 }
100 
101 void KFPTrackVector::Set(KFPTrackVector& v, int vSize, int offset)
102 {
109  for(int iV=0; iV<vSize; iV++)
110  {
111  for(int i=0; i<6; i++)
112  fP[i][offset+iV] = v.fP[i][iV];
113  for(int i=0; i<21; i++)
114  fC[i][offset+iV] = v.fC[i][iV];
115 #ifdef NonhomogeneousField
116  for(int i=0; i<10; i++)
117  fField[i][offset+iV] = v.fField[i][iV];
118 #endif
119 // fChi2[offset+iV] = v.fChi2[iV];
120 // fNDF[offset+iV] = v.fNDF[iV];
121  fId[offset+iV] = v.fId[iV];
122  fPDG[offset+iV] = v.fPDG[iV];
123  fQ[offset+iV] = v.fQ[iV];
124  fPVIndex[offset+iV] = v.fPVIndex[iV];
125  fNPixelHits[offset+iV] = v.fNPixelHits[iV];
126  }
127 }
128 
129 void KFPTrackVector::SetTracks(const KFPTrackVector& track, const kfvector_uint& trackIndex, const int nIndexes)
130 {
136  if(nIndexes == 0) return;
137 
138  Resize(nIndexes);
139 
140  for(int iP=0; iP<6; iP++)
141  {
142  int iElement = 0;
143  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
144  {
145  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
146  float_v& vec = reinterpret_cast<float_v&>(fP[iP][iElement]);
147  vec.gather(&(track.fP[iP][0]), index);
148  }
149  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
150  float_v& vec = reinterpret_cast<float_v&>(fP[iP][iElement]);
151  vec.gather(&(track.fP[iP][0]), index, simd_cast<float_m>(iElement+uint_v::IndexesFromZero()<nIndexes));
152 
153  }
154  for(int iC=0; iC<21; iC++)
155  {
156  int iElement=0;
157  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
158  {
159  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
160  float_v& vec = reinterpret_cast<float_v&>(fC[iC][iElement]);
161  vec.gather(&(track.fC[iC][0]), index);
162  }
163  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
164  float_v& vec = reinterpret_cast<float_v&>(fC[iC][iElement]);
165  vec.gather(&(track.fC[iC][0]), index, simd_cast<float_m>(iElement+uint_v::IndexesFromZero()<nIndexes));
166  }
167 #ifdef NonhomogeneousField
168  for(int iP=0; iP<10; iP++)
169  {
170  int iElement = 0;
171  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
172  {
173  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
174  float_v& vec = reinterpret_cast<float_v&>(fField[iP][iElement]);
175  vec.gather(&(track.fField[iP][0]), index);
176  }
177  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
178  float_v& vec = reinterpret_cast<float_v&>(fField[iP][iElement]);
179  vec.gather(&(track.fField[iP][0]), index, simd_cast<float_m>(iElement+uint_v::IndexesFromZero()<nIndexes));
180  }
181 #endif
182  {
183  int iElement=0;
184  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
185  {
186  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
187  int_v& vec = reinterpret_cast<int_v&>(fId[iElement]);
188  vec.gather(&(track.fId[0]), index);
189  }
190  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
191  int_v& vec = reinterpret_cast<int_v&>(fId[iElement]);
192  vec.gather(&(track.fId[0]), index, int_m(iElement+uint_v::IndexesFromZero()<nIndexes));
193  }
194  {
195  int iElement=0;
196  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
197  {
198  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
199  int_v& vec = reinterpret_cast<int_v&>(fPDG[iElement]);
200  vec.gather(&(track.fPDG[0]), index);
201  }
202  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
203  int_v& vec = reinterpret_cast<int_v&>(fPDG[iElement]);
204  vec.gather(&(track.fPDG[0]), index, int_m(iElement+uint_v::IndexesFromZero()<nIndexes));
205  }
206  {
207  int iElement=0;
208  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
209  {
210  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
211  int_v& vec = reinterpret_cast<int_v&>(fQ[iElement]);
212  vec.gather(&(track.fQ[0]), index);
213  }
214  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
215  int_v& vec = reinterpret_cast<int_v&>(fQ[iElement]);
216  vec.gather(&(track.fQ[0]), index, int_m(iElement+uint_v::IndexesFromZero()<nIndexes));
217  }
218  {
219  int iElement=0;
220  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
221  {
222  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
223  int_v& vec = reinterpret_cast<int_v&>(fPVIndex[iElement]);
224  vec.gather(&(track.fPVIndex[0]), index);
225  }
226  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
227  int_v& vec = reinterpret_cast<int_v&>(fPVIndex[iElement]);
228  vec.gather(&(track.fPVIndex[0]), index, int_m(iElement+uint_v::IndexesFromZero()<nIndexes));
229  }
230  {
231  int iElement=0;
232  for(iElement=0; iElement<nIndexes-float_vLen; iElement += float_vLen)
233  {
234  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
235  int_v& vec = reinterpret_cast<int_v&>(fNPixelHits[iElement]);
236  vec.gather(&(track.fNPixelHits[0]), index);
237  }
238  const uint_v& index = reinterpret_cast<const uint_v&>(trackIndex[iElement]);
239  int_v& vec = reinterpret_cast<int_v&>(fNPixelHits[iElement]);
240  vec.gather(&(track.fNPixelHits[0]), index, int_m(iElement+uint_v::IndexesFromZero()<nIndexes));
241  }
242 }
243 
244 void KFPTrackVector::GetTrack(KFPTrack& track, const int n)
245 {
250  track.SetParameters(fP[0][n],fP[1][n],fP[2][n],fP[3][n],fP[4][n],fP[5][n]);
251  for(int i=0; i<21; i++)
252  track.SetCovariance(i,fC[i][n]);
253 // track.SetChi2(fChi2[n]);
254 // track.SetNDF(fNDF[n]);
255  track.SetId(fId[n]);
256  track.SetCharge(fQ[n]);
257 
258 #ifdef NonhomogeneousField
259  for(int i=0; i<10; i++)
260  track.SetFieldCoeff( fField[i][n], i);
261 #endif
262 }
263 
264 void KFPTrackVector::RotateXY( float_v alpha, int firstElement )
265 {
280  const float_v cA = KFPMath::Cos( alpha );
281  const float_v sA = KFPMath::Sin( alpha );
282 
283  const float_v xInit = reinterpret_cast<const float_v&>(fP[0][firstElement]);
284  const float_v yInit = reinterpret_cast<const float_v&>(fP[1][firstElement]);
285 
286  float_v& x = reinterpret_cast<float_v&>(fP[0][firstElement]);
287  float_v& y = reinterpret_cast<float_v&>(fP[1][firstElement]);
288 
289  x = -(xInit*sA + yInit*cA);
290  y = xInit*cA - yInit*sA;
291 
292  const float_v pxInit = reinterpret_cast<const float_v&>(fP[3][firstElement]);
293  const float_v pyInit = reinterpret_cast<const float_v&>(fP[4][firstElement]);
294 
295  float_v& px = reinterpret_cast<float_v&>(fP[3][firstElement]);
296  float_v& py = reinterpret_cast<float_v&>(fP[4][firstElement]);
297 
298  px = -(pxInit*sA + pyInit*cA);
299  py = pxInit*cA - pyInit*sA;
300 
301  float_v cov[21];
302  for(int iC=0; iC<21; iC++)
303  cov[iC] = reinterpret_cast<const float_v&>(fC[iC][firstElement]);
304 
305  reinterpret_cast<float_v&>(fC[0][firstElement]) = cA*cA* cov[2] + 2* cA* cov[1]* sA + cov[0]*sA* sA;
306 
307  reinterpret_cast<float_v&>(fC[1][firstElement]) = -(cA*cA * cov[1]) + cA* (-cov[0] + cov[2])* sA + cov[1]*sA* sA;
308  reinterpret_cast<float_v&>(fC[2][firstElement]) = cA*cA* cov[0] - 2* cA* cov[1]* sA + cov[2]*sA* sA;
309 
310  reinterpret_cast<float_v&>(fC[3][firstElement]) = -(cA* cov[4]) - cov[3]* sA;
311  reinterpret_cast<float_v&>(fC[4][firstElement]) = cA* cov[3] - cov[4]* sA;
312 // reinterpret_cast<float_v&>(fC[5][firstElement]) = cov[5];
313 
314  reinterpret_cast<float_v&>(fC[6][firstElement]) = cA*cA* cov[11] + cA *(cov[10] + cov[7])* sA + cov[6]*sA* sA;
315  reinterpret_cast<float_v&>(fC[7][firstElement]) = -(cA*cA * cov[10]) + cA* (cov[11] - cov[6])* sA + cov[7] *sA*sA;
316  reinterpret_cast<float_v&>(fC[8][firstElement]) = -(cA *cov[12]) - cov[8] *sA;
317  reinterpret_cast<float_v&>(fC[9][firstElement]) = cA*cA* cov[14] + 2 *cA* cov[13]* sA + cov[9]* sA*sA;
318 
319  reinterpret_cast<float_v&>(fC[10][firstElement]) = -(cA*cA* cov[7]) + cA* (cov[11] - cov[6])* sA + cov[10]*sA* sA;
320  reinterpret_cast<float_v&>(fC[11][firstElement]) = cA*cA* cov[6] - cA* (cov[10] + cov[7]) *sA + cov[11]*sA* sA;
321  reinterpret_cast<float_v&>(fC[12][firstElement]) = cA* cov[8] - cov[12]* sA;
322  reinterpret_cast<float_v&>(fC[13][firstElement]) = -(cA*cA* cov[13]) + cA* (cov[14] - cov[9])* sA + cov[13]* sA*sA;
323  reinterpret_cast<float_v&>(fC[14][firstElement]) = cA*cA* cov[9] - 2* cA* cov[13]* sA + cov[14]* sA*sA;
324 
325  reinterpret_cast<float_v&>(fC[15][firstElement]) = -(cA* cov[16]) - cov[15]* sA;
326  reinterpret_cast<float_v&>(fC[16][firstElement]) = cA* cov[15] - cov[16]* sA;
327 // reinterpret_cast<float_v&>(fC[17][firstElement]) = cov[17];
328  reinterpret_cast<float_v&>(fC[18][firstElement]) = -(cA* cov[19]) - cov[18]* sA;
329  reinterpret_cast<float_v&>(fC[19][firstElement]) = cA* cov[18] - cov[19]* sA;
330 // reinterpret_cast<float_v&>(fC[20][firstElement]) = cov[20];
331 } // RotateXY
332 
333 
334 void KFPTrackVector::PrintTrack(int n)
335 {
339  for(int i=0; i<6; i++)
340  std::cout << fP[i][n] << " ";
341  std::cout << std::endl;
342 
343  for(int i=0; i<21; i++)
344  std::cout << fC[i][n] << " ";
345  std::cout << std::endl;
346 
347  std::cout << fId[n] << " " << fPDG[n] << " " << fQ[n] << " " << fPVIndex[n] << " " << fNPixelHits[n] << std::endl;
348 }
349 
350 void KFPTrackVector::Print()
351 {
354  std::cout << "NTracks " << Size() << std::endl;
355  if( Size()==0 ) return;
356 
357  std::cout << "Parameters: " << std::endl;
358  for(int iP=0; iP<6; iP++)
359  {
360  std::cout << " iP " << iP << ": ";
361  for(int iTr=0; iTr<Size(); iTr++)
362  std::cout << Parameter(iP)[iTr]<< " ";
363  std::cout << std::endl;
364  }
365 
366  std::cout << "Cov matrix: " << std::endl;
367  for(int iC=0; iC<21; iC++)
368  {
369  std::cout << " iC " << iC << ": ";
370  for(int iTr=0; iTr<Size(); iTr++)
371  std::cout << Covariance(iC)[iTr]<< " ";
372  std::cout << std::endl;
373  }
374 
375  std::cout << "Id: " << std::endl;
376  for(int iTr=0; iTr<Size(); iTr++)
377  std::cout << Id()[iTr] << " ";
378  std::cout << std::endl;
379 
380  std::cout << "Pdg: " << std::endl;
381  for(int iTr=0; iTr<Size(); iTr++)
382  std::cout << PDG()[iTr] << " ";
383  std::cout << std::endl;
384 
385  std::cout << "Q: " << std::endl;
386  for(int iTr=0; iTr<Size(); iTr++)
387  std::cout << Q()[iTr] << " ";
388  std::cout << std::endl;
389 
390  std::cout << "PV index: " << std::endl;
391  for(int iTr=0; iTr<Size(); iTr++)
392  std::cout << PVIndex()[iTr] << " ";
393  std::cout << std::endl;
394 
395  std::cout << "fNPixelHits: " << std::endl;
396  for(int iTr=0; iTr<Size(); iTr++)
397  std::cout << NPixelHits()[iTr] << " ";
398  std::cout << std::endl;
399 
400 #ifdef NonhomogeneousField
401  std::cout << "Field: " << std::endl;
402  for(int iF=0; iF<6; iF++)
403  {
404  std::cout << " iF " << iF << ": ";
405  for(int iTr=0; iTr<Size(); iTr++)
406  std::cout << FieldCoefficient(iF)[iTr]<< " ";
407  std::cout << std::endl;
408  }
409 #endif
410 
411  std::cout << "Last particle index: " << std::endl;
412  std::cout << LastElectron() << " "
413  << LastMuon() << " "
414  << LastPion() << " "
415  << LastKaon() << " "
416  << LastProton() << " "
417  << LastDeuteron() << " "
418  << LastTritium() << " "
419  << LastHe3() << " "
420  << LastHe4() << std::endl;
421 }
422