CbmRoot
CbmKFParticleFinder.cxx
Go to the documentation of this file.
1 //-----------------------------------------------------------
2 //-----------------------------------------------------------
3 
4 // Cbm Headers ----------------------
5 #include "CbmKFParticleFinder.h"
6 #include "CbmEvent.h"
7 #include "CbmKFParticleFinderPID.h"
8 #include "CbmKFVertex.h"
9 #include "CbmL1PFFitter.h"
10 #include "CbmMCDataArray.h"
11 #include "CbmMCDataManager.h"
12 #include "CbmMCEventList.h"
13 #include "CbmMCTrack.h"
14 #include "FairRunAna.h"
15 #include "L1Field.h"
16 
17 //KF Particle headers
18 #include "KFPTrackVector.h"
19 #include "KFParticleTopoReconstructor.h"
20 
21 #include "FairLogger.h"
22 
23 //ROOT headers
24 #include "TClonesArray.h" //to get arrays from the FairRootManager
25 #include "TMath.h" //to calculate Prob function
26 #include "TStopwatch.h" //to measure the time
27 
28 //c++ and std headers
29 #include <cmath>
30 #include <iostream>
31 #include <vector>
32 using std::vector;
33 
34 CbmKFParticleFinder::CbmKFParticleFinder(const char* name, Int_t iVerbose)
35  : FairTask(name, iVerbose)
36  , fStsTrackBranchName("StsTrack")
37  , fTrackArray(0)
38  , fEvents(0)
39  , fTopoReconstructor(0)
40  , fPVFindMode(2)
41  , fPID(0)
42  , fSuperEventAnalysis(0)
43  , fSETracks(0)
44  , fSEField(0)
45  , fSEpdg(0)
46  , fSETrackId(0)
47  , fSEChiPrim(0)
48  , fTimeSliceMode(0) {
49  fTopoReconstructor = new KFParticleTopoReconstructor;
50 
51  // set default cuts
52  SetPrimaryProbCut(0.0001); // 0.01% to consider primary track as a secondary;
53 }
54 
55 CbmKFParticleFinder::~CbmKFParticleFinder() {
56  if (fTopoReconstructor) delete fTopoReconstructor;
57 }
58 
59 InitStatus CbmKFParticleFinder::Init() {
60  //Get ROOT Manager
61  FairRootManager* ioman = FairRootManager::Instance();
62 
63  if (ioman == 0) {
64  Error("CbmKFParticleFinder::Init", "RootManager not instantiated!");
65  return kERROR;
66  }
67 
68  //check the mode
69  fTimeSliceMode = 0;
70  if (ioman->CheckBranch("Event")) {
71  fTimeSliceMode = 1;
72  LOG(info) << GetName() << ": Running in the timeslice mode.";
73  } else
74  LOG(info) << GetName() << ": Running in the event by event mode.";
75 
76  // Get reconstructed events
77  if (fTimeSliceMode) {
78  fEvents = (TClonesArray*) ioman->GetObject("Event");
79  if (fEvents == 0)
80  Fatal("CbmKFParticleFinder::Init",
81  "No events available. Running in the event-by-event mode.");
82  }
83 
84  // Get input collection
85  fTrackArray = (TClonesArray*) ioman->GetObject(fStsTrackBranchName);
86  if (fTrackArray == 0) {
87  Error("CbmKFParticleFinder::Init", "track-array not found!");
88  return kERROR;
89  }
90 
91  //In case of reconstruction of pure signal no PV is defined. The MC PV is used.
92  if (fPVFindMode == 0) {
93  if (fTimeSliceMode) {
94  CbmMCDataManager* mcManager =
95  (CbmMCDataManager*) ioman->GetObject("MCDataManager");
96  if (mcManager == 0)
97  Error("CbmKFParticleFinder::Init", "MC Data Manager not found!");
98 
99  fMCTrackArray = mcManager->InitBranch("MCTrack");
100 
101  if (fMCTrackArray == 0) {
102  Error("CbmKFParticleFinder::Init", "MC track array not found!");
103  return kERROR;
104  }
105 
106  fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
107  if (fEventList == 0) {
108  Error("CbmKFParticleFinder::Init", "MC Event List not found!");
109  return kERROR;
110  }
111  } else {
112  fMCTrackArrayEvent = (TClonesArray*) ioman->GetObject("MCTrack");
113  if (fMCTrackArrayEvent == 0) {
114  Error("CbmKFParticleFinder::Init", "MC track array not found!");
115  return kERROR;
116  }
117  }
118  }
119 
120  fCbmPrimVertex = (CbmVertex*) ioman->GetObject("PrimaryVertex.");
121 
122  return kSUCCESS;
123 }
124 
125 void CbmKFParticleFinder::Exec(Option_t* /*opt*/) {
126  fTopoReconstructor->Clear();
127 
128  int nEvents = 1;
129  if (fTimeSliceMode) nEvents = fEvents->GetEntriesFast();
130 
131  vector<KFParticleTopoReconstructor> eventTopoReconstructor;
132  eventTopoReconstructor.resize(nEvents);
133 
134  for (int iEvent = 0; iEvent < nEvents; iEvent++) {
135  CbmEvent* event = 0;
136  if (fTimeSliceMode) event = (CbmEvent*) fEvents->At(iEvent);
137  eventTopoReconstructor[iEvent].SetChi2PrimaryCut(
138  InversedChi2Prob(0.0001, 2));
139  eventTopoReconstructor[iEvent].CopyCuts(fTopoReconstructor);
140  eventTopoReconstructor[iEvent].GetKFParticleFinder()->SetReconstructionList(
141  fTopoReconstructor->GetKFParticleFinder()->GetReconstructionList());
142 
143  int nTracksEvent = 0;
144  if (fTimeSliceMode)
145  nTracksEvent = event->GetNofStsTracks();
146  else
147  nTracksEvent = fTrackArray->GetEntriesFast();
148 
149  Int_t ntracks = 0; //fTrackArray->GetEntriesFast();
150 
151  //calculate number of d-He4 candidates
152  int nCandidatesDHe4 = 0;
153  if (fPID)
154  for (int iTr = 0; iTr < nTracksEvent; iTr++)
155  if (TMath::Abs(fPID->GetPID()[iTr]) == 1000010029) nCandidatesDHe4++;
156 
157  vector<CbmStsTrack> vRTracks(nTracksEvent + nCandidatesDHe4);
158  vector<int> pdg(nTracksEvent + nCandidatesDHe4, -1);
159  vector<int> trackId(nTracksEvent + nCandidatesDHe4, -1);
160 
161  for (int iTr = 0; iTr < nTracksEvent; iTr++) {
162  int stsTrackIndex = 0;
163  if (fTimeSliceMode)
164  stsTrackIndex = event->GetStsTrackIndex(iTr);
165  else
166  stsTrackIndex = iTr;
167 
168  CbmStsTrack* stsTrack = ((CbmStsTrack*) fTrackArray->At(stsTrackIndex));
169  const FairTrackParam* parameters = stsTrack->GetParamFirst();
170 
171  Double_t V[15] = {0.f};
172  for (Int_t i = 0, iCov = 0; i < 5; i++)
173  for (Int_t j = 0; j <= i; j++, iCov++)
174  V[iCov] = parameters->GetCovariance(i, j);
175 
176  if (stsTrack->GetNofHits() < 3) continue;
177 
178  bool ok = 1;
179  ok = ok && finite(parameters->GetX());
180  ok = ok && finite(parameters->GetY());
181  ok = ok && finite(parameters->GetZ());
182  ok = ok && finite(parameters->GetTx());
183  ok = ok && finite(parameters->GetTy());
184  ok = ok && finite(parameters->GetQp());
185 
186  for (unsigned short iC = 0; iC < 15; iC++)
187  ok = ok && finite(V[iC]);
188  ok = ok && (V[0] < 1. && V[0] > 0.) && (V[2] < 1. && V[2] > 0.)
189  && (V[5] < 1. && V[5] > 0.) && (V[9] < 1. && V[9] > 0.)
190  && (V[14] < 1. && V[14] > 0.);
191  ok = ok && stsTrack->GetChiSq() < 10 * stsTrack->GetNDF();
192  if (!ok) continue;
193 
194  if (fPID) {
195  if (fPID->GetPID()[stsTrackIndex] == -2) continue;
196 
197  //not clear separation between d and He4
198  if (TMath::Abs(fPID->GetPID()[stsTrackIndex]) == 1000010029) {
199  int sgn = fPID->GetPID()[stsTrackIndex]
200  / TMath::Abs(fPID->GetPID()[stsTrackIndex]);
201  pdg[ntracks] = sgn * 1000010020;
202  vRTracks[ntracks] = *stsTrack;
203  trackId[ntracks] = stsTrackIndex;
204  ntracks++;
205 
206  pdg[ntracks] = sgn * 1000020040;
207  } else
208  pdg[ntracks] = fPID->GetPID()[stsTrackIndex];
209  } else
210  pdg[ntracks] = -1;
211 
212  vRTracks[ntracks] = *stsTrack;
213  trackId[ntracks] = stsTrackIndex;
214 
215  ntracks++;
216  }
217 
218  vRTracks.resize(ntracks);
219  pdg.resize(ntracks);
220  trackId.resize(ntracks);
221 
222  CbmL1PFFitter fitter;
223  vector<float> vChiToPrimVtx;
224  CbmKFVertex kfVertex;
225  if (fPVFindMode == 0) {
226  int nMCEvents = 1;
227  if (fTimeSliceMode) nMCEvents = fEventList->GetNofEvents();
228 
229  float mcpv[3] = {0.f, 0.f, 0.f};
230  bool isMCPVFound = false;
231  for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) {
232  int nMCTracks = 0;
233  if (fTimeSliceMode)
234  nMCTracks = fMCTrackArray->Size(0, iMCEvent);
235  else
236  nMCTracks = fMCTrackArrayEvent->GetEntriesFast();
237 
238  for (Int_t iMC = 0; iMC < nMCTracks; iMC++) {
239  CbmMCTrack* cbmMCTrack;
240  if (fTimeSliceMode)
241  cbmMCTrack = (CbmMCTrack*) fMCTrackArray->Get(0, iMCEvent, iMC);
242  else
243  cbmMCTrack = (CbmMCTrack*) fMCTrackArrayEvent->At(iMC);
244 
245  if (cbmMCTrack->GetMotherId() < 0) {
246  mcpv[0] = cbmMCTrack->GetStartX();
247  mcpv[1] = cbmMCTrack->GetStartY();
248  mcpv[2] = cbmMCTrack->GetStartZ();
249  isMCPVFound = true;
250  break;
251  }
252  }
253  if (isMCPVFound) break;
254  }
255 
256  kfVertex.GetRefX() = mcpv[0];
257  kfVertex.GetRefY() = mcpv[1];
258  kfVertex.GetRefZ() = mcpv[2];
259  }
260 
261  if (fPVFindMode == 3) {
262  kfVertex.GetRefX() = fCbmPrimVertex->GetX();
263  kfVertex.GetRefY() = fCbmPrimVertex->GetY();
264  kfVertex.GetRefZ() = fCbmPrimVertex->GetZ();
265  }
266 
267  vector<L1FieldRegion> vField, vFieldAtLastPoint;
268  fitter.Fit(vRTracks, pdg);
269  fitter.GetChiToVertex(vRTracks, vField, vChiToPrimVtx, kfVertex, 3);
270  fitter.CalculateFieldRegionAtLastPoint(vRTracks, vFieldAtLastPoint);
271  vector<KFFieldVector> vFieldVector(ntracks),
272  vFieldVectorAtLastPoint(ntracks);
273  for (Int_t iTr = 0; iTr < ntracks; iTr++) {
274  int entrSIMD = iTr % fvecLen;
275  int entrVec = iTr / fvecLen;
276  vFieldVector[iTr].fField[0] = vField[entrVec].cx0[entrSIMD];
277  vFieldVector[iTr].fField[1] = vField[entrVec].cx1[entrSIMD];
278  vFieldVector[iTr].fField[2] = vField[entrVec].cx2[entrSIMD];
279  vFieldVector[iTr].fField[3] = vField[entrVec].cy0[entrSIMD];
280  vFieldVector[iTr].fField[4] = vField[entrVec].cy1[entrSIMD];
281  vFieldVector[iTr].fField[5] = vField[entrVec].cy2[entrSIMD];
282  vFieldVector[iTr].fField[6] = vField[entrVec].cz0[entrSIMD];
283  vFieldVector[iTr].fField[7] = vField[entrVec].cz1[entrSIMD];
284  vFieldVector[iTr].fField[8] = vField[entrVec].cz2[entrSIMD];
285  vFieldVector[iTr].fField[9] = vField[entrVec].z0[entrSIMD];
286  }
287  for (Int_t iTr = 0; iTr < ntracks; iTr++) {
288  int entrSIMD = iTr % fvecLen;
289  int entrVec = iTr / fvecLen;
290  vFieldVectorAtLastPoint[iTr].fField[0] =
291  vFieldAtLastPoint[entrVec].cx0[entrSIMD];
292  vFieldVectorAtLastPoint[iTr].fField[1] =
293  vFieldAtLastPoint[entrVec].cx1[entrSIMD];
294  vFieldVectorAtLastPoint[iTr].fField[2] =
295  vFieldAtLastPoint[entrVec].cx2[entrSIMD];
296  vFieldVectorAtLastPoint[iTr].fField[3] =
297  vFieldAtLastPoint[entrVec].cy0[entrSIMD];
298  vFieldVectorAtLastPoint[iTr].fField[4] =
299  vFieldAtLastPoint[entrVec].cy1[entrSIMD];
300  vFieldVectorAtLastPoint[iTr].fField[5] =
301  vFieldAtLastPoint[entrVec].cy2[entrSIMD];
302  vFieldVectorAtLastPoint[iTr].fField[6] =
303  vFieldAtLastPoint[entrVec].cz0[entrSIMD];
304  vFieldVectorAtLastPoint[iTr].fField[7] =
305  vFieldAtLastPoint[entrVec].cz1[entrSIMD];
306  vFieldVectorAtLastPoint[iTr].fField[8] =
307  vFieldAtLastPoint[entrVec].cz2[entrSIMD];
308  vFieldVectorAtLastPoint[iTr].fField[9] =
309  vFieldAtLastPoint[entrVec].z0[entrSIMD];
310  }
311 
312  if (!fSuperEventAnalysis) {
313  KFPTrackVector tracks;
314  FillKFPTrackVector(
315  &tracks, vRTracks, vFieldVector, pdg, trackId, vChiToPrimVtx);
316  KFPTrackVector tracksAtLastPoint;
317  FillKFPTrackVector(&tracksAtLastPoint,
318  vRTracks,
319  vFieldVectorAtLastPoint,
320  pdg,
321  trackId,
322  vChiToPrimVtx,
323  0);
324 
325  TStopwatch timer;
326  timer.Start();
327 
328  eventTopoReconstructor[iEvent].Init(tracks, tracksAtLastPoint);
329 
330  if (fPVFindMode == 0 || fPVFindMode == 3) {
331  KFPVertex primVtx_tmp;
332  primVtx_tmp.SetXYZ(
333  kfVertex.GetRefX(), kfVertex.GetRefY(), kfVertex.GetRefZ());
334  primVtx_tmp.SetCovarianceMatrix(0, 0, 0, 0, 0, 0);
335  primVtx_tmp.SetNContributors(0);
336  primVtx_tmp.SetChi2(-100);
337 
338  vector<int> pvTrackIds;
339  KFVertex pv(primVtx_tmp);
340  eventTopoReconstructor[iEvent].AddPV(pv, pvTrackIds);
341  } else if (fPVFindMode == 1)
342  eventTopoReconstructor[iEvent].ReconstructPrimVertex();
343  else if (fPVFindMode == 2)
344  eventTopoReconstructor[iEvent].ReconstructPrimVertex(0);
345 
346  eventTopoReconstructor[iEvent].SortTracks();
347  eventTopoReconstructor[iEvent].ReconstructParticles();
348 
349  timer.Stop();
350  eventTopoReconstructor[iEvent].SetTime(timer.RealTime());
351  } else {
352  for (int iTr = 0; iTr < ntracks; iTr++) {
353  const FairTrackParam* parameters = vRTracks[iTr].GetParamFirst();
354  float a = parameters->GetTx(), b = parameters->GetTy(),
355  qp = parameters->GetQp();
356  Int_t q = 0;
357  if (qp > 0.f) q = 1;
358  if (qp < 0.f) q = -1;
359  float c2 = 1.f / (1.f + a * a + b * b);
360  float pq = 1.f / qp * TMath::Abs(q);
361  float p2 = pq * pq;
362  float pz = sqrt(p2 * c2);
363  float px = a * pz;
364  float py = b * pz;
365  float pt = sqrt(px * px + py * py);
366 
367  bool save = 0;
368 
369  if (vChiToPrimVtx[iTr] < 3) {
370  if ((fabs(pdg[iTr]) == 11 && pt > 0.2f) || (fabs(pdg[iTr]) == 13)
371  || (fabs(pdg[iTr]) == 19))
372  save = 1;
373  }
374 
375  if (vChiToPrimVtx[iTr] > 3) {
376  if ((fabs(pdg[iTr]) == 211 || fabs(pdg[iTr]) == 321
377  || fabs(pdg[iTr]) == 2212 || pdg[iTr] == -1)
378  && pt > 0.2f)
379  save = 1;
380  }
381 
382  if (save) {
383  fSETracks.push_back(vRTracks[iTr]);
384  fSEField.push_back(vFieldVector[iTr]);
385  fSEpdg.push_back(pdg[iTr]);
386  fSETrackId.push_back(fSETrackId.size());
387  fSEChiPrim.push_back(vChiToPrimVtx[iTr]);
388  }
389  }
390  }
391  }
392 
393 
394  vector<int> PVTrackIndexArray;
395  int indexAdd = 0;
396 
397  for (int iEvent = 0; iEvent < nEvents; iEvent++) {
398  const KFParticleTopoReconstructor* eventTR =
399  &eventTopoReconstructor[iEvent];
400 
401  for (int iPV = 0; iPV < eventTR->NPrimaryVertices(); iPV++) {
402  PVTrackIndexArray = eventTR->GetPVTrackIndexArray(iPV);
403  for (unsigned int iTr = 0; iTr < PVTrackIndexArray.size(); iTr++)
404  PVTrackIndexArray[iTr] = PVTrackIndexArray[iTr] + indexAdd;
405 
406  fTopoReconstructor->AddPV(eventTR->GetPrimKFVertex(iPV),
407  PVTrackIndexArray);
408  PVTrackIndexArray.clear();
409  }
410  for (unsigned int iP = 0; iP < eventTR->GetParticles().size(); iP++) {
411  KFParticle particle;
412  const KFParticle& particleEvent = eventTR->GetParticles()[iP];
413  particle = particleEvent;
414  particle.CleanDaughtersId();
415  int idP = particleEvent.Id() + indexAdd;
416  int idDaughter = 0;
417  for (int nD = 0; nD < particleEvent.NDaughters(); nD++) {
418  if (particleEvent.NDaughters() == 1)
419  idDaughter = particleEvent.DaughterIds()[nD];
420  if (particleEvent.NDaughters() > 1)
421  idDaughter = particleEvent.DaughterIds()[nD] + indexAdd;
422  particle.AddDaughterId(idDaughter);
423  }
424  particle.SetId(idP);
425  fTopoReconstructor->AddParticle(particle);
426  }
427  indexAdd += eventTR->GetParticles().size();
428  }
429 }
430 
431 void CbmKFParticleFinder::Finish() {
432  if (fSuperEventAnalysis) {
433  KFPTrackVector tracks;
434  FillKFPTrackVector(
435  &tracks, fSETracks, fSEField, fSEpdg, fSETrackId, fSEChiPrim);
436  KFPTrackVector tracksAtLastPoint;
437 
438  LOG(info) << "CbmKFParticleFinder: Start SE analysis";
439  TStopwatch timer;
440  timer.Start();
441 
442  fTopoReconstructor->Init(tracks, tracksAtLastPoint);
443 
444  KFPVertex primVtx_tmp;
445  primVtx_tmp.SetXYZ(0, 0, 0);
446  primVtx_tmp.SetCovarianceMatrix(0, 0, 0, 0, 0, 0);
447  primVtx_tmp.SetNContributors(0);
448  primVtx_tmp.SetChi2(-100);
449  vector<int> pvTrackIds;
450  KFVertex pv(primVtx_tmp);
451  fTopoReconstructor->AddPV(pv, pvTrackIds);
452 
453  fTopoReconstructor->SortTracks();
454  fTopoReconstructor->ReconstructParticles();
455 
456  timer.Stop();
457  fTopoReconstructor->SetTime(timer.RealTime());
458  LOG(info) << "CbmKFParticleFinder: Finish SE analysis" << timer.RealTime();
459  }
460 }
461 
462 void CbmKFParticleFinder::FillKFPTrackVector(
463  KFPTrackVector* tracks,
464  const vector<CbmStsTrack>& vRTracks,
465  const vector<KFFieldVector>& vField,
466  const vector<int>& pdg,
467  const vector<int>& trackId,
468  const vector<float>& vChiToPrimVtx,
469  bool atFirstPoint) const {
470  int ntracks = vRTracks.size();
471  tracks->Resize(ntracks);
472  //fill vector with tracks
473  for (Int_t iTr = 0; iTr < ntracks; iTr++) {
474  const FairTrackParam* parameters;
475  if (atFirstPoint)
476  parameters = vRTracks[iTr].GetParamFirst();
477  else
478  parameters = vRTracks[iTr].GetParamLast();
479 
480  double par[6] = {0.f};
481 
482  double tx = parameters->GetTx(), ty = parameters->GetTy(),
483  qp = parameters->GetQp();
484 
485  Int_t q = 0;
486  if (qp > 0.f) q = 1;
487  if (qp < 0.f) q = -1;
488  if (TMath::Abs(pdg[iTr]) == 1000020030
489  || TMath::Abs(pdg[iTr]) == 1000020040)
490  q *= 2;
491 
492 
493  double c2 = 1.f / (1.f + tx * tx + ty * ty);
494  double pq = 1.f / qp * TMath::Abs(q);
495  double p2 = pq * pq;
496  double pz = sqrt(p2 * c2);
497  double px = tx * pz;
498  double py = ty * pz;
499 
500  par[0] = parameters->GetX();
501  par[1] = parameters->GetY();
502  par[2] = parameters->GetZ();
503  par[3] = px;
504  par[4] = py;
505  par[5] = pz;
506 
507  //calculate covariance matrix
508  double t = sqrt(1.f + tx * tx + ty * ty);
509  double t3 = t * t * t;
510  double dpxdtx = q / qp * (1.f + ty * ty) / t3;
511  double dpxdty = -q / qp * tx * ty / t3;
512  double dpxdqp = -q / (qp * qp) * tx / t;
513  double dpydtx = -q / qp * tx * ty / t3;
514  double dpydty = q / qp * (1.f + tx * tx) / t3;
515  double dpydqp = -q / (qp * qp) * ty / t;
516  double dpzdtx = -q / qp * tx / t3;
517  double dpzdty = -q / qp * ty / t3;
518  double dpzdqp = -q / (qp * qp * t3);
519 
520  double F[6][5] = {{1.f, 0.f, 0.f, 0.f, 0.f},
521  {0.f, 1.f, 0.f, 0.f, 0.f},
522  {0.f, 0.f, 0.f, 0.f, 0.f},
523  {0.f, 0.f, dpxdtx, dpxdty, dpxdqp},
524  {0.f, 0.f, dpydtx, dpydty, dpydqp},
525  {0.f, 0.f, dpzdtx, dpzdty, dpzdqp}};
526 
527  double VFT[5][6];
528  for (int i = 0; i < 5; i++)
529  for (int j = 0; j < 6; j++) {
530  VFT[i][j] = 0;
531  for (int k = 0; k < 5; k++) {
532  VFT[i][j] += parameters->GetCovariance(i, k) * F[j][k];
533  }
534  }
535 
536  double cov[21];
537  for (int i = 0, l = 0; i < 6; i++)
538  for (int j = 0; j <= i; j++, l++) {
539  cov[l] = 0;
540  for (int k = 0; k < 5; k++) {
541  cov[l] += F[i][k] * VFT[k][j];
542  }
543  }
544 
545  for (Int_t iP = 0; iP < 6; iP++)
546  tracks->SetParameter(par[iP], iP, iTr);
547  for (Int_t iC = 0; iC < 21; iC++)
548  tracks->SetCovariance(cov[iC], iC, iTr);
549  for (Int_t iF = 0; iF < 10; iF++)
550  tracks->SetFieldCoefficient(vField[iTr].fField[iF], iF, iTr);
551  tracks->SetId(trackId[iTr], iTr);
552  tracks->SetPDG(pdg[iTr], iTr);
553  tracks->SetQ(q, iTr);
554  tracks->SetNPixelHits(vRTracks[iTr].GetNofMvdHits(), iTr);
555 
556  if (fPVFindMode == 0 || fPVFindMode == 3) {
557  if (vChiToPrimVtx[iTr] < 3)
558  tracks->SetPVIndex(0, iTr);
559  else
560  tracks->SetPVIndex(-1, iTr);
561  } else
562  tracks->SetPVIndex(-1, iTr);
563  }
564 }
565 
566 double CbmKFParticleFinder::InversedChi2Prob(double p, int ndf) const {
567  double epsilon = 1.e-14;
568  double chi2Left = 0.f;
569  double chi2Right = 10000.f;
570 
571  double probLeft = p - TMath::Prob(chi2Left, ndf);
572 
573  double chi2Centr = (chi2Left + chi2Right) / 2.f;
574  double probCentr = p - TMath::Prob(chi2Centr, ndf);
575 
576  while (TMath::Abs(chi2Right - chi2Centr) / chi2Centr > epsilon) {
577  if (probCentr * probLeft > 0.f) {
578  chi2Left = chi2Centr;
579  probLeft = probCentr;
580  } else {
581  chi2Right = chi2Centr;
582  }
583 
584  chi2Centr = (chi2Left + chi2Right) / 2.f;
585  probCentr = p - TMath::Prob(chi2Centr, ndf);
586  }
587 
588  return chi2Centr;
589 }
590 
591 void CbmKFParticleFinder::SetPrimaryProbCut(float prob) {
592  fTopoReconstructor->SetChi2PrimaryCut(InversedChi2Prob(prob, 2));
593 }
594 
595 void CbmKFParticleFinder::SetSuperEventAnalysis() {
596  fSuperEventAnalysis = 1;
597  fPVFindMode = 0;
598  fTopoReconstructor->SetMixedEventAnalysis();
599 }
600 
601 KFParticleFinder* CbmKFParticleFinder::GetKFParticleFinder() {
602  return fTopoReconstructor->GetKFParticleFinder();
603 }
604 void CbmKFParticleFinder::SetMaxDistanceBetweenParticlesCut(float cut) {
605  GetKFParticleFinder()->SetMaxDistanceBetweenParticlesCut(cut);
606 }
607 void CbmKFParticleFinder::SetLCut(float cut) {
608  GetKFParticleFinder()->SetLCut(cut);
609 }
610 void CbmKFParticleFinder::SetChiPrimaryCut2D(float cut) {
611  GetKFParticleFinder()->SetChiPrimaryCut2D(cut);
612 }
613 void CbmKFParticleFinder::SetChi2Cut2D(float cut) {
614  GetKFParticleFinder()->SetChi2Cut2D(cut);
615 }
616 void CbmKFParticleFinder::SetLdLCut2D(float cut) {
617  GetKFParticleFinder()->SetLdLCut2D(cut);
618 }
619 void CbmKFParticleFinder::SetLdLCutXiOmega(float cut) {
620  GetKFParticleFinder()->SetLdLCutXiOmega(cut);
621 }
622 void CbmKFParticleFinder::SetChi2TopoCutXiOmega(float cut) {
623  GetKFParticleFinder()->SetChi2TopoCutXiOmega(cut);
624 }
625 void CbmKFParticleFinder::SetChi2CutXiOmega(float cut) {
626  GetKFParticleFinder()->SetChi2CutXiOmega(cut);
627 }
628 void CbmKFParticleFinder::SetChi2TopoCutResonances(float cut) {
629  GetKFParticleFinder()->SetChi2TopoCutResonances(cut);
630 }
631 void CbmKFParticleFinder::SetChi2CutResonances(float cut) {
632  GetKFParticleFinder()->SetChi2CutResonances(cut);
633 }
634 void CbmKFParticleFinder::SetPtCutLMVM(float cut) {
635  GetKFParticleFinder()->SetPtCutLMVM(cut);
636 }
637 void CbmKFParticleFinder::SetPCutLMVM(float cut) {
638  GetKFParticleFinder()->SetPCutLMVM(cut);
639 }
640 void CbmKFParticleFinder::SetPtCutJPsi(float cut) {
641  GetKFParticleFinder()->SetPtCutJPsi(cut);
642 }
643 void CbmKFParticleFinder::SetPtCutCharm(float cut) {
644  GetKFParticleFinder()->SetPtCutCharm(cut);
645 }
646 void CbmKFParticleFinder::SetChiPrimaryCutCharm(float cut) {
647  GetKFParticleFinder()->SetChiPrimaryCutCharm(cut);
648 }
649 void CbmKFParticleFinder::SetLdLCutCharmManybodyDecays(float cut) {
650  GetKFParticleFinder()->SetLdLCutCharmManybodyDecays(cut);
651 }
652 void CbmKFParticleFinder::SetChi2TopoCutCharmManybodyDecays(float cut) {
653  GetKFParticleFinder()->SetChi2TopoCutCharmManybodyDecays(cut);
654 }
655 void CbmKFParticleFinder::SetChi2CutCharmManybodyDecays(float cut) {
656  GetKFParticleFinder()->SetChi2CutCharmManybodyDecays(cut);
657 }
658 void CbmKFParticleFinder::SetLdLCutCharm2D(float cut) {
659  GetKFParticleFinder()->SetLdLCutCharm2D(cut);
660 }
661 void CbmKFParticleFinder::SetChi2TopoCutCharm2D(float cut) {
662  GetKFParticleFinder()->SetChi2TopoCutCharm2D(cut);
663 }
664 void CbmKFParticleFinder::SetChi2CutCharm2D(float cut) {
665  GetKFParticleFinder()->SetChi2CutCharm2D(cut);
666 }
667 void CbmKFParticleFinder::AddDecayToReconstructionList(int pdg) {
668  GetKFParticleFinder()->AddDecayToReconstructionList(pdg);
669 }
670 
671 ClassImp(CbmKFParticleFinder);
CbmKFParticleFinder::fMCTrackArray
CbmMCDataArray * fMCTrackArray
Definition: CbmKFParticleFinder.h:106
CbmKFParticleFinder::SetPrimaryProbCut
void SetPrimaryProbCut(float prob)
Definition: CbmKFParticleFinder.cxx:591
CbmKFParticleFinder::fMCTrackArrayEvent
TClonesArray * fMCTrackArrayEvent
Definition: CbmKFParticleFinder.h:107
CbmKFParticleFinder::Exec
virtual void Exec(Option_t *opt)
Definition: CbmKFParticleFinder.cxx:125
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmKFParticleFinder::SetChi2CutResonances
void SetChi2CutResonances(float cut)
Definition: CbmKFParticleFinder.cxx:631
CbmTrack::GetChiSq
Double_t GetChiSq() const
Definition: CbmTrack.h:58
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition: CbmMCTrack.h:75
CbmMCDataManager.h
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
CbmKFParticleFinder::fStsTrackBranchName
TString fStsTrackBranchName
Definition: CbmKFParticleFinder.h:101
CbmKFParticleFinder::SetLCut
void SetLCut(float cut)
Definition: CbmKFParticleFinder.cxx:607
CbmKFParticleFinder.h
CbmKFParticleFinder::InversedChi2Prob
double InversedChi2Prob(double p, int ndf) const
Definition: CbmKFParticleFinder.cxx:566
CbmKFParticleFinder::AddDecayToReconstructionList
void AddDecayToReconstructionList(int pdg)
Definition: CbmKFParticleFinder.cxx:667
CbmKFParticleFinder
Definition: CbmKFParticleFinder.h:26
CbmKFParticleFinderPID::GetPID
const std::vector< int > & GetPID() const
Definition: CbmKFParticleFinderPID.h:80
CbmKFParticleFinder::SetChi2TopoCutCharm2D
void SetChi2TopoCutCharm2D(float cut)
Definition: CbmKFParticleFinder.cxx:661
L1Field.h
CbmL1PFFitter
Definition: CbmL1PFFitter.h:31
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, std::vector< int > &pidHypo)
Definition: CbmL1PFFitter.cxx:81
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition: CbmMCDataArray.cxx:133
CbmKFParticleFinder::SetChiPrimaryCutCharm
void SetChiPrimaryCutCharm(float cut)
Definition: CbmKFParticleFinder.cxx:646
CbmKFParticleFinderPID.h
CbmKFParticleFinder::SetChi2Cut2D
void SetChi2Cut2D(float cut)
Definition: CbmKFParticleFinder.cxx:613
CbmKFParticleFinder::SetMaxDistanceBetweenParticlesCut
void SetMaxDistanceBetweenParticlesCut(float cut)
Definition: CbmKFParticleFinder.cxx:604
CbmKFVertex::GetRefX
Double_t & GetRefX()
Definition: CbmKFVertex.h:23
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
ClassImp
ClassImp(CbmKFParticleFinder)
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmKFParticleFinder::SetSuperEventAnalysis
void SetSuperEventAnalysis()
Definition: CbmKFParticleFinder.cxx:595
CbmMCDataArray.h
CbmKFParticleFinder::fSEChiPrim
std::vector< float > fSEChiPrim
Definition: CbmKFParticleFinder.h:125
CbmKFParticleFinder::SetLdLCutXiOmega
void SetLdLCutXiOmega(float cut)
Definition: CbmKFParticleFinder.cxx:619
CbmKFParticleFinder::Init
virtual InitStatus Init()
Definition: CbmKFParticleFinder.cxx:59
CbmVertex::GetX
Double_t GetX() const
Definition: CbmVertex.h:68
CbmStsTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmStsTrack.h:76
CbmKFParticleFinder::fCbmPrimVertex
CbmVertex * fCbmPrimVertex
Definition: CbmKFParticleFinder.h:109
CbmKFParticleFinder::SetPtCutCharm
void SetPtCutCharm(float cut)
Definition: CbmKFParticleFinder.cxx:643
CbmKFParticleFinder::SetChi2TopoCutXiOmega
void SetChi2TopoCutXiOmega(float cut)
Definition: CbmKFParticleFinder.cxx:622
CbmL1PFFitter::CalculateFieldRegionAtLastPoint
void CalculateFieldRegionAtLastPoint(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field)
Definition: CbmL1PFFitter.cxx:668
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition: CbmMCTrack.h:77
CbmKFParticleFinder::SetChi2CutCharm2D
void SetChi2CutCharm2D(float cut)
Definition: CbmKFParticleFinder.cxx:664
CbmEvent.h
CbmKFParticleFinder::SetPtCutJPsi
void SetPtCutJPsi(float cut)
Definition: CbmKFParticleFinder.cxx:640
sgn
friend F32vec4 sgn(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:65
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition: CbmMCEventList.h:90
finite
T finite(T x)
Definition: CbmL1Def.h:21
fvecLen
const int fvecLen
Definition: L1/vectors/P4_F32vec4.h:251
CbmKFParticleFinder::fEventList
CbmMCEventList * fEventList
Definition: CbmKFParticleFinder.h:108
CbmKFParticleFinder::fSEpdg
std::vector< int > fSEpdg
Definition: CbmKFParticleFinder.h:123
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition: CbmKFVertex.h:25
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition: CbmL1PFFitter.cxx:403
tracks
TClonesArray * tracks
Definition: Analyze_matching.h:17
CbmKFParticleFinder::fPVFindMode
int fPVFindMode
Definition: CbmKFParticleFinder.h:114
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmKFParticleFinder::GetKFParticleFinder
KFParticleFinder * GetKFParticleFinder()
Definition: CbmKFParticleFinder.cxx:601
CbmKFParticleFinder::fSETracks
std::vector< CbmStsTrack > fSETracks
Definition: CbmKFParticleFinder.h:121
CbmKFParticleFinder::SetPCutLMVM
void SetPCutLMVM(float cut)
Definition: CbmKFParticleFinder.cxx:637
CbmKFParticleFinder::SetChi2TopoCutCharmManybodyDecays
void SetChi2TopoCutCharmManybodyDecays(float cut)
Definition: CbmKFParticleFinder.cxx:652
CbmVertex
Definition: CbmVertex.h:26
CbmKFParticleFinder::fPID
CbmKFParticleFinderPID * fPID
Definition: CbmKFParticleFinder.h:117
CbmKFParticleFinder::SetChi2TopoCutResonances
void SetChi2TopoCutResonances(float cut)
Definition: CbmKFParticleFinder.cxx:628
CbmL1PFFitter.h
CbmKFParticleFinder::CbmKFParticleFinder
CbmKFParticleFinder(const char *name="CbmKFParticleFinder", Int_t iVerbose=0)
Definition: CbmKFParticleFinder.cxx:34
nMCTracks
Int_t nMCTracks
Definition: CbmHadronAnalysis.cxx:63
CbmTrack::GetNDF
Int_t GetNDF() const
Definition: CbmTrack.h:59
CbmVertex::GetZ
Double_t GetZ() const
Definition: CbmVertex.h:70
CbmKFParticleFinder::SetChi2CutCharmManybodyDecays
void SetChi2CutCharmManybodyDecays(float cut)
Definition: CbmKFParticleFinder.cxx:655
CbmMCEventList
Container class for MC events with number, file and start time.
Definition: CbmMCEventList.h:38
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmKFParticleFinder::fSuperEventAnalysis
bool fSuperEventAnalysis
Definition: CbmKFParticleFinder.h:120
CbmKFParticleFinder::fTopoReconstructor
KFParticleTopoReconstructor * fTopoReconstructor
Definition: CbmKFParticleFinder.h:112
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
CbmKFParticleFinder::FillKFPTrackVector
void FillKFPTrackVector(KFPTrackVector *tracks, const std::vector< CbmStsTrack > &vRTracks, const std::vector< KFFieldVector > &vField, const std::vector< int > &pdg, const std::vector< int > &trackId, const std::vector< float > &vChiToPrimVtx, bool atFirstPoint=1) const
Definition: CbmKFParticleFinder.cxx:462
CbmMCEventList.h
CbmKFParticleFinder::Finish
virtual void Finish()
Definition: CbmKFParticleFinder.cxx:431
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition: CbmMCTrack.h:76
CbmKFParticleFinder::SetLdLCutCharm2D
void SetLdLCutCharm2D(float cut)
Definition: CbmKFParticleFinder.cxx:658
CbmMCTrack.h
CbmKFParticleFinder::fTrackArray
TClonesArray * fTrackArray
Name of the input TCA with reco tracks.
Definition: CbmKFParticleFinder.h:104
CbmVertex::GetY
Double_t GetY() const
Definition: CbmVertex.h:69
CbmMCTrack
Definition: CbmMCTrack.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition: CbmEvent.h:30
CbmKFParticleFinder::SetPtCutLMVM
void SetPtCutLMVM(float cut)
Definition: CbmKFParticleFinder.cxx:634
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition: CbmKFVertex.h:24
CbmKFParticleFinder::SetLdLCut2D
void SetLdLCut2D(float cut)
Definition: CbmKFParticleFinder.cxx:616
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmKFParticleFinder::fSEField
std::vector< KFFieldVector > fSEField
Definition: CbmKFParticleFinder.h:122
CbmKFParticleFinder::SetLdLCutCharmManybodyDecays
void SetLdLCutCharmManybodyDecays(float cut)
Definition: CbmKFParticleFinder.cxx:649
CbmKFParticleFinder::fSETrackId
std::vector< int > fSETrackId
Definition: CbmKFParticleFinder.h:124
CbmKFParticleFinder::~CbmKFParticleFinder
~CbmKFParticleFinder()
Definition: CbmKFParticleFinder.cxx:55
CbmKFVertex.h
CbmKFParticleFinder::fEvents
TClonesArray * fEvents
Definition: CbmKFParticleFinder.h:105
CbmKFParticleFinder::fTimeSliceMode
bool fTimeSliceMode
Definition: CbmKFParticleFinder.h:127
CbmKFParticleFinder::SetChiPrimaryCut2D
void SetChiPrimaryCut2D(float cut)
Definition: CbmKFParticleFinder.cxx:610
CbmKFVertex
Definition: CbmKFVertex.h:6
CbmKFParticleFinder::SetChi2CutXiOmega
void SetChi2CutXiOmega(float cut)
Definition: CbmKFParticleFinder.cxx:625