FairRoot/PandaRoot
KFParticleFinder.cxx
Go to the documentation of this file.
1 #include "KFParticleFinder.h"
2 
3 #define cnst static const fvec
4 
5 //for particle finding
6 #include "KFPTrack.h"
7 #include "PndFTSCADef.h"
8 #include "PndFTSCAMath.h"
9 #include <map>
10 using std::map;
11 using std::vector;
12 
13 
14 const float KFParticleFinder::DefaultCuts[2][3] = {{3.,3.,-100.},{3.,3.,-100.}};
15 
16 KFParticleFinder::KFParticleFinder()
17 {
18 }
19 
20 void KFParticleFinder::FindParticles(vector<KFPTrack> &vRTracks, vector<float>& ChiToPrimVtx,
21 #ifdef NonhomogeniousField
22  vector<L1FieldRegion>& vField,
23 #endif
24  vector<KFParticle>& Particles,
25  KFParticleSIMD& PrimVtx, const vector<int>& vTrackPDG, const float cuts[2][3])
26 {
27  //* Finds particles (K0s and Lambda) from a given set of tracks
28 
29  static const int NTrackTypes = 8;
30 
31  int pdgPos[NTrackTypes]={-11,-13, 211, 321, 2212, 211, 321, 2212};
32  int pdgNeg[NTrackTypes]={ 11, 13,-211,-321,-2212, -211,-321,-2212};
33 
34  vector<short> idPosSec[NTrackTypes];
35  vector<short> idNegSec[NTrackTypes];
36 
37  vector<short> idPosPrim[NTrackTypes];
38  vector<short> idNegPrim[NTrackTypes];
39 
40  for(unsigned short iTr=0; iTr < vRTracks.size(); iTr++)
41  {
42  KFPTrack &kfTrack = vRTracks[iTr];
43  bool ok = 1;
44  for(unsigned short iT=0; iT<6; iT++)
45  ok = ok && finite(kfTrack.GetTrack()[iT]);
46  for(unsigned short iC=0; iC<21; iC++)
47  ok = ok && finite(kfTrack.GetCovMatrix()[iC]);
48  ok = ok && (kfTrack.GetCovMatrix()[0] < 100. && kfTrack.GetCovMatrix()[0] > 0.)
49  && (kfTrack.GetCovMatrix()[2] < 100. && kfTrack.GetCovMatrix()[2] > 0.)
50  && (kfTrack.GetCovMatrix()[5] < 100. && kfTrack.GetCovMatrix()[5] > 0.)
51  && (kfTrack.GetCovMatrix()[9] < 1. && kfTrack.GetCovMatrix()[9] > 0.)
52  && (kfTrack.GetCovMatrix()[14] < 1. && kfTrack.GetCovMatrix()[14] > 0.)
53  && (kfTrack.GetCovMatrix()[20] < 1. && kfTrack.GetCovMatrix()[20] > 0.);
54  ok = ok && kfTrack.GetChi2() < 10*kfTrack.GetNDF();
55 // if(!ok) continue;
56  const int pdg = abs(vTrackPDG[iTr]);
57 
58  short pdgIndex = -1;
59  switch (pdg)
60  {
61  case 11: pdgIndex = 0; break;
62  case 13: pdgIndex = 1; break;
63  case 211: pdgIndex = 2; break;
64  case 321: pdgIndex = 3; break;
65  case 2212: pdgIndex = 4; break;
66  }
67 
68  short incr = 3;
69  short pdgIndexMax = pdgIndex+incr;
70 
71  if(pdgIndex<2)
72  {
73  incr = 1;
74  pdgIndexMax = pdgIndex;
75  }
76 
77  if(pdgIndex < 0)
78  {
79  pdgIndex = 5;
80  pdgIndexMax = 7;
81  incr = 1;
82  }
83 
84  if( ChiToPrimVtx[iTr] < cuts[0][0] )
85  {
86  if(kfTrack.Charge() >= 0.)
87  {
88  for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
89  idPosPrim[ipdg].push_back(iTr);
90  }
91  if(kfTrack.Charge() < 0.)
92  for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
93  idNegPrim[ipdg].push_back(iTr);
94  }
95  else
96  {
97  if(kfTrack.Charge() >= 0.)
98  for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
99  idPosSec[ipdg].push_back(iTr);
100  if(kfTrack.Charge() < 0.)
101  for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
102  idNegSec[ipdg].push_back(iTr);
103  }
104  }
105 
106  const int nPart = idPosSec[5].size() * idNegSec[5].size()+
107  idPosSec[5].size() * idNegSec[7].size()+
108  idPosSec[7].size() * idNegSec[5].size()+
109  idPosPrim[2].size() * idNegPrim[3].size() +
110  idPosPrim[3].size() * idNegPrim[2].size() +
111  idPosPrim[3].size() * idNegPrim[3].size() +
112  idPosPrim[4].size() * idNegPrim[3].size() +
113  idPosPrim[3].size() * idNegPrim[4].size() +
114  idPosPrim[0].size() * idNegPrim[0].size() +
115  idPosPrim[1].size() * idNegPrim[1].size();
116 
117 //std::cout << "NPart estim " << nPart << std::endl;
118  Particles.reserve(vRTracks.size() + nPart);
119 
120  const float massLambdaPDG = 1.115683;
121  const float massK0sPDG = 0.497614;
122  const float massXiPDG = 1.32171;
123 
124 #if defined(PANDA_STT) || defined(PANDA_FTS)
125  const float massLambdaSigma=3.7e-3;
126  const float massK0sSigma=2.2e-3; //TODO tune
127  const float massXiSigma=2.e-3; //TODO tune
128 #endif
129 #ifdef ALICE_ITS
130  const float massLambdaSigma=5.9e-3;
131  const float massK0sSigma=17.7e-3;
132  const float massXiSigma=7.3e-3;
133 #endif
134 #ifdef STAR_HFT
135  const float massLambdaSigma=5.9e-3;
136  const float massK0sSigma=17.7e-3;
137  const float massXiSigma=7.3e-3;
138 #endif
139 
140  for(unsigned short iTr=0; iTr < vRTracks.size(); iTr++) {
141  KFPTrack& kfTrack = vRTracks[iTr] ;
142 // kfTrack.SetId(iTr);
143  KFParticle tmp(kfTrack, vTrackPDG[iTr]);
144  tmp.SetPDG(211);//vTrackPDG[iTr]);
145  tmp.SetId(Particles.size());
146  tmp.SetNDaughters(1);
147  tmp.AddDaughterId( kfTrack.Id() );
148  Particles.push_back(tmp);
149  }
150 
151  vector<float> vLambdaTopoChi2Ndf;
152  vector<KFParticle> vLambdaSec;
153  vector<KFParticle> vLambdaPrim;
154 
155  vector<float> vLambdaBarTopoChi2Ndf;
156  vector<KFParticle> vLambdaBarSec;
157  vector<KFParticle> vLambdaBarPrim;
158 
159  vector<float> vK0sTopoChi2Ndf;
160  vector<KFParticle> vK0sPrim;
161 
162  vector<KFParticle> vXiPrim;
163  vector<KFParticle> vXiSec;
164  vector<KFParticle> vXiBarPrim;
165 
166  vector<KFParticle> vXiStarPrim;
167  vector<KFParticle> vXiStarBarPrim;
168 
169  const float SecCuts[3] = {3.f,5.f,5.f};
170 // const float SecCuts[3] = {3.f,5.f,4.f};
171 
172 
173  //K0s -> pi+ pi-
174  Find2DaughterDecay(vRTracks,
175 #ifdef NonhomogeniousField
176  vField,
177 #endif
178  Particles, pdgNeg[5], pdgPos[5], 310,
179  idNegSec[5], idPosSec[5], PrimVtx, cuts[0], 0, &vK0sTopoChi2Ndf,
180  SecCuts, massK0sPDG, massK0sSigma, &vK0sPrim, 0);
181  //Lambda -> p pi-
182  Find2DaughterDecay(vRTracks,
183 #ifdef NonhomogeniousField
184  vField,
185 #endif
186  Particles, pdgNeg[5], pdgPos[7], 3122,
187  idNegSec[5], idPosSec[7], PrimVtx, cuts[1], 0,&vLambdaTopoChi2Ndf,
188  SecCuts, massLambdaPDG, massLambdaSigma, &vLambdaPrim, &vLambdaSec);
189  //Lambda_bar -> pi+ p-
190  Find2DaughterDecay(vRTracks,
191 #ifdef NonhomogeniousField
192  vField,
193 #endif
194  Particles, pdgPos[5], pdgNeg[4], -3122,
195  idPosSec[5], idNegSec[4], PrimVtx, cuts[1], 0, &vLambdaBarTopoChi2Ndf,
196  SecCuts, massLambdaPDG, massLambdaSigma, &vLambdaBarPrim, &vLambdaBarSec);
197  //K*0 -> K+ pi-
198  Find2DaughterDecay(vRTracks,
199 #ifdef NonhomogeniousField
200  vField,
201 #endif
202  Particles, pdgPos[3], pdgNeg[2], 313,
203  idPosPrim[3], idNegPrim[2], PrimVtx, cuts[1], 1);
204  //K*0_bar -> pi+ K-
205  Find2DaughterDecay(vRTracks,
206 #ifdef NonhomogeniousField
207  vField,
208 #endif
209  Particles, pdgNeg[3], pdgPos[2], -313,
210  idNegPrim[3], idPosPrim[2], PrimVtx, cuts[1], 1);
211  //Lambda* -> p K-
212  Find2DaughterDecay(vRTracks,
213 #ifdef NonhomogeniousField
214  vField,
215 #endif
216  Particles, pdgNeg[3], pdgPos[4], 3124,
217  idNegPrim[3], idPosPrim[4], PrimVtx, cuts[1], 1);
218  //Lambda*_bar -> K+ p-
219  Find2DaughterDecay(vRTracks,
220 #ifdef NonhomogeniousField
221  vField,
222 #endif
223  Particles, pdgNeg[4], pdgPos[3], -3124,
224  idNegPrim[4], idPosPrim[3], PrimVtx, cuts[1], 1);
225  //phi -> K+ K-
226  Find2DaughterDecay(vRTracks,
227 #ifdef NonhomogeniousField
228  vField,
229 #endif
230  Particles, pdgNeg[3], pdgPos[3], 333,
231  idNegPrim[3], idPosPrim[3], PrimVtx, cuts[1], 1);
232 // //rho -> pi+ pi-
233 // Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[2], pdgPos[2], 113,
234 // idNegPrim[2], idPosPrim[2],
235 // PrimVtx, cuts[1]);
236  //gamma -> e+ e-
237  Find2DaughterDecay(vRTracks,
238 #ifdef NonhomogeniousField
239  vField,
240 #endif
241  Particles, pdgNeg[0], pdgPos[0], 22,
242  idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1);
243  Find2DaughterDecay(vRTracks,
244 #ifdef NonhomogeniousField
245  vField,
246 #endif
247  Particles, pdgNeg[0], pdgPos[0], 22,
248  idNegSec[0], idPosSec[0], PrimVtx, cuts[1], 0);
249  Find2DaughterDecay(vRTracks,
250 #ifdef NonhomogeniousField
251  vField,
252 #endif
253  Particles, pdgNeg[0], pdgPos[0], 22,
254  idNegSec[0], idPosPrim[0], PrimVtx, cuts[1], 0);
255  Find2DaughterDecay(vRTracks,
256 #ifdef NonhomogeniousField
257  vField,
258 #endif
259  Particles, pdgNeg[0], pdgPos[0], 22,
260  idNegPrim[0], idPosSec[0], PrimVtx, cuts[1], 0);
261  //JPsi-> e+ e-
262  Find2DaughterDecay(vRTracks,
263 #ifdef NonhomogeniousField
264  vField,
265 #endif
266  Particles, pdgNeg[0], pdgPos[0], 443,
267  idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1, 1.f);
268  //JPsi-> mu+ mu-
269  Find2DaughterDecay(vRTracks,
270 #ifdef NonhomogeniousField
271  vField,
272 #endif
273  Particles, pdgNeg[1], pdgPos[1], 100443,
274  idNegPrim[1], idPosPrim[1], PrimVtx, cuts[1], 1, 1.f);
275  //rho, omega, phi -> e+ e-
276  Find2DaughterDecay(vRTracks,
277 #ifdef NonhomogeniousField
278  vField,
279 #endif
280  Particles, pdgNeg[0], pdgPos[0], 100113,
281  idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1, 0.2f);
282  //rho, omega, phi -> mu+ mu-
283  const float PCut = 1.f;
284  Find2DaughterDecay(vRTracks,
285 #ifdef NonhomogeniousField
286  vField,
287 #endif
288  Particles, pdgNeg[1], pdgPos[1], 200113,
289  idNegPrim[1], idPosPrim[1], PrimVtx, cuts[1], 1, 0.2f, -100, 0, &PCut);
290 
291  ExtrapolateToPV(vLambdaPrim,PrimVtx);
292  ExtrapolateToPV(vLambdaBarPrim,PrimVtx);
293  ExtrapolateToPV(vK0sPrim,PrimVtx);
294 
295  // Find Xi-
296  float cutXi[3] = {10.,5.,6.};
297 // float cutXi[3] = {-300.,10.,10.};
298  FindTrackV0Decay(3312, Particles, vLambdaSec, vRTracks,
299 #ifdef NonhomogeniousField
300  vField,
301 #endif
302  pdgNeg[5], idNegSec[5],
303  PrimVtx, cutXi, 0, 0, &vXiPrim, massXiPDG, massXiSigma );
304 
305 
306  float cutLL[3] = {10.,10000000.,3.};
307  float cutLL2[3] = {10.,3.,3.};
308  vector<KFParticle> vLL;
309  FindTrackV0Decay(3002, vLL, vLambdaSec, vRTracks,
310 #ifdef NonhomogeniousField
311  vField,
312 #endif
313  pdgNeg[5], idNegSec[5],
314  PrimVtx, cutLL, 0, &ChiToPrimVtx);
315  // Find H0->Lambda p pi-
316  //Find Omega*-
317  FindTrackV0Decay(3001, Particles, vLL, vRTracks,
318 #ifdef NonhomogeniousField
319  vField,
320 #endif
321  pdgPos[4], idPosSec[4],
322  PrimVtx, cutLL2, 0, &ChiToPrimVtx);
323  // Find Xi+
324  float cutXiPlus[3] = {10.,5.,6.};
325  FindTrackV0Decay(-3312, Particles, vLambdaBarSec, vRTracks,
326 #ifdef NonhomogeniousField
327  vField,
328 #endif
329  pdgPos[5], idPosSec[5],
330  PrimVtx, cutXiPlus, 0, 0, &vXiBarPrim, massXiPDG, massXiSigma);
331  //Find Omega-
332  float cutOmega[3] = {10.,3.,3.};
333  FindTrackV0Decay(3334, Particles, vLambdaSec, vRTracks,
334 #ifdef NonhomogeniousField
335  vField,
336 #endif
337  pdgNeg[6], idNegSec[6],
338  PrimVtx, cutOmega, 0, &ChiToPrimVtx);
339  //Find Omega+
340  float cutOmegaPlus[3] = {10.,3.,3.};
341  FindTrackV0Decay(-3334, Particles, vLambdaBarSec, vRTracks,
342 #ifdef NonhomogeniousField
343  vField,
344 #endif
345  pdgPos[6], idPosSec[6],
346  PrimVtx, cutOmegaPlus, 0, &ChiToPrimVtx);
347  //Find Xi*-
348  float cutXiStarMinus[3] = {-100.,10000.,3.};
349  FindTrackV0Decay(1003314, Particles, vLambdaPrim, vRTracks,
350 #ifdef NonhomogeniousField
351  vField,
352 #endif
353  pdgNeg[3], idNegPrim[3],
354  PrimVtx, cutXiStarMinus, 1);
355  //Find Xi*+
356  float cutXiStarPlus[3] = {-100.,10000.,3.};
357  FindTrackV0Decay(-1003314, Particles, vLambdaBarPrim, vRTracks,
358 #ifdef NonhomogeniousField
359  vField,
360 #endif
361  pdgPos[3], idPosPrim[3],
362  PrimVtx, cutXiStarPlus, 1);
363 
364  ExtrapolateToPV(vXiPrim,PrimVtx);
365  ExtrapolateToPV(vXiBarPrim,PrimVtx);
366 
367  //Find Xi*0
368  float cutXiStar0[3] = {-100.,10000.,3.};
369  FindTrackV0Decay(3324, Particles, vXiPrim, vRTracks,
370 #ifdef NonhomogeniousField
371  vField,
372 #endif
373  pdgPos[5], idPosPrim[5],
374  PrimVtx, cutXiStar0, 1, 0, &vXiStarPrim);
375  //Find Xi*0 bar
376  float cutXiBarStar0[3] = {-100.,10000.,3.};
377  FindTrackV0Decay(-3324, Particles, vXiBarPrim, vRTracks,
378 #ifdef NonhomogeniousField
379  vField,
380 #endif
381  pdgNeg[5], idNegPrim[5],
382  PrimVtx, cutXiBarStar0, 1, 0, &vXiStarBarPrim);
383  //Find Omega*-
384  const float cutOmegaStar[2] = {-100., 3.};
385  for(unsigned int iPart=0; iPart<vXiStarPrim.size(); iPart++)
386  CombineTrackPart(vRTracks,
387 #ifdef NonhomogeniousField
388  vField,
389 #endif
390  Particles, vXiStarPrim[iPart], pdgNeg[3],
391  1003334, idNegPrim[3], cutOmegaStar, 0);
392  //Find Omega*+
393  for(unsigned int iPart=0; iPart<vXiStarBarPrim.size(); iPart++)
394  CombineTrackPart(vRTracks,
395 #ifdef NonhomogeniousField
396  vField,
397 #endif
398  Particles, vXiStarBarPrim[iPart], pdgPos[3],
399  -1003334, idPosPrim[3], cutOmegaStar, 0);
400  // Find K*+
401  float cutKStarPlus[3] = {-100.,10000.,3.};
402  FindTrackV0Decay(323, Particles, vK0sPrim, vRTracks,
403 #ifdef NonhomogeniousField
404  vField,
405 #endif
406  pdgPos[5], idPosPrim[5],
407  PrimVtx, cutKStarPlus, 1, 0);
408  // Find K*-
409  float cutKStarMinus[3] = {-100.,10000.,3.};
410  FindTrackV0Decay(-323, Particles, vK0sPrim, vRTracks,
411 #ifdef NonhomogeniousField
412  vField,
413 #endif
414  pdgNeg[5], idNegPrim[5],
415  PrimVtx, cutKStarMinus, 1, 0);
416  // Find Sigma*+
417  float cutSigmaStarPlus[3] = {-100.,10000.,3.};
418  FindTrackV0Decay(3224, Particles, vLambdaPrim, vRTracks,
419 #ifdef NonhomogeniousField
420  vField,
421 #endif
422  pdgPos[5], idPosPrim[5],
423  PrimVtx, cutSigmaStarPlus, 1, 0);
424  // Find Sigma*+ bar
425  float cutSigmaStarPlusBar[3] = {-100.,10000.,3.};
426  FindTrackV0Decay(-3114, Particles, vLambdaBarPrim, vRTracks,
427 #ifdef NonhomogeniousField
428  vField,
429 #endif
430  pdgPos[5], idPosPrim[5],
431  PrimVtx, cutSigmaStarPlusBar, 1, 0);
432  // Find Sigma*-
433  float cutSigmaStarMinus[3] = {-100.,10000.,3.};
434  FindTrackV0Decay(3114, Particles, vLambdaPrim, vRTracks,
435 #ifdef NonhomogeniousField
436  vField,
437 #endif
438  pdgNeg[5], idNegPrim[5],
439  PrimVtx, cutSigmaStarMinus, 1, 0);
440  // Find Sigma*- bar
441  float cutSigmaStarMinusBar[3] = {-100.,10000.,3.};
442  FindTrackV0Decay(-3224, Particles, vLambdaBarPrim, vRTracks,
443 #ifdef NonhomogeniousField
444  vField,
445 #endif
446  pdgNeg[5], idNegPrim[5],
447  PrimVtx, cutSigmaStarMinusBar, 1, 0);
448 
449 
450  // Find H-dibarion
451  vector<KFParticle> vHdibarion;
452  float cutHdb[3] = {3.,3.,3.};
453  for(unsigned short iL=0; iL < vLambdaSec.size(); iL++)
454  {
455  KFParticleSIMD vDaughters[2] = {KFParticleSIMD(),KFParticleSIMD(vLambdaSec[iL])};
456 
457  vector<int> daughterIds;
458  for(unsigned int iD=0; iD<vLambdaSec[iL].DaughterIds().size(); iD++)
459  daughterIds.push_back(vLambdaSec[iL].DaughterIds()[iD]);
460  FindHyperons(3000, vDaughters, daughterIds, vLambdaSec, vHdibarion, PrimVtx, cutHdb, iL+1);
461  }
462 
463  for(unsigned int iH=0; iH<vHdibarion.size(); iH++)
464  {
465  vHdibarion[iH].SetId(Particles.size());
466  Particles.push_back(vHdibarion[iH]);
467  }
468  // chi2_prim chi2_geo z pt chi2_topo
469  const float cutsD[8][8] = {{ 6., 3., 0.04, 0.3, 3.}, //D0 -> pi+ K-
470  { 6., 3., 0.04, 0.3, 3.}, //D+ -> K- pi+ pi+
471  { 6., 3., 0.04, 0.3, 3.}, //D0 -> pi+ pi+ pi- K-
472  { 6., 3., 0.04, 0.3, 3.}, //Ds+ -> K- K+ pi+
473  { 6., 3., 0.04, 0.3, 3.}, //Lambdac -> pi+ K- p
474  { 6., 3., -100., 0.3, -100.}, //D*0 -> D+ pi-
475  { 6., 3., -100., 0.3, -100.}, //D*+ -> D0 pi+
476  { 6., 3., -100., 0.3, -100.}}; //D*+4 -> D04 pi+
477 
478  const int DMesLambdcDaughterPDG[5] = { 211, -321, -211, 321, 2212 };
479  const int DMesLambdcMotherPDG[8] = { 421, 411, 100421, 431, 4122, 10421, 10411, 20411 };
480  vector<short>* DMesLambdcIdTrack[5] = {&idPosSec[5],
481  &idNegSec[3],
482  &idNegSec[5],
483  &idPosSec[3],
484  &idPosSec[4]};
485  FindDMesLambdac(vRTracks,
486 #ifdef NonhomogeniousField
487  vField,
488 #endif
489  Particles, DMesLambdcDaughterPDG, DMesLambdcMotherPDG,
490  DMesLambdcIdTrack, PrimVtx, cutsD, ChiToPrimVtx);
491 
492  const int DMesLambdcBarDaughterPDG[5] = { -211, 321, 211, -321, -2212 };
493  const int DMesLambdcBarMotherPDG[8] = { -421, -411, -100421, -431, -4122, -10421, -10411, -20411 };
494  vector<short>* DMesLambdcBarIdTrack[5] = {&idNegSec[5],
495  &idPosSec[3],
496  &idPosSec[5],
497  &idNegSec[3],
498  &idNegSec[4]};
499  FindDMesLambdac(vRTracks,
500 #ifdef NonhomogeniousField
501  vField,
502 #endif
503  Particles, DMesLambdcBarDaughterPDG, DMesLambdcBarMotherPDG,
504  DMesLambdcBarIdTrack, PrimVtx, cutsD, ChiToPrimVtx);
505 // std::cout << "NPart " << Particles.size() << std::endl;
506 }
507 
508 void KFParticleFinder::ExtrapolateToPV(vector<KFParticle>& vParticles, KFParticleSIMD& PrimVtx)
509 {
510  for(unsigned int iL=0; iL<vParticles.size(); iL += fvecLen)
511  {
512  KFParticle* parts[fvecLen];
513  unsigned int nPart = vParticles.size();
514 
515  unsigned int nEntries = (iL + fvecLen < nPart) ? fvecLen : (nPart - iL);
516 
517  for(unsigned short iv=0; iv<nEntries; iv++)
518  parts[iv] = &vParticles[iL+iv];
519 
520 
521  KFParticleSIMD tmp(parts,nEntries);
522  tmp.TransportToPoint(PrimVtx.Parameters());
523 
524  for(unsigned int iv=0; iv<nEntries; iv++)
525  {
526  tmp.GetKFParticle(vParticles[iL+iv], iv);
527  }
528  }
529 }
530 
531 fvec KFParticleFinder::GetChi2BetweenParticles(KFParticleSIMD &p1, KFParticleSIMD &p2)
532 {
533  const fvec& x1 = p1.GetX();
534  const fvec& y1 = p1.GetY();
535  const fvec& z1 = p1.GetZ();
536 
537  const fvec& x2 = p2.GetX();
538  const fvec& y2 = p2.GetY();
539  const fvec& z2 = p2.GetZ();
540 
541  const fvec dx = x1 - x2;
542  const fvec dy = y1 - y2;
543  const fvec dz = z1 - z2;
544 
545  const fvec& c0 = p1.GetCovariance(0) + p2.GetCovariance(0);
546  const fvec& c1 = p1.GetCovariance(1) + p2.GetCovariance(1);
547  const fvec& c2 = p1.GetCovariance(2) + p2.GetCovariance(2);
548  const fvec& c3 = p1.GetCovariance(3) + p2.GetCovariance(3);
549  const fvec& c4 = p1.GetCovariance(4) + p2.GetCovariance(4);
550  const fvec& c5 = p1.GetCovariance(5) + p2.GetCovariance(5);
551 
552  const fvec r2 = dx*dx + dy*dy + dz*dz;
553  const fvec err2 = c0*dx*dx + c2*dy*dy + c5*dz*dz + 2*( c1*dx*dy + c3*dx*dz + c4*dy*dz );
554 
555  return (r2*r2/err2);
556 }
557 
558 void KFParticleFinder::Find2DaughterDecay(vector<KFPTrack>& vTracks,
559 #ifdef NonhomogeniousField
560  const vector<L1FieldRegion>& vField,
561 #endif
562  vector<KFParticle>& Particles,
563  const int DaughterNegPDG,
564  const int DaughterPosPDG,
565  const int MotherPDG,
566  vector<short>& idNeg,
567  vector<short>& idPos,
568  KFParticleSIMD& PrimVtx,
569  const float* cuts,
570  bool isPrimary,
571  vector<float>* vMotherTopoChi2Ndf,
572  const float* secCuts,
573  const float massMotherPDG,
574  const float massMotherPDGSigma,
575  vector<KFParticle>* vMotherPrim,
576  vector<KFParticle>* vMotherSec )
577 {
578  const int NPositive = idPos.size();
579  KFParticle mother_temp;
580 
581  KFParticleSIMD mother;
582  mother.SetPDG( MotherPDG );
583  KFParticleSIMD motherTopo;
584 
585  const short NPosVect = NPositive%fvecLen == 0 ? NPositive/fvecLen : NPositive/fvecLen+1;
586  vector<KFParticleSIMD> posPart(NPosVect);
587 
588  // create particles from tracks
589  KFPTrack* vvPos[fvecLen];
590 
591  int iPosVect = 0;
592  for(unsigned short iTrP=0; iTrP < NPositive; iTrP += fvecLen)
593  {
594  const unsigned short NTracks = (iTrP + fvecLen < NPositive) ? fvecLen : (NPositive - iTrP);
595 
596 #ifdef NonhomogeniousField
597  L1FieldRegion posField;
598 #endif
599  for(unsigned short iv=0; iv<NTracks; iv++)
600  {
601  vvPos[iv] = &vTracks[idPos[iTrP+iv]];
602 #ifdef NonhomogeniousField
603  int entrSIMDPos = idPos[iTrP+iv] % fvecLen;
604  int entrVecPos = idPos[iTrP+iv] / fvecLen;
605  posField.SetOneEntry(iv,vField[entrVecPos],entrSIMDPos);
606 #endif
607  }
608 
609  posPart[iPosVect].Create(vvPos,NTracks,0,&DaughterPosPDG);
610 #ifdef NonhomogeniousField
611  posPart[iPosVect].SetField(posField);
612 #endif
613  fvec posId;
614  for(int iv=0; iv<NTracks; iv++)
615  posId[iv] = idPos[iTrP+iv];
616  posPart[iPosVect].SetId(posId);
617  iPosVect++;
618  }
619 
620  for(unsigned short iTrN=0; iTrN < idNeg.size(); iTrN++)
621  {
622  KFPTrack &kfTrackNeg = vTracks[idNeg[iTrN]];
623  KFParticleSIMD vDaughters[2] = {KFParticleSIMD(kfTrackNeg,0,&DaughterNegPDG), KFParticleSIMD()};
624 #ifdef NonhomogeniousField
625  int entrSIMD = idNeg[iTrN] % fvecLen;
626  int entrVec = idNeg[iTrN] / fvecLen;
627  vDaughters[0].SetField(vField[entrVec],1,entrSIMD);
628 #endif
629  vDaughters[0].SetId(idNeg[iTrN]);
630 
631  for(unsigned short iTrP=0; iTrP < NPositive; iTrP += fvecLen)
632  {
633  const unsigned short NTracks = (iTrP + fvecLen < NPositive) ? fvecLen : (NPositive - iTrP);
634 
635  vDaughters[1] = posPart[iTrP/fvecLen];
636 
637  if(isPrimary)
638  {
639  fvec errGuess[3] = {100*sqrt(PrimVtx.CovarianceMatrix()[0]),
640  100*sqrt(PrimVtx.CovarianceMatrix()[2]),
641  100*sqrt(PrimVtx.CovarianceMatrix()[5])};
642  mother.SetVtxGuess(PrimVtx.X(), PrimVtx.Y(), PrimVtx.Z());
643  mother.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
644  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
645  mother.Construct(vDaughtersPointer, 2, 0);
646  }
647  else
648  {
649  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
650  mother.Construct(vDaughtersPointer, 2, 0);
651  }
652 
653  if(vMotherTopoChi2Ndf)
654  {
655  motherTopo = mother;
656  motherTopo.SetProductionVertex(PrimVtx);
657  }
658 
659 // if( (fabs(DaughterNegPDG) == 211) && (fabs(DaughterPosPDG) == 211) )
660 // {
661 // fvec pl0[3] = {vDaughters[0].GetPx(),vDaughters[0].GetPy(),vDaughters[0].GetPz()};
662 // fvec pl1[3] = {vDaughters[1].GetPx(),vDaughters[1].GetPy(),vDaughters[1].GetPz()};
663 //
664 // fvec pcm[3] = {pl0[0]+pl1[0], pl0[1]+pl1[1], pl0[2]+pl1[2]};
665 // fvec ecm = vDaughters[0].GetE()+vDaughters[1].GetE();
666 //
667 // fvec gamma = sqrt(fvec(1.) - (pcm[0]*pcm[0]+pcm[1]*pcm[1]+pcm[2]*pcm[2])/(ecm*ecm));
668 //
669 // fvec e0 = gamma*(vDaughters[0].GetE() - (pcm[0]*pl0[0]+pcm[1]*pl0[1]+pcm[2]*pl0[2])/ecm);
670 // fvec e1 = gamma*(vDaughters[1].GetE() - (pcm[0]*pl1[0]+pcm[1]*pl1[1]+pcm[2]*pl1[2])/ecm);
671 // if( fabs(e0[0]-e1[0]) > 1.e-7 )
672 // std::cout << e0 << " " << e1 << std::endl;
673 // }
674 
675  // KFParticleSIMD mother_temp[fvecLen];
676  // mother.GetKFParticle(mother_temp, NTracks);
677 
678  for(int iv=0; iv<NTracks; iv++)
679  {
680  if(!finite(mother.GetChi2()[iv])) continue;
681  if(!(mother.GetChi2()[iv] > 0.0f)) continue;
682  if(!(mother.GetChi2()[iv]==mother.GetChi2()[iv])) continue;
683  fvec l, dl, isParticleFromVertex;
684  mother.GetDistanceToVertexLine(PrimVtx, l, dl, &isParticleFromVertex);
685  if((l/dl)[iv] < cuts[2]) continue;
686  if(isPrimary && ((l/dl)[iv]>3) ) continue;
687  if(!(isParticleFromVertex[iv])) continue;
688  if( mother.GetChi2()[iv]/mother.GetNDF()[iv] < cuts[1] )
689  {
690  mother.GetKFParticle(mother_temp, iv);
691  mother_temp.SetId(Particles.size());
692  Particles.push_back(mother_temp);
693  float motherTopoChi2Ndf(0);
694  if(vMotherTopoChi2Ndf)
695  motherTopoChi2Ndf = motherTopo.GetChi2()[iv]/motherTopo.GetNDF()[iv];
696 
697  if(vMotherPrim)
698  {
699  Double_t mass, errMass;
700 
701  mother_temp.GetMass(mass, errMass);
702  if( (fabs(mass - massMotherPDG)/massMotherPDGSigma) > secCuts[0] ) continue;
703  if((l/dl)[iv] < secCuts[2]) continue;
704 
705  mother_temp.SetNonlinearMassConstraint(massMotherPDG);
706 
707  if( motherTopoChi2Ndf < secCuts[1] )
708  {
709  vMotherPrim->push_back(mother_temp);
710  }
711  else if(vMotherSec)
712  {
713  //if((l/dl)[iv] < secCuts[2]) continue;
714  vMotherSec->push_back(mother_temp);
715  }
716  }
717  }
718  }
719  }
720  }
721 }
722 
723 void KFParticleFinder::Find2DaughterDecay(vector<KFPTrack>& vTracks,
724 #ifdef NonhomogeniousField
725  const vector<L1FieldRegion>& vField,
726 #endif
727  vector<KFParticle>& Particles,
728  const int DaughterNegPDG,
729  const int DaughterPosPDG,
730  const int MotherPDG,
731  vector<short>& idNeg,
732  vector<short>& idPos,
733  KFParticleSIMD& PrimVtx,
734  const float* cuts,
735  bool isPrimary,
736  const float PtCut,
737  const float Chi2PrimCut,
738  vector<float>* ChiToPrimVtx,
739  const float* PCut)
740 {
741  vector<short> idPosPt;
742  vector<short> idNegPt;
743 
744  for(unsigned int iEl=0; iEl<idPos.size(); iEl++)
745  {
746  KFPTrack& kfTrack = vTracks[idPos[iEl]];
747  float pt = kfTrack.GetPt();
748  float p = kfTrack.GetP();
749  if(PCut)
750  if(p < *PCut) continue;
751  if(pt<PtCut) continue;
752  if(Chi2PrimCut > -1)
753  if(ChiToPrimVtx->at(idPos[iEl]) < Chi2PrimCut) continue;
754  idPosPt.push_back(idPos[iEl]);
755  }
756 
757  for(unsigned int iEl=0; iEl<idNeg.size(); iEl++)
758  {
759  KFPTrack& kfTrack = vTracks[idNeg[iEl]];
760  float pt = kfTrack.GetPt();
761  float p = kfTrack.GetP();
762  if(PCut)
763  if(p < *PCut) continue;
764  if(pt<PtCut) continue;
765  if(Chi2PrimCut > -1)
766  if(ChiToPrimVtx->at(idNeg[iEl]) < Chi2PrimCut) continue;
767  idNegPt.push_back(idNeg[iEl]);
768  }
769 
770  Find2DaughterDecay(vTracks,
771 #ifdef NonhomogeniousField
772  vField,
773 #endif
774  Particles, DaughterNegPDG, DaughterPosPDG, MotherPDG,
775  idNegPt, idPosPt,
776  PrimVtx, cuts, isPrimary);
777 }
778 
779 void KFParticleFinder::FindTrackV0Decay(const int MotherPDG,
780  vector<KFParticle>& Particles,
781  vector<KFParticle>& vV0,
782  vector<KFPTrack>& vTracks,
783 #ifdef NonhomogeniousField
784  const vector<L1FieldRegion>& vField,
785 #endif
786  const int DaughterPDG,
787  vector<short>& idTrack,
788  KFParticleSIMD& PrimVtx,
789  const float* cuts,
790  bool isPrimary,
791  vector<float>* ChiToPrimVtx,
792  vector<KFParticle>* vHyperonPrim,
793  float hyperonPrimMass,
794  float hyperonPrimMassErr,
795  vector<KFParticle>* vHyperonSec)
796 {
797  KFParticle hyperon_temp;
798  KFParticleSIMD hyperon;
799  hyperon.SetPDG( MotherPDG );
800 
801  for(unsigned short iV0=0; iV0 < vV0.size(); iV0++)
802  {
803  unsigned short nElements = 0;
804  KFParticleSIMD vDaughters[2]= {KFParticleSIMD(vV0[iV0]),KFParticleSIMD()};
805 
806  KFPTrack* vvTr[fvecLen];
807 #ifdef NonhomogeniousField
808  L1FieldRegion field;
809 #endif
810  fvec trId;
811 
812  for(unsigned short iTr=0; iTr < idTrack.size(); iTr++)
813  {
814  bool ok = 1;
815  if(ChiToPrimVtx)
816  if( (ChiToPrimVtx->at(idTrack[iTr]) < 7) ) ok=0; //TODO 7 for Omega
817 
818  if(ok)
819  {
820  trId[nElements] = idTrack[iTr];
821  vvTr[nElements] = &vTracks[idTrack[iTr]];
822 #ifdef NonhomogeniousField
823  int entrSIMD = idTrack[iTr] % fvecLen;
824  int entrVec = idTrack[iTr] / fvecLen;
825  field.SetOneEntry(nElements,vField[entrVec],entrSIMD);
826 #endif
827  nElements++;
828  }
829  else if( (iTr != idTrack.size()-1) ) continue;
830 
831  if( (nElements == fvecLen) || ((iTr == idTrack.size()-1)&&(nElements>0)) )
832  {
833  vDaughters[1].Create(vvTr,nElements,0,&DaughterPDG);
834 #ifdef NonhomogeniousField
835  vDaughters[1].SetField(field);
836 #endif
837  vDaughters[1].SetId(trId);
838 
839  if(isPrimary)
840  {
841  fvec errGuess[3] = {100*sqrt(PrimVtx.CovarianceMatrix()[0]),
842  100*sqrt(PrimVtx.CovarianceMatrix()[2]),
843  100*sqrt(PrimVtx.CovarianceMatrix()[5])};
844  hyperon.SetVtxGuess(PrimVtx.X(), PrimVtx.Y(), PrimVtx.Z());
845  hyperon.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
846  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
847  hyperon.Construct(vDaughtersPointer, 2, 0, -1, 0, 1);
848  }
849  else
850  {
851  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
852  hyperon.Construct(vDaughtersPointer, 2, 0);
853  }
854 
855  KFParticleSIMD hyperonTopo(hyperon);
856  hyperonTopo.SetProductionVertex(PrimVtx);
857 
858  for(unsigned int iv=0; iv<nElements; iv++)
859  {
860  bool isSameTrack = 0;
861  for(unsigned short iD=0; iD<vV0[iV0].DaughterIds().size(); iD++)
862  if(vV0[iV0].DaughterIds()[iD] == trId[iv]) isSameTrack=1;
863 
864  if(isSameTrack) continue;
865  if(!finite(hyperon.GetChi2()[iv])) continue;
866  if(!(hyperon.GetChi2()[iv] > 0.0f)) continue;
867  if(!(hyperon.GetChi2()[iv]==hyperon.GetChi2()[iv])) continue;
868 
869  fvec l, dl, isParticleFromVertex;
870  hyperon.GetDistanceToVertexLine(PrimVtx, l, dl, &isParticleFromVertex);
871  if(!(isParticleFromVertex[iv])) continue;
872 
873 
874  if(((l/dl)[iv] < cuts[0]) ) continue;
875  if(!isPrimary)
876  {
877  fvec l1, dl1;
878  vDaughters[0].GetDistanceToVertexLine(hyperon, l1, dl1, &isParticleFromVertex);
879  if(!(isParticleFromVertex[iv])) continue;
880  }
881 
882  if(hyperonTopo.GetChi2()[iv]/hyperonTopo.GetNDF()[iv] > cuts[1] ) continue;
883 
884  if( hyperon.GetChi2()[iv]/hyperon.GetNDF()[iv] > cuts[2] ) continue;
885  hyperon.GetKFParticle(hyperon_temp, iv);
886 
887  hyperon_temp.SetId(Particles.size());
888  Particles.push_back(hyperon_temp);
889 
890  if(vHyperonPrim)
891  {
892  Double_t mass, errMass;
893 
894  hyperon_temp.GetMass(mass, errMass);
895  hyperon_temp.SetNonlinearMassConstraint(hyperonPrimMass);
896 
897  if( (fabs(mass - hyperonPrimMass)/hyperonPrimMassErr) <= 3 )
898  vHyperonPrim->push_back(hyperon_temp);
899  else
900  if( vHyperonSec )
901  vHyperonSec->push_back(hyperon_temp);
902  }
903  }
904  nElements=0;
905  }
906  }
907  }
908 }
909 
910 void KFParticleFinder::FindHyperons(int PDG,
911  KFParticleSIMD vDaughters[2],
912  vector<int>& daughterIds,
913  vector<KFParticle>& vLambdaSec,
914  vector<KFParticle>& vHyperon,
915  KFParticleSIMD& PrimVtx,
916  const float *cuts,
917  int startIndex)
918 {
919  KFParticle* lambdas[fvecLen];
920  int nLambdasSec = vLambdaSec.size();
921 
922  for(unsigned short iL=startIndex; iL < vLambdaSec.size(); iL += fvecLen)
923  {
924  unsigned int nEntries = (iL + fvecLen < nLambdasSec) ? fvecLen : (nLambdasSec - iL);
925 
926  for(unsigned short iv=0; iv<nEntries; iv++)
927  lambdas[iv] = &vLambdaSec[iL+iv];
928 
929  vDaughters[0] = KFParticleSIMD(lambdas,nEntries);
930 
931  KFParticleSIMD hyperon;
932 
933  hyperon.SetPDG( PDG );
934  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
935  hyperon.Construct(vDaughtersPointer, 2, 0);
936 
937  KFParticleSIMD hyperonTopo(hyperon);
938  hyperonTopo.SetProductionVertex(PrimVtx);
939 
940  for(unsigned int iv=0; iv<nEntries; iv++)
941  {
942  bool isSameTrack = 0;
943  for(unsigned short iD=0; iD<lambdas[iv]->DaughterIds().size(); iD++)
944  for(unsigned short iD0=0; iD0<daughterIds.size(); iD0++)
945  if(lambdas[iv]->DaughterIds()[iD] == daughterIds[iD0]) isSameTrack=1;
946 
947  if(isSameTrack) continue;
948  if(!finite(hyperon.GetChi2()[iv])) continue;
949  if(!(hyperon.GetChi2()[iv] > 0.0f)) continue;
950  if(!(hyperon.GetChi2()[iv]==hyperon.GetChi2()[iv])) continue;
951 
952  fvec l, dl, isParticleFromVertex;
953  fvec l1, dl1, l2, dl2;
954  hyperon.GetDistanceToVertexLine(PrimVtx, l, dl, &isParticleFromVertex);
955  vDaughters[0].GetDistanceToVertexLine(hyperon, l1, dl1);
956  vDaughters[1].GetDistanceToVertexLine(hyperon, l2, dl2);
957  if(!(isParticleFromVertex[iv])) continue;
958 
959  if(((l/dl)[iv] < cuts[0]) ) continue;
960  if(((l1 - l)[iv] < 0) || ((l2 - l)[iv] < 0)) continue;
961 
962  if(hyperonTopo.GetChi2()[iv]/hyperonTopo.GetNDF()[iv] > cuts[1] ) continue;
963 
964  if( hyperon.GetChi2()[iv]/hyperon.GetNDF()[iv] > cuts[2] ) continue;
965  KFParticle hyperon_temp;
966  hyperon.GetKFParticle(hyperon_temp, iv);
967  vHyperon.push_back(hyperon_temp);
968  }
969  }
970 }
971 
972 void KFParticleFinder::FindDMesLambdac(vector<KFPTrack>& vTracks,
973 #ifdef NonhomogeniousField
974  const vector<L1FieldRegion>& vField,
975 #endif
976  vector<KFParticle>& Particles,
977  const int DaughterPDG[5], //pi, K_b, pi_b, K, p
978  const int MotherPDG[8], //D0, D+, D0_4, Ds+, Lc
979  vector<short>* idTrack[5], //pi, K_b, pi_b, K, p
980  KFParticleSIMD& PrimVtx,
981  const float cuts[8][8],
982  vector<float> ChiToPrimVtx)
983 {
984  vector<short> id[5]; //pi, K_b, pi_b, K, p
985  for(int iId=0; iId<5; iId++)
986  {
987  for(unsigned int iTr=0; iTr<idTrack[iId]->size(); iTr++)
988  {
989  KFPTrack& kfTrack = vTracks[idTrack[iId]->at(iTr)];
990  float pt = kfTrack.GetPt();
991  if(pt<cuts[0][3]) continue;
992  if(ChiToPrimVtx[idTrack[iId]->at(iTr)] < cuts[0][0]) continue;
993  id[iId].push_back(idTrack[iId]->at(iTr));
994  }
995  }
996 
997  vector<KFParticle> kpi;
998  vector<KFParticle> kpipi;
999  vector<KFParticle> kpipipi;
1000  vector<KFParticle> kpik;
1001  vector<KFParticle> kpip;
1002 
1003  const float cutskpi[3] = {3., 3., -100.};
1004  Find2DaughterDecay(vTracks,
1005 #ifdef NonhomogeniousField
1006  vField,
1007 #endif
1008  kpi, DaughterPDG[0], DaughterPDG[1], MotherPDG[0],
1009  id[0], id[1], PrimVtx, cutskpi, 0);
1010 
1011  for(unsigned int iKPi=0; iKPi<kpi.size(); iKPi++)
1012  {
1013  unsigned short startPiIndex = kpi[iKPi].DaughterIds()[0]+1;
1014 
1015  CombineTrackPart(vTracks,
1016 #ifdef NonhomogeniousField
1017  vField,
1018 #endif
1019  kpipi,kpi[iKPi],DaughterPDG[0], MotherPDG[1], id[0], cuts[1], startPiIndex, 1);
1020  CombineTrackPart(vTracks,
1021 #ifdef NonhomogeniousField
1022  vField,
1023 #endif
1024  kpik ,kpi[iKPi],DaughterPDG[3], MotherPDG[3], id[3], cuts[3], 0, 1);
1025  CombineTrackPart(vTracks,
1026 #ifdef NonhomogeniousField
1027  vField,
1028 #endif
1029  kpip ,kpi[iKPi],DaughterPDG[4], MotherPDG[4], id[4], cuts[4], 0, 1);
1030  }
1031  for(unsigned int iKPiPi=0; iKPiPi<kpipi.size(); iKPiPi++)
1032  CombineTrackPart(vTracks,
1033 #ifdef NonhomogeniousField
1034  vField,
1035 #endif
1036  kpipipi,kpipi[iKPiPi],DaughterPDG[2], MotherPDG[2], id[2], cuts[2], 0, 1);
1037 
1038  SelectParticleCandidates(Particles, kpi, PrimVtx, cuts[0]);
1039  SelectParticleCandidates(Particles, kpipi, PrimVtx, cuts[1]);
1040  SelectParticleCandidates(Particles, kpipipi, PrimVtx, cuts[2]);
1041  SelectParticleCandidates(Particles, kpik, PrimVtx, cuts[3]);
1042  SelectParticleCandidates(Particles, kpip, PrimVtx, cuts[4]);
1043 
1044  vector<KFParticle> d0pi;
1045  vector<KFParticle> d2pi;
1046  vector<KFParticle> d4pi;
1047 
1048  for(unsigned int iKPiPi=0; iKPiPi<kpipi.size(); iKPiPi++)
1049  CombineTrackPart(vTracks,
1050 #ifdef NonhomogeniousField
1051  vField,
1052 #endif
1053  d0pi,kpipi[iKPiPi],DaughterPDG[2], MotherPDG[5], id[2], cuts[5], 0, 0);
1054 
1055  for(unsigned int iKPi=0; iKPi<kpi.size(); iKPi++)
1056  CombineTrackPart(vTracks,
1057 #ifdef NonhomogeniousField
1058  vField,
1059 #endif
1060  d2pi,kpi[iKPi],DaughterPDG[0], MotherPDG[6], id[0], cuts[6], 0, 0);
1061 
1062  for(unsigned int iKPiPiPi=0; iKPiPiPi<kpipipi.size(); iKPiPiPi++)
1063  CombineTrackPart(vTracks,
1064 #ifdef NonhomogeniousField
1065  vField,
1066 #endif
1067  d4pi,kpipipi[iKPiPiPi],DaughterPDG[0], MotherPDG[7], id[0], cuts[7], 0, 0);
1068 
1069  SelectParticleCandidates(Particles, d0pi, PrimVtx, cuts[5]);
1070  SelectParticleCandidates(Particles, d2pi, PrimVtx, cuts[6]);
1071  SelectParticleCandidates(Particles, d4pi, PrimVtx, cuts[7]);
1072 }
1073 
1074 void KFParticleFinder::CombineTrackPart(vector<KFPTrack>& vTracks,
1075 #ifdef NonhomogeniousField
1076  const vector<L1FieldRegion>& vField,
1077 #endif
1078  vector<KFParticle>& Particles,
1079  KFParticle& part,
1080  const int DaughterPDG,
1081  const int MotherPDG,
1082  vector<short>& id,
1083  const float* cuts,
1084  const unsigned short startIndex,
1085  const bool IsSamePart)
1086 {
1087  KFParticleSIMD vDaughters[2] = {KFParticleSIMD(part),KFParticleSIMD()};
1088 
1089  KFPTrack* vTr[fvecLen];
1090 
1091  fvec vtxGuess[3] = {part.GetX(), part.GetY(), part.GetZ()};
1092  fvec errGuess[3] = {fvec(10.f*sqrt(part.CovarianceMatrix()[0])),
1093  fvec(10.f*sqrt(part.CovarianceMatrix()[2])),
1094  fvec(10.f*sqrt(part.CovarianceMatrix()[5]))};
1095 
1096  const unsigned short NTr = id.size();
1097  for(unsigned short iTr=startIndex; iTr < NTr; iTr += fvecLen)
1098  {
1099  const unsigned short NTracks = (iTr + fvecLen < NTr) ? fvecLen : (NTr - iTr);
1100 
1101 #ifdef NonhomogeniousField
1102  L1FieldRegion trField;
1103 #endif
1104  for(unsigned short iv=0; iv<NTracks; iv++)
1105  {
1106  vTr[iv] = &vTracks[id[iTr+iv]];
1107 #ifdef NonhomogeniousField
1108  int entrSIMDPos = id[iTr+iv] % fvecLen;
1109  int entrVecPos = id[iTr+iv] / fvecLen;
1110  trField.SetOneEntry(iv,vField[entrVecPos],entrSIMDPos);
1111 #endif
1112  }
1113 
1114  vDaughters[1].Create(vTr,NTracks,0,&DaughterPDG);
1115 #ifdef NonhomogeniousField
1116  vDaughters[1].SetField(trField);
1117 #endif
1118  fvec trId;
1119  for(int iv=0; iv<NTracks; iv++)
1120  trId[iv] = id[iTr+iv];
1121  vDaughters[1].SetId(trId);
1122 
1123  KFParticleSIMD mother;
1124  mother.SetPDG( MotherPDG );
1125  if(IsSamePart)
1126  {
1127  mother.SetVtxGuess(vtxGuess[0], vtxGuess[1], vtxGuess[2]);
1128  mother.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
1129  }
1130  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
1131  mother.Construct(vDaughtersPointer, 2, 0);
1132 
1133  for(int iv=0; iv<NTracks; iv++)
1134  {
1135  if(!finite(mother.GetChi2()[iv])) continue;
1136  if(!(mother.GetChi2()[iv] > 0.0f)) continue;
1137  if(!(mother.GetChi2()[iv]==mother.GetChi2()[iv])) continue;
1138 
1139  if( mother.GetChi2()[iv]/mother.GetNDF()[iv] > cuts[1] ) continue;
1140 
1141  bool isSameTrack = 0;
1142  for(unsigned short iD=0; iD<part.DaughterIds().size(); iD++)
1143  if(part.DaughterIds()[iD] == id[iTr+iv]) isSameTrack=1;
1144  if(isSameTrack) continue;
1145 
1146  KFParticle mother_temp;
1147  mother.GetKFParticle(mother_temp, iv);
1148  mother_temp.CleanDaughtersId();
1149  mother_temp.AddDaughter(id[iTr+iv]);
1150  for(unsigned short iD=0; iD<part.DaughterIds().size(); iD++)
1151  mother_temp.AddDaughter(part.DaughterIds()[iD]);
1152  Particles.push_back(mother_temp);
1153  }
1154  }
1155 }
1156 
1157 void KFParticleFinder::SelectParticleCandidates(vector<KFParticle>& Particles,
1158  vector<KFParticle>& vCandidates,
1159  KFParticleSIMD& PrimVtx,
1160  const float cuts[5])
1161 {
1162  KFParticle* cand[fvecLen];
1163  int nCand = vCandidates.size();
1164 
1165  for(unsigned short iC=0; iC < nCand; iC += fvecLen)
1166  {
1167  unsigned int nEntries = (iC + fvecLen < nCand) ? fvecLen : (nCand - iC);
1168 
1169  for(unsigned short iv=0; iv<nEntries; iv++)
1170  cand[iv] = &vCandidates[iC+iv];
1171 
1172  KFParticleSIMD candTopo(cand,nEntries);
1173 
1174  candTopo.SetProductionVertex(PrimVtx);
1175 
1176  for(unsigned int iv=0; iv<nEntries; iv++)
1177  {
1178  if(!finite(candTopo.GetChi2()[iv])) continue;
1179  if(!(candTopo.GetChi2()[iv] > 0.0f)) continue;
1180  if(!(candTopo.GetChi2()[iv]==candTopo.GetChi2()[iv])) continue;
1181 
1182  fvec l, dl, isParticleFromVertex;
1183  candTopo.GetDistanceToVertexLine(PrimVtx, l, dl, &isParticleFromVertex);
1184  if(!(isParticleFromVertex[iv])) continue;
1185  if(((l/dl)[iv] < cuts[2]) ) continue;
1186 
1187  if(candTopo.GetChi2()[iv]/candTopo.GetNDF()[iv] > cuts[4] ) continue;
1188 
1189  Particles.push_back(vCandidates[iC+iv]);
1190  }
1191  }
1192 }
1193 
1194 
1195 
1196 #undef cnst
c5
c5
Definition: plot_dirc.C:75
KFParticleSIMD::X
const fvec & X() const
Definition: KFParticleSIMD.h:132
p
Double_t p
Definition: anasim.C:58
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
KFPTrack::Id
int Id() const
Definition: KFPTrack.h:108
KFParticleSIMD::SetField
static void SetField(fvec Bz)
Definition: KFParticleSIMD.h:410
KFParticleFinder::ExtrapolateToPV
static void ExtrapolateToPV(std::vector< KFParticle > &vParticles, KFParticleSIMD &PrimVtx)
Definition: KFParticleFinder.cxx:508
KFParticleSIMD::GetChi2
fvec GetChi2() const
Definition: KFParticleSIMD.h:496
KFParticleSIMD::Construct
void Construct(const KFParticleSIMD *vDaughters[], int nDaughters, const KFParticleSIMD *ProdVtx=0, Float_t Mass=-1, Bool_t IsConstrained=0, Bool_t isAtVtxGuess=0)
Definition: KFParticleSIMD.h:880
c4
c4
Definition: plot_dirc.C:71
KFParticleFinder::FindTrackV0Decay
static void FindTrackV0Decay(const int MotherPDG, std::vector< KFParticle > &Particles, std::vector< KFParticle > &vV0, std::vector< KFPTrack > &vTracks, const int DaughterPDG, std::vector< short > &idTrack, KFParticleSIMD &PrimVtx, const float *cuts=0, bool isPrimary=0, std::vector< float > *ChiToPrimVtx=0, std::vector< KFParticle > *vHyperonPrim=0, float hyperonPrimMass=0, float hyperonPrimMassErr=0, std::vector< KFParticle > *vHyperonSec=0)
Definition: KFParticleFinder.cxx:779
KFParticleBase::SetNonlinearMassConstraint
void SetNonlinearMassConstraint(Double_t Mass)
Definition: KFParticleBase.cxx:1393
KFPTrack::GetPt
float GetPt() const
Definition: KFPTrack.h:51
KFParticle::GetX
Double_t GetX() const
Definition: KFParticle.h:447
map
PndTransMap * map
Definition: sim_emc_apd.C:99
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
KFParticleSIMD::Y
const fvec & Y() const
Definition: KFParticleSIMD.h:133
KFParticleBaseSIMD::SetVtxErrGuess
void SetVtxErrGuess(fvec &x, fvec &y, fvec &z)
Definition: KFParticleBaseSIMD.cxx:115
push_back
labels push_back("electron")
KFParticle::CovarianceMatrix
Double_t * CovarianceMatrix()
Definition: KFParticle.h:824
KFParticleBaseSIMD::GetDistanceToVertexLine
void GetDistanceToVertexLine(const KFParticleBaseSIMD &Vertex, fvec &l, fvec &dl, fvec *isParticleFromVertex=0) const
Definition: KFParticleBaseSIMD.cxx:1632
KFParticleFinder::FindHyperons
static void FindHyperons(int PDG, KFParticleSIMD vDaughters[2], std::vector< int > &daughterIds, std::vector< KFParticle > &vLambdaSec, std::vector< KFParticle > &vHyperon, KFParticleSIMD &PrimVtx, const float *cuts=0, int startIndex=0)
Definition: KFParticleFinder.cxx:910
KFParticleSIMD::SetVtxGuess
void SetVtxGuess(fvec x, fvec y, fvec z)
Definition: KFParticleSIMD.h:446
c2
c2
Definition: plot_dirc.C:39
cuts
TString cuts[MAX]
Definition: autocutx.C:35
KFPTrack::Charge
int Charge() const
Definition: KFPTrack.h:62
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
KFParticleFinder::SelectParticleCandidates
static void SelectParticleCandidates(std::vector< KFParticle > &Particles, std::vector< KFParticle > &vCandidates, KFParticleSIMD &PrimVtx, const float cuts[5])
Definition: KFParticleFinder.cxx:1157
KFParticleFinder::GetChi2BetweenParticles
static fvec GetChi2BetweenParticles(KFParticleSIMD &p1, KFParticleSIMD &p2)
Definition: KFParticleFinder.cxx:531
KFParticleSIMD::GetNDF
fvec GetNDF() const
Definition: KFParticleSIMD.h:501
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
KFPTrack::GetTrack
const float * GetTrack() const
Definition: KFPTrack.h:67
KFParticleFinder::FindDMesLambdac
static void FindDMesLambdac(std::vector< KFPTrack > &vTracks, std::vector< KFParticle > &Particles, const int DaughterPDG[5], const int MotherPDG[8], std::vector< short > *idTrack[5], KFParticleSIMD &PrimVtx, const float cuts[8][8], std::vector< float > ChiToPrimVtx)
Definition: KFParticleFinder.cxx:972
KFParticleSIMD::CovarianceMatrix
fvec * CovarianceMatrix()
Definition: KFParticleSIMD.h:848
KFParticleSIMD::GetCovariance
fvec GetCovariance(int i) const
Definition: KFParticleSIMD.h:511
KFPTrack::GetP
float GetP() const
Definition: KFPTrack.h:52
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
ok
float_m ok
Definition: PndFTSCAGBTracker.cxx:820
KFParticleSIMD::TransportToPoint
void TransportToPoint(const fvec xyz[])
Definition: KFParticleSIMD.h:920
f
TFile * f
Definition: bump_analys.C:12
KFParticleSIMD::Z
const fvec & Z() const
Definition: KFParticleSIMD.h:134
fvec
F32vec4 fvec
Definition: P4_F32vec4.h:218
KFParticle::SetNDaughters
void SetNDaughters(int n)
Definition: KFParticle.h:91
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
c1
c1
Definition: plot_dirc.C:35
KFParticleBase::SetPDG
void SetPDG(int pdg)
Definition: KFParticleBase.h:276
c3
c3
Definition: plot_dirc.C:50
p2
TPad * p2
Definition: hist-t7.C:117
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
KFParticleFinder::DefaultCuts
static const float DefaultCuts[2][3]
Definition: KFParticleFinder.h:141
KFParticleBase::AddDaughterId
void AddDaughterId(int id)
Definition: KFParticleBase.h:274
fvecLen
const int fvecLen
Definition: P4_F32vec4.h:220
KFParticleFinder::FindParticles
static void FindParticles(std::vector< KFPTrack > &vRTracks, std::vector< float > &ChiToPrimVtx, std::vector< KFParticle > &Particles, KFParticleSIMD &PrimVtx, const std::vector< int > &vTrackPDG, const float cuts[2][3]=DefaultCuts)
Definition: KFParticleFinder.cxx:20
KFParticleSIMD::GetZ
fvec GetZ() const
Definition: KFParticleSIMD.h:461
lambdas
TCut lambdas
Definition: plotDiaBnew.C:15
KFParticleFinder::Find2DaughterDecay
static void Find2DaughterDecay(std::vector< KFPTrack > &vTracks, std::vector< KFParticle > &Particles, const int DaughterNegPDG, const int DaughterPosPDG, const int MotherPDG, std::vector< short > &idNeg, std::vector< short > &idPos, KFParticleSIMD &PrimVtx, const float *cuts=0, bool isPrimary=0, std::vector< float > *vMotherTopoChi2Ndf=0, const float *secCuts=0, const float massMotherPDG=0, const float massMotherPDGSigma=0, std::vector< KFParticle > *vMotherPrim=0, std::vector< KFParticle > *vMotherSec=0)
Definition: KFParticleFinder.cxx:558
KFParticleFinder::CombineTrackPart
static void CombineTrackPart(std::vector< KFPTrack > &vTracks, std::vector< KFParticle > &Particles, KFParticle &part, const int DaughterPDG, const int MotherPDG, std::vector< short > &id, const float *cuts, const unsigned short startIndex=0, const bool IsSamePart=0)
Definition: KFParticleFinder.cxx:1074
KFParticleSIMD::GetX
fvec GetX() const
Definition: KFParticleSIMD.h:451
KFParticleBase::DaughterIds
const std::vector< int > & DaughterIds() const
Definition: KFParticleBase.h:271
p1
TPad * p1
Definition: hist-t7.C:116
KFParticle::GetZ
Double_t GetZ() const
Definition: KFParticle.h:457
KFPTrack::GetCovMatrix
const float * GetCovMatrix() const
Definition: KFPTrack.h:68
KFParticleSIMD::GetKFParticle
void GetKFParticle(KFParticle &Part, int iPart=0)
Definition: KFParticleSIMD.cxx:819
KFParticle::GetY
Double_t GetY() const
Definition: KFParticle.h:452
KFParticle::GetMass
Double_t GetMass() const
Definition: KFParticle.h:553
KFParticle::AddDaughter
void AddDaughter(int id)
Definition: KFParticle.h:92
KFParticleSIMD::Parameters
fvec * Parameters()
Definition: KFParticleSIMD.h:843
KFPTrack::GetNDF
int GetNDF() const
Definition: KFPTrack.h:65
KFParticleBase::SetId
void SetId(int id)
Definition: KFParticleBase.h:273
KFParticleSIMD::Create
void Create(const fvec Param[], const fvec Cov[], fvec Charge, fvec mass)
Definition: KFParticleSIMD.cxx:50
KFPTrack::GetChi2
float GetChi2() const
Definition: KFPTrack.h:64
r2
double r2
Definition: RiemannTrackFinderTest.C:19
KFParticle::CleanDaughtersId
void CleanDaughtersId()
Definition: KFParticle.h:90
L1FieldRegion::SetOneEntry
void SetOneEntry(const int i0, const L1FieldRegion &f1, const int i1)
Definition: L1Field.h:211
KFParticleSIMD::GetY
fvec GetY() const
Definition: KFParticleSIMD.h:456
KFParticleSIMD::SetProductionVertex
void SetProductionVertex(const KFParticleSIMD &Vtx)
Definition: KFParticleSIMD.h:865
Generated on Wed Apr 3 2019 10:02:47 for FairRoot/PandaRoot by   doxygen 1.8.5