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