Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFParticleFinder.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFParticleFinder.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 "KFParticleFinder.h"
23 
24 using std::map;
25 using std::vector;
26 
27 #include "KFParticleDatabase.h"
28 #include "KFPEmcCluster.h"
29 
30 KFParticleFinder::KFParticleFinder():
31  fNPV(-1),fNThreads(1),fDistanceCut(1.f),fLCut(-5.f),fCutCharmPt(0.2f),fCutCharmChiPrim(85.f),fCutLVMPt(0.0f),fCutLVMP(0.0f),fCutJPsiPt(1.0f),
32  fD0(0), fD0bar(0), fD04(0), fD04bar(0), fD0KK(0), fD0pipi(0), fDPlus(0), fDMinus(0),
33  fDPlus3Pi(0), fDMinus3Pi(0), fDsPlusK2Pi(0), fDsMinusK2Pi(0), fLcPlusP2Pi(0), fLcMinusP2Pi(0),
34  fLPi(0), fLPiPIndex(0), fHe3Pi(0), fHe3PiBar(0), fHe4Pi(0), fHe4PiBar(0),
35  fHe4L(0), fHe5L(0), fLLn(0), fH5LL(0),
36  fSecCandidates(), fPrimCandidates(), fPrimCandidatesTopo(),fPrimCandidatesTopoMass(),
37  fEmcClusters(0), fMixedEventAnalysis(0), fDecayReconstructionList()
38 {
40  //Cuts
41  //track + track
42  //chi2_prim chi2_geo ldl
43  fCuts2D[0] = 3.f; fCuts2D[1] = 3.f; fCuts2D[2] = 5.f;
44  //cuts to select primary and secondary particles
45  //mass chi2_topo ldl
46 #ifdef PANDA_STT
47  fSecCuts[0] = 3.f; fSecCuts[1] = -3.f; fSecCuts[2] = 10.f;
48 #else
49  fSecCuts[0] = 3.f; fSecCuts[1] = 5.f; fSecCuts[2] = 10.f;
50 #endif
51 
52 #ifdef __ROOT__
53  fCutCharmChiPrim = 8;
54 #endif
55 
56  //track + particle
57  // ldl chi2_topo chi2_geo
58  fCutsTrackV0[0][0] = 5; fCutsTrackV0[0][1] = 5; fCutsTrackV0[0][2] = 6; //Xi, Omega
59  fCutsTrackV0[1][0] = 5; fCutsTrackV0[1][1] = 5; fCutsTrackV0[1][2] = 6; //Charm, H0, Sigma+
60  fCutsTrackV0[2][0] = -100.; fCutsTrackV0[2][1] = 10000.; fCutsTrackV0[2][2] = 3; //resonances
61 
62  //charm
63  //chi2 l/dl chi2_topo
64  fCutsCharm[0] = 3.f; fCutsCharm[1] = 10.f; fCutsCharm[2] = 3.f; //D0 -> pi+ K-
65 
66  //cuts on particles reconstructed from short-lived particles
67  //ldl, chi2_topo chi2_geo
68  //H0 -> Lambda Lambda, Xi0 -> Lambda pi0
69  fCutsPartPart[0][0] = 10; fCutsPartPart[0][1] = 3; fCutsPartPart[0][2] = 3;
70  //Sigma0 -> Lambda Gamma, pi0 -> Gamma Gamma, K* -> K pi0, Sigma*0 -> Lambda pi0, Xi* -> Xi pi0
71  fCutsPartPart[1][0] = -10; fCutsPartPart[1][1] = 3; fCutsPartPart[1][2] = 3;
72 }
73 
74 //________________________________________________________________________________
75 void KFParticleFinder::Init(int nPV)
76 {
82  fNPV = nPV;
83 // Particles.reserve(vRTracks.size() + nPart);
84 
85  fD0.clear();
86  fD0bar.clear();
87  fD04.clear();
88  fD04bar.clear();
89  fD0KK.clear();
90  fD0pipi.clear();
91  fDPlus.clear();
92  fDMinus.clear();
93  fDPlus3Pi.clear();
94  fDMinus3Pi.clear();
95  fDsPlusK2Pi.clear();
96  fDsMinusK2Pi.clear();
97  fLcPlusP2Pi.clear();
98  fLcMinusP2Pi.clear();
99  fLPi.clear();
100  fLPiPIndex.clear();
101  fHe3Pi.clear();
102  fHe3PiBar.clear();
103  fHe4Pi.clear();
104  fHe4PiBar.clear();
105  fHe4L.clear();
106  fHe5L.clear();
107  fLLn.clear();
108  fH5LL.clear();
109 
110  for(int iCandidates=0; iCandidates<fNSecCandidatesSets; iCandidates++)
111  fSecCandidates[iCandidates].clear();
112 
113  for(int iCandidates=0; iCandidates<fNPrimCandidatesSets; iCandidates++)
114  {
115  fPrimCandidates[iCandidates].clear();
116  fPrimCandidates[iCandidates].resize(fNPV);
117  }
118 
119  for(int iCandidates=0; iCandidates<fNPrimCandidatesTopoSets; iCandidates++)
120  {
121  fPrimCandidatesTopo[iCandidates].clear();
122  fPrimCandidatesTopo[iCandidates].resize(fNPV);
123 
124  fPrimCandidatesTopoMass[iCandidates].clear();
125  fPrimCandidatesTopoMass[iCandidates].resize(fNPV);
126  }
127 
128 }
129 //________________________________________________________________________________
130 
131 void KFParticleFinder::FindParticles(KFPTrackVector* vRTracks, kfvector_float* ChiToPrimVtx, std::vector<KFParticle>& Particles,
132  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx, int nPV)
133 {
159  Init(nPV);
160  const int nPartPrim = vRTracks[2].NPions() * vRTracks[3].NKaons() +
161  vRTracks[3].NPions() * vRTracks[2].NKaons() +
162  vRTracks[2].NKaons() * vRTracks[3].NKaons() +
163  vRTracks[2].NKaons() * vRTracks[3].NProtons() +
164  vRTracks[3].NKaons() * vRTracks[2].NProtons() +
165  vRTracks[2].NElectrons() * vRTracks[3].NElectrons() +
166  vRTracks[2].NMuons() * vRTracks[3].NMuons();
167 
168  const int nPart = vRTracks[0].NPions() * vRTracks[1].NPions() +
169  vRTracks[0].NPions() * vRTracks[1].NProtons() +
170  vRTracks[1].NPions() * vRTracks[0].NProtons() + nPartPrim;
171  int nEmcClusters = 0;
172  if(fEmcClusters)
173  nEmcClusters = fEmcClusters->Size();
174  vector<KFParticle> vGammaPrimEmc;
175 
176  int nPartEstimation = nPart+vRTracks[0].Size()+vRTracks[1].Size()+vRTracks[2].Size()+vRTracks[3].Size() + nEmcClusters;
177 
178  if(nPartEstimation < 100000)
179  Particles.reserve(nPartEstimation);
180  //* Finds particles (K0s and Lambda) from a given set of tracks
181  {
182  KFPTrack kfTrack;
183  for(int iV=0; iV<4; iV++)
184  {
185  for(int iTr=0; iTr < vRTracks[iV].Size(); iTr++)
186  {
187  vRTracks[iV].GetTrack(kfTrack, iTr);
188  int pdg = vRTracks[iV].PDG()[iTr];
189  if( pdg == 19 ) pdg = 13;
190  if( pdg ==-19 ) pdg = -13;
191  KFParticle tmp(kfTrack, pdg);
192  tmp.SetPDG(pdg);
193  tmp.SetId(Particles.size());
194  vRTracks[iV].SetId(Particles.size(),iTr);
195  if(vRTracks[iV+4].Size() > 0)
196  vRTracks[iV+4].SetId(Particles.size(),iTr);
197  tmp.AddDaughterId( kfTrack.Id() );
198 #ifdef NonhomogeneousField
199  for(int iF=0; iF<10; iF++)
200  tmp.SetFieldCoeff( vRTracks[iV].FieldCoefficient(iF)[iTr], iF);
201 #endif
202  Particles.push_back(tmp);
203  }
204  }
205 
206  if(fEmcClusters)
207  {
208  KFParticleSIMD tmpGammaSIMD;
209  KFParticle tmpGamma;
210 
211  for(int iEmc=0; iEmc < fEmcClusters->Size(); iEmc += float_vLen)
212  {
213  const int NClustersVec = (iEmc + float_vLen < fEmcClusters->Size()) ? float_vLen : (fEmcClusters->Size() - iEmc);
214  tmpGammaSIMD.Load(*fEmcClusters, iEmc, PrimVtx[0]);
215  for(int iV=0; iV<NClustersVec; iV++)
216  {
217  tmpGammaSIMD.GetKFParticle(tmpGamma, iV);
218  tmpGamma.SetPDG(22); //gamma pdg
219  tmpGamma.SetId(Particles.size());
220  tmpGamma.CleanDaughtersId();
221  tmpGamma.AddDaughterId(fEmcClusters->Id()[iEmc+iV]);
222  Particles.push_back(tmpGamma);
223  vGammaPrimEmc.push_back(tmpGamma);
224  }
225  }
226  }
227  }
228 
229 
230 
231  Find2DaughterDecay(vRTracks, ChiToPrimVtx,
232  Particles, PrimVtx, fCuts2D,
233  fSecCuts, fPrimCandidates, fSecCandidates);
234 
235  if(!fMixedEventAnalysis)
236  {
237  //Construct two-particle background from positive primary tracks for subtraction from the resonance spectra
238  ConstructPrimaryBG(vRTracks, Particles, PrimVtx, fCuts2D, fSecCuts, fPrimCandidates, fSecCandidates);
239 
240  for(int iPV=0; iPV<fNPV; iPV++ )
241  {
242  ExtrapolateToPV(fPrimCandidates[0][iPV],PrimVtx[iPV]);
243  ExtrapolateToPV(fPrimCandidates[1][iPV],PrimVtx[iPV]);
244  ExtrapolateToPV(fPrimCandidates[2][iPV],PrimVtx[iPV]);
245  ExtrapolateToPV(fPrimCandidates[3][iPV],PrimVtx[iPV]);
246  }
247 
248  NeutralDaughterDecay(vRTracks, Particles);
249 
250  //Xi- -> Lambda pi-, Omega- -> Lambda K-
251  FindTrackV0Decay(fSecCandidates[1], 3122, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastKaon(),
252  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]), &fPrimCandidates[5]);
253  //Xi+ -> Lambda pi+, Omega+ -> Lambda K+
254  FindTrackV0Decay(fSecCandidates[2], -3122, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastKaon(),
255  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]), &fPrimCandidates[6]);
256 
257  for(int iPV=0; iPV<fNPV; iPV++ )
258  {
259  ExtrapolateToPV(fPrimCandidates[5][iPV],PrimVtx[iPV]);
260  ExtrapolateToPV(fPrimCandidates[6][iPV],PrimVtx[iPV]);
261 
262  ExtrapolateToPV(fPrimCandidates[7][iPV],PrimVtx[iPV]);
263  ExtrapolateToPV(fPrimCandidates[8][iPV],PrimVtx[iPV]);
264  }
265 
266  //K*+ -> K0 pi+
267  for(int iPV=0; iPV<fNPV; iPV++)
268  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
269  Particles, PrimVtx, iPV, 0);
270  //K*- -> K0 pi-
271  for(int iPV=0; iPV<fNPV; iPV++)
272  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
273  Particles, PrimVtx, iPV, 0);
274  //Sigma*+ -> Lambda pi+
275  for(int iPV=0; iPV<fNPV; iPV++)
276  FindTrackV0Decay(fPrimCandidates[1][iPV], 3122, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
277  Particles, PrimVtx, iPV, 0);
278  //Sigma*- -> Lambda pi-, Xi*- -> Lambda K-
279  for(int iPV=0; iPV<fNPV; iPV++)
280  FindTrackV0Decay(fPrimCandidates[1][iPV], 3122, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastKaon(),
281  Particles, PrimVtx, iPV, 0);
282  //Sigma*+_bar -> Lambda_bar pi-
283  for(int iPV=0; iPV<fNPV; iPV++)
284  FindTrackV0Decay(fPrimCandidates[2][iPV], -3122, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
285  Particles, PrimVtx, iPV, 0);
286  //Sigma*-_bar -> Lambda_bar pi+, Xi*+ -> Lambda_bar + K+
287  for(int iPV=0; iPV<fNPV; iPV++)
288  FindTrackV0Decay(fPrimCandidates[2][iPV], -3122, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastKaon(),
289  Particles, PrimVtx, iPV, 0);
290  //Xi*0 -> Xi- pi+
291  for(int iPV=0; iPV<fNPV; iPV++)
292  FindTrackV0Decay(fPrimCandidates[5][iPV], 3312, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
293  Particles, PrimVtx, iPV, 0, &fPrimCandidates[9]);
294  //Xi*0_bar -> Xi+ pi-
295  for(int iPV=0; iPV<fNPV; iPV++)
296  FindTrackV0Decay(fPrimCandidates[6][iPV], -3312, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
297  Particles, PrimVtx, iPV, 0, &fPrimCandidates[10]);
298  //Omega*- -> Xi- pi+ K-
299  for(int iPV=0; iPV<fNPV; iPV++)
300  FindTrackV0Decay(fPrimCandidates[9][iPV], 3324, vRTracks[3], -1, vRTracks[3].FirstKaon(), vRTracks[3].LastKaon(),
301  Particles, PrimVtx, iPV, 0);
302  //Omega*+ -> Xi+ pi- K+
303  for(int iPV=0; iPV<fNPV; iPV++)
304  FindTrackV0Decay(fPrimCandidates[10][iPV], -3324, vRTracks[2], 1, vRTracks[2].FirstKaon(), vRTracks[2].LastKaon(),
305  Particles, PrimVtx, iPV, 0);
306  //Hypernuclei
307  //He4L -> He3 p pi-
308  FindTrackV0Decay(fHe3Pi , 3004, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(), Particles, PrimVtx, -1, 0, 0, &fHe4L);
309  //LLn -> H3L pi-
310  FindTrackV0Decay(fHe3Pi , 3004, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(), Particles, PrimVtx, -1, 0 );
311  //He4L_bar -> He3_bar p- pi+
312  FindTrackV0Decay(fHe3PiBar,-3004, vRTracks[1], -1, vRTracks[1].FirstProton(), vRTracks[1].LastProton(), Particles, PrimVtx, -1, 0);
313  //He5L -> He4 p pi-
314  FindTrackV0Decay(fHe4Pi , 3005, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(), Particles, PrimVtx, -1, 0, 0, &fHe5L);
315  //He5L_bar -> He4_bar p- pi+
316  FindTrackV0Decay(fHe4PiBar,-3005, vRTracks[1], -1, vRTracks[1].FirstProton(), vRTracks[1].LastProton(), Particles, PrimVtx, -1, 0);
317  //H4LL -> He4L pi-
318  FindTrackV0Decay(fHe4L , 3006, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(), Particles, PrimVtx, -1, 0 );
319  //H5LL -> He5L pi-
320  FindTrackV0Decay(fHe5L , 3007, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(), Particles, PrimVtx, -1, 0 );
321  //H4LL -> H3L p pi-
322  FindTrackV0Decay(fLLn , 3203, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(), Particles, PrimVtx, -1, 0 );
323  //He6LL -> He5L p pi-
324  FindTrackV0Decay(fH5LL , 3010, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(), Particles, PrimVtx, -1, 0 );
325  // Charm
326  //LambdaC -> pi+ K- p, Ds+ -> pi+ K- K+, D+ -> pi+ K- pi+
327  FindTrackV0Decay(fD0, 421, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastProton(),
328  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]));
329  //LambdaC_bar -> pi- K+ p-, Ds- -> pi- K+ K-, D- -> pi- K+ pi-
330  FindTrackV0Decay(fD0bar, -421, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastProton(),
331  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]));
332  //D0->pi+ K- pi+ pi-
333  FindTrackV0Decay(fDPlus, 411, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(),
334  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]));
335  //D0_bar->pi- K+ pi- pi+
336  FindTrackV0Decay(fDMinus, -411, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastPion(),
337  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]));
338  //B+ -> D0_bar pi+, B+ -> D0_bar K+
339  FindTrackV0Decay(fD0bar, -421, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastKaon(),
340  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]));
341  //B- -> D0 pi-, B- -> D0 K-
342  FindTrackV0Decay(fD0, 421, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastKaon(),
343  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]));
344  //B0 -> D- pi+, B0 -> D- K+
345  FindTrackV0Decay(fDMinus, -419, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastKaon(),
346  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]));
347  //B0_bar -> D+ pi-, B0_bar -> D+ K-
348  FindTrackV0Decay(fDPlus, 419, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastKaon(),
349  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]));
350  //D0 -> pi+ K-
351  SelectParticles(Particles,fD0,PrimVtx,fCutsCharm[2],fCutsCharm[1],
352  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
353  //D0_bar -> pi+ K-
354  SelectParticles(Particles,fD0bar,PrimVtx,fCutsCharm[2],fCutsCharm[1],
355  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
356  //D*+->D0 pi+
357  for(int iPV=0; iPV<fNPV; iPV++)
358  FindTrackV0Decay(fD0, 421, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
359  Particles, PrimVtx, iPV, 0);
360  //D*- -> D0_bar pi-
361  for(int iPV=0; iPV<fNPV; iPV++)
362  FindTrackV0Decay(fD0, -421, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
363  Particles, PrimVtx, iPV, 0);
364  //D0 -> pi+ K- pi+ pi-
365  SelectParticles(Particles,fD04,PrimVtx,fCutsCharm[2],fCutsCharm[1],
366  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
367  //D0_bar -> pi- K+ pi- pi+
368  SelectParticles(Particles,fD04bar,PrimVtx,fCutsCharm[2],fCutsCharm[1],
369  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
370  //D*+->D0 pi+
371  for(int iPV=0; iPV<fNPV; iPV++)
372  FindTrackV0Decay(fD04, 429, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
373  Particles, PrimVtx, iPV, 0);
374  //D0*- -> D0_bar pi-
375  for(int iPV=0; iPV<fNPV; iPV++)
376  FindTrackV0Decay(fD04bar, -429, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
377  Particles, PrimVtx, iPV, 0);
378  //D+
379  SelectParticles(Particles,fDPlus,PrimVtx,fCutsCharm[2],fCutsCharm[1],
380  KFParticleDatabase::Instance()->GetDPlusMass(), KFParticleDatabase::Instance()->GetDPlusMassSigma(), fSecCuts[0]);
381  //D-
382  SelectParticles(Particles,fDMinus,PrimVtx,fCutsCharm[2],fCutsCharm[1],
383  KFParticleDatabase::Instance()->GetDPlusMass(), KFParticleDatabase::Instance()->GetDPlusMassSigma(), fSecCuts[0]);
384  //D*0->D+ pi-
385  for(int iPV=0; iPV<fNPV; iPV++)
386  FindTrackV0Decay(fDPlus, 411, vRTracks[3], -1, vRTracks[3].FirstPion(), vRTracks[3].LastPion(),
387  Particles, PrimVtx, iPV, 0);
388  //D*0_bar->D- pi+
389  for(int iPV=0; iPV<fNPV; iPV++)
390  FindTrackV0Decay(fDMinus, -411, vRTracks[2], 1, vRTracks[2].FirstPion(), vRTracks[2].LastPion(),
391  Particles, PrimVtx, iPV, 0);
392 
393  float cutsD0[3] = {fCutsCharm[1], fCutsCharm[2], fCutsCharm[0]};
394 // float cutsD0[3] = {-5, 1e10, 1e10};
395  //D0 -> K0 pi+ pi-
396  for(int iPV=0; iPV<fNPV; iPV++)
397  CombinePartPart(fD0pipi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, 425, 0, 1);
398  //D0 -> K0 K+ K-
399  for(int iPV=0; iPV<fNPV; iPV++)
400  CombinePartPart(fD0KK, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, 427, 0, 1);
401 
402 
403  //LambdaC -> p pi+ pi-, Ds+ -> K+ pi+ pi-, D+ -> pi+ pi+ pi-
404  FindTrackV0Decay(fD0pipi, 420, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastProton(),
405  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]));
406  //LambdaC_bar -> p_bar pi+ pi-, Ds- -> K- pi+ pi-, D- -> pi+ pi- pi-
407  FindTrackV0Decay(fD0pipi, 420, vRTracks[1],-1, vRTracks[1].FirstPion(), vRTracks[1].LastProton(),
408  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]));
409 
410  //D+ -> K0 pi+ pi+ pi-
411  for(int iPV=0; iPV<fNPV; iPV++)
412  CombinePartPart(fDPlus3Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, 200411, 0, 1);
413  //D- -> K0 pi+ pi- pi-
414  for(int iPV=0; iPV<fNPV; iPV++)
415  CombinePartPart(fDMinus3Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, -200411, 0, 1);
416  //Lc+ -> Lambda pi+ pi+ pi-
417  for(int iPV=0; iPV<fNPV; iPV++)
418  CombinePartPart(fDPlus3Pi, fPrimCandidates[1][iPV], Particles, PrimVtx, cutsD0, -1, 404122, 0, 1);
419  //Lc- -> Lambda_bar pi+ pi- pi-
420  for(int iPV=0; iPV<fNPV; iPV++)
421  CombinePartPart(fDMinus3Pi, fPrimCandidates[2][iPV], Particles, PrimVtx, cutsD0, -1, -404122, 0, 1);
422  //Xic0 -> Xi- pi+ pi+ pi-
423  for(int iPV=0; iPV<fNPV; iPV++)
424  CombinePartPart(fDPlus3Pi, fPrimCandidates[5][iPV], Particles, PrimVtx, cutsD0, -1, 4132, 0, 1);
425  //Xic0_bar -> Xi+ pi+ pi- pi-
426  for(int iPV=0; iPV<fNPV; iPV++)
427  CombinePartPart(fDMinus3Pi, fPrimCandidates[6][iPV], Particles, PrimVtx, cutsD0, -1, -4132, 0, 1);
428 
429  //Ds+ -> K0 K+ pi+ pi-
430  for(int iPV=0; iPV<fNPV; iPV++)
431  CombinePartPart(fDsPlusK2Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, 300431, 0, 1);
432  //Ds- -> K0 K- pi+ pi-
433  for(int iPV=0; iPV<fNPV; iPV++)
434  CombinePartPart(fDsMinusK2Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, -300431, 0, 1);
435 
436  //Lc+ -> p K0 pi+ pi-
437  for(int iPV=0; iPV<fNPV; iPV++)
438  CombinePartPart(fLcPlusP2Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, 204122, 0, 1);
439  //Lc- -> p- K0 pi+ pi-
440  for(int iPV=0; iPV<fNPV; iPV++)
441  CombinePartPart(fLcMinusP2Pi, fPrimCandidates[0][iPV], Particles, PrimVtx, cutsD0, -1, -204122, 0, 1);
442 
443  //D0 -> pi+ pi-
444  SelectParticles(Particles,fD0pipi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
445  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
446  //D0 -> K+ K-
447  SelectParticles(Particles,fD0KK,PrimVtx,fCutsCharm[2],fCutsCharm[1],
448  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
449  //D+ -> pi+ pi+ pi-
450  SelectParticles(Particles,fDPlus3Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
451  KFParticleDatabase::Instance()->GetDPlusMass(), KFParticleDatabase::Instance()->GetDPlusMassSigma(), fSecCuts[0]);
452  //D- -> pi+ pi- pi-
453  SelectParticles(Particles,fDMinus3Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
454  KFParticleDatabase::Instance()->GetDPlusMass(), KFParticleDatabase::Instance()->GetDPlusMassSigma(), fSecCuts[0]);
455  //Ds+ -> K+ pi+ pi-,
456  SelectParticles(Particles,fDsPlusK2Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
457  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
458  //Ds- -> K- pi+ pi-
459  SelectParticles(Particles,fDsMinusK2Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
460  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
461  //LambdaC -> p pi+ pi-
462  SelectParticles(Particles,fLcPlusP2Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
463  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
464  //LambdaC_bar -> p_bar pi+ pi-
465  SelectParticles(Particles,fLcMinusP2Pi,PrimVtx,fCutsCharm[2],fCutsCharm[1],
466  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
467 
468 
469  //D+ -> K0 pi+
470  for(int iPV=0; iPV<fNPV; iPV++)
471  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastPion(),
472  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]) );
473  //D- -> K0 pi-
474  for(int iPV=0; iPV<fNPV; iPV++)
475  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(),
476  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]) );
477  //Ds+ -> K0 K+
478  for(int iPV=0; iPV<fNPV; iPV++)
479  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[0], 1, vRTracks[0].FirstKaon(), vRTracks[0].LastKaon(),
480  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]) );
481  //Ds- -> K0 K-
482  for(int iPV=0; iPV<fNPV; iPV++)
483  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[1], -1, vRTracks[1].FirstKaon(), vRTracks[1].LastKaon(),
484  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]) );
485  //Lambdac+ -> Lambda pi+
486  for(int iPV=0; iPV<fNPV; iPV++)
487  FindTrackV0Decay(fPrimCandidates[1][iPV], 3122, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastPion(),
488  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]) );
489  //Lambdac_bar- -> Lambda_bar pi-
490  for(int iPV=0; iPV<fNPV; iPV++)
491  FindTrackV0Decay(fPrimCandidates[2][iPV], -3122, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(),
492  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]) );
493  //Lambdac+ -> p K0s
494  for(int iPV=0; iPV<fNPV; iPV++)
495  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(),
496  Particles, PrimVtx, -1, &(ChiToPrimVtx[0]) );
497  //Lambdac_bar- -> p_bar K0s
498  for(int iPV=0; iPV<fNPV; iPV++)
499  FindTrackV0Decay(fPrimCandidates[0][iPV], 310, vRTracks[1], -1, vRTracks[1].FirstProton(), vRTracks[1].LastProton(),
500  Particles, PrimVtx, -1, &(ChiToPrimVtx[1]) );
501 
502  //H0 -> Lambda Lambda
503  CombinePartPart(fSecCandidates[1], fSecCandidates[1], Particles, PrimVtx, fCutsPartPart[0], -1, 3000, 1, 1);
504  //H0 -> Lambda p pi-
505  FindTrackV0Decay(fLPi, 3002, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(),
506  Particles, PrimVtx, -1);
507  //Sigma0 -> Lambda Gamma
508  for(int iPV=0; iPV<fNPV; iPV++)
509  CombinePartPart(fPrimCandidates[3][iPV], fPrimCandidates[1][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 3212);
510  //Sigma0_bar -> Lambda_bar Gamma
511  for(int iPV=0; iPV<fNPV; iPV++)
512  CombinePartPart(fPrimCandidates[3][iPV], fPrimCandidates[2][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, -3212);
513  //pi0 -> gamma gamma
514  const float& mPi0 = KFParticleDatabase::Instance()->GetPi0Mass();
515  const float& mPi0Sigma = KFParticleDatabase::Instance()->GetPi0MassSigma();
516  CombinePartPart(fSecCandidates[3], fSecCandidates[3], Particles, PrimVtx, fCutsPartPart[1], -1, 111, 1, 0, &fPrimCandidates[4], &fSecCandidates[4], mPi0, mPi0Sigma);
517  for(int iPV=0; iPV<fNPV; iPV++)
518  {
519  CombinePartPart(fPrimCandidates[3][iPV], fPrimCandidates[3][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 111, 1, 0, &fPrimCandidates[4], &fSecCandidates[4], mPi0, mPi0Sigma);
520  CombinePartPart(fSecCandidates[3], fPrimCandidates[3][iPV], Particles, PrimVtx, fCutsPartPart[1], -1, 111, 0, 0, &fPrimCandidates[4], &fSecCandidates[4], mPi0, mPi0Sigma);
521  }
522  for(int iPV=0; iPV<fNPV; iPV++ )
523  ExtrapolateToPV(fPrimCandidates[4][iPV],PrimVtx[iPV]);
524  //eta -> pi0 pi0 pi0
525  //TODO implement this
526  //Sigma+ -> p pi0
527  FindTrackV0Decay(fSecCandidates[4], 111, vRTracks[0], 1, vRTracks[0].FirstProton(), vRTracks[0].LastProton(),
528  Particles, PrimVtx, -1);
529  //Sigma+_bar -> p- pi0
530  FindTrackV0Decay(fSecCandidates[4], 111, vRTracks[1], -1, vRTracks[1].FirstProton(), vRTracks[1].LastProton(),
531  Particles, PrimVtx, -1);
532  //Xi0 -> Lambda pi0
533  CombinePartPart(fSecCandidates[4], fSecCandidates[1], Particles, PrimVtx, fCutsPartPart[0], -1, 3322);
534  //Xi0_bar -> Lambda_bar pi0
535  CombinePartPart(fSecCandidates[4], fSecCandidates[2], Particles, PrimVtx, fCutsPartPart[0], -1, -3322);
536  //K*+ -> K+ pi0
537  for(int iPV=0; iPV<fNPV; iPV++)
538  FindTrackV0Decay(fPrimCandidates[4][iPV], 111, vRTracks[2], 1, vRTracks[2].FirstKaon(), vRTracks[2].LastKaon(),
539  Particles, PrimVtx, -1);
540  //K*- -> K- pi0
541  for(int iPV=0; iPV<fNPV; iPV++)
542  FindTrackV0Decay(fPrimCandidates[4][iPV], 111, vRTracks[3], -1, vRTracks[3].FirstKaon(), vRTracks[3].LastKaon(),
543  Particles, PrimVtx, -1);
544  //K*0 -> K0 pi0
545  for(int iPV=0; iPV<fNPV; iPV++)
546  CombinePartPart(fPrimCandidates[4][iPV], fPrimCandidates[0][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 100313, 0, 1);
547  //Sigma*0 -> Lambda pi0
548  for(int iPV=0; iPV<fNPV; iPV++)
549  CombinePartPart(fPrimCandidates[4][iPV], fPrimCandidates[1][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 3214, 0, 1);
550  //Sigma*0_bar -> Lambda_bar pi0
551  for(int iPV=0; iPV<fNPV; iPV++)
552  CombinePartPart(fPrimCandidates[4][iPV], fPrimCandidates[2][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, -3214, 0, 1);
553  //Xi*- -> Xi- pi0
554  for(int iPV=0; iPV<fNPV; iPV++)
555  CombinePartPart(fPrimCandidates[4][iPV], fPrimCandidates[5][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 3314, 0, 1);
556  //Xi*+ -> Xi+ pi0
557  for(int iPV=0; iPV<fNPV; iPV++)
558  CombinePartPart(fPrimCandidates[4][iPV], fPrimCandidates[6][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, -3314, 0, 1);
559  //JPsi -> Lambda Lambda_bar
560  for(int iPV=0; iPV<fNPV; iPV++)
561  CombinePartPart(fPrimCandidates[1][iPV], fPrimCandidates[2][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 300443, 0, 1);
562  //JPsi -> Xi- Xi+
563  for(int iPV=0; iPV<fNPV; iPV++)
564  CombinePartPart(fPrimCandidates[5][iPV], fPrimCandidates[6][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 400443, 0, 1);
565  //Psi -> Omega- Omega+
566  for(int iPV=0; iPV<fNPV; iPV++)
567  CombinePartPart(fPrimCandidates[7][iPV], fPrimCandidates[8][iPV], Particles, PrimVtx, fCutsPartPart[1], iPV, 500443, 0, 1);
568 
569  //reconstruct particles with daughters in ElectroMagnetic Calorimeter
570 // if(fEmcClusters)
571 // {
572 // //pi0 -> gamma gamma, EMC
573 // vector< vector<KFParticle> > vPi0PrimEmc(1);
574 // vector<KFParticle> vPi0SecEmc;
575 // vector< vector<KFParticle> > vD0PrimEmc(1);
576 // CombinePartPart(vGammaPrimEmc, vGammaPrimEmc, Particles, PrimVtx, fCutsPartPart[1], 0, 111, 1, 0, &vPi0PrimEmc, &vPi0SecEmc, mPi0, mPi0Sigma);
577 //
578 // //D+ -> K0 pi+
579 // FindTrackV0Decay(fSecCandidates[0], 310, vRTracks[0], 1, vRTracks[0].FirstPion(), vRTracks[0].LastPion(), Particles, PrimVtx, -1/*, &(ChiToPrimVtx[0])*/);
580 // //D0 -> K0 pi+ pi-
581 // FindTrackV0Decay(fK0PiPlus, 100411, vRTracks[1], -1, vRTracks[1].FirstPion(), vRTracks[1].LastPion(), Particles, PrimVtx, -1/*, &(ChiToPrimVtx[0])*/);
582 // //D0 -> K0 pi+ pi- pi0
583 // CombinePartPart(fK0PiPi, vPi0PrimEmc[0], Particles, PrimVtx, fCutsPartPart[1], -1, 428, 0, 0, &vD0PrimEmc, 0,
584 // KFParticleDatabase::Instance()->GetD0Mass(), 0.025);
585 //
586 // for(int iPV=0; iPV<1; iPV++ )
587 // {
588 // ExtrapolateToPV(vPi0PrimEmc[iPV],PrimVtx[iPV]);
589 // ExtrapolateToPV(vD0PrimEmc[iPV],PrimVtx[iPV]);
590 // }
591 // //D0* -> D0 pi0
592 // CombinePartPart(vD0PrimEmc[0], vPi0PrimEmc[0], Particles, PrimVtx, fCutsPartPart[1], 0, 10428);
593 // }
594  }
595  else
596  {
597  //D0 -> pi+ K-
598  SelectParticles(Particles,fD0,PrimVtx,fCutsCharm[2],fCutsCharm[1],
599  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
600  //D0_bar -> pi+ K-
601  SelectParticles(Particles,fD0bar,PrimVtx,fCutsCharm[2],fCutsCharm[1],
602  KFParticleDatabase::Instance()->GetD0Mass(), KFParticleDatabase::Instance()->GetD0MassSigma(), fSecCuts[0]);
603  }
604 }
605 
606 void KFParticleFinder::ExtrapolateToPV(vector<KFParticle>& vParticles, KFParticleSIMD& PrimVtx)
607 {
612  KFParticle* parts[float_vLen];
613  KFParticle tmpPart[float_vLen];
614 
615  for(int iv=0; iv<float_vLen; iv++)
616  parts[iv] = &tmpPart[iv];
617 
618  for(unsigned int iL=0; iL<vParticles.size(); iL += float_vLen)
619  {
620 
621  unsigned int nPart = vParticles.size();
622  unsigned int nEntries = (iL + float_vLen < nPart) ? float_vLen : (nPart - iL);
623 
624 
625  for(unsigned int iv=0; iv<nEntries; iv++)
626  tmpPart[iv] = vParticles[iL+iv];
627 
628  KFParticleSIMD tmp(parts,nEntries);
629 
630  tmp.TransportToPoint(PrimVtx.Parameters());
631 
632  for(unsigned int iv=0; iv<nEntries; iv++)
633  {
634  tmp.GetKFParticle(vParticles[iL+iv], iv);
635  }
636  }
637 }
638 
639 inline void KFParticleFinder::ConstructV0(KFPTrackVector* vTracks,
640  int iTrTypePos,
641  int iTrTypeNeg,
642  uint_v& idPosDaughters,
643  uint_v& idNegDaughters,
644  int_v& daughterPosPDG,
645  int_v& daughterNegPDG,
646  KFParticleSIMD& mother,
647  KFParticle& mother_temp,
648  const unsigned short NTracks,
649  kfvector_floatv& l,
650  kfvector_floatv& dl,
651  vector<KFParticle>& Particles,
652  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
653  const float* cuts,
654  const int_v& pvIndex,
655  const float* secCuts,
656  const float_v& massMotherPDG,
657  const float_v& massMotherPDGSigma,
658  KFParticleSIMD& motherPrimSecCand,
659  int& nPrimSecCand,
660  vector< vector<KFParticle> >* vMotherPrim,
661  vector<KFParticle>* vMotherSec
662  )
663 {
698  float_m isPrimary(simd_cast<float_m>(pvIndex>-1));
699  int_v trackId;
700  KFParticleSIMD posDaughter(vTracks[iTrTypePos],idPosDaughters, daughterPosPDG);
701  trackId.gather( &(vTracks[iTrTypePos].Id()[0]), idPosDaughters );
702  posDaughter.SetId(trackId);
703 
704  KFParticleSIMD negDaughter(vTracks[iTrTypeNeg],idNegDaughters, daughterNegPDG);
705  trackId.gather( &(vTracks[iTrTypeNeg].Id()[0]), idNegDaughters );
706  negDaughter.SetId(trackId);
707 #ifdef CBM
708  float_v ds[2] = {0.f,0.f};
709  float_v dsdr[4][6];
710  negDaughter.GetDStoParticle( posDaughter, ds, dsdr );
711  negDaughter.TransportToDS(ds[0], dsdr[0]);
712  posDaughter.TransportToDS(ds[1], dsdr[3]);
713 #endif
714  const KFParticleSIMD* vDaughtersPointer[2] = {&negDaughter, &posDaughter};
715  mother.Construct(vDaughtersPointer, 2, 0);
716 
717  float_m saveParticle(simd_cast<float_m>(int_v::IndexesFromZero() < int(NTracks)));
718  float_v chi2Cut = cuts[1];
719  float_v ldlCut = cuts[2];
720  if( !(simd_cast<float_m>(abs(mother.PDG()) == 421 || abs(mother.PDG()) == 426 || abs(mother.PDG()) == 420)).isEmpty() )
721  {
722  chi2Cut( simd_cast<float_m>(abs(mother.PDG()) == 421 || abs(mother.PDG()) == 426 || abs(mother.PDG()) == 420) ) = fCutsCharm[0];
723  ldlCut( simd_cast<float_m>(abs(mother.PDG()) == 421 || abs(mother.PDG()) == 426 || abs(mother.PDG()) == 420) ) = -1;//fCutsCharm[1];
724  }
725 
726  saveParticle &= (mother.Chi2()/simd_cast<float_v>(mother.NDF()) < chi2Cut );
727  saveParticle &= KFPMath::Finite(mother.GetChi2());
728  saveParticle &= (mother.GetChi2() > 0.0f);
729  saveParticle &= (mother.GetChi2() == mother.GetChi2());
730 
731  if( saveParticle.isEmpty() ) return;
732 
733  float_v lMin(1.e8f);
734  float_v ldlMin(1.e8f);
735  float_m isParticleFromVertex(false);
736 
737  for(int iP=0; iP<fNPV; iP++)
738  {
739  float_m isParticleFromVertexLocal;
740  mother.GetDistanceToVertexLine(PrimVtx[iP], l[iP], dl[iP], &isParticleFromVertexLocal);
741  isParticleFromVertex |= isParticleFromVertexLocal;
742  float_v ldl = (l[iP]/dl[iP]);
743  lMin( (l[iP] < lMin) && saveParticle) = l[iP];
744  ldlMin( (ldl < ldlMin) && saveParticle) = ldl;
745  }
746 
747  saveParticle &= (lMin < 200.f);
748 #ifdef NonhomogeneousField
749  KFParticleSIMD motherTopo;
750  ldlMin = 1.e8f;
751  for(int iP=0; iP<fNPV; iP++)
752  {
753  motherTopo = mother;
754  motherTopo.SetProductionVertex(PrimVtx[iP]);
755  motherTopo.GetDecayLength(l[iP], dl[iP]);
756  float_v ldl = (l[iP]/dl[iP]);
757  ldlMin( (ldl < ldlMin) && saveParticle) = ldl;
758  }
759 #endif
760  saveParticle &= ( (float_m(!isPrimary) && ldlMin > ldlCut) || float_m(isPrimary) );
761 
762  saveParticle &= (float_m(!isPrimary) && isParticleFromVertex) || isPrimary;
763  if( saveParticle.isEmpty() ) return;
764 
765  float_m isK0 = saveParticle && simd_cast<float_m>(mother.PDG() == int_v(310));
766  float_m isLambda = saveParticle && simd_cast<float_m>(abs(mother.PDG()) == int_v(3122));
767  float_m isGamma = saveParticle && simd_cast<float_m>(mother.PDG() == int_v(22));
768  float_m isHyperNuclei = saveParticle && simd_cast<float_m>(abs(mother.PDG()) > 3000 && abs(mother.PDG()) < 3104);
769 
770  saveParticle &= ( ((isK0 || isLambda || isHyperNuclei) && lMin > float_v(fLCut)) || !(isK0 || isLambda || isHyperNuclei) );
771 
772  float_m saveMother(false);
773 
774  if( !(isK0.isEmpty()) || !(isLambda.isEmpty()) || !(isGamma.isEmpty()))
775  {
776  float_v mass, errMass;
777 
778  mother.GetMass(mass, errMass);
779  saveMother = saveParticle;
780  saveMother &= (abs(mass - massMotherPDG)/massMotherPDGSigma) < secCuts[0];
781  saveMother &= ((ldlMin > secCuts[2]) && !isGamma) || isGamma;
782  saveMother &= (isK0 || isLambda || isGamma);
783  }
784 
785  for(int iv=0; iv<NTracks; iv++)
786  {
787  if(!saveParticle[iv]) continue;
788 
789  mother.GetKFParticle(mother_temp, iv);
790  int motherId = Particles.size();
791  mother_temp.SetId(Particles.size());
792  if( mother.PDG()[iv] == 421 )
793  {
794  fD0.push_back(mother_temp);
795  continue;
796  }
797  if( mother.PDG()[iv] == -421 )
798  {
799  fD0bar.push_back(mother_temp);
800  continue;
801  }
802  if( mother.PDG()[iv] == 426 )
803  {
804  fD0KK.push_back(mother_temp);
805  continue;
806  }
807  if( mother.PDG()[iv] == 420 )
808  {
809  fD0pipi.push_back(mother_temp);
810  continue;
811  }
812  if( mother.PDG()[iv] == 3004)
813  fHe3Pi.push_back(mother_temp);
814  if( mother.PDG()[iv] == -3004)
815  fHe3PiBar.push_back(mother_temp);
816  if( mother.PDG()[iv] == 3005)
817  fHe4Pi.push_back(mother_temp);
818  if( mother.PDG()[iv] == -3005)
819  fHe4PiBar.push_back(mother_temp);
820 
821  Particles.push_back(mother_temp);
822 
823  if( mother.PDG()[iv] == 22 && isPrimary[iv] )
824  {
825  float negPt2 = negDaughter.Px()[iv]*negDaughter.Px()[iv] + negDaughter.Py()[iv]*negDaughter.Py()[iv];
826  float posPt2 = posDaughter.Px()[iv]*posDaughter.Px()[iv] + posDaughter.Py()[iv]*posDaughter.Py()[iv];
827 
828  if( (negPt2 >fCutLVMPt*fCutLVMPt) && (posPt2 >fCutLVMPt*fCutLVMPt) )
829  {
830  mother_temp.SetPDG(100113);
831  mother_temp.SetId(Particles.size());
832  Particles.push_back(mother_temp);
833 
834  if( (negPt2 >fCutJPsiPt*fCutJPsiPt) && (posPt2 >fCutJPsiPt*fCutJPsiPt) )
835  {
836  mother_temp.SetPDG(443);
837  mother_temp.SetId(Particles.size());
838  Particles.push_back(mother_temp);
839  }
840  }
841  }
842 
843  if( mother.PDG()[iv] == 200113 )
844  {
845  float negPt2 = negDaughter.Px()[iv]*negDaughter.Px()[iv] + negDaughter.Py()[iv]*negDaughter.Py()[iv];
846  float posPt2 = posDaughter.Px()[iv]*posDaughter.Px()[iv] + posDaughter.Py()[iv]*posDaughter.Py()[iv];
847 
848  if( (negPt2 >fCutJPsiPt*fCutJPsiPt) && (posPt2 >fCutJPsiPt*fCutJPsiPt) && (abs(daughterPosPDG[iv]) == 13) && (abs(daughterNegPDG[iv]) == 13))
849  {
850  mother_temp.SetPDG(100443);
851  mother_temp.SetId(Particles.size());
852  Particles.push_back(mother_temp);
853  }
854  }
855 
856  if(saveMother[iv])
857  {
858  mother.SetId(motherId);
859  motherPrimSecCand.SetOneEntry(nPrimSecCand,mother,iv);
860 
861  nPrimSecCand++;
862  if(nPrimSecCand==float_vLen)
863  {
864  SaveV0PrimSecCand(motherPrimSecCand,nPrimSecCand,mother_temp,PrimVtx,secCuts,vMotherPrim,vMotherSec);
865  nPrimSecCand = 0;
866  }
867  }
868  }
869 }
870 
871 inline void KFParticleFinder::SaveV0PrimSecCand(KFParticleSIMD& mother,
872  int& NParticles,
873  KFParticle& mother_temp,
874  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
875  const float* secCuts,
876  vector< vector<KFParticle> >* vMotherPrim,
877  vector<KFParticle>* vMotherSec)
878 {
890  KFParticleSIMD motherTopo;
891  float_v massMotherPDG, massMotherPDGSigma;
892 
893  float_m isSec(false);
894  float_m isPrim(false);
895  vector<int> iPrimVert[float_vLen];
896 
897  KFParticleDatabase::Instance()->GetMotherMass(mother.PDG(),massMotherPDG,massMotherPDGSigma);
898 
899  float_m isK0 = simd_cast<float_m>(mother.PDG() == int_v(310));
900  float_m isLambda = simd_cast<float_m>(abs(mother.PDG()) == int_v(3122));
901  float_m isGamma = simd_cast<float_m>(mother.PDG() == int_v(22));
902 
903  int_v arrayIndex(-1); //for saving primary candidates;
904 
905  arrayIndex(mother.PDG() == int_v(310)) = 0;
906  arrayIndex(mother.PDG() == int_v(3122)) = 1;
907  arrayIndex(mother.PDG() == int_v(-3122)) = 2;
908  arrayIndex(mother.PDG() == int_v(22)) = 3;
909 
910  float_m isPrimaryPart(false);
911 
912  float_v chi2TopoMin = 1.e4f;
913 
914  for(int iP=0; iP< fNPV; iP++)
915  {
916  motherTopo = mother;
917  motherTopo.SetProductionVertex(PrimVtx[iP]);
918 
919  const float_v& motherTopoChi2Ndf = motherTopo.GetChi2()/simd_cast<float_v>(motherTopo.GetNDF());
920  chi2TopoMin(motherTopoChi2Ndf < chi2TopoMin) = motherTopoChi2Ndf;
921  const float_m isPrimaryPartLocal = ( motherTopoChi2Ndf < secCuts[1] );
922  if(isPrimaryPartLocal.isEmpty()) continue;
923  isPrimaryPart |= isPrimaryPartLocal;
924  for(int iV=0; iV<NParticles; iV++)
925  {
926  if(isPrimaryPartLocal[iV])
927  {
928  motherTopo.GetKFParticle(mother_temp, iV);
929  fPrimCandidatesTopo[arrayIndex[iV]][iP].push_back(mother_temp);
930  iPrimVert[iV].push_back(iP);
931  }
932  }
933 
934  motherTopo.SetNonlinearMassConstraint(massMotherPDG);
935  for(int iV=0; iV<NParticles; iV++)
936  {
937  if(isPrimaryPartLocal[iV])
938  {
939  motherTopo.GetKFParticle(mother_temp, iV);
940  fPrimCandidatesTopoMass[arrayIndex[iV]][iP].push_back(mother_temp);
941  }
942  }
943  }
944 
945  isPrim |= ( ( isPrimaryPart ) && (isK0 || isLambda || isGamma) );
946 #ifdef __ROOT__
947  isSec |= ( (!isPrimaryPart ) && (isK0 || isLambda || isGamma) && (chi2TopoMin < float_v(500.f)) );
948 #else
949  isSec |= ( (!isPrimaryPart ) && (isK0 || isLambda || isGamma) );
950 #endif
951 
952  mother.SetNonlinearMassConstraint(massMotherPDG);
953 
954  for(int iv=0; iv<NParticles; iv++)
955  {
956  if(isPrim[iv] || isSec[iv])
957  {
958  mother.GetKFParticle(mother_temp, iv);
959 
960  if(isPrim[iv] )
961  {
962  for(unsigned int iP = 0; iP<iPrimVert[iv].size(); iP++)
963  vMotherPrim[arrayIndex[iv]][iPrimVert[iv][iP]].push_back(mother_temp);
964  }
965 
966  if(isSec[iv] )
967  vMotherSec[arrayIndex[iv]].push_back(mother_temp);
968  }
969  }
970 }
971 
972 void KFParticleFinder::Find2DaughterDecay(KFPTrackVector* vTracks, kfvector_float* ChiToPrimVtx,
973  vector<KFParticle>& Particles,
974  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
975  const float* cuts,
976  const float* secCuts,
977  vector< vector<KFParticle> >* vMotherPrim,
978  vector<KFParticle>* vMotherSec )
979 {
998  KFParticle mother_temp;
999  KFParticleSIMD mother;
1000  kfvector_floatv l(fNPV), dl(fNPV);
1001 
1002  KFParticleSIMD daughterNeg, daughterPos;
1003 
1004  // for secondary V0
1005  unsigned int nBufEntry = 0;
1006  uint_v idNegDaughters;
1007  uint_v idPosDaughters;
1008  int_v daughterPosPDG(-1);
1009  int_v daughterNegPDG(-1);
1010 
1011  int_v pvIndexMother(-1);
1012 
1013  float_v massMotherPDG(Vc::Zero), massMotherPDGSigma(Vc::Zero);
1014  int_v V0PDG(Vc::Zero);
1015 
1016  KFParticleSIMD motherPrimSecCand;
1017  int nPrimSecCand =0;
1018 
1019  int trTypeIndexPos[2] = {0,2};
1020  int trTypeIndexNeg[2] = {1,3};
1021 
1022  for( int iTrTypeNeg = 0; iTrTypeNeg<2; iTrTypeNeg++)
1023  {
1024  KFPTrackVector& negTracks = vTracks[ trTypeIndexNeg[iTrTypeNeg] ];
1025 
1026  for(int iTrTypePos=0; iTrTypePos<2; iTrTypePos++)
1027  {
1028  KFPTrackVector& posTracks = vTracks[ trTypeIndexPos[iTrTypePos] ];
1029  int_v negTracksSize = negTracks.Size();
1030  int nPositiveTracks = posTracks.Size();
1031 
1032  //track categories
1033  int nTC = 5;
1034  int startTCPos[5] = {0};
1035  int endTCPos[5] = {0};
1036  int startTCNeg[5] = {0};
1037  int endTCNeg[5] = {0};
1038 
1039  if((iTrTypeNeg == 0) && (iTrTypePos == 0))
1040  {
1041  // Secondary particles
1042  nTC = 5;
1043  // e-
1044  startTCPos[0] = 0; endTCPos[0] = nPositiveTracks; //posTracks.LastElectron();
1045  startTCNeg[0] = 0; endTCNeg[0] = negTracksSize[0]; //negTracks.LastElectron();
1046  //mu-
1047  startTCPos[1] = 0; endTCPos[1] = 0;
1048  startTCNeg[1] = 0; endTCNeg[1] = 0;
1049  //pi- + ghosts
1050  startTCPos[2] = posTracks.FirstPion(); endTCPos[2] = nPositiveTracks;
1051  startTCNeg[2] = negTracks.FirstPion(); endTCNeg[2] = negTracks.LastPion();
1052  //K-
1053  startTCPos[3] = posTracks.FirstPion(); endTCPos[3] = posTracks.LastKaon();
1054  startTCNeg[3] = negTracks.FirstKaon(); endTCNeg[3] = negTracks.LastKaon();
1055  //p-, d-, t-, he3-, he4-
1056  startTCPos[4] = posTracks.FirstPion(); endTCPos[4] = posTracks.LastPion();
1057  startTCNeg[4] = negTracks.FirstProton(); endTCNeg[4] = negTracksSize[0];
1058  }
1059 
1060  if( iTrTypeNeg != iTrTypePos )
1061  {
1062  //Mixed particles - only gamma -> e+ e-
1063  nTC = 1;
1064  startTCPos[0] = 0; endTCPos[0] = nPositiveTracks; //posTracks.LastElectron();
1065  startTCNeg[0] = 0; endTCNeg[0] = negTracksSize[0]; //negTracks.LastElectron();
1066  }
1067 
1068  if((iTrTypeNeg == 1) && (iTrTypePos == 1))
1069  {
1070  //primary particles
1071  nTC = 5;
1072  // e-
1073  startTCPos[0] = 0; endTCPos[0] = nPositiveTracks; //posTracks.LastElectron();
1074  startTCNeg[0] = 0; endTCNeg[0] = negTracksSize[0]; //negTracks.LastElectron();
1075  //mu-
1076  startTCPos[1] = posTracks.FirstMuon(); endTCPos[1] = posTracks.LastMuon();
1077  startTCNeg[1] = negTracks.FirstMuon(); endTCNeg[1] = negTracks.LastMuon();
1078  //pi- + ghosts
1079  startTCPos[2] = posTracks.FirstPion(); endTCPos[2] = posTracks.LastProton();
1080  startTCNeg[2] = negTracks.FirstPion(); endTCNeg[2] = negTracks.LastPion();
1081  //K-
1082  startTCPos[3] = posTracks.FirstPion(); endTCPos[3] = nPositiveTracks;
1083  startTCNeg[3] = negTracks.FirstKaon(); endTCNeg[3] = negTracks.LastKaon();
1084  //p-
1085  startTCPos[4] = posTracks.FirstPion(); endTCPos[4] = posTracks.LastProton();
1086  startTCNeg[4] = negTracks.FirstProton(); endTCNeg[4] = negTracks.LastProton();
1087  }
1088 
1089  for(int iTC=0; iTC<nTC; iTC++)
1090  {
1091  for(int iTrN=startTCNeg[iTC]; iTrN < endTCNeg[iTC]; iTrN += float_vLen)
1092  {
1093  const int NTracksNeg = (iTrN + float_vLen < negTracks.Size()) ? float_vLen : (negTracks.Size() - iTrN);
1094 
1095  int_v negInd = int_v::IndexesFromZero() + int(iTrN);
1096 
1097  int_v negPDG = reinterpret_cast<const int_v&>(negTracks.PDG()[iTrN]);
1098  int_v negPVIndex = reinterpret_cast<const int_v&>(negTracks.PVIndex()[iTrN]);
1099  int_v negNPixelHits = reinterpret_cast<const int_v&>(negTracks.NPixelHits()[iTrN]);
1100 
1101  int_v trackPdgNeg = negPDG;
1102  int_m activeNeg = (negPDG != -1);
1103 #ifdef CBM
1104  if( !((negPDG == -1).isEmpty()) )
1105  {
1106  trackPdgNeg(negPVIndex<0 && (negPDG == -1) ) = -211;
1107 
1108  activeNeg |= int_m(negPVIndex < 0) && int_m(negPDG == -1) ;
1109  }
1110 #endif
1111  activeNeg &= (int_v::IndexesFromZero() < int(NTracksNeg));
1112 
1113  daughterNeg.Load(negTracks, iTrN, negPDG);
1114 
1115  float_v chiPrimNeg(Vc::Zero);
1116  float_v chiPrimPos(Vc::Zero);
1117 
1118  if( (iTrTypeNeg == 0) && (iTrTypePos == 0) )
1119  chiPrimNeg = reinterpret_cast<const float_v&>( ChiToPrimVtx[trTypeIndexNeg[iTrTypeNeg]][iTrN]);
1120 
1121  for(int iTrP=startTCPos[iTC]; iTrP < endTCPos[iTC]; iTrP += float_vLen)
1122  {
1123  const int NTracks = (iTrP + float_vLen < nPositiveTracks) ? float_vLen : (nPositiveTracks - iTrP);
1124 
1125  const int_v& posPDG = reinterpret_cast<const int_v&>(posTracks.PDG()[iTrP]);
1126  const int_v& posPVIndex = reinterpret_cast<const int_v&>(posTracks.PVIndex()[iTrP]);
1127  const int_v& posNPixelHits = reinterpret_cast<const int_v&>(posTracks.NPixelHits()[iTrP]);
1128  const int_m& isPosSecondary = (posPVIndex < 0);
1129 
1130  daughterPos.Load(posTracks, iTrP, posPDG);
1131 
1132  if( (iTrTypeNeg == 0) && (iTrTypePos == 0) )
1133  chiPrimPos = reinterpret_cast<const float_v&>( ChiToPrimVtx[trTypeIndexPos[iTrTypePos]][iTrP]);
1134 
1135  for(int iRot = 0; iRot<float_vLen; iRot++)
1136  {
1137 // if(iRot>0)
1138  {
1139  negPDG = negPDG.rotated(1);
1140  negPVIndex = negPVIndex.rotated(1);
1141  negNPixelHits = negNPixelHits.rotated(1);
1142  negInd = negInd.rotated(1);
1143  trackPdgNeg = trackPdgNeg.rotated(1);
1144 
1145  daughterNeg.Rotate();
1146  chiPrimNeg = chiPrimNeg.rotated(1);
1147 
1148  activeNeg = ( (negPDG != -1) || ( (negPVIndex < 0) && (negPDG == -1) ) ) && (negInd < negTracksSize);
1149  }
1150  const int_m& isSecondary = int_m( negPVIndex < 0 ) && isPosSecondary;
1151  const int_m& isPrimary = int_m( negPVIndex >= 0 ) && (!isPosSecondary);
1152 
1153  float_m closeDaughters = simd_cast<float_m>(activeNeg && (int_v::IndexesFromZero() < int_v(NTracks)));
1154 
1155  if(closeDaughters.isEmpty() && (iTC != 0)) continue;
1156 
1157 
1158  int_v trackPdgPos[2];
1159  int_m active[2];
1160 
1161  active[0] = (posPDG != -1);
1162  active[0] &= ((isPrimary && (posPVIndex == negPVIndex)) || !(isPrimary));
1163 
1164  active[1] = int_m(false);
1165 
1166  trackPdgPos[0] = posPDG;
1167 #ifdef CBM
1168  int nPDGPos = 2;
1169  if( (posPDG == -1).isEmpty() && (posPDG > 1000000000).isEmpty() && (posPDG == 211).isEmpty() )
1170  {
1171  nPDGPos = 1;
1172  }
1173  else
1174  {
1175  trackPdgPos[0](isSecondary && posPDG == -1) = 211;
1176  trackPdgPos[1] = 2212;
1177 
1178  active[0] |= isSecondary && int_m(posPDG == -1);
1179  active[1] = isSecondary && (int_m(posPDG == -1) || (posPDG > 1000000000) || (posPDG == 211));
1180  }
1181 #else
1182  int nPDGPos = 1;
1183 #endif
1184  active[0] &= simd_cast<int_m>(closeDaughters);
1185  active[1] &= simd_cast<int_m>(closeDaughters);
1186 
1187  if(iTC==0)
1188  {
1189  nPDGPos = 1;
1190  active[0] = (negInd < negTracksSize) && (int_v::IndexesFromZero() < int_v(NTracks));
1191  }
1192 
1193  for(int iPDGPos=0; iPDGPos<nPDGPos; iPDGPos++)
1194  {
1195  if(active[iPDGPos].isEmpty()) continue;
1196 
1197  //detetrmine a pdg code of the mother particle
1198 
1199  int_v motherPDG(-1);
1200 
1201  if(!fMixedEventAnalysis)
1202  {
1203  if(iTC==0)
1204  {
1205  motherPDG( (abs(trackPdgPos[iPDGPos]) == 11) || int_m(abs(trackPdgNeg) == 11) || isSecondary ) = 22; //gamma -> e+ e-
1206  }
1207  else if(iTC==1)
1208  {
1209  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 13 || abs(trackPdgPos[iPDGPos])==19)
1210  && (int_m(abs(trackPdgNeg) == 13) || int_m(abs(trackPdgNeg) == 19)) ) = 200113; //rho -> mu+ mu-
1211  }
1212  else if(iTC==2)
1213  {
1214  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 211) ) = 310; //K0 -> pi+ pi-
1215  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 2212) && int_m(abs(trackPdgNeg) == 211) ) = 3122; //Lambda -> p+ pi-
1216  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 1000010020) && int_m(abs(trackPdgNeg) == 211) ) = 3003; //LambdaN -> d+ pi-
1217  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 1000010030) && int_m(abs(trackPdgNeg) == 211) ) = 3103; //LambdaNN -> t+ pi-
1218  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 1000020030) && int_m(abs(trackPdgNeg) == 211) ) = 3004; //H3Lambda -> He3+ pi-
1219  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 1000020040) && int_m(abs(trackPdgNeg) == 211) ) = 3005; //H4Lambda -> He4+ pi-
1220  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 211) ) = -421; //D0_bar -> pi- K+
1221  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 211) ) = 313; //K*0 -> K+ pi-
1222  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 211) ) = 113; //rho -> pi+ pi-
1223  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 2212) && int_m(abs(trackPdgNeg) == 211) ) = 2114; //Delta0 -> p pi-
1224  }
1225  else if(iTC==3)
1226  {
1227  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 321) ) = 421; //D0 -> pi+ K-
1228  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 321) ) = 426; //D0 -> K+ K-
1229  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 321) ) = -313; //K*0_bar -> K- pi+
1230  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 2212) && int_m(abs(trackPdgNeg) == 321) ) = 3124; //Lambda* -> p K-
1231  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 321) ) = 333; //phi -> K+ K-
1232  }
1233  else if(iTC==4)
1234  {
1235  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 2212) ) = -3122; //Lambda_bar -> p- pi+
1236  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 1000010020) ) = -3003; //LambdaN_bar -> d- pi+
1237  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 1000010030) ) = -3103; //LambdaNN_bar -> t- pi+
1238  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 1000020030) ) = -3004; //H3Lambda_bar -> He3- pi+
1239  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 1000020040) ) = -3005; //H4Lambda_bar -> He4- pi+
1240  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 2212) ) = -3124; //Lambda*_bar -> p- K+
1241  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 2212) && int_m(abs(trackPdgNeg) == 2212) ) = 200443; //JPsi -> p- p
1242  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 2212) ) = -2114; //Delta0_bar -> p- pi+
1243  }
1244  }
1245  else
1246  {
1247  if(iTC==0)
1248  motherPDG( (abs(trackPdgPos[iPDGPos])== 11) && int_m(abs(trackPdgNeg) == 11) ) = 22; //gamma -> e+ e-
1249  else if(iTC==1)
1250  motherPDG( isPrimary && (abs(trackPdgPos[iPDGPos])== 13 || abs(trackPdgPos[iPDGPos])==19)
1251  && (int_m(abs(trackPdgNeg) == 13) || int_m(abs(trackPdgNeg) == 19)) ) = 200113; //rho -> mu+ mu-
1252  else if(iTC==2)
1253  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 321) && int_m(abs(trackPdgNeg) == 211) ) = -421; //D0_bar -> pi- K+
1254  else if(iTC==3)
1255  motherPDG( isSecondary && (abs(trackPdgPos[iPDGPos])== 211) && int_m(abs(trackPdgNeg) == 321) ) = 421; //D0 -> pi+ K-
1256  }
1257 
1258  if( (iTrTypeNeg == 0) && (iTrTypePos == 0) )
1259  {
1260  float_v chiprimCut = fCuts2D[0];
1261  chiprimCut( simd_cast<float_m>(abs(motherPDG) == 421 || abs(motherPDG) == 426) ) = fCutCharmChiPrim;
1262  active[iPDGPos] &= simd_cast<int_m>(chiPrimNeg > chiprimCut && chiPrimPos > chiprimCut);
1263  }
1264 
1265  active[iPDGPos] &= (motherPDG != -1);
1266  if(!(fDecayReconstructionList.empty()))
1267  {
1268  for(int iV=0; iV<float_vLen; iV++)
1269  {
1270  if(!(active[iPDGPos][iV])) continue;
1271  if(fDecayReconstructionList.find(motherPDG[iV]) == fDecayReconstructionList.end())
1272  motherPDG[iV] = -1;
1273  }
1274  active[iPDGPos] &= (motherPDG != -1);
1275  }
1276  if(active[iPDGPos].isEmpty()) continue;
1277 
1278  if(!( (iTrTypePos == 1) && (iTrTypeNeg == 1) ) )
1279  {
1280  float_v dS[2];
1281  daughterNeg.GetDStoParticleFast( daughterPos, dS );
1282  float_v negParameters[8], posParameters[8];
1283  daughterNeg.TransportFast( dS[0], negParameters );
1284  daughterPos.TransportFast( dS[1], posParameters );
1285  float_v dx = negParameters[0]-posParameters[0];
1286  float_v dy = negParameters[1]-posParameters[1];
1287  float_v dz = negParameters[2]-posParameters[2];
1288  float_v dr = sqrt(dx*dx+dy*dy+dz*dz);
1289 
1290  active[iPDGPos] &= simd_cast<int_m>(dr < float_v(fDistanceCut));
1291  if(active[iPDGPos].isEmpty()) continue;
1292 
1293  float_v p1p2 = posParameters[3]*negParameters[3] + posParameters[4]*negParameters[4] + posParameters[5]*negParameters[5];
1294  float_v p12 = posParameters[3]*posParameters[3] + posParameters[4]*posParameters[4] + posParameters[5]*posParameters[5];
1295  float_v p22 = negParameters[3]*negParameters[3] + negParameters[4]*negParameters[4] + negParameters[5]*negParameters[5];
1296  active[iPDGPos] &= simd_cast<int_m>(p1p2 > -p12);
1297  active[iPDGPos] &= simd_cast<int_m>(p1p2 > -p22);
1298  }
1299 
1300  const float_v& ptNeg2 = daughterNeg.Px()*daughterNeg.Px() + daughterNeg.Py()*daughterNeg.Py();
1301  const float_v& ptPos2 = daughterPos.Px()*daughterPos.Px() + daughterPos.Py()*daughterPos.Py();
1302  if( !((abs(motherPDG) == 421 || abs(motherPDG) == 426).isEmpty()) )
1303  {
1304  active[iPDGPos] &= ( (abs(motherPDG) == 421 || abs(motherPDG) == 426) &&
1305  simd_cast<int_m>(ptNeg2 >= fCutCharmPt*fCutCharmPt) &&
1306  simd_cast<int_m>(ptPos2 >= fCutCharmPt*fCutCharmPt) &&
1307  simd_cast<int_m>(chiPrimNeg > fCutCharmChiPrim) && simd_cast<int_m>(chiPrimPos > fCutCharmChiPrim) &&
1308  int_m(negNPixelHits >= int_v(3)) && int_m(posNPixelHits >= int_v(3)) )
1309  || (!(abs(motherPDG) == 421 || abs(motherPDG) == 426));
1310  }
1311 
1312  if(active[iPDGPos].isEmpty()) continue;
1313 
1314  for(int iV=0; iV<float_vLen; iV++)
1315  {
1316  if(!(active[iPDGPos][iV])) continue;
1317 
1318 
1319  idPosDaughters[nBufEntry] = iTrP+iV;
1320  idNegDaughters[nBufEntry] = negInd[iV];
1321 
1322  daughterPosPDG[nBufEntry] = trackPdgPos[iPDGPos][iV];
1323  daughterNegPDG[nBufEntry] = trackPdgNeg[iV];
1324 
1325  if(motherPDG[iV] == 22)
1326  {
1327  daughterPosPDG[nBufEntry] = -11;
1328  daughterNegPDG[nBufEntry] = 11;
1329  }
1330 
1331  pvIndexMother[nBufEntry] = isPrimary[iV] ? negPVIndex[iV] : -1;
1332 
1333  if( iTrTypeNeg != iTrTypePos ) pvIndexMother[nBufEntry] = 0;
1334 
1335  V0PDG[nBufEntry] = motherPDG[iV];
1336 
1337  nBufEntry++;
1338 
1339  if(int(nBufEntry) == float_vLen)
1340  {
1341  KFParticleDatabase::Instance()->GetMotherMass(V0PDG,massMotherPDG,massMotherPDGSigma);
1342  mother.SetPDG( V0PDG );
1343  ConstructV0(vTracks, trTypeIndexPos[iTrTypePos], trTypeIndexNeg[iTrTypeNeg],
1344  idPosDaughters, idNegDaughters, daughterPosPDG, daughterNegPDG,
1345  mother, mother_temp,
1346  nBufEntry, l, dl, Particles, PrimVtx,
1347  cuts, pvIndexMother, secCuts, massMotherPDG,
1348  massMotherPDGSigma, motherPrimSecCand, nPrimSecCand, vMotherPrim, vMotherSec);
1349  nBufEntry = 0;
1350  }
1351 
1352  //TODO optimize this part of code for D-mesons
1353  if(motherPDG[iV] == 310 &&
1354  (fDecayReconstructionList.empty() ||
1355  (!(fDecayReconstructionList.empty()) && !(fDecayReconstructionList.find(420) == fDecayReconstructionList.end()) ) ) &&
1356  negNPixelHits[iV] >= 3 && posNPixelHits[iV] >= 3 &&
1357  chiPrimNeg[iV] > fCutCharmChiPrim && chiPrimPos[iV] > fCutCharmChiPrim &&
1358  ptNeg2[iV] >= fCutCharmPt*fCutCharmPt && ptPos2[iV] >= fCutCharmPt*fCutCharmPt )
1359  {
1360  idPosDaughters[nBufEntry] = iTrP+iV;
1361  idNegDaughters[nBufEntry] = negInd[iV];
1362 
1363  daughterPosPDG[nBufEntry] = trackPdgPos[iPDGPos][iV];
1364  daughterNegPDG[nBufEntry] = trackPdgNeg[iV];
1365 
1366  pvIndexMother[nBufEntry] = isPrimary[iV] ? negPVIndex[iV] : -1;
1367 
1368  V0PDG[nBufEntry] = 420;
1369 
1370  nBufEntry++;
1371 
1372  if(int(nBufEntry) == float_vLen)
1373  {
1374  KFParticleDatabase::Instance()->GetMotherMass(V0PDG,massMotherPDG,massMotherPDGSigma);
1375  mother.SetPDG( V0PDG );
1376  ConstructV0(vTracks, trTypeIndexPos[iTrTypePos], trTypeIndexNeg[iTrTypeNeg],
1377  idPosDaughters, idNegDaughters, daughterPosPDG, daughterNegPDG,
1378  mother, mother_temp,
1379  nBufEntry, l, dl, Particles, PrimVtx,
1380  cuts, pvIndexMother, secCuts, massMotherPDG,
1381  massMotherPDGSigma, motherPrimSecCand, nPrimSecCand, vMotherPrim, vMotherSec);
1382  nBufEntry = 0;
1383  }
1384  }
1385 
1386  }//iV
1387  }//iPDGPos
1388  }//iRot
1389  }//iTrP
1390  }//iTrN
1391 
1392  if( nBufEntry>0 )
1393  {
1394  for(int iV=nBufEntry; iV<float_vLen; iV++)
1395  {
1396  idPosDaughters[iV] = idPosDaughters[0];
1397  idNegDaughters[iV] = idNegDaughters[0];
1398  }
1399 
1400  KFParticleDatabase::Instance()->GetMotherMass(V0PDG,massMotherPDG,massMotherPDGSigma);
1401  mother.SetPDG( V0PDG );
1402 
1403  ConstructV0(vTracks, trTypeIndexPos[iTrTypePos], trTypeIndexNeg[iTrTypeNeg],
1404  idPosDaughters, idNegDaughters, daughterPosPDG, daughterNegPDG,
1405  mother, mother_temp,
1406  nBufEntry, l, dl, Particles, PrimVtx,
1407  cuts, pvIndexMother, secCuts, massMotherPDG,
1408  massMotherPDGSigma, motherPrimSecCand, nPrimSecCand, vMotherPrim, vMotherSec);
1409  nBufEntry = 0;
1410  }
1411 
1412  if(nPrimSecCand>0)
1413  {
1414  SaveV0PrimSecCand(motherPrimSecCand,nPrimSecCand,mother_temp,PrimVtx,secCuts,vMotherPrim,vMotherSec);
1415  nPrimSecCand = 0;
1416  }
1417  }//iTC
1418  }//iTrTypeNeg
1419  }//iTrTypePos
1420 }
1421 
1422 void KFParticleFinder::ConstructPrimaryBG(KFPTrackVector* vTracks,
1423  vector<KFParticle>& Particles,
1424  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
1425  const float* cuts,
1426  const float* secCuts,
1427  vector< vector<KFParticle> >* vMotherPrim,
1428  vector<KFParticle>* vMotherSec )
1429 {
1447  KFParticle mother_temp;
1448  KFParticleSIMD mother;
1449  kfvector_floatv l(fNPV), dl(fNPV);
1450 
1451  KFParticleSIMD daughterNeg, daughterPos;
1452 
1453  // for secondary V0
1454  unsigned int nBufEntry = 0;
1455  float_v dS;
1456  uint_v idNegDaughters;
1457  uint_v idPosDaughters;
1458  int_v daughterPosPDG(-1);
1459  int_v daughterNegPDG(-1);
1460 
1461  int_v pvIndexMother(-1);
1462 
1463  float_v massMotherPDG(Vc::Zero), massMotherPDGSigma(Vc::Zero);
1464  int_v V0PDG(Vc::Zero);
1465 
1466  KFParticleSIMD motherPrimSecCand;
1467  int nPrimSecCand =0;
1468 
1469  for(int iSet=2; iSet<4; iSet++)
1470  {
1471  int signPDG = 1;
1472  if(iSet == 3)
1473  signPDG = -1;
1474 
1475  KFPTrackVector& positivePrimaryTracks = vTracks[iSet];
1476  int nPositiveTracks = positivePrimaryTracks.Size();
1477 
1478  for(int iTrN = positivePrimaryTracks.FirstPion(); iTrN < positivePrimaryTracks.LastProton(); iTrN += float_vLen)
1479  {
1480  int_v negPDG = positivePrimaryTracks.PDG()[iTrN];
1481  int_v negPVIndex = positivePrimaryTracks.PVIndex()[iTrN];
1482 
1483  int_m activeNeg = (negPDG != -1);
1484 
1485  for(int iTrP = iTrN+1; iTrP < positivePrimaryTracks.LastProton(); iTrP += float_vLen)
1486  {
1487  const int NTracks = (iTrP + float_vLen < nPositiveTracks) ? float_vLen : (nPositiveTracks - iTrP);
1488 
1489  int_v posPDG(0); // = reinterpret_cast<const int_v&>(positivePrimaryTracks.PDG()[iTrP]);
1490  int_v posPVIndex(0); // = reinterpret_cast<const int_v&>(positivePrimaryTracks.PVIndex()[iTrP]);
1491  for(int iV=0; iV<NTracks; iV++)
1492  {
1493  posPDG[iV] = positivePrimaryTracks.PDG()[iTrP+iV];
1494  posPVIndex[iV] = positivePrimaryTracks.PVIndex()[iTrP+iV];
1495  }
1496 
1497  int_m active = (activeNeg && (int_v::IndexesFromZero() < int_v(NTracks)));
1498  if(active.isEmpty()) continue;
1499 
1500  active &= (posPDG != int_v(-1));
1501  active &= (posPVIndex == negPVIndex);
1502 
1503  if(active.isEmpty()) continue;
1504 
1505  //detetrmine a pdg code of the mother particle
1506 
1507  int_v motherPDG(-1);
1508 
1509  motherPDG( (abs(posPDG)== 211) && int_m(abs(negPDG) == 211) ) = signPDG*9001; //pi+pi+
1510  motherPDG( (abs(posPDG)== 321) && int_m(abs(negPDG) == 211) ) = signPDG*9002; //pi+K+
1511  motherPDG( (abs(posPDG)== 211) && int_m(abs(negPDG) == 321) ) = signPDG*9002; //pi+K+
1512  motherPDG( (abs(posPDG)== 211) && int_m(abs(negPDG) == 2212) ) = signPDG*2224; //pi+p
1513  motherPDG( (abs(posPDG)== 2212) && int_m(abs(negPDG) == 211) ) = signPDG*2224; //pi+p
1514  motherPDG( (abs(posPDG)== 321) && int_m(abs(negPDG) == 321) ) = signPDG*9003; //K+K+
1515  motherPDG( (abs(posPDG)== 321) && int_m(abs(negPDG) == 2212) ) = signPDG*9004; //K+p
1516  motherPDG( (abs(posPDG)== 2212) && int_m(abs(negPDG) == 321) ) = signPDG*9004; //K+p
1517 
1518  active &= (motherPDG != -1);
1519  if(!(fDecayReconstructionList.empty()))
1520  {
1521  for(int iV=0; iV<float_vLen; iV++)
1522  {
1523  if(!(active[iV])) continue;
1524  if(fDecayReconstructionList.find(motherPDG[iV]) == fDecayReconstructionList.end())
1525  motherPDG[iV] = -1;
1526  }
1527  active &= (motherPDG != -1);
1528  }
1529  if(active.isEmpty()) continue;
1530 
1531 
1532  if(active.isEmpty()) continue;
1533 
1534  for(int iV=0; iV<NTracks; iV++)
1535  {
1536  if(!(active[iV])) continue;
1537 
1538  idPosDaughters[nBufEntry] = iTrP+iV;
1539  idNegDaughters[nBufEntry] = iTrN;
1540 
1541  daughterPosPDG[nBufEntry] = posPDG[iV];
1542  daughterNegPDG[nBufEntry] = negPDG[iV];
1543 
1544  pvIndexMother[nBufEntry] = negPVIndex[iV];
1545 
1546  V0PDG[nBufEntry] = motherPDG[iV];
1547 
1548  nBufEntry++;
1549 
1550  if(int(nBufEntry) == float_vLen)
1551  {
1552  KFParticleDatabase::Instance()->GetMotherMass(V0PDG,massMotherPDG,massMotherPDGSigma);
1553  mother.SetPDG( V0PDG );
1554  ConstructV0(vTracks, iSet, iSet,
1555  idPosDaughters, idNegDaughters, daughterPosPDG, daughterNegPDG,
1556  mother, mother_temp,
1557  nBufEntry, l, dl, Particles, PrimVtx,
1558  cuts, pvIndexMother, secCuts, massMotherPDG,
1559  massMotherPDGSigma, motherPrimSecCand, nPrimSecCand, vMotherPrim, vMotherSec);
1560  nBufEntry = 0;
1561  }
1562  }//iV
1563  }//iTrP
1564 
1565  if( nBufEntry>0 )
1566  {
1567  for(int iV=nBufEntry; iV<float_vLen; iV++)
1568  {
1569  idPosDaughters[iV] = idPosDaughters[0];
1570  idNegDaughters[iV] = idNegDaughters[0];
1571  }
1572 
1573  KFParticleDatabase::Instance()->GetMotherMass(V0PDG,massMotherPDG,massMotherPDGSigma);
1574  mother.SetPDG( V0PDG );
1575 
1576  ConstructV0(vTracks, iSet, iSet,
1577  idPosDaughters, idNegDaughters, daughterPosPDG, daughterNegPDG,
1578  mother, mother_temp,
1579  nBufEntry, l, dl, Particles, PrimVtx,
1580  cuts, pvIndexMother, secCuts, massMotherPDG,
1581  massMotherPDGSigma, motherPrimSecCand, nPrimSecCand, vMotherPrim, vMotherSec);
1582  nBufEntry = 0;
1583  }
1584 
1585  if(nPrimSecCand>0)
1586  {
1587  SaveV0PrimSecCand(motherPrimSecCand,nPrimSecCand,mother_temp,PrimVtx,secCuts,vMotherPrim,vMotherSec);
1588  nPrimSecCand = 0;
1589  }
1590  }
1591  }
1592 }
1593 
1594 void KFParticleFinder::ConstructTrackV0Cand(KFPTrackVector& vTracks,
1595  uint_v& idTracks,
1596  int_v& trackPDG,
1597  KFParticle* vV0[],
1598  KFParticleSIMD& mother,
1599  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& motherTopo,
1600  KFParticle& mother_temp,
1601  const unsigned short nElements,
1602  kfvector_floatv& l,
1603  kfvector_floatv& dl,
1604  std::vector<KFParticle>& Particles,
1605  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
1606  const float_v* cuts,
1607  const int_v& pvIndex,
1608  const float_v& massMotherPDG,
1609  const float_v& massMotherPDGSigma,
1610  std::vector< std::vector<KFParticle> >* vMotherPrim,
1611  std::vector<KFParticle>* vMotherSec)
1612 {
1635  float_m isPrimary(simd_cast<float_m>(pvIndex>-1));
1636 
1637  int_v trackId( &(vTracks.Id()[0]), idTracks );
1638  KFParticleSIMD V0(vV0,nElements);
1639  KFParticleSIMD track(vTracks, idTracks, trackPDG);
1640  track.SetId(trackId);
1641 
1642  float_m isSameParticle = simd_cast<float_m>((abs(mother.PDG()) == int_v(4122)) ||
1643  (abs(mother.PDG()) == int_v(114122)) ||
1644  (abs(mother.PDG()) == int_v(204122)) ||
1645  (abs(mother.PDG()) == int_v(504122)) ||
1646  (abs(mother.PDG()) == int_v(404122)) ||
1647  (abs(mother.PDG()) == int_v(425)) ||
1648  (abs(mother.PDG()) == int_v(427)) ||
1649  (abs(mother.PDG()) == int_v(200411)) ||
1650  (abs(mother.PDG()) == int_v(300411)) ||
1651  (abs(mother.PDG()) == int_v(300431)) ||
1652  (abs(mother.PDG()) == int_v(400431)) ||
1653  (abs(mother.PDG()) == int_v(411)) ||
1654  (abs(mother.PDG()) == int_v(431)) ||
1655  (abs(mother.PDG()) == int_v(429)) ||
1656  (abs(mother.PDG()) == int_v(1003334)) ||
1657  (abs(mother.PDG()) == int_v(3001)) ||
1658  (abs(mother.PDG()) == int_v(3006)) ||
1659  (abs(mother.PDG()) == int_v(3007)) ||
1660  (abs(mother.PDG()) == int_v(3009)) ||
1661  (abs(mother.PDG()) == int_v(3011)) );
1662  if( isSameParticle.isEmpty() )
1663  {
1664 #ifdef CBM
1665  float_v ds[2] = {0.f,0.f};
1666  float_v dsdr[4][6];
1667  track.GetDStoParticle( V0, ds, dsdr );
1668  track.TransportToDS(ds[0], dsdr[0]);
1669  V0.TransportToDS(ds[1], dsdr[3]);
1670 #endif
1671  const KFParticleSIMD* vDaughtersPointer[2] = {&track, &V0};
1672  mother.Construct(vDaughtersPointer, 2, 0);
1673  }
1674  else
1675  {
1676  int_v motherPDG = mother.PDG();
1677  mother = V0;
1678  mother.SetPDG(motherPDG);
1679  track.TransportToPoint(V0.Parameters());
1680  mother += track;
1681  }
1682 
1683  float_m active = simd_cast<float_m>(int_v::IndexesFromZero() < int(nElements));
1684 
1685  float_m saveParticle(active);
1686  saveParticle &= (mother.Chi2()/simd_cast<float_v>(mother.NDF()) < cuts[2] );
1687  saveParticle &= KFPMath::Finite(mother.GetChi2());
1688  saveParticle &= (mother.GetChi2() > 0.0f);
1689  saveParticle &= (mother.GetChi2() == mother.GetChi2());
1690 
1691  if( saveParticle.isEmpty() ) { return; }
1692 
1693  int_m isSameTrack(false);
1694  for(unsigned int iD=0; iD<V0.DaughterIds().size(); iD++)
1695  isSameTrack |= ( int_v(V0.DaughterIds()[iD]) == int_v(trackId) );
1696 
1697  saveParticle &= ( !simd_cast<float_m>(isSameTrack));
1698  if( saveParticle.isEmpty() ) { return; }
1699 
1700  float_v lMin(1.e8f);
1701  float_v ldlMin(1.e8f);
1702  float_m isParticleFromVertex(false);
1703 
1704  for(int iP=0; iP<fNPV; iP++)
1705  {
1706  float_m isParticleFromVertexLocal;
1707  mother.GetDistanceToVertexLine(PrimVtx[iP], l[iP], dl[iP], &isParticleFromVertexLocal);
1708  isParticleFromVertex |= isParticleFromVertexLocal;
1709  float_v ldl = (l[iP]/dl[iP]);
1710  lMin( (l[iP] < lMin) && active) = l[iP];
1711  ldlMin( (ldl < ldlMin) && active) = ldl;
1712  }
1713  saveParticle &= (lMin < 200.f);
1714  saveParticle &= ((float_m(!isPrimary) && isParticleFromVertex) || float_m(isPrimary) );
1715  if( saveParticle.isEmpty() ) { return; }
1716 
1717  isSameParticle = isSameParticle || isPrimary;
1718  if(!((isSameParticle).isFull()))
1719  {
1720  float_m isParticleFromVertexLocal;
1721  float_v l1, dl1;
1722  V0.GetDistanceToVertexLine(mother, l1, dl1, &isParticleFromVertexLocal);
1723 
1724  saveParticle &= ( isSameParticle || ((!isSameParticle) && isParticleFromVertexLocal));
1725  if( saveParticle.isEmpty() ) { return; }
1726  }
1727 
1728  saveParticle &= ( (float_m(!isPrimary) && ldlMin > cuts[0]) || float_m(isPrimary) );
1729 
1730  int_m setLCut = abs(mother.PDG()) == 3312 || abs(mother.PDG()) == 3334 || abs(mother.PDG()) == 3001;
1731  saveParticle &= ( (simd_cast<float_m>(setLCut) && lMin > float_v(fLCut)) || simd_cast<float_m>(!setLCut) );
1732 
1733  ldlMin = 1.e8f;
1734  for(int iP=0; iP<fNPV; iP++)
1735  {
1736  motherTopo[iP] = mother;
1737  motherTopo[iP].SetProductionVertex(PrimVtx[iP]);
1738  motherTopo[iP].GetDecayLength(l[iP], dl[iP]);
1739  float_v ldl = (l[iP]/dl[iP]);
1740  ldlMin( (ldl < ldlMin) && active) = ldl;
1741  }
1742 
1743  vector<int> iPrimVert[float_vLen];
1744  float_m isPrimaryPart(false);
1745 
1746  for(int iP=0; iP<fNPV; iP++)
1747  {
1748  const float_v& motherTopoChi2Ndf = motherTopo[iP].GetChi2()/simd_cast<float_v>(motherTopo[iP].GetNDF());
1749  const float_m isPrimaryPartLocal = ( motherTopoChi2Ndf < cuts[1] );
1750  isPrimaryPart |= isPrimaryPartLocal;
1751  for(int iV=0; iV<float_vLen; iV++)
1752  {
1753  if(isPrimaryPartLocal[iV])
1754  iPrimVert[iV].push_back(iP);
1755  }
1756  }
1757 
1758  for(unsigned int iv=0; iv<nElements; iv++)
1759  {
1760  if(!saveParticle[iv]) continue;
1761 
1762  mother.GetKFParticle(mother_temp, iv);
1763  if( mother.PDG()[iv] == 3312 )
1764  {
1765  fLPi.push_back(mother_temp);
1766  fLPiPIndex.push_back( V0.DaughterIds()[1][iv] );
1767  }
1768 
1769 // if( mother.PDG()[iv] == 100411 )
1770 // {
1771 // fK0PiPlus.push_back(mother_temp);
1772 // fK0PiMinusIndex.push_back( V0.DaughterIds()[0][iv] );
1773 // continue;
1774 // }
1775 
1776  //TODO check if needed with current implementation
1777  // reset daughter ids for 3- and 4-particle decays
1778  if( (abs(mother.PDG()[iv]) == 411) ||
1779  (abs(mother.PDG()[iv]) == 429) ||
1780  (abs(mother.PDG()[iv]) == 431) ||
1781  (abs(mother.PDG()[iv]) == 4122) ||
1782  (abs(mother.PDG()[iv]) == 114122) ||
1783  (abs(mother.PDG()[iv]) == 204122) ||
1784  (abs(mother.PDG()[iv]) == 314122) ||
1785  (abs(mother.PDG()[iv]) == 404122) ||
1786  (abs(mother.PDG()[iv]) == 504122) ||
1787  (abs(mother.PDG()[iv]) == 425) ||
1788  (abs(mother.PDG()[iv]) == 427) ||
1789  (abs(mother.PDG()[iv]) == 200411) ||
1790  (abs(mother.PDG()[iv]) == 300411) ||
1791  (abs(mother.PDG()[iv]) == 300431) ||
1792  (abs(mother.PDG()[iv]) == 400431) ||
1793  (abs(mother.PDG()[iv]) == 1003334) ||
1794  (abs(mother.PDG()[iv]) == 3001) )
1795  {
1796  mother_temp.CleanDaughtersId();
1797  for(int iD=0; iD < vV0[iv]->NDaughters(); iD++)
1798  mother_temp.AddDaughterId( vV0[iv]->DaughterIds()[iD] );
1799  mother_temp.AddDaughterId(trackId[iv]);
1800  }
1801 
1802  if( mother.PDG()[iv] == 411 )
1803  {
1804  fDPlus.push_back(mother_temp);
1805  continue;
1806  }
1807  if( mother.PDG()[iv] == -411 )
1808  {
1809  fDMinus.push_back(mother_temp);
1810  continue;
1811  }
1812  if( mother.PDG()[iv] == 300411 )
1813  {
1814  fDPlus3Pi.push_back(mother_temp);
1815  continue;
1816  }
1817  if( mother.PDG()[iv] == -300411 )
1818  {
1819  fDMinus3Pi.push_back(mother_temp);
1820  continue;
1821  }
1822  if( mother.PDG()[iv] == 400431 )
1823  {
1824  fDsPlusK2Pi.push_back(mother_temp);
1825  continue;
1826  }
1827  if( mother.PDG()[iv] == -400431 )
1828  {
1829  fDsMinusK2Pi.push_back(mother_temp);
1830  continue;
1831  }
1832  if( mother.PDG()[iv] == 504122 )
1833  {
1834  fLcPlusP2Pi.push_back(mother_temp);
1835  continue;
1836  }
1837  if( mother.PDG()[iv] == -504122 )
1838  {
1839  fLcMinusP2Pi.push_back(mother_temp);
1840  continue;
1841  }
1842 // if( mother.PDG()[iv] == 427 )
1843 // {
1844 // fK0PiPi.push_back(mother_temp);
1845 // continue;
1846 // }
1847 
1848  if( mother.PDG()[iv] == 429 )
1849  {
1850  fD04.push_back(mother_temp);
1851  continue;
1852  }
1853  if( mother.PDG()[iv] == -429 )
1854  {
1855  fD04bar.push_back(mother_temp);
1856  continue;
1857  }
1858 
1859  if( mother.PDG()[iv] == 3203 )
1860  fLLn.push_back(mother_temp);
1861  if( mother.PDG()[iv] == 3010 )
1862  fH5LL.push_back(mother_temp);
1863 
1864  mother_temp.SetId(Particles.size());
1865 
1866  if(!(isPrimaryPart[iv]))
1867  {
1868  if( vMotherSec )
1869  {
1870  float mass, errMass;
1871  mother_temp.GetMass(mass, errMass);
1872  if(abs(mother.PDG()[iv]) == 3324)
1873  {
1874  vMotherSec->push_back(mother_temp);
1875  }
1876  else
1877  {
1878  if( (fabs(mass - massMotherPDG[iv])/massMotherPDGSigma[iv]) <= 3 )
1879  {
1880  KFParticle mother_sec = mother_temp;
1881  mother_sec.SetNonlinearMassConstraint(massMotherPDG[iv]);
1882  vMotherSec->push_back(mother_temp);
1883  }
1884  }
1885  }
1886  if(!(mother.PDG()[iv] == 3006 || mother.PDG()[iv] == 3007))
1887  continue;
1888  }
1889 
1890  //check Ds+ and Lc+ candidates not to be D+
1891 // if(abs(mother_temp.GetPDG())==431 || abs(mother_temp.GetPDG())==4122)
1892 // {
1893 // KFPTrack dPionTrack;
1894 // vTracks.GetTrack(dPionTrack, idTracks[iv]);
1895 // KFParticle dPion(dPionTrack, 211);
1896 //
1897 // KFParticle dMeson = *vV0[iv];
1898 // dMeson += dPion;
1899 // float dMass, dMassError;
1900 // dMeson.GetMass(dMass, dMassError);
1901 // if(fabs(dMass - KFParticleDatabase::Instance()->GetDPlusMass())/KFParticleDatabase::Instance()->GetDPlusMassSigma() < 3) continue;
1902 // }
1903 
1904  if(abs(mother.GetPDG()[iv]) == 521 || abs(mother.GetPDG()[iv]) == 529 || abs(mother.GetPDG()[iv]) == 511 || abs(mother.GetPDG()[iv]) == 519)
1905  {
1906  KFParticle daughter_temp = *vV0[iv];
1907  float massPDG = KFParticleDatabase::Instance()->GetDPlusMass();
1908  float massSigmaPDG = KFParticleDatabase::Instance()->GetDPlusMassSigma();
1909  if(abs(mother.GetPDG()[iv]) == 521 || abs(mother.GetPDG()[iv]) == 529)
1910  {
1911  massPDG = KFParticleDatabase::Instance()->GetD0Mass();
1912  massSigmaPDG = KFParticleDatabase::Instance()->GetD0MassSigma();
1913  }
1914  float mass, dm;
1915  daughter_temp.GetMass(mass,dm);
1916  if( (fabs(mass - massPDG)/massSigmaPDG) > 3 ) continue;
1917 
1918 // KFParticleSIMD daughter_tempSIMD(daughter_temp);
1919 // daughter_tempSIMD.SetProductionVertex(PrimVtx[0]);
1920 // if(daughter_tempSIMD.GetChi2()[0]/daughter_tempSIMD.GetNDF()[0] < 3. ) continue;
1921 
1922  daughter_temp.SetId(Particles.size());
1923  daughter_temp.SetPDG(-1);
1924  mother_temp.SetId(Particles.size()+1);
1925  mother_temp.CleanDaughtersId();
1926  mother_temp.AddDaughterId(Particles.size());
1927  mother_temp.AddDaughterId(trackId[iv]);
1928  Particles.push_back(daughter_temp);
1929  }
1930  Particles.push_back(mother_temp);
1931 
1932  if( abs(mother.GetPDG()[iv]) == 3334 ) //Omega-
1933  {
1934  float mass, errMass;
1935  mother_temp.GetMass(mass, errMass);
1936 
1937  vector< vector<KFParticle> >* motherVector = &fPrimCandidates[7];
1938  if( mother.GetPDG()[iv] == 3334 )
1939  motherVector = &fPrimCandidates[8];
1940 
1941  mother_temp.SetNonlinearMassConstraint(massMotherPDG[iv]);
1942 
1943  if( (fabs(mass - massMotherPDG[iv])/massMotherPDGSigma[iv]) <= 3 )
1944  for(unsigned int iP=0; iP<iPrimVert[iv].size(); iP++)
1945  (*motherVector)[iPrimVert[iv][iP]].push_back(mother_temp);
1946  }
1947 
1948  if(vMotherPrim)
1949  {
1950  if( !((abs(mother.GetPDG()[iv]) == 3312) || (abs(mother.GetPDG()[iv]) == 3324))) continue;
1951  float mass, errMass;
1952 
1953  mother_temp.GetMass(mass, errMass);
1954  if(abs(mother.PDG()[iv]) == 3324)
1955  {
1956  for(unsigned int iP=0; iP<iPrimVert[iv].size(); iP++)
1957  (*vMotherPrim)[iPrimVert[iv][iP]].push_back(mother_temp);
1958  }
1959  else
1960  {
1961  mother_temp.SetNonlinearMassConstraint(massMotherPDG[iv]);
1962 
1963  if( (fabs(mass - massMotherPDG[iv])/massMotherPDGSigma[iv]) <= 3 )
1964  for(unsigned int iP=0; iP<iPrimVert[iv].size(); iP++)
1965  (*vMotherPrim)[iPrimVert[iv][iP]].push_back(mother_temp);
1966  }
1967  }
1968  }
1969 }
1970 
1971 void KFParticleFinder::FindTrackV0Decay(vector<KFParticle>& vV0,
1972  const int V0PDG,
1973  KFPTrackVector& vTracks,
1974  const int q,
1975  const int firstTrack,
1976  const int lastTrack,
1977  vector<KFParticle>& Particles,
1978  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
1979  int v0PVIndex,
1980  kfvector_float* ChiToPrimVtx,
1981  vector< vector<KFParticle> >* vMotherPrim,
1982  vector<KFParticle>* vMotherSec)
1983 {
1999  if( (vV0.size() < 1) || ((lastTrack-firstTrack) < 1) ) return;
2000  KFParticle mother_temp;
2001 
2002  KFParticle* v0Pointer[float_v::Size];
2003 
2004  KFParticleSIMD mother, track;
2005  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> > motherTopo(fNPV);
2006 
2007  kfvector_floatv l(fNPV), dl(fNPV);
2008 
2009  float_v cuts[3];
2010 
2011  // for secondary V0
2012  unsigned int nBufEntry = 0;
2013  float_v dS;
2014  uint_v idTrack;
2015  int_v trackPDGMother(-1);
2016 
2017  int_v pvIndexMother(-1);
2018 
2019  float_v massMotherPDG(Vc::Zero), massMotherPDGSigma(Vc::Zero);
2020  int_v motherParticlePDG(Vc::Zero);
2021 // Particles.reserve(Particles.size() + vV0.size());
2022 
2023  bool isCharm = ((abs(V0PDG) == 421) || (abs(V0PDG) == 411) || (abs(V0PDG) == 429) || (abs(V0PDG) == 420) || (abs(V0PDG) == 419)) && (v0PVIndex<0);
2024 
2025  for(unsigned int iV0=0; iV0 < vV0.size(); iV0++)
2026  {
2027  int iNegDaughter = vV0[iV0].DaughterIds()[0];
2028  int iPosDaughter = vV0[iV0].DaughterIds()[1];
2029 
2030  for(int iTr=firstTrack; iTr < lastTrack; iTr += float_vLen)
2031  {
2032  const int NTracks = (iTr + float_vLen < lastTrack) ? float_vLen : (lastTrack - iTr);
2033 
2034  const int_v& trackPDG = reinterpret_cast<const int_v&>(vTracks.PDG()[iTr]);
2035  const int_v& trackPVIndex = reinterpret_cast<const int_v&>(vTracks.PVIndex()[iTr]);
2036 
2037  const int_m& isTrackSecondary = (trackPVIndex < 0);
2038  const int_m& isSecondary = int_m( v0PVIndex < 0 ) && isTrackSecondary;
2039  const int_m& isPrimary = int_m( v0PVIndex >= 0 ) && (!isTrackSecondary);
2040  const int_m& isSamePV = (isPrimary && (v0PVIndex == trackPVIndex)) || !(isPrimary);
2041 
2042  float_m closeDaughters = simd_cast<float_m>(isSamePV) && simd_cast<float_m>(int_v::IndexesFromZero() < int(NTracks));
2043 
2044 // if(v0PVIndex < 0)
2045 // {
2046 // KFParticleSIMD v0(vV0[iV0]);
2047 // track.Load(vTracks, iTr, trackPDG);
2048 
2049 // float_v dsV0, dsTrack;
2050 // float_v dsdrV0[6] = {0.f,0.f,0.f,0.f,0.f,0.f};
2051 // float_v dsdrTrack[6] = {0.f,0.f,0.f,0.f,0.f,0.f};
2052 // float_v par1[8], cov1[36], par2[8], cov2[36];
2053 // v0.GetDStoParticle(track, dsV0, dsTrack);
2054 // v0.Transport(dsV0, dsdrV0, par1, cov1);
2055 // track.Transport(dsTrack, dsdrTrack, par2, cov2);
2056 //
2057 // const float_v& dx = par1[0] - par2[0];
2058 // const float_v& dy = par1[1] - par2[1];
2059 // const float_v& dz = par1[2] - par2[2];
2060 // const float_v& r2 = dx*dx + dy*dy + dz*dz;
2061 //
2062 // const float_v vtx[3] = {(par1[0] + par2[0])/2.f,
2063 // (par1[1] + par2[1])/2.f,
2064 // (par1[2] + par2[2])/2.f, };
2065 //
2066 // v0.CorrectErrorsOnS(par1, vtx, cov1);
2067 // track.CorrectErrorsOnS(par2, vtx, cov2);
2068 //
2069 // const float_v cov[6] = {cov1[0]+cov2[0],
2070 // cov1[1]+cov2[1],
2071 // cov1[2]+cov2[2],
2072 // cov1[3]+cov2[3],
2073 // cov1[4]+cov2[4],
2074 // cov1[5]+cov2[5] };
2075 // const float_v& err2 = cov[0]*dx*dx + cov[2]*dy*dy + cov[5]*dz*dz + 2.f*( cov[1]*dx*dy + cov[3]*dx*dz + cov[4]*dy*dz );
2076 //
2077 // closeDaughters &= ( (r2 < float_v(1.f)) && (r2*r2/err2) < float_v(3.f) && isSecondary);
2078 // closeDaughters &= v0.GetDeviationFromParticle(track) < float_v(10.f);
2079 // }
2080 
2081  if(v0PVIndex < 0)
2082  {
2083  KFParticleSIMD v0(vV0[iV0]);
2084  track.Load(vTracks, iTr, trackPDG);
2085  closeDaughters &= v0.GetDistanceFromParticle(track) < float_v(fDistanceCut);
2086  if(closeDaughters.isEmpty()) continue;
2087  }
2088 
2089  int_v trackPdgPos[2];
2090  int_m active[2];
2091 
2092  int nPDGPos = 2;
2093 
2094  active[0] = simd_cast<int_m>(closeDaughters);
2095  active[1] = (trackPDG == -1) && isSecondary && simd_cast<int_m>(closeDaughters);
2096 
2097  trackPdgPos[0] = trackPDG;
2098 
2099  if( (trackPDG == -1).isEmpty() || (abs(V0PDG) == 421) || (abs(V0PDG) == 411) )
2100  {
2101  nPDGPos = 1;
2102  }
2103  else
2104  {
2105  trackPdgPos[0](trackPDG == -1) = q*211;
2106  nPDGPos = 1;//TODO
2107  trackPdgPos[1](isSecondary) = q*321;
2108  }
2109 
2110  for(int iPDGPos=0; iPDGPos<nPDGPos; iPDGPos++)
2111  {
2112 
2113  if(active[iPDGPos].isEmpty()) continue;
2114 
2115  //detetrmine a pdg code of the mother particle
2116 
2117  int_v motherPDG(-1);
2118 
2119  if( V0PDG == 3122 )
2120  {
2121  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 3312;
2122  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 304122;
2123  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 321) ) = 3334;
2124  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = 3224;
2125  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = 3114;
2126  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -321) ) = 1003314;
2127  }
2128  else if( V0PDG == -3122 )
2129  {
2130  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = -3312;
2131  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -304122;
2132  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 321) ) = -3334;
2133  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = -3224;
2134  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = -3114;
2135  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 321) ) = -1003314;
2136  }
2137  else if( V0PDG == 310)
2138  {
2139  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = 323;
2140  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = -323;
2141  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 100411;
2142  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -100411;
2143  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 321) ) = 100431;
2144  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = -100431;
2145  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 104122;
2146  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -2212) ) = -104122;
2147  }
2148  else if( V0PDG == 3312 )
2149  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = 3324;
2150  else if( V0PDG == -3312)
2151  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = -3324;
2152  else if( V0PDG == 3324 )
2153  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -321) ) = 1003334;
2154  else if( V0PDG == -3324 )
2155  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 321) ) = -1003334;
2156  else if(V0PDG == 421)
2157  {
2158  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 411;
2159  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 321) ) = 431;
2160  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 4122;
2161  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2162  int_m isDMeson = isSecondary && int_m(trackPdgPos[iPDGPos] == 211);
2163  active[iPDGPos] &= (!(isDMeson)) || (isDMeson && ( id > iPosDaughter) );
2164  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -521;
2165  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = -529;
2166  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = 10411;
2167  }
2168  else if(V0PDG == -421)
2169  {
2170  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -411;
2171  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = -431;
2172  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] ==-2212) ) = -4122;
2173  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2174  int_m isDMeson = isSecondary && int_m(trackPdgPos[iPDGPos] == -211);
2175  active[iPDGPos] &= (!(isDMeson)) || (isDMeson && ( id > iNegDaughter) );
2176  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 521;
2177  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 321) ) = 529;
2178  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = -10411;
2179  }
2180  else if(V0PDG == 420 && q>0)
2181  {
2182  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 211) ) = 300411;
2183  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 321) ) = 400431;
2184  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 2212) ) = 504122;
2185  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2186  int_m isDMeson = isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 211);
2187  active[iPDGPos] &= (!(isDMeson)) || (isDMeson && ( id > iPosDaughter) );
2188  }
2189  else if(V0PDG == 420 && q<0)
2190  {
2191  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 211) ) = -300411;
2192  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 321) ) = -400431;
2193  motherPDG( isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 2212) ) = -504122;
2194  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2195  int_m isDMeson = isSecondary && int_m(abs(trackPdgPos[iPDGPos]) == 211);
2196  active[iPDGPos] &= (!(isDMeson)) || (isDMeson && ( id > iNegDaughter) );
2197  }
2198  else if(V0PDG == 411)
2199  {
2200  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 429;
2201  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = 10421;
2202  }
2203  else if(V0PDG == -411)
2204  {
2205  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = -429;
2206  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = -10421;
2207  }
2208  else if(V0PDG == 419)
2209  {
2210  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -511;
2211  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = -519;
2212  }
2213  else if(V0PDG == -419)
2214  {
2215  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 511;
2216  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 321) ) = 519;
2217  }
2218  else if(V0PDG == 429)
2219  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 211) ) = 20411;
2220  else if(V0PDG == -429)
2221  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -211) ) = -20411;
2222  else if( V0PDG == 3002 )
2223  {
2224  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2225  int_m isSameProton = (id == fLPiPIndex[iV0]);
2226  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) && (!isSameProton)) = 3001;
2227  }
2228  else if( V0PDG == 100411 )
2229  {
2230  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 425;
2231  }
2232  else if( V0PDG == 100431 )
2233  {
2234  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = 427;
2235  }
2236  else if( V0PDG == 425)
2237  {
2238  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 200411;
2239  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -200411;
2240  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 321) ) = 300431;
2241  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -321) ) = -300431;
2242  }
2243  else if( V0PDG == 111 )
2244  {
2245  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 3222;
2246  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -2212) ) = -3222;
2247  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == 321) ) = 100323;
2248  motherPDG( isPrimary && int_m(trackPdgPos[iPDGPos] == -321) ) = -100323;
2249  }
2250  else if( V0PDG == 3004 )
2251  {
2252  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 3006;
2253  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 3203;
2254  }
2255  else if( V0PDG == -3004 )
2256  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -2212) ) = -3006;
2257  else if( V0PDG == 3005 )
2258  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 3007;
2259  else if( V0PDG == -3005 )
2260  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -2212) ) = -3007;
2261  else if( V0PDG == 3006 )
2262  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 3008;
2263  else if( V0PDG == 3007 )
2264  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 3010;
2265  else if( V0PDG == 3203 )
2266  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) ) = 3009;
2267  else if( V0PDG == 3010 )
2268  {
2269  const int_v& id = reinterpret_cast<const int_v&>(vTracks.Id()[iTr]);
2270  int_m isSameProton = (id == Particles[vV0[iV0].DaughterIds()[1]].DaughterIds()[2]);
2271  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 2212) && (!isSameProton)) = 3011;
2272  }
2273  else if(V0PDG == 304122)
2274  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 314122;
2275  else if(V0PDG == -304122)
2276  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -314122;
2277  else if(V0PDG == 314122)
2278  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 404122;
2279  else if(V0PDG == -314122)
2280  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = -404122;
2281  else if(V0PDG == 104122)
2282  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = 114122;
2283  else if(V0PDG == -104122)
2284  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = -114122;
2285  else if(V0PDG == 114122)
2286  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == -211) ) = 204122;
2287  else if(V0PDG == -114122)
2288  motherPDG( isSecondary && int_m(trackPdgPos[iPDGPos] == 211) ) = -204122;
2289 
2290  active[iPDGPos] &= (motherPDG != -1);
2291  if(!(fDecayReconstructionList.empty()))
2292  {
2293  for(int iV=0; iV<float_vLen; iV++)
2294  {
2295  if(!(active[iPDGPos][iV])) continue;
2296  if(fDecayReconstructionList.find(motherPDG[iV]) == fDecayReconstructionList.end())
2297  motherPDG[iV] = -1;
2298  }
2299  active[iPDGPos] &= (motherPDG != -1);
2300  }
2301  if(ChiToPrimVtx)
2302  active[iPDGPos] &= ( !( (abs(motherPDG) == 3334 || abs(motherPDG) == 3312 ) ) ||
2303  ( (abs(motherPDG) == 3334 || abs(motherPDG) == 3312 ) && simd_cast<int_m>(reinterpret_cast<const float_v&>((*ChiToPrimVtx)[iTr]) > float_v(fCuts2D[0])) ) );
2304 
2305  if(active[iPDGPos].isEmpty()) continue;
2306 
2307  if(isCharm)
2308  {
2309  track.Load(vTracks, iTr, trackPDG);
2310  const float_v& trackPt = track.Px()*track.Px() + track.Py()*track.Py();
2311  const int_v& nPixelHits = reinterpret_cast<const int_v&>(vTracks.NPixelHits()[iTr]);
2312 
2313  active[iPDGPos] &= simd_cast<int_m>(trackPt >= fCutCharmPt*fCutCharmPt) && simd_cast<int_m>(reinterpret_cast<const float_v&>((*ChiToPrimVtx)[iTr]) > fCutCharmChiPrim ) && int_m(nPixelHits >= int_v(3));
2314  }
2315  {
2316  int_m isCharmParticle = (abs(motherPDG) == 104122) ||
2317  (abs(motherPDG) == 204122) ||
2318  (abs(motherPDG) == 304122) ||
2319  (abs(motherPDG) == 404122) ||
2320  (abs(motherPDG) == 425) ||
2321  (abs(motherPDG) == 426) ||
2322  (abs(motherPDG) == 427) ||
2323  (abs(motherPDG) == 100411) ||
2324  (abs(motherPDG) == 200411) ||
2325  (abs(motherPDG) == 100431) ||
2326  (abs(motherPDG) == 300431) ;
2327 
2328  if(!(isCharmParticle.isEmpty()))
2329  {
2330  track.Load(vTracks, iTr, trackPDG);
2331  const float_v& trackPt = track.Px()*track.Px() + track.Py()*track.Py();
2332  const int_v& nPixelHits = reinterpret_cast<const int_v&>(vTracks.NPixelHits()[iTr]);
2333 
2334  active[iPDGPos] &= ( (simd_cast<int_m>(trackPt >= fCutCharmPt*fCutCharmPt) && simd_cast<int_m>(reinterpret_cast<const float_v&>((*ChiToPrimVtx)[iTr]) > fCutCharmChiPrim ) && (nPixelHits >= int_v(3)) ) && isCharmParticle ) || (!isCharmParticle);
2335  }
2336  }
2337 
2338  for(int iV=0; iV<NTracks; iV++)
2339  {
2340  if(!(active[iPDGPos][iV])) continue;
2341 
2342 
2343  idTrack[nBufEntry] = iTr+iV;
2344  v0Pointer[nBufEntry] = &vV0[iV0];
2345 
2346  trackPDGMother[nBufEntry] = trackPdgPos[iPDGPos][iV];
2347 
2348  pvIndexMother[nBufEntry] = v0PVIndex;
2349 
2350  float massMother, massMotherSigma;
2351  KFParticleDatabase::Instance()->GetMotherMass(motherPDG[iV],massMother,massMotherSigma);
2352 
2353  massMotherPDG[nBufEntry] = massMother;
2354  massMotherPDGSigma[nBufEntry] = massMotherSigma;
2355  motherParticlePDG[nBufEntry] = motherPDG[iV];
2356 
2357  int motherType = 0;
2358 
2359  switch (abs(motherPDG[iV]))
2360  {
2361  case 3312: motherType = 0; break; //Xi
2362  case 3334: motherType = 0; break; //Omega
2363  case 4122: motherType = 1; break; //LambdaC
2364  case 104122: motherType = 1; break; //LambdaC
2365  case 304122: motherType = 1; break; //LambdaC
2366  case 504122: motherType = 1; break; //LambdaC
2367  case 425: motherType = 1; break; //D0
2368  case 426: motherType = 1; break; //D0
2369  case 427: motherType = 1; break; //D0
2370  case 100411: motherType = 1; break; //D+
2371  case 300411: motherType = 1; break; //D+
2372  case 100431: motherType = 1; break; //Ds+
2373  case 400431: motherType = 1; break; //Ds+
2374  case 431: motherType = 1; break; //Ds+-
2375  case 411: motherType = 1; break; //D+-
2376  case 428: motherType = 1; break; //D0
2377  case 429: motherType = 1; break; //D0
2378  case 521: motherType = 1; break; //B+
2379  case 529: motherType = 1; break; //B+
2380  case 511: motherType = 1; break; //B0
2381  case 519: motherType = 1; break; //B0
2382  case 3001: motherType = 1; break; //H0
2383  case 3222: motherType = 1; break; //Sigma+
2384  case 3006: motherType = 1; break; //He4L
2385  case 3007: motherType = 1; break; //He5L
2386  case 3008: motherType = 1; break; //H4LL
2387  case 3009: motherType = 1; break; //H4LL
2388  case 3011: motherType = 1; break; //He6LL
2389  default: motherType = 2; break; //resonances
2390  }
2391  for(int iCut=0; iCut<3; iCut++)
2392  cuts[iCut][nBufEntry] = fCutsTrackV0[motherType][iCut];
2393 
2394  nBufEntry++;
2395 
2396  if(int(nBufEntry) == float_vLen)
2397  {
2398  mother.SetPDG( motherParticlePDG );
2399  ConstructTrackV0Cand(vTracks,
2400  idTrack, trackPDGMother, v0Pointer,
2401  mother, motherTopo, mother_temp,
2402  nBufEntry, l, dl, Particles, PrimVtx,
2403  cuts, pvIndexMother, massMotherPDG,
2404  massMotherPDGSigma, vMotherPrim, vMotherSec);
2405  nBufEntry = 0;
2406  }
2407  }//iV
2408  }//iPDGPos
2409  }//iTr
2410  }
2411 
2412  if(nBufEntry > 0)
2413  {
2414  for(int iV=nBufEntry; iV<float_vLen; iV++)
2415  idTrack[iV] = idTrack[0];
2416 
2417  mother.SetPDG( motherParticlePDG );
2418  ConstructTrackV0Cand(vTracks,
2419  idTrack, trackPDGMother, v0Pointer,
2420  mother, motherTopo, mother_temp,
2421  nBufEntry, l, dl, Particles, PrimVtx,
2422  cuts, pvIndexMother, massMotherPDG,
2423  massMotherPDGSigma, vMotherPrim, vMotherSec);
2424  nBufEntry = 0;
2425  }
2426 }
2427 
2428 void KFParticleFinder::SelectParticles(vector<KFParticle>& Particles,
2429  vector<KFParticle>& vCandidates,
2430  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
2431  const float& cutChi2Topo,
2432  const float& cutLdL,
2433  const float& mass,
2434  const float& massErr,
2435  const float& massCut)
2436 {
2450  KFParticle* cand[float_vLen];
2451  int nCand = vCandidates.size();
2452 
2453  vector<KFParticle> newCandidates;
2454  kfvector_floatv l(fNPV), dl(fNPV);
2455 
2456  for(int iC=0; iC < nCand; iC += float_vLen)
2457  {
2458  int nEntries = (iC + float_vLen < nCand) ? float_vLen : (nCand - iC);
2459 
2460  for(int iv=0; iv<nEntries; iv++)
2461  cand[iv] = &vCandidates[iC+iv];
2462 
2463  KFParticleSIMD mother(cand,nEntries);
2464 
2465  float_m saveParticle(simd_cast<float_m>(int_v::IndexesFromZero() < int(nEntries)));
2466 
2467  float_v lMin(1.e8f);
2468  float_v ldlMin(1.e8f);
2469  float_m isParticleFromVertex(false);
2470 
2471  for(int iP=0; iP<fNPV; iP++)
2472  {
2473  float_m isParticleFromVertexLocal;
2474  mother.GetDistanceToVertexLine(PrimVtx[iP], l[iP], dl[iP], &isParticleFromVertexLocal);
2475  isParticleFromVertex |= isParticleFromVertexLocal;
2476  float_v ldl = (l[iP]/dl[iP]);
2477  lMin( (l[iP] < lMin) && saveParticle) = l[iP];
2478  ldlMin( (ldl < ldlMin) && saveParticle) = ldl;
2479  }
2480  saveParticle &= ldlMin > cutLdL;
2481  saveParticle &= (lMin < 200.f);
2482  saveParticle &= isParticleFromVertex;
2483  if( saveParticle.isEmpty() ) continue;
2484 
2485  KFParticleSIMD* candTopo = new KFParticleSIMD[fNPV];
2486 
2487  for(int iP=0; iP<fNPV; iP++)
2488  {
2489  candTopo[iP] = mother;
2490  candTopo[iP].SetProductionVertex(PrimVtx[iP]);
2491  }
2492 
2493  for(int iv=0; iv<nEntries; iv++)
2494  {
2495  if(!saveParticle[iv]) continue;
2496 
2497  bool isPrimary = 0;
2498  for(int iP=0; iP<fNPV; iP++)
2499  {
2500  if( !(KFPMath::Finite(candTopo[iP].GetChi2())[iv]) ) continue;
2501  if(!(candTopo[iP].GetChi2()[iv] > 0.0f)) continue;
2502  if(!(candTopo[iP].GetChi2()[iv]==candTopo[iP].GetChi2()[iv])) continue;
2503 
2504  if(float(candTopo[iP].GetChi2()[iv])/float(candTopo[iP].GetNDF()[iv]) <= cutChi2Topo )
2505  isPrimary = 1;
2506  }
2507  if(!isPrimary)
2508  continue;
2509 
2510  vCandidates[iC+iv].SetId(Particles.size());
2511  Particles.push_back(vCandidates[iC+iv]);
2512 
2513  float m, dm;
2514  vCandidates[iC+iv].GetMass(m,dm);
2515  if( (fabs(m - mass)/massErr) > massCut ) continue;
2516 
2517  vCandidates[iC+iv].SetNonlinearMassConstraint(mass);
2518  newCandidates.push_back(vCandidates[iC+iv]);
2519  }
2520  if(candTopo) delete [] candTopo;
2521  }
2522 
2523  vCandidates = newCandidates;
2524 }
2525 
2526 void KFParticleFinder::CombinePartPart(vector<KFParticle>& particles1,
2527  vector<KFParticle>& particles2,
2528  vector<KFParticle>& Particles,
2529  std::vector<KFParticleSIMD, KFPSimdAllocator<KFParticleSIMD> >& PrimVtx,
2530  const float* cuts,
2531  int iPV,
2532  const int MotherPDG,
2533  bool isSameInputPart,
2534  bool saveOnlyPrimary,
2535  vector< vector<KFParticle> >* vMotherPrim,
2536  vector<KFParticle>* vMotherSec,
2537  float massMotherPDG,
2538  float massMotherPDGSigma)
2539 {
2555  if( (particles1.size() == 0) || (particles2.size() == 0) ) return;
2556  if(!(fDecayReconstructionList.empty()) && (fDecayReconstructionList.find(MotherPDG) == fDecayReconstructionList.end())) return;
2557 
2558  KFParticle mother_temp;
2559  KFParticleSIMD mother;
2560  KFParticleSIMD *motherTopo = new KFParticleSIMD[fNPV];
2561  mother.SetPDG( MotherPDG );
2562 
2563  kfvector_floatv l(fNPV), dl(fNPV);
2564 
2565  KFParticle* tmpPart2[float_vLen];
2566  int nPart2 = particles2.size();
2567 
2568  bool isPrimary = (iPV >= 0);
2569  bool isCharm = (MotherPDG == 425) ||
2570  (MotherPDG == 427) ||
2571  (abs(MotherPDG) == 200411) ||
2572  (abs(MotherPDG) == 404122) ||
2573  (abs(MotherPDG) == 4132) ||
2574  (abs(MotherPDG) == 300431) ||
2575  (abs(MotherPDG) == 204122);
2576 
2577  for(unsigned int iP1=0; iP1 < particles1.size(); iP1++)
2578  {
2579  KFParticleSIMD vDaughters[2] = {KFParticleSIMD(particles1[iP1]), KFParticleSIMD()};
2580 
2581  unsigned int startIndex=0;
2582  if(isSameInputPart) startIndex=iP1+1;
2583  for(int iP2=startIndex; iP2 < nPart2; iP2 += float_vLen)
2584  {
2585  int nElements = (iP2 + float_vLen < nPart2) ? float_vLen : (nPart2 - iP2);
2586  float_m active(simd_cast<float_m>(int_v::IndexesFromZero() < int(nElements)));
2587 
2588  for(int iv=0; iv<nElements; iv++)
2589  tmpPart2[iv] = &particles2[iP2+iv];
2590 
2591  vDaughters[1] = KFParticleSIMD(tmpPart2,nElements);
2592 
2593 // if( reconstructPi0 )
2594 // {
2595 // int indexOffset = fEmcClusters->Id()[0];
2596 // uint_v gammaIndex1( (unsigned int)0);
2597 // uint_v gammaIndex2( (unsigned int)0);
2598 // for(int iv=0; iv<nElements; iv++)
2599 // {
2600 // gammaIndex1[iv] = Particles[ particles2[iP2+iv].DaughterIds()[0] ].DaughterIds()[0] - indexOffset;
2601 // gammaIndex2[iv] = Particles[ particles2[iP2+iv].DaughterIds()[1] ].DaughterIds()[0] - indexOffset;
2602 // }
2603 //
2604 // KFParticleSIMD gamma1(*fEmcClusters, gammaIndex1, vDaughters[0]);
2605 // KFParticleSIMD gamma2(*fEmcClusters, gammaIndex2, vDaughters[0]);
2606 // const KFParticleSIMD* pi0Daughters[2] = {&gamma1, &gamma2};
2607 //
2608 // int_v gammaId = vDaughters[1].Id();
2609 // vDaughters[1].SetVtxGuess(vDaughters[0].X(), vDaughters[0].Y(), vDaughters[0].Z());
2610 // vDaughters[1].Construct(pi0Daughters, 2);
2611 // vDaughters[1].SetId(gammaId);
2612 //
2613 // float_v mass, dm;
2614 // vDaughters[1].GetMass(mass,dm);
2615 // const float& mPi0 = KFParticleDatabase::Instance()->GetPi0Mass();
2616 // const float& mPi0Sigma = KFParticleDatabase::Instance()->GetPi0MassSigma();
2617 // active &= (abs(mass - mPi0)/mPi0Sigma) < 3.f;
2618 // vDaughters[1].SetNonlinearMassConstraint(mPi0);
2619 // if(active.isEmpty()) continue;
2620 // }
2621 
2622  if(isCharm)
2623  {
2624  mother = vDaughters[0];
2625  mother += vDaughters[1];
2626  mother.SetPDG( MotherPDG );
2627  }
2628  else
2629  {
2630  const KFParticleSIMD* vDaughtersPointer[2] = {&vDaughters[0], &vDaughters[1]};
2631  mother.Construct(vDaughtersPointer, 2, 0);
2632  }
2633 
2634  float_m saveParticle(active);
2635  saveParticle &= (mother.Chi2()/simd_cast<float_v>(mother.NDF()) < cuts[2] );
2636  saveParticle &= KFPMath::Finite(mother.GetChi2());
2637  saveParticle &= (mother.GetChi2() >= 0.0f);
2638  saveParticle &= (mother.GetChi2() == mother.GetChi2());
2639 
2640  if( saveParticle.isEmpty() ) { continue; }
2641 
2642  int_m isSameTrack(false);
2643  for(unsigned int iD=0; iD<vDaughters[0].DaughterIds().size(); iD++)
2644  for(unsigned int iD1=0; iD1<vDaughters[1].DaughterIds().size(); iD1++)
2645  isSameTrack |= ( int_v(vDaughters[0].DaughterIds()[iD]) == int_v(vDaughters[1].DaughterIds()[iD1]) );
2646  saveParticle &= ( !simd_cast<float_m>(isSameTrack));
2647  if( saveParticle.isEmpty() ) { continue; }
2648 
2649  float_v lMin(1.e8f);
2650  float_v ldlMin(1.e8f);
2651  float_m isParticleFromVertex(false);
2652 
2653  for(int iP=0; iP<fNPV; iP++)
2654  {
2655  if( (iPV > -1) && (iP !=iPV) ) continue;
2656  float_m isParticleFromVertexLocal;
2657  mother.GetDistanceToVertexLine(PrimVtx[iP], l[iP], dl[iP], &isParticleFromVertexLocal);
2658  isParticleFromVertex |= isParticleFromVertexLocal;
2659  float_v ldl = (l[iP]/dl[iP]);
2660  lMin( (l[iP] < lMin) && active) = l[iP];
2661  ldlMin( (ldl < ldlMin) && active) = ldl;
2662  }
2663  saveParticle &= ( (float_m(!isPrimary) && ldlMin > cuts[0]) || float_m(isPrimary) );
2664  saveParticle &= (lMin < 200.f);
2665 
2666  int_m setLCut = abs(mother.PDG()) == 3000;
2667  saveParticle &= ( (simd_cast<float_m>(setLCut) && lMin > float_v(fLCut)) || simd_cast<float_m>(!setLCut) );
2668 
2669 // if(isPrimary && (float(ldlMin > 3) )) continue;
2670  saveParticle &= ((float_m(!isPrimary) && isParticleFromVertex) || float_m(isPrimary) );
2671  if( saveParticle.isEmpty() ) { continue; }
2672 
2673  float_m isSameParticle(isPrimary || isCharm);
2674  if(!((isSameParticle).isFull()))
2675  {
2676  float_m isParticleFromVertexLocal;
2677  float_v l1, dl1;
2678  vDaughters[0].GetDistanceToVertexLine(mother, l1, dl1, &isParticleFromVertexLocal);
2679 
2680  saveParticle &= ( isSameParticle || ((!isSameParticle) && isParticleFromVertexLocal));
2681  if( saveParticle.isEmpty() ) { continue; }
2682  }
2683 
2684  for(int iP=0; iP<fNPV; iP++)
2685  {
2686  if( (iPV > -1) && (iP !=iPV) ) continue;
2687  motherTopo[iP] = mother;
2688  motherTopo[iP].SetProductionVertex(PrimVtx[iP]);
2689  }
2690 
2691  vector<int> iPrimVert[float_vLen];
2692  float_m isPrimaryPart(false);
2693 
2694  for(int iP=0; iP<fNPV; iP++)
2695  {
2696  if( (iPV > -1) && (iP !=iPV) ) continue;
2697  const float_v& motherTopoChi2Ndf = motherTopo[iP].GetChi2()/simd_cast<float_v>(motherTopo[iP].GetNDF());
2698  const float_m isPrimaryPartLocal = ( motherTopoChi2Ndf < float_v(cuts[1]) );
2699  isPrimaryPart |= isPrimaryPartLocal;
2700  for(int iV=0; iV<float_vLen; iV++)
2701  {
2702  if(isPrimaryPartLocal[iV])
2703  iPrimVert[iV].push_back(iP);
2704  }
2705  }
2706 
2707  for(int iv=0; iv<nElements; iv++)
2708  {
2709  if(!saveParticle[iv]) continue;
2710 
2711  mother.GetKFParticle(mother_temp, iv);
2712 
2713  // reset daughter ids for 3- and 4-particle decays
2714  if( (abs(mother.PDG()[iv]) == 428))
2715  {
2716  mother_temp.CleanDaughtersId();
2717  for(int iD=0; iD < particles1[iP1].NDaughters(); iD++)
2718  mother_temp.AddDaughterId( particles1[iP1].DaughterIds()[iD] );
2719  mother_temp.AddDaughterId(tmpPart2[iv]->Id());
2720  }
2721 
2722  if(saveOnlyPrimary)
2723  {
2724  if(isPrimaryPart[iv])
2725  {
2726  mother_temp.SetId(Particles.size());
2727  Particles.push_back(mother_temp);
2728  }
2729  }
2730  else
2731  {
2732  mother_temp.SetId(Particles.size());
2733  Particles.push_back(mother_temp);
2734  }
2735 
2736  if(vMotherPrim || vMotherSec)
2737  {
2738  float mass, errMass;
2739  mother_temp.GetMass(mass, errMass);
2740  if( (fabs(mass - massMotherPDG)/massMotherPDGSigma) > 3.f ) continue;
2741  mother_temp.SetNonlinearMassConstraint(massMotherPDG);
2742 
2743  if(MotherPDG == 428)
2744  {
2745  mother_temp.CleanDaughtersId();
2746  for(int iD=0; iD < tmpPart2[iv]->NDaughters(); iD++)
2747  mother_temp.AddDaughterId( tmpPart2[iv]->DaughterIds()[iD] );
2748  for(int iD=0; iD < particles1[iP1].NDaughters(); iD++)
2749  mother_temp.AddDaughterId( particles1[iP1].DaughterIds()[iD] );
2750  }
2751 
2752  if(vMotherSec && (!(isPrimaryPart[iv])) )
2753  vMotherSec->push_back(mother_temp);
2754  if(vMotherPrim)
2755  for(unsigned int iP=0; iP<iPrimVert[iv].size(); iP++)
2756  (*vMotherPrim)[iPrimVert[iv][iP]].push_back(mother_temp);
2757  }
2758  }
2759  }
2760  }
2761 
2762  if(motherTopo) delete [] motherTopo;
2763 }
2764 
2765 
2766 void KFParticleFinder::NeutralDaughterDecay(KFPTrackVector* vTracks,
2767  vector<KFParticle>& Particles)
2768 {
2781  KFParticle mother_temp;
2782  KFParticleSIMD ChargedDaughter, MotherTrack;
2783 
2784  uint_v idMotherTrack;
2785  uint_v idChargedDaughter;
2786  int_v ChargedDaughterPDG(-1);
2787 
2788  int_v pvIndexMother(-1);
2789 
2790  int outNeutralDaughterPDG[4][5]; //[iTC][iHypothesis]
2791  int outMotherPDG[4][5];
2792 
2793  int trTypeIndexMother[2] = {6,7};
2794  int trTypeIndexDaughter[2] = {0,1};
2795 
2796  for( int iTrTypeDaughter = 0; iTrTypeDaughter<2; iTrTypeDaughter++)
2797  {
2798  KFPTrackVector& DaughterTracks = vTracks[ trTypeIndexDaughter[iTrTypeDaughter] ];
2799  KFPTrackVector& MotherTracks = vTracks[ trTypeIndexMother[iTrTypeDaughter] ];
2800 
2801  int_v DaughterTracksSize = DaughterTracks.Size();
2802  int MotherTracksSize = MotherTracks.Size();
2803 
2804  //track categories
2805  int nTC = 4;
2806  int startTCMother[4] = {0,0,0,0};
2807  int endTCMother[4] = {0,0,0,0};
2808  int startTCDaughter[4] = {0,0,0,0};
2809  int endTCDaughter[4] = {0,0,0,0};
2810 
2811  nTC = 4;
2812  vector<int> nMotherHypothesis(nTC,0);
2813  vector< vector<int> > motherPDGHypothesis(nTC);
2814  vector< vector<float> > neutralDaughterMassHypothesis(nTC);
2815 
2816 
2817  //mu+, mu-
2818  startTCMother[0] = 0; endTCMother[0] = MotherTracksSize;
2819  startTCDaughter[0] = DaughterTracks.FirstMuon(); endTCDaughter[0] = DaughterTracks.LastMuon();
2820 
2821  nMotherHypothesis[0] = 2;
2822 
2823 
2824  motherPDGHypothesis[0].push_back(211);
2825  motherPDGHypothesis[0].push_back(321);
2826 
2827  neutralDaughterMassHypothesis[0].push_back(0.);
2828  neutralDaughterMassHypothesis[0].push_back(0.);
2829 
2830  outNeutralDaughterPDG[0][0]=-7000014;
2831  outNeutralDaughterPDG[0][1]=-8000014;
2832 
2833  outMotherPDG[0][0]=-7000211;
2834  outMotherPDG[0][1]=-7000321;
2835 
2836  //Pi+, Pi-
2837  startTCMother[1] = 0; endTCMother[1] = MotherTracksSize;
2838  startTCDaughter[1] = DaughterTracks.FirstPion(); endTCDaughter[1] = DaughterTracks.LastPion();
2839 
2840  nMotherHypothesis[1] = 5;
2841 
2842  motherPDGHypothesis[1].push_back(3112);
2843  motherPDGHypothesis[1].push_back(3222);
2844  motherPDGHypothesis[1].push_back(3312);
2845  motherPDGHypothesis[1].push_back(3334);
2846  motherPDGHypothesis[1].push_back(321);
2847 
2848  neutralDaughterMassHypothesis[1].push_back(0.939565);
2849  neutralDaughterMassHypothesis[1].push_back(0.939565);
2850  neutralDaughterMassHypothesis[1].push_back(1.115683);
2851  neutralDaughterMassHypothesis[1].push_back(1.31486);
2852  neutralDaughterMassHypothesis[1].push_back(0.1349766);
2853 
2854  outNeutralDaughterPDG[1][0]= 7002112;
2855  outNeutralDaughterPDG[1][1]=-8002112;
2856  outNeutralDaughterPDG[1][2]= 7003122;
2857  outNeutralDaughterPDG[1][3]= 7003322;
2858  outNeutralDaughterPDG[1][4]=-9000111;
2859 
2860  outMotherPDG[1][0]= 7003112;
2861  outMotherPDG[1][1]=-7003222;
2862  outMotherPDG[1][2]= 7003312;
2863  outMotherPDG[1][3]= 7003334;
2864  outMotherPDG[1][4]=-9000321;
2865 
2866  //K+, K-
2867  startTCMother[2] = 0; endTCMother[2] = MotherTracksSize;
2868  startTCDaughter[2] = DaughterTracks.FirstKaon(); endTCDaughter[2] = DaughterTracks.LastKaon();
2869 
2870  nMotherHypothesis[2] = 1;
2871 
2872  motherPDGHypothesis[2].push_back(3334);
2873 
2874  neutralDaughterMassHypothesis[2].push_back(1.115683);
2875 
2876  outNeutralDaughterPDG[2][0]= 8003122;
2877 
2878  outMotherPDG[2][0]= 8003334;
2879 
2880  //p+, p-
2881  startTCMother[3] = 0; endTCMother[3] = MotherTracksSize;
2882  startTCDaughter[3] = DaughterTracks.FirstProton(); endTCDaughter[3] = DaughterTracks.LastProton();
2883 
2884  nMotherHypothesis[3] = 1;
2885 
2886  motherPDGHypothesis[3].push_back(3222);
2887 
2888  neutralDaughterMassHypothesis[3].push_back(0.1349766);
2889 
2890  outNeutralDaughterPDG[3][0]=-8000111;
2891 
2892  outMotherPDG[3][0]=-8003222;
2893 
2894 
2895 
2896  for(int iTC=0; iTC<nTC; iTC++)
2897  {
2898  for(unsigned short iTrD=startTCDaughter[iTC]; iTrD < endTCDaughter[iTC]; iTrD += float_vLen)
2899  {
2900  const unsigned short NTracksDaughter = (iTrD + float_vLen < DaughterTracks.Size()) ? float_vLen : (DaughterTracks.Size() - iTrD);
2901 
2902  int_v DaughterInd = int_v::IndexesFromZero() + int(iTrD);
2903 
2904  int_v DaughterPDG = reinterpret_cast<const int_v&>(DaughterTracks.PDG()[iTrD]);
2905  int_v DaughterPVIndex = reinterpret_cast<const int_v&>(DaughterTracks.PVIndex()[iTrD]);
2906  int_v daughterId = reinterpret_cast<const int_v&>(DaughterTracks.Id()[iTrD]);
2907 
2908  int_v trackPdgDaughter = DaughterPDG;
2909  int_m activeDaughter = (DaughterPDG != -1);
2910 
2911  if( !((DaughterPDG == -1).isEmpty()) )
2912  {
2913  trackPdgDaughter(DaughterPVIndex<0 && (DaughterPDG == -1) ) = 211;
2914 
2915 // activeDaughter |= int_m(DaughterPVIndex < 0) && int_m(DaughterPDG == -1) ;
2916  }
2917 
2918  activeDaughter = (int_v::IndexesFromZero() < int(NTracksDaughter));
2919 
2920  ChargedDaughter.Load(DaughterTracks, iTrD, DaughterPDG);
2921  ChargedDaughter.SetId(daughterId);
2922 
2923  for(unsigned short iTrM=startTCMother[iTC]; iTrM < endTCMother[iTC]; iTrM += float_vLen)
2924  {
2925  const unsigned short NTracks = (iTrM + float_vLen < MotherTracksSize) ? float_vLen : (MotherTracksSize - iTrM);
2926 
2927  //const int_v& MotherPDG = reinterpret_cast<const int_v&>(MotherTracks.PDG()[iTrM]);
2928  //const int_v& MotherPVIndex = reinterpret_cast<const int_v&>(MotherTracks.PVIndex()[iTrM]);
2929  const int_v& motherTrackId = reinterpret_cast<const int_v&>(MotherTracks.Id()[iTrM]);
2930 
2931  for(int iRot = 0; iRot<float_vLen; iRot++)
2932  {
2933  if(iRot>0)
2934  {
2935  DaughterPDG = DaughterPDG.rotated(1);
2936  DaughterPVIndex = DaughterPVIndex.rotated(1);
2937  DaughterInd = DaughterInd.rotated(1);
2938  trackPdgDaughter = trackPdgDaughter.rotated(1);
2939 
2940  ChargedDaughter.Rotate();
2941 
2942  activeDaughter = /*( (DaughterPDG != -1) || ( (DaughterPVIndex < 0) && (DaughterPDG == -1) ) ) &&*/ (DaughterInd < DaughterTracksSize);
2943  }
2944 
2945  int_v trackPdgMother;
2946 
2947  if(iTC==0)
2948  activeDaughter &= abs(DaughterPDG)==13;
2949  if(iTC==1)
2950  activeDaughter &= abs(DaughterPDG)==211;
2951  if(iTC==2)
2952  activeDaughter &= abs(DaughterPDG)==321;
2953  if(iTC==3)
2954  activeDaughter &= abs(DaughterPDG)==2212;
2955  if (activeDaughter.isEmpty()) continue;
2956 
2957 
2958  for(int iHypothesis=0; iHypothesis<nMotherHypothesis[iTC]; iHypothesis++)
2959  {
2960  int motherKFPDG = outMotherPDG[iTC][iHypothesis];
2961  if(iTrTypeDaughter==0) motherKFPDG = -outMotherPDG[iTC][iHypothesis];
2962  if(!(fDecayReconstructionList.empty()) && (fDecayReconstructionList.find(motherKFPDG) == fDecayReconstructionList.end())) continue;
2963 
2964  int_m active = activeDaughter && (int_v::IndexesFromZero() < int(NTracks));
2965 
2966  MotherTrack.Load(MotherTracks, iTrM, motherPDGHypothesis[iTC][iHypothesis]);
2967 
2968  float_v zMother = MotherTrack.Z();
2969  float_v zCD = ChargedDaughter.Z();
2970 
2971  //daughter particle should start after the last hit of a mother track
2972  active &= simd_cast<int_m>(zCD >= (zMother - float_v(0.5f)));
2973  if( active.isEmpty() ) continue;
2974 
2975  KFParticleSIMD neutralDaughter = MotherTrack;
2976  //energy of the mother particle should be greater then of the daughter particle
2977  active &= simd_cast<int_m>(neutralDaughter.E() > ChargedDaughter.E());
2978  if( active.isEmpty() ) continue;
2979 
2980  neutralDaughter.AddDaughterId(motherTrackId);
2981  neutralDaughter.NDF() = -1;
2982  neutralDaughter.Chi2() = 0.f;
2983  neutralDaughter.SubtractDaughter(ChargedDaughter);
2984 
2985  //decay point shoud be between mother and daughter tracks
2986  active &= simd_cast<int_m>(neutralDaughter.Z() >= zMother - float_v(10.0f));
2987  active &= simd_cast<int_m>(neutralDaughter.Z() <= zCD + float_v(10.0f));
2988  //set cut on chi2 of the fit of the neutral daughter
2989  active &= simd_cast<int_m>(neutralDaughter.NDF() >= int_v(Vc::Zero));
2990  active &= simd_cast<int_m>(neutralDaughter.Chi2()/simd_cast<float_v>(neutralDaughter.NDF()) <= fCuts2D[1]);
2991  //fit should converge
2992  active &= simd_cast<int_m>(neutralDaughter.Chi2() >= float_v(Vc::Zero));
2993  active &= simd_cast<int_m>(neutralDaughter.Chi2() == neutralDaughter.Chi2());
2994  if( active.isEmpty() ) continue;
2995 
2996  //kill particle-candidates produced by clones
2997  active &= simd_cast<int_m>( neutralDaughter.GetRapidity()<6.f && neutralDaughter.GetRapidity()>0.f);
2998  if ((iTC==1 && iHypothesis<4) || iTC==2)
2999  active &= simd_cast<int_m>( !( (neutralDaughter.GetPt())<0.5f && neutralDaughter.GetRapidity()<0.5f ) );
3000  if (iTC==3)
3001  active &= simd_cast<int_m>( !( (neutralDaughter.GetPt())<0.2f && neutralDaughter.GetRapidity()<1.f ) );
3002  if( active.isEmpty() ) continue;
3003 
3004  KFParticleSIMD neutralDaughterUnconstr = neutralDaughter;
3005  neutralDaughter.SetNonlinearMassConstraint(neutralDaughterMassHypothesis[iTC][iHypothesis]);
3006 
3007  const KFParticleSIMD* daughters[2] = {&neutralDaughter, &ChargedDaughter};
3008  KFParticleSIMD mother;
3009  mother.Construct(daughters, 2);
3010 
3011  //decay point shoud be between mother and daughter tracks
3012  active &= simd_cast<int_m>(mother.Z() >= zMother);
3013  active &= simd_cast<int_m>(mother.Z() <= zCD);
3014  //set cut on chi2 of the fit of the mother particle
3015  active &= simd_cast<int_m>(mother.NDF() >= int_v(Vc::Zero));
3016  active &= simd_cast<int_m>(mother.Chi2()/simd_cast<float_v>(mother.NDF()) <= fCuts2D[1]);
3017  //fit should converge
3018  active &= simd_cast<int_m>(mother.Chi2() >= float_v(Vc::Zero));
3019  active &= simd_cast<int_m>(mother.Chi2() == mother.Chi2());
3020  if( active.isEmpty() ) continue;
3021 
3022  for(int iV=0; iV<NTracks; iV++)
3023  {
3024  if(!active[iV]) continue;
3025 
3026  neutralDaughterUnconstr.GetKFParticle(mother_temp, iV);
3027  int neutralId = Particles.size();
3028  mother_temp.SetId(neutralId);
3029  if (iTrTypeDaughter==0)
3030  mother_temp.SetPDG(-outNeutralDaughterPDG[iTC][iHypothesis]);
3031  else
3032  mother_temp.SetPDG(outNeutralDaughterPDG[iTC][iHypothesis]);
3033  Particles.push_back(mother_temp);
3034 
3035  mother.GetKFParticle(mother_temp, iV);
3036  mother_temp.SetId(Particles.size());
3037  mother_temp.CleanDaughtersId();
3038  mother_temp.AddDaughterId(ChargedDaughter.Id()[iV]);
3039  mother_temp.AddDaughterId(neutralId);
3040 
3041  if (iTrTypeDaughter==0)
3042  mother_temp.SetPDG(-outMotherPDG[iTC][iHypothesis]);
3043  else
3044  mother_temp.SetPDG(outMotherPDG[iTC][iHypothesis]);
3045  Particles.push_back(mother_temp);
3046  }
3047  }
3048  }//iRot
3049  }//iTrM
3050  }//iTrD
3051  }//iTC
3052  }//iTrTypeDaughter
3053 }
3054 
3055 void KFParticleFinder::AddCandidate(const KFParticle& candidate, int iPV)
3056 {
3071  //0 Ks, 1 Lambda,2 LambdaBar, 3 gamma, 4 pi0, 5 Xi, 6 XiBar, 7 Omega, 8 OmegaBar, 9 XiStar
3072  int iSet = -1;
3073  if(candidate.GetPDG() == 310) iSet = 0;
3074  if(candidate.GetPDG() == 3122) iSet = 1;
3075  if(candidate.GetPDG() == -3122) iSet = 2;
3076  if(candidate.GetPDG() == 22) iSet = 3;
3077  if(candidate.GetPDG() == 111) iSet = 4;
3078  if(candidate.GetPDG() == 3312) iSet = 5;
3079  if(candidate.GetPDG() == -3312) iSet = 6;
3080  if(candidate.GetPDG() == 3334) iSet = 7;
3081  if(candidate.GetPDG() == -3334) iSet = 8;
3082 
3083  if(iSet > -1)
3084  {
3085  if(iPV >= 0 && iPV<fNPV)
3086  {
3087  if(candidate.NDF() == 2)
3088  fPrimCandidates[iSet][iPV].push_back(candidate);
3089  if(candidate.NDF() == 3)
3090  fPrimCandidatesTopo[iSet][iPV].push_back(candidate);
3091  if(candidate.NDF() == 4)
3092  fPrimCandidatesTopoMass[iSet][iPV].push_back(candidate);
3093  }
3094  else if(iPV < 0)
3095  {
3096  fSecCandidates[iSet].push_back(candidate);
3097  }
3098  }
3099 }
3100 
3101 void KFParticleFinder::SetNPV(int nPV)
3102 {
3108  fNPV = nPV;
3109 
3110  for(int iCandidates=0; iCandidates<fNPrimCandidatesSets; iCandidates++)
3111  {
3112  fPrimCandidates[iCandidates].clear();
3113  fPrimCandidates[iCandidates].resize(fNPV);
3114  }
3115 
3116  for(int iCandidates=0; iCandidates<fNPrimCandidatesTopoSets; iCandidates++)
3117  {
3118  fPrimCandidatesTopo[iCandidates].clear();
3119  fPrimCandidatesTopo[iCandidates].resize(fNPV);
3120 
3121  fPrimCandidatesTopoMass[iCandidates].clear();
3122  fPrimCandidatesTopoMass[iCandidates].resize(fNPV);
3123  }
3124 }