Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFTopoPerformance.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFTopoPerformance.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 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
23 
24 #include "KFTopoPerformance.h"
25 
26 #ifdef KFPWITHTRACKER
27 #ifdef MAIN_DRAW
28 #include "AliHLTTPCCADisplay.h"
29 #endif //DRAW
30 #include "AliHLTTPCCATracker.h"
31 #include "AliHLTTPCCAPerformance.h"
32 #endif
33 
34 #include "KFParticleTopoReconstructor.h"
35 #include "KFParticleSIMD.h"
36 #include "KFPHistogram/KFPHistogram.h"
37 #include "TParticlePDG.h"
38 #include "TDatabasePDG.h"
39 
40 #include "TMath.h"
41 #include "TH1.h"
42 #include "TH2.h"
43 #include "TH3.h"
44 #include "TProfile.h"
45 #include "TProfile2D.h"
46 
47 #include <map>
48 #include <algorithm>
49 using std::sort;
50 using std::vector;
51 
52 KFTopoPerformance::KFTopoPerformance():KFParticlePerformanceBase(),fTopoReconstructor(0),fPrimVertices(0), fMCTrackToMCPVMatch(0),
53  fPVPurity(0), fNCorrectPVTracks(0), fTrackMatch(0), vMCTracks(0), vMCParticles(0), fNeutralIndex(0), MCtoRParticleId(0), RtoMCParticleId(0),
54  MCtoRPVId(0), RtoMCPVId(0), fPrintEffFrequency(1), fCentralityBin(-1), fCentralityWeight(0.f)
55 {
56 }
57 
58 KFTopoPerformance::~KFTopoPerformance()
59 {
60 }
61 
62 #ifdef KFPWITHTRACKER
63 void KFTopoPerformance::SetNewEvent(
64  const AliHLTTPCCAGBTracker * const tracker,
65  AliHLTResizableArray<AliHLTTPCCAHitLabel> *hitLabels,
66  AliHLTResizableArray<AliHLTTPCCAMCTrack> *mcTracks,
67  AliHLTResizableArray<AliHLTTPCCALocalMCPoint> *localMCPoints)
68 {
69  vMCTracks.resize(mcTracks->Size());
70  for(int iTr=0; iTr<vMCTracks.size(); iTr++)
71  {
72  for(int iP=0; iP<3; iP++)
73  vMCTracks[iTr].SetPar(iP, (*mcTracks)[iTr].Par(iP));
74  for(int iP=3; iP<6; iP++)
75  vMCTracks[iTr].SetPar(iP, (*mcTracks)[iTr].Par(iP) * (*mcTracks)[iTr].P());
76  vMCTracks[iTr].SetPar(6, (*mcTracks)[iTr].Par(6));
77 
78  vMCTracks[iTr].SetPDG( (*mcTracks)[iTr].PDG() );
79  vMCTracks[iTr].SetMotherId( (*mcTracks)[iTr].MotherId() );
80  vMCTracks[iTr].SetNMCPoints( (*mcTracks)[iTr].NMCPoints() );
81  }
82 } // void KFTopoPerformance::SetNewEvent
83 #endif
84 
85 void KFTopoPerformance::SetTopoReconstructor( const KFParticleTopoReconstructor * const TopoReconstructor)
86 {
88  fTopoReconstructor = TopoReconstructor;
89 } // void KFTopoPerformance::SetTopoReconstructor
90 
91 void KFTopoPerformance::CheckMCTracks()
92 {
94  fMCTrackToMCPVMatch.clear();
95  fPrimVertices.clear();
96 
97  // find prim vertex
98  if (vMCTracks.size() <= 0) return;
99 
100  fMCTrackToMCPVMatch.resize(vMCTracks.size(),-1);
101 
102  vector<int> pvIndex;
103  for(unsigned int iTr=0; iTr<vMCTracks.size(); iTr++)
104  {
105  KFMCTrack &mctr = vMCTracks[iTr];
106  int motherID = mctr.MotherId();
107  if(motherID <0 )
108  {
109  bool newPV = 1;
110  for(unsigned int iPV=0; iPV<pvIndex.size(); iPV++)
111  {
112  if(motherID == pvIndex[iPV])
113  {
114  fPrimVertices[iPV].AddDaughterTrack(iTr);
115  fMCTrackToMCPVMatch[iTr] = iPV;
116  newPV = 0;
117  }
118  }
119  if(newPV)
120  {
121  KFMCVertex primVertex;
122  primVertex.SetX( mctr.X() );
123  primVertex.SetY( mctr.Y() );
124  primVertex.SetZ( mctr.Z() );
125  primVertex.AddDaughterTrack(iTr);
126  if(motherID == -1)
127  primVertex.SetTriggerPV();
128  fPrimVertices.push_back(primVertex);
129  fMCTrackToMCPVMatch[iTr] = pvIndex.size();
130  pvIndex.push_back(motherID);
131  }
132  }
133  }
134 } // void KFTopoPerformance::CheckMCTracks()
135 
136 void KFTopoPerformance::GetMCParticles()
137 {
140  vMCParticles.clear();
141  vMCParticles.reserve(vMCTracks.size());
142  // all MC tracks are copied into KF MC Particles
143  for(unsigned int iMC=0; iMC < vMCTracks.size(); iMC++)
144  {
145  KFMCTrack &mtra = vMCTracks[iMC];
146  KFMCParticle part;
147  part.SetMCTrackID( iMC );
148  part.SetMotherId ( mtra.MotherId() );
149  part.SetPDG ( mtra.PDG() );
150  vMCParticles.push_back( part );
151  }
152  // find relations between mother and daughter MC particles
153  const int nMCParticles = vMCParticles.size();
154  for (int iP = 0; iP < nMCParticles; iP++ ) {
155  KFMCParticle &part = vMCParticles[iP];
156  int motherId = part.GetMotherId();
157  if(motherId < 0) continue;
158  if(motherId >= nMCParticles)
159  {
160  std::cout << "ERROR!!!!! KF Particle Performance: MC track mother Id is out of range." << std::endl;
161  exit(1);
162  }
163  KFMCParticle &motherPart = vMCParticles[motherId];
164  motherPart.AddDaughter(iP);
165  }
166 
167  fNeutralIndex.clear();
168  fNeutralIndex.resize(nMCParticles, -1);
169 
170  for(unsigned int iMC=0; iMC < vMCParticles.size(); iMC++)
171  {
172  KFMCParticle &part = vMCParticles[iMC];
173  part.SetMCTrackID( iMC );
174  }
175 
176  const int NmmPDG = 7;
177  vector<int> mmMotherPDG(NmmPDG); //PDG for particles found by the missing mass method
178  vector<int> mmChargedDaughterPDG(NmmPDG);
179  vector<int> mmNeutralDaughterPDG(NmmPDG);
180  vector<int> newMotherPDG(NmmPDG);
181  vector<int> newNeutralPDG(NmmPDG);
182 
183  mmMotherPDG[ 0] = 3112; mmChargedDaughterPDG[ 0] = 211; mmNeutralDaughterPDG[ 0] = 2112;
184  mmMotherPDG[ 1] = 3222; mmChargedDaughterPDG[ 1] = 211; mmNeutralDaughterPDG[ 1] = 2112;
185  mmMotherPDG[ 2] = 3312; mmChargedDaughterPDG[ 2] = 211; mmNeutralDaughterPDG[ 2] = 3122;
186  mmMotherPDG[ 3] = 3334; mmChargedDaughterPDG[ 3] = 211; mmNeutralDaughterPDG[ 3] = 3322;
187  mmMotherPDG[ 4] = 321; mmChargedDaughterPDG[ 4] = 211; mmNeutralDaughterPDG[ 4] = 111;
188  mmMotherPDG[ 5] = 3334; mmChargedDaughterPDG[ 5] = 321; mmNeutralDaughterPDG[ 5] = 3122;
189  mmMotherPDG[ 6] = 3222; mmChargedDaughterPDG[ 6] = 2212; mmNeutralDaughterPDG[ 6] = 111;
190 
191  newMotherPDG[ 0] = 7003112; newNeutralPDG[ 0] = 7002112;
192  newMotherPDG[ 1] = 7003222; newNeutralPDG[ 1] = 8002112;
193  newMotherPDG[ 2] = 7003312; newNeutralPDG[ 2] = 7003122;
194  newMotherPDG[ 3] = 7003334; newNeutralPDG[ 3] = 7003322;
195  newMotherPDG[ 4] = 9000321; newNeutralPDG[ 4] = 9000111;
196  newMotherPDG[ 5] = 8003334; newNeutralPDG[ 5] = 8003122;
197  newMotherPDG[ 6] = 8003222; newNeutralPDG[ 6] = 8000111;
198 
199  //add neutrinos, if they are not saved
200  for(int iMC=0; iMC<nMCParticles; iMC++)
201  {
202  if( abs(vMCParticles[iMC].GetPDG()) == 211 || abs(vMCParticles[iMC].GetPDG()) == 321 )
203  {
204  int muonIndex = -1;
205  for(int iD=0; iD<vMCParticles[iMC].NDaughters(); iD++)
206  if( abs(vMCParticles[vMCParticles[iMC].GetDaughterIds()[iD]].GetPDG()) == 13 )
207  muonIndex = vMCParticles[iMC].GetDaughterIds()[iD];
208 
209  if(muonIndex > -1)
210  {
211  int newPDG = 0;
212  if(vMCParticles[iMC].GetPDG() >0)
213  newPDG = vMCParticles[iMC].GetPDG() + 7000000;
214  else
215  newPDG = vMCParticles[iMC].GetPDG() - 7000000;
216  KFMCParticle motherPart = vMCParticles[iMC];
217  KFMCTrack motherTrack = vMCTracks[motherPart.GetMCTrackID()];
218  motherTrack.SetPDG(newPDG);
219  motherTrack.SetNotReconstructed();
220  int newMotherIndex = vMCTracks.size();
221  motherPart.SetPDG(newPDG);
222  motherPart.SetMCTrackID(newMotherIndex);
223 
224  const KFMCParticle& daughterPart = vMCParticles[muonIndex];
225  const KFMCTrack& daughterTrack = vMCTracks[daughterPart.GetMCTrackID()];
226 
227  int neutrinoPDG = 7000014;
228  if(vMCParticles[iMC].GetPDG() == -211) neutrinoPDG = -7000014;
229  if(vMCParticles[iMC].GetPDG() == 321) neutrinoPDG = 8000014;
230  if(vMCParticles[iMC].GetPDG() == -321) neutrinoPDG = -8000014;
231 
232  int neutrinoIndex = vMCTracks.size()+1;
233  vMCParticles[iMC].AddDaughter(neutrinoIndex);
234 
235  fNeutralIndex[iMC] = neutrinoIndex;
236 
237  KFMCParticle neutrinoPart;
238  KFMCTrack neutrinoTrack;
239  neutrinoTrack.SetX(daughterTrack.X());
240  neutrinoTrack.SetY(daughterTrack.Y());
241  neutrinoTrack.SetZ(daughterTrack.Z());
242  neutrinoTrack.SetPx(motherTrack.Px() - daughterTrack.Px());
243  neutrinoTrack.SetPy(motherTrack.Py() - daughterTrack.Py());
244  neutrinoTrack.SetPz(motherTrack.Pz() - daughterTrack.Pz());
245  neutrinoTrack.SetQP(0);
246  neutrinoTrack.SetMotherId(newMotherIndex);
247  neutrinoTrack.SetPDG(neutrinoPDG);
248 
249  motherPart.CleanDaughters();
250  motherPart.AddDaughter(muonIndex);
251  motherPart.AddDaughter(neutrinoIndex);
252  motherPart.SetInitialParticleId(iMC);
253  vMCTracks.push_back(motherTrack);
254  vMCParticles.push_back(motherPart);
255  fNeutralIndex.push_back(-1);
256 
257  neutrinoPart.SetMCTrackID(neutrinoIndex);
258  neutrinoPart.SetMotherId(newMotherIndex);
259  neutrinoPart.SetPDG(neutrinoPDG);
260  neutrinoPart.AddDaughter(iMC);
261  neutrinoPart.AddDaughter(muonIndex);
262  vMCTracks.push_back(neutrinoTrack);
263  vMCParticles.push_back(neutrinoPart);
264  fNeutralIndex.push_back(-1);
265 
266  vMCParticles[iMC].SetAsReconstructable(4);
267  vMCParticles[muonIndex].SetAsReconstructable(4);
268  }
269  }
270  // add sigmas, omegas, xis ...
271  if( vMCParticles[iMC].NDaughters() >= 2 &&
272  (abs(vMCParticles[iMC].GetPDG()) == 3112 ||
273  abs(vMCParticles[iMC].GetPDG()) == 3222 ||
274  abs(vMCParticles[iMC].GetPDG()) == 3312 ||
275  abs(vMCParticles[iMC].GetPDG()) == 3334 ||
276  abs(vMCParticles[iMC].GetPDG()) == 321) )
277  {
278  int neutralDaughterId = -1, chargedDaughterId = -1;
279 
280  int newPDG = 0;
281  int neutralPDG = 0;
282 
283  for(int iPDG=0; iPDG<NmmPDG; iPDG++)
284  {
285  if(abs(vMCParticles[iMC].GetPDG()) == mmMotherPDG[iPDG])
286  {
287  bool isDaughter[2] = {0,0};
288 
289  vector<float> xDaughter;
290  vector<float> yDaughter;
291  vector<float> zDaughter;
292  vector< vector<int> > nDaughtersAtPoint;
293 
294  for(int iMCDaughter=0; iMCDaughter<vMCParticles[iMC].NDaughters(); iMCDaughter++)
295  {
296  bool isNewDecayPoint = 1;
297 
298  for(unsigned int iPoint=0; iPoint<xDaughter.size(); iPoint++)
299  {
300  float dx = fabs(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].X() - xDaughter[iPoint]);
301  float dy = fabs(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].Y() - yDaughter[iPoint]);
302  float dz = fabs(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].Z() - zDaughter[iPoint]);
303 
304  bool isSamePoint = (dx < 1.e-5 && dy < 1.e-5 && dz < 1.e-5);
305  if(isSamePoint)
306  nDaughtersAtPoint[iPoint].push_back(iMCDaughter);
307 
308  isNewDecayPoint &= !isSamePoint;
309  }
310 
311  if(isNewDecayPoint)
312  {
313  xDaughter.push_back(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].X());
314  yDaughter.push_back(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].Y());
315  zDaughter.push_back(vMCTracks[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].Z());
316  vector<int> newPointIndex;
317  newPointIndex.push_back(iMCDaughter);
318  nDaughtersAtPoint.push_back(newPointIndex);
319  }
320  }
321 
322  for(unsigned int iPoint = 0; iPoint<nDaughtersAtPoint.size(); iPoint++)
323  {
324  if(nDaughtersAtPoint[iPoint].size() == 2)
325  {
326  for(unsigned int iDaughter=0; iDaughter<nDaughtersAtPoint[iPoint].size(); iDaughter++)
327  {
328  int iMCDaughter = nDaughtersAtPoint[iPoint][iDaughter];
329  if(abs(vMCParticles[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].GetPDG()) == mmChargedDaughterPDG[iPDG])
330  {
331  isDaughter[0] = 1;
332  chargedDaughterId = vMCParticles[iMC].GetDaughterIds()[iMCDaughter];
333  }
334  if(abs(vMCParticles[vMCParticles[iMC].GetDaughterIds()[iMCDaughter]].GetPDG()) == mmNeutralDaughterPDG[iPDG])
335  {
336  isDaughter[1] = 1;
337  neutralDaughterId = vMCParticles[iMC].GetDaughterIds()[iMCDaughter];
338  }
339  }
340 
341  if(isDaughter[0] && isDaughter[1])
342  {
343  int signPDG = vMCParticles[iMC].GetPDG()/abs(vMCParticles[iMC].GetPDG());
344  newPDG = signPDG * newMotherPDG[iPDG];
345  neutralPDG = signPDG * newNeutralPDG[iPDG];
346  }
347  }
348  }
349  }
350  }
351 
352  if(newPDG != 0)
353  {
354  KFMCParticle motherPart = vMCParticles[iMC];
355  KFMCTrack motherTrack = vMCTracks[motherPart.GetMCTrackID()];
356  motherTrack.SetPDG(newPDG);
357  motherTrack.SetNotReconstructed();
358  int newMotherIndex = vMCTracks.size();
359  motherPart.SetPDG(newPDG);
360  motherPart.SetMCTrackID(newMotherIndex);
361 
362  int neutrinoIndex = vMCTracks.size()+1;
363  fNeutralIndex[iMC] = neutrinoIndex;
364 
365  KFMCTrack neutralTrack = vMCTracks[neutralDaughterId];
366  neutralTrack.SetMotherId(newMotherIndex);
367  neutralTrack.SetPDG(neutralPDG);
368 
369  motherPart.CleanDaughters();
370  motherPart.AddDaughter(chargedDaughterId);
371  motherPart.AddDaughter(neutrinoIndex);
372  motherPart.SetInitialParticleId(iMC);
373  vMCTracks.push_back(motherTrack);
374  vMCParticles.push_back(motherPart);
375  fNeutralIndex.push_back(-1);
376 
377  KFMCParticle neutralPart;
378  neutralPart.SetMCTrackID(neutrinoIndex);
379  neutralPart.SetMotherId(newMotherIndex);
380  neutralPart.SetPDG(neutralPDG);
381  neutralPart.AddDaughter(iMC);
382  neutralPart.AddDaughter(chargedDaughterId);
383  vMCTracks.push_back(neutralTrack);
384  vMCParticles.push_back(neutralPart);
385  fNeutralIndex.push_back(-1);
386 
387  vMCParticles[iMC].SetAsReconstructable(4);
388  vMCParticles[chargedDaughterId].SetAsReconstructable(4);
389  }
390  }
391  }
392 
393  //clean Lambda c daughters
394  for(unsigned int iMC=0; iMC < vMCParticles.size(); iMC++)
395  {
396  KFMCParticle &part = vMCParticles[iMC];
397 
398 // if(abs(part.GetPDG()) == 4122)
399  {
400  //add daughters into one pool
401  vector<int> newDaughters;
402  for(unsigned int iD=0; iD<part.GetDaughterIds().size(); iD++)
403  {
404  KFMCParticle &d = vMCParticles[part.GetDaughterIds()[iD]];
405  if( abs(d.GetPDG())==3224 ||
406  abs(d.GetPDG())==3114 ||
407  abs(d.GetPDG())==113 ||
408  abs(d.GetPDG())==313 ||
409  abs(d.GetPDG())==323 ||
410  abs(d.GetPDG())==2224 ||
411  abs(d.GetPDG())==2214 ||
412  abs(d.GetPDG())==2114 ||
413  abs(d.GetPDG())==1114
414  )
415  {
416  for(unsigned int iDaughter=0; iDaughter<d.GetDaughterIds().size(); iDaughter++)
417  {
418  newDaughters.push_back(d.GetDaughterIds()[iDaughter]);
419  vMCParticles[d.GetDaughterIds()[iDaughter]].SetMotherId(iMC);
420  }
421  }
422  else
423 // if(d.GetDaughterIds().size() == 0 || abs(d.GetPDG())==310 || abs(d.GetPDG())==3122)
424  newDaughters.push_back(part.GetDaughterIds()[iD]);
425  }
426  part.CleanDaughters();
427  for(unsigned int iDaughter=0; iDaughter<newDaughters.size(); iDaughter++)
428  part.AddDaughter(newDaughters[iDaughter]);
429 
430  //change PDG to separate channels
431  int indexPDG = fParteff.GetParticleIndex(part.GetPDG());
432  for(int iPDG=indexPDG; iPDG<indexPDG+10; iPDG++)
433  {
434  const int nDaughters = part.GetDaughterIds().size();
435 
436  if(int(fParteff.partDaughterPdg[iPDG].size()) != nDaughters)
437  continue;
438 
439  vector<bool> isDaughterFound(nDaughters);
440  vector<bool> isDaughterUsed(nDaughters);
441  for(int iDMC=0; iDMC<nDaughters; iDMC++)
442  {
443  isDaughterFound[iDMC] = 0;
444  isDaughterUsed[iDMC] = 0;
445  }
446 
447  bool isCorrectPDG = 1;
448  for(int iDMC=0; iDMC<nDaughters; iDMC++)
449  for(int iD=0; iD<nDaughters; iD++)
450  {
451  if(isDaughterUsed[iD]) continue;
452  if(vMCParticles[part.GetDaughterIds()[iD]].GetPDG() == fParteff.partDaughterPdg[iPDG][iDMC])
453  {
454  isDaughterUsed[iD] = 1;
455  isDaughterFound[iDMC] = 1;
456  break;
457  }
458  }
459 
460  for(int iDMC=0; iDMC<nDaughters; iDMC++)
461  isCorrectPDG &= isDaughterFound[iDMC];
462 
463  if(isCorrectPDG)
464  {
465  part.SetPDG(fParteff.partPDG[iPDG]);
466  break;
467  }
468  }
469  }
470  }
471 }
472 
473 void KFTopoPerformance::FindReconstructableMCParticles()
474 {
476  const unsigned int nMCParticles = vMCParticles.size();
477 
478  for ( unsigned int iP = 0; iP < nMCParticles; iP++ ) {
479  KFMCParticle &part = vMCParticles[iP];
480  CheckMCParticleIsReconstructable(part);
481  }
482 }
483 
484 void KFTopoPerformance::CheckMCParticleIsReconstructable(KFMCParticle &part)
485 {
487  if ( part.IsReconstructable(0) ) return;
488  if ( vMCTracks[part.GetMCTrackID()].IsOutOfDetector() ) return;
489 
490  if( abs(part.GetPDG()) == 211 ||
491  abs(part.GetPDG()) == 2212 ||
492  abs(part.GetPDG()) == 321 ||
493  abs(part.GetPDG()) == 13 ||
494  abs(part.GetPDG()) == 3112 ||
495  abs(part.GetPDG()) == 3222 ||
496  abs(part.GetPDG()) == 3312 ||
497  abs(part.GetPDG()) == 3334 )
498  {
499  int iMCTrack = part.GetMCTrackID();
500  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
501 
502  // reconstructable in 4pi is defined in GetMCParticles, when decay is found
503  if(mcTrack.IsReconstructed())
504  part.SetAsReconstructable(3);
505  }
506  // tracks
507  if ( abs(part.GetPDG()) == 211 ||
508  abs(part.GetPDG()) == 2212 ||
509  abs(part.GetPDG()) == 321 ||
510  abs(part.GetPDG()) == 11 ||
511  abs(part.GetPDG()) == 13 ||
512  abs(part.GetPDG()) == 1000010020 ||
513  abs(part.GetPDG()) == 1000010030 ||
514  abs(part.GetPDG()) == 1000020030 ||
515  abs(part.GetPDG()) == 1000020040 ||
516  ( (part.GetPDG() == 22) && (vMCTracks[part.GetMCTrackID()].IsReconstructed()) ) )
517  {
518  int iMCTrack = part.GetMCTrackID();
519  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
520 
521  part.SetAsReconstructable(0);
522 
523  if(mcTrack.NMCPoints() >= 15)
524  part.SetAsReconstructable(1);
525 // if(mc.IsReconstructable())
526 // part.SetAsReconstructable2();
527 
528  if(mcTrack.IsReconstructed())
529  part.SetAsReconstructable(2);
530  }
531  // mother particles
532  else
533  {
534  //Check if the particle is V0
535 
536  if(part.NDaughters() >= 2)
537  {
538  bool isPositiveDaughter[5] = {0,0,0,0,0};
539  bool isNegativeDaughter[5] = {0,0,0,0,0};
540 
541  int nRecoDaughters[5] = {0,0,0,0,0};
542 
543  for(int iD=0; iD < part.NDaughters(); iD++)
544  {
545  KFMCParticle &daughter = vMCParticles[part.GetDaughterIds()[iD]];
546  CheckMCParticleIsReconstructable(daughter);
547 
548  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(daughter.GetPDG());
549  Double_t charge = (particlePDG) ? particlePDG->Charge()/3 : 0;
550 
551  for(int iEff=0; iEff<5; iEff++)
552  {
553  if(charge > 0)
554  isPositiveDaughter[iEff] |= daughter.IsReconstructable(iEff);
555  else
556  isNegativeDaughter[iEff] |= daughter.IsReconstructable(iEff);
557 
558  if(daughter.IsReconstructable(iEff))
559  nRecoDaughters[iEff]++;
560  }
561  }
562 
563 // for(int iEff=0; iEff<3; iEff++)
564 // if(isPositiveDaughter[iEff] && isNegativeDaughter[iEff])
565 // part.SetAsReconstructableV0(iEff);
566  for(int iEff=0; iEff<3; iEff++)
567  if(nRecoDaughters[iEff] > 1)
568  part.SetAsReconstructableV0(iEff);
569  }
570 
571  for(int iPart=0; iPart<fParteff.nParticles; iPart++)
572  {
573  if(part.GetPDG() == fParteff.partPDG[iPart])
574  {
575  const unsigned int nDaughters = fParteff.partDaughterPdg[iPart].size();
576  if( part.GetDaughterIds().size() != nDaughters ) return;
577  vector<int> pdg(nDaughters);
578 
579  for(unsigned int iD=0; iD<nDaughters; iD++)
580  pdg[iD] = vMCParticles[part.GetDaughterIds()[iD]].GetPDG();
581 
582  vector<bool> isDaughterFound(nDaughters);
583  for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
584  isDaughterFound[iDMC] = 0;
585 
586  bool isReco = 1;
587  for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
588  for(unsigned int iD=0; iD<nDaughters; iD++)
589  if(pdg[iD] == fParteff.partDaughterPdg[iPart][iDMC]) isDaughterFound[iDMC] = 1;
590 
591  for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
592  isReco = isReco && isDaughterFound[iDMC];
593 
594  if(!isReco) return;
595  }
596  }
597 
598  const vector<int>& dIds = part.GetDaughterIds();
599  const unsigned int nD = dIds.size();
600  if(nD == 0) return; //TODO optimize for all species
601 
602  bool reco1 = 1;
603  bool reco2 = 1;
604  bool reco3 = 1;
605  bool reco4 = 1;
606  bool reco5 = 1;
607 
608  for ( unsigned int iD = 0; iD < nD && (reco1 || reco2 || reco3 || reco4 || reco5); iD++ ) {
609  KFMCParticle &dp = vMCParticles[dIds[iD]];
610  CheckMCParticleIsReconstructable(dp);
611  reco1 &= dp.IsReconstructable(0);
612  reco2 &= dp.IsReconstructable(1);
613  reco3 &= dp.IsReconstructable(2);
614  reco4 &= dp.IsReconstructable(3);
615  reco5 &= dp.IsReconstructable(4);
616  }
617 
618  if (reco1) part.SetAsReconstructable(0);
619  if (reco2) part.SetAsReconstructable(1);
620  if (reco3) part.SetAsReconstructable(2);
621  int iParticle = fParteff.GetParticleIndex(part.GetPDG());
622  if (reco4 && iParticle>=KFPartEfficiencies::fFirstMissingMassParticleIndex &&
623  iParticle<=KFPartEfficiencies::fLastMissingMassParticleIndex ) part.SetAsReconstructable(3);
624  if (reco5 && iParticle>=KFPartEfficiencies::fFirstMissingMassParticleIndex &&
625  iParticle<=KFPartEfficiencies::fLastMissingMassParticleIndex ) part.SetAsReconstructable(4);
626  }
627 }
628 
629 
630 void KFTopoPerformance::FindReconstructableMCVertices()
631 {
633  const unsigned int nMCVertices = fPrimVertices.size();
634 
635  for ( unsigned int iV = 0; iV < nMCVertices; iV++ ) {
636  KFMCVertex &vert = fPrimVertices[iV];
637 
638  int nReconstructableDaughters = 0;
639  int nMCReconstructableDaughters = 0;
640 
641  for(int iP=0; iP<vert.NDaughterTracks(); iP++)
642  {
643  int idDaughter = vert.DaughterTrack(iP);
644  KFMCParticle &part = vMCParticles[idDaughter];
645 
646  if(part.IsReconstructable(2)) nReconstructableDaughters++;
647  if(part.IsReconstructable(1)) nMCReconstructableDaughters++;
648  }
649 
650  if(nReconstructableDaughters >= 2) vert.SetReconstructable();
651  else vert.SetUnReconstructable();
652 
653  if(nMCReconstructableDaughters >= 2) vert.SetMCReconstructable();
654  else vert.SetMCUnReconstructable();
655  }
656 }
657 
658 void KFTopoPerformance::MatchParticles()
659 {
661  MCtoRParticleId.clear();
662  RtoMCParticleId.clear();
663  MCtoRParticleId.resize(vMCParticles.size());
664  RtoMCParticleId.resize(fTopoReconstructor->GetParticles().size() );
665 
666  // match tracks ( particles which are direct copies of tracks )
667  for( unsigned int iRP = 0; iRP < fTopoReconstructor->GetParticles().size(); iRP++ )
668  {
669  const KFParticle &rPart = fTopoReconstructor->GetParticles()[iRP];
670 
671  if (rPart.NDaughters() != 1) continue;
672 
673  const int rTrackId = rPart.DaughterIds()[0];
674  const int mcTrackId = fTrackMatch[rTrackId];
675 
676  if(mcTrackId < 0) continue;
677 
678  KFMCParticle &mPart = vMCParticles[mcTrackId];
679  if( mPart.GetPDG() == rPart.GetPDG() )
680  {
681  MCtoRParticleId[mcTrackId].ids.push_back(iRP);
682  RtoMCParticleId[iRP].ids.push_back(mcTrackId);
683  }
684  else {
685  MCtoRParticleId[mcTrackId].idsMI.push_back(iRP);
686  RtoMCParticleId[iRP].idsMI.push_back(mcTrackId);
687  }
688  }
689 
690  // match created mother particles
691  for( unsigned int iRP = 0; iRP < fTopoReconstructor->GetParticles().size(); iRP++ ) {
692  const KFParticle &rPart = fTopoReconstructor->GetParticles()[iRP];
693  const unsigned int NRDaughters = rPart.NDaughters();
694 
695  if (NRDaughters < 2) continue;
696 
697  bool isMissingMass = ((abs(rPart.GetPDG()) == 7000211)||(abs(rPart.GetPDG()) == 7000321)||(abs(rPart.GetPDG()) == 7003112) || (abs(rPart.GetPDG()) == 7003222)|| (abs(rPart.GetPDG()) == 7003312)|| (abs(rPart.GetPDG()) == 7003334)|| (abs(rPart.GetPDG()) == 9000321)|| (abs(rPart.GetPDG()) == 8003334)||(abs(rPart.GetPDG()) == 8003222));
698 
699  //missing mass method
700  if ( (abs(rPart.GetPDG()) == 7000014 ) || (abs(rPart.GetPDG()) == 8000014 ) || (abs(rPart.GetPDG()) == 7002112 ) ||
701  (abs(rPart.GetPDG()) == 8002112 ) || (abs(rPart.GetPDG()) == 7003122 ) || (abs(rPart.GetPDG()) == 7003322 ) ||
702  (abs(rPart.GetPDG()) == 9000111 ) || (abs(rPart.GetPDG()) == 8003122 ) || (abs(rPart.GetPDG()) == 8000111 ) )
703  {
704  //During the reconstruction 1st daughter - mother particle, 2nd daughter - charged daughter
705 
706  int mcNeutralDaughterId = -1;
707  const int recoMotherId = rPart.DaughterIds()[0];
708  if ( !RtoMCParticleId[recoMotherId].IsMatched() ) continue;
709 
710  const int mcMotherId = RtoMCParticleId[recoMotherId].GetBestMatch();
711 
712  const int recoChargedDaughterId = rPart.DaughterIds()[1];
713  if ( !RtoMCParticleId[recoChargedDaughterId].IsMatched() ) continue;
714 
715  const int mcChargedDaughterId = RtoMCParticleId[recoChargedDaughterId].GetBestMatch();
716  const KFMCParticle& chargedDaughter = vMCParticles[mcChargedDaughterId];
717 
718  if(chargedDaughter.GetMotherId() != mcMotherId) continue;
719  const KFMCParticle& mother = vMCParticles[mcMotherId];
720 
721  if(fNeutralIndex[mcMotherId] > -1)
722  mcNeutralDaughterId = fNeutralIndex[mcMotherId];
723 
724  if(mcNeutralDaughterId > -1)
725  {
726  KFMCParticle &neutralDaughter = vMCParticles[mcNeutralDaughterId];
727 
728  int iParticle = fParteff.GetParticleIndex(rPart.GetPDG());
729 
730  bool allCorrectDaughters = mother.GetPDG() == fParteff.partDaughterPdg[iParticle][0] &&
731  chargedDaughter.GetPDG() == fParteff.partDaughterPdg[iParticle][1];
732 
733  if( neutralDaughter.GetPDG() == rPart.GetPDG() &&
734  neutralDaughter.NDaughters() == rPart.NDaughters() &&
735  allCorrectDaughters) {
736  MCtoRParticleId[mcNeutralDaughterId].ids.push_back(iRP);
737  RtoMCParticleId[iRP].ids.push_back(mcNeutralDaughterId);
738  }
739  else {
740  MCtoRParticleId[mcNeutralDaughterId].idsMI.push_back(iRP);
741  RtoMCParticleId[iRP].idsMI.push_back(mcNeutralDaughterId);
742  }
743  }
744  }
745 
746 
747 
748  //normal decays
749  else
750  {
751  unsigned int iD = 0;
752  vector<int> mcDaughterIds;
753  int mmId = -2; // MC id for rPart
754  {
755  const int rdId = rPart.DaughterIds()[iD];
756  if ( !RtoMCParticleId[rdId].IsMatched() ) continue;
757  const int mdId = RtoMCParticleId[rdId].GetBestMatch();
758  mcDaughterIds.push_back(mdId);
759  mmId = vMCParticles[mdId].GetMotherId();
760  }
761 
762  iD++;
763  for ( ; iD < NRDaughters; iD++ ) {
764  const int rdId = rPart.DaughterIds()[iD];
765  if ( !RtoMCParticleId[rdId].IsMatched() ) break;
766  const int mdId = RtoMCParticleId[rdId].GetBestMatch();
767  mcDaughterIds.push_back(mdId);
768 
769  if(isMissingMass)
770  {
771  const KFMCParticle &neutralDaughter = vMCParticles[mdId];
772  if(mmId != vMCParticles[neutralDaughter.GetMotherId()].InitialParticleId()) break;
773  mmId = neutralDaughter.GetMotherId();
774  }
775 
776  if( !(isMissingMass) && (vMCParticles[mdId].GetMotherId() != mmId) ) break;
777  }
778 
779  int nClones = 0;
780  sort(mcDaughterIds.begin(), mcDaughterIds.end());
781  for(unsigned int ie=1; ie<mcDaughterIds.size(); ie++)
782  {
783  if(mcDaughterIds[ie] == mcDaughterIds[ie-1])
784  nClones++;
785  }
786 
787  if(nClones > 0) continue;
788 
789  if ( iD == NRDaughters && mmId > -1 ) { // match is found and it is not primary vertex
790  KFMCParticle &mmPart = vMCParticles[mmId];
791 
792  if( mmPart.GetPDG() == rPart.GetPDG() &&
793  mmPart.NDaughters() == rPart.NDaughters() ) {
794  MCtoRParticleId[mmId].ids.push_back(iRP);
795  RtoMCParticleId[iRP].ids.push_back(mmId);
796  }
797  else {
798  MCtoRParticleId[mmId].idsMI.push_back(iRP);
799  RtoMCParticleId[iRP].idsMI.push_back(mmId);
800  }
801  }
802  }
803  }
804 }
805 
806 void KFTopoPerformance::MatchPV()
807 {
809  MCtoRPVId.clear();
810  RtoMCPVId.clear();
811  MCtoRPVId.resize(fPrimVertices.size());
812  RtoMCPVId.resize(fTopoReconstructor->NPrimaryVertices() );
813 
814  fPVPurity.clear();
815  fPVPurity.resize(fTopoReconstructor->NPrimaryVertices(), 0.);
816  fNCorrectPVTracks.clear();
817  fNCorrectPVTracks.resize(fTopoReconstructor->NPrimaryVertices(), 0);
818 
819  for(int iE = 0; iE<4; iE++)
820  {
821  fPVTracksRate[iE].clear();
822  fPVTracksRate[iE].resize(fTopoReconstructor->NPrimaryVertices(),0);
823  }
824 
825  for( unsigned int iMCPV=0; iMCPV<fPrimVertices.size(); iMCPV++)
826  {
827  KFMCVertex &mcPV = fPrimVertices[iMCPV];
828  int nReconstructedDaughters = 0;
829  for(int iD=0; iD<mcPV.NDaughterTracks(); iD++)
830  if(MCtoRParticleId[ mcPV.DaughterTrack(iD) ].IsMatched())
831  nReconstructedDaughters++;
832 
833  mcPV.SetNReconstructedDaughters(nReconstructedDaughters);
834  }
835 
836  for( int iPV = 0; iPV < fTopoReconstructor->NPrimaryVertices(); iPV++ ) {
837 
838  vector<int> &tracks = fTopoReconstructor->GetPVTrackIndexArray(iPV);
839 
840  int nPVTracks = tracks.size()>0 ? tracks.size() : -1;//tracks.size();
841 
842  vector<short int> nTracksFromMCPV(fPrimVertices.size() + vMCParticles.size(), 0);
843  vector< vector<int> > mcIDs(fPrimVertices.size() + vMCParticles.size());
844 
845  int nGhostTracks = 0;
846  int nTriggerTracks = 0;
847  int nPileupTracks = 0;
848  int nBGTracks = 0;
849 
850  for(int iRP=0; iRP < nPVTracks; iRP++)
851  {
852  const int rTrackId = tracks[iRP];
853  if ( !RtoMCParticleId[rTrackId].IsMatched() )
854  {
855  nGhostTracks++;
856  continue;
857  }
858 
859  int iMCPart = RtoMCParticleId[rTrackId].GetBestMatch();
860  KFMCParticle &mcPart = vMCParticles[iMCPart];
861  int iMCTrack = mcPart.GetMCTrackID();
862  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
863 
864  int motherId = mcTrack.MotherId();
865 
866  if(motherId < 0) //real PV
867  {
868  if(motherId == -1)
869  nTriggerTracks++;
870  if(motherId < -1)
871  nPileupTracks++;
872 
873  int iMCPV = fMCTrackToMCPVMatch[iMCTrack];
874  if(iMCPV < 0)
875  {
876  std::cout << "Error!!! iMCPV < 0" << std::endl;
877  continue;
878  }
879 
880  nTracksFromMCPV[iMCPV]++;
881  mcIDs[iMCPV].push_back(iMCTrack);
882  }
883  else // pchysics background
884  {
885  if(motherId >-1)
886  {
887  nBGTracks++;
888 
889  int iMCPV = motherId + fPrimVertices.size();
890 
891  nTracksFromMCPV[iMCPV]++;
892  mcIDs[iMCPV].push_back(iMCTrack);
893  }
894  }
895  }
896 
897 
898 
899  for(unsigned int iMCPV=0; iMCPV<nTracksFromMCPV.size(); iMCPV++)
900  {
901  int nClones = 0;
902  sort(mcIDs[iMCPV].begin(), mcIDs[iMCPV].end());
903  for(unsigned int ie=1; ie<mcIDs[iMCPV].size(); ie++)
904  {
905  if(mcIDs[iMCPV][ie] == mcIDs[iMCPV][ie-1])
906  nClones++;
907  }
908  nTracksFromMCPV[iMCPV] = nTracksFromMCPV[iMCPV] - nClones;
909  nPVTracks -= nClones;
910  }
911 
912  // calculate rate of each type of tracks in reconstructed PV
913  fPVTracksRate[0][iPV] = double(nGhostTracks)/double(nPVTracks);
914  fPVTracksRate[1][iPV] = double(nTriggerTracks)/double(nPVTracks);
915  fPVTracksRate[2][iPV] = double(nPileupTracks)/double(nPVTracks);
916  fPVTracksRate[3][iPV] = double(nBGTracks)/double(nPVTracks);
917 
918  int iBestMCPV=-1;
919  int nTracksBestMCPV = 1;
920  for(unsigned int iMCPV=0; iMCPV<nTracksFromMCPV.size(); iMCPV++ )
921  {
922  if(nTracksFromMCPV[iMCPV] > nTracksBestMCPV )
923  {
924  nTracksBestMCPV = nTracksFromMCPV[iMCPV];
925  iBestMCPV = iMCPV;
926  }
927  }
928 
929  if( (iBestMCPV > -1) && (iBestMCPV<int(fPrimVertices.size())) )
930  {
931  fNCorrectPVTracks[iPV] = nTracksFromMCPV[iBestMCPV];
932  }
933  else
934  fNCorrectPVTracks[iPV] = 0;
935 
936  double purity = double(nTracksBestMCPV)/double(nPVTracks);
937 
938  fPVPurity[iPV] = purity;
939 // if(purity < 0.7) continue;
940 
941  if(iBestMCPV < 0) continue;
942 
943  if(iBestMCPV < int(fPrimVertices.size()))
944  {
945  fPrimVertices[iBestMCPV].SetReconstructed();
946 
947  MCtoRPVId[iBestMCPV].ids.push_back(iPV);
948  RtoMCPVId[iPV].ids.push_back(iBestMCPV);
949  }
950  else
951  {
952  RtoMCPVId[iPV].idsMI.push_back(iBestMCPV);
953  }
954  }
955 }
956 
957 void KFTopoPerformance::MatchTracks()
958 {
960 #ifdef KFPWITHTRACKER
961  for(int iTr=0; iTr<vMCTracks.size(); iTr++)
962  {
963  if( (AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance")->GetMCData())[iTr].IsReconstructed() )
964  vMCTracks[iTr].SetReconstructed();
965  }
966  fTrackMatch.resize(AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance")->GetRecoData().size());
967  for(int iTr=0; iTr<fTrackMatch.size(); iTr++)
968  {
969  const AliHLTTPCCAPerformanceRecoTrackData& matchInfo = (AliHLTTPCCAPerformance::Instance().GetSubPerformance("Global Performance")->GetRecoData())[iTr];
970  if( matchInfo.IsGhost(PParameters::MinTrackPurity) || matchInfo.GetMCTrackId()<0 )
971  fTrackMatch[iTr] = -1;
972  else
973  fTrackMatch[iTr] = matchInfo.GetMCTrackId();
974  }
975 #endif
976 
977  GetMCParticles();
978  FindReconstructableMCParticles();
979  MatchParticles();
980  CalculateEfficiency();
981 
982  FindReconstructableMCVertices();
983  MatchPV();
984  CalculatePVEfficiency();
985 } // void KFTopoPerformance::MatchTracks()
986 
987 void KFTopoPerformance::CalculateEfficiency()
988 {
990  KFPartEfficiencies partEff; // efficiencies for current event
991 
992  const int NRP = fTopoReconstructor->GetParticles().size();
993  for ( int iP = 0; iP < NRP; ++iP ) {
994 // const CbmKFParticle &part = fPF->GetParticles()[iP];
995  const KFParticle &part = fTopoReconstructor->GetParticles()[iP];
996  const int pdg = part.GetPDG();
997 
998  const bool isBG = RtoMCParticleId[iP].idsMI.size() != 0;
999  const bool isGhost = !RtoMCParticleId[iP].IsMatched();
1000 
1001  for(int iPart=0; iPart<fParteff.nParticles; iPart++)
1002  if ( pdg == fParteff.partPDG[iPart] )
1003  partEff.IncReco(isGhost, isBG, fParteff.partName[iPart].data());
1004 
1005  // Calculate the gost level for V0
1006  if(abs(pdg) == 310 /*||
1007  CAMath::Abs(pdg) == 3122 ||
1008  CAMath::Abs(pdg) == 421 ||
1009  CAMath::Abs(pdg) == 22 */)
1010  {
1011  partEff.IncReco(isGhost, 0, fParteff.partName[fParteff.nParticles - 1].data());
1012  }
1013  }
1014 
1015  const int NMP = vMCParticles.size();
1016  for ( int iP = 0; iP < NMP; ++iP ) {
1017  const KFMCParticle &part = vMCParticles[iP];
1018  const int pdg = part.GetPDG();
1019  const int mId = part.GetMotherId();
1020 
1021  vector<bool> isReco;
1022  vector<int> nClones;
1023 
1024  vector<int> iParticle;
1025  iParticle.push_back(fParteff.GetParticleIndex(pdg));
1026  vector< vector<bool> > isReconstructable;
1027  vector<bool> isRecPart;
1028 
1029  const std::map<int,bool>& decays = fTopoReconstructor->GetKFParticleFinder()->GetReconstructionList();
1030  if(!(decays.empty()) && (iParticle[0] < fParteff.fFirstStableParticleIndex || iParticle[0] > fParteff.fLastStableParticleIndex))
1031  if(decays.find(pdg) == decays.end()) continue;
1032 
1033  if( fParteff.GetParticleIndex(pdg)>=KFPartEfficiencies::fFirstMissingMassParticleIndex &&
1034  fParteff.GetParticleIndex(pdg)<=KFPartEfficiencies::fLastMissingMassParticleIndex )
1035  {
1036  isRecPart.push_back(part.IsReconstructable(4));
1037  isRecPart.push_back(part.IsReconstructable(1));
1038  isRecPart.push_back(part.IsReconstructable(3));
1039  }
1040  else
1041  for(int iEff = 0; iEff < 3; iEff++)
1042  isRecPart.push_back(part.IsReconstructable(iEff));
1043 
1044  isReconstructable.push_back(isRecPart);
1045  isReco.push_back( MCtoRParticleId[iP].ids.size() != 0 );
1046  nClones.push_back( MCtoRParticleId[iP].ids.size() - 1 );
1047 
1048  if(decays.empty() && (part.IsReconstructableV0(0) || part.IsReconstructableV0(1) || part.IsReconstructableV0(2)) )
1049  {
1050  iParticle.push_back(fParteff.nParticles - 1);
1051  vector<bool> isRecV0;
1052  for(int iEff = 0; iEff < 3; iEff++)
1053  isRecV0.push_back(part.IsReconstructableV0(iEff));
1054  isReconstructable.push_back(isRecV0);
1055  isReco.push_back( (MCtoRParticleId[iP].ids.size() != 0) || (MCtoRParticleId[iP].idsMI.size() != 0) );
1056 
1057  int nClonesV0 = MCtoRParticleId[iP].ids.size() + MCtoRParticleId[iP].idsMI.size() - 1;
1058  nClones.push_back( nClonesV0 );
1059  }
1060 
1061  {
1062  for(unsigned int iPType=0; iPType<iParticle.size(); iPType++)
1063  {
1064  int iPart = iParticle[iPType];
1065  if(iPart<0) continue;
1066 
1067  partEff.Inc(isReco[iPType], nClones[iPType], isReconstructable[iPType][0], isReconstructable[iPType][1], isReconstructable[iPType][2], fParteff.partName[iPart].data());
1068  if ( mId == -1 )
1069  partEff.Inc(isReco[iPType], nClones[iPType], isReconstructable[iPType][0], isReconstructable[iPType][1], isReconstructable[iPType][2], (fParteff.partName[iPart]+"_prim").data());
1070  else
1071  partEff.Inc(isReco[iPType], nClones[iPType], isReconstructable[iPType][0], isReconstructable[iPType][1], isReconstructable[iPType][2], (fParteff.partName[iPart]+"_sec").data());
1072 
1073  for(int iEff=0; iEff<3; iEff++)
1074  {
1075  if(!isReconstructable[iPType][iEff]) continue;
1076 
1077  int iMCTrack = part.GetMCTrackID();
1078  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
1079 
1080  Double_t massMC = fParteff.partMass[iPart];
1081  Double_t E = sqrt(mcTrack.P()*mcTrack.P() + massMC*massMC);
1082  Double_t Y = 0.5*log((E + mcTrack.Pz())/(E - mcTrack.Pz()));
1083  Double_t Z = mcTrack.Z();
1084  Double_t R = -1, L=-1;
1085  Double_t Mt_mc = sqrt(mcTrack.Pt()*mcTrack.Pt()+massMC*massMC)-massMC;
1086  Double_t cT = -1.e10;
1087  Double_t decayLength = -1.e10;
1088 
1089  if(part.NDaughters() > 0)
1090  {
1091  int mcDaughterId = part.GetDaughterIds()[0];
1092  KFMCTrack &mcDaughter = vMCTracks[mcDaughterId];
1093  R = sqrt(mcDaughter.X()*mcDaughter.X() + mcDaughter.Y()*mcDaughter.Y());
1094  L = sqrt(mcDaughter.X()*mcDaughter.X() + mcDaughter.Y()*mcDaughter.Y());
1095  Z = mcDaughter.Z();
1096 
1097  if(mcTrack.MotherId() < 0)
1098  {
1099  KFParticle motherKFParticle;
1100  float decayPoint[3] = { mcDaughter.X(), mcDaughter.Y(), mcDaughter.Z() };
1101  for(int iP=0; iP<6; iP++)
1102  motherKFParticle.Parameter(iP) = mcTrack.Par()[iP];
1103 
1104  float dsdr[6];
1105  double s = motherKFParticle.GetDStoPoint(decayPoint, dsdr);
1106  int jParticlePDG = fParteff.GetParticleIndex(mcTrack.PDG());
1107  Double_t massMC = (jParticlePDG>=0) ? fParteff.partMass[jParticlePDG] :0.13957;
1108 
1109  cT = s*massMC;
1110  decayLength = s*mcTrack.P();
1111  }
1112  }
1113 
1114  if(fStoreMCHistograms)
1115  {
1116  hPartEfficiency[iPart][iEff][0]->Fill( mcTrack.P(), isReco[iPType] );
1117  hPartEfficiency[iPart][iEff][1]->Fill( mcTrack.Pt(), isReco[iPType] );
1118  hPartEfficiency[iPart][iEff][2]->Fill( Y, isReco[iPType] );
1119  hPartEfficiency[iPart][iEff][3]->Fill( Z, isReco[iPType] );
1120  if(cT > -1.e10) hPartEfficiency[iPart][iEff][4]->Fill( cT, isReco[iPType] );
1121  if(decayLength > -1.e10) hPartEfficiency[iPart][iEff][5]->Fill( decayLength, isReco[iPType] );
1122  hPartEfficiency[iPart][iEff][3]->Fill( Z, isReco[iPType] );
1123  hPartEfficiency[iPart][iEff][6]->Fill( L, isReco[iPType] );
1124  hPartEfficiency[iPart][iEff][7]->Fill( R, isReco[iPType] );
1125  hPartEfficiency[iPart][iEff][8]->Fill( Mt_mc, isReco[iPType] );
1126 
1127  hPartEfficiency2D[iPart][iEff][0]->Fill( Y, mcTrack.Pt(), isReco[iPType] );
1128  hPartEfficiency2D[iPart][iEff][1]->Fill( Y, Mt_mc, isReco[iPType] );
1129  }
1130  }
1131  }
1132  }
1133  }
1134 
1135  fNEvents++;
1136 
1137  fParteff += partEff;
1138 
1139  partEff.CalcEff();
1140  fParteff.CalcEff();
1141 
1142  if(fNEvents%fPrintEffFrequency == 0)
1143  {
1144  std::cout << " ---- KF Particle finder --- " << std::endl;
1145  std::cout << "ACCUMULATED STAT : " << fNEvents << " EVENTS " << std::endl << std::endl;
1146  fParteff.PrintEff();
1147  std::cout<<std::endl;
1148  }
1149 }
1150 
1151 void KFTopoPerformance::CalculatePVEfficiency()
1152 {
1154  KFPVEfficiencies pvEff; // efficiencies for current event
1155  KFPVEfficiencies pvEffMCReconstructable;
1156  int nTracks = 0;
1157  //calculate N reco tracks
1158  for( unsigned int iRP = 0; iRP < fTopoReconstructor->GetParticles().size(); iRP++ ) {
1159 // CbmKFParticle &rPart = fTopoReconstructor->GetParticles()[iRP];
1160  const KFParticle &rPart = fTopoReconstructor->GetParticles()[iRP];
1161 
1162  if (rPart.NDaughters() != 1) continue;
1163 
1164  nTracks++;
1165  }
1166 
1167  const int NRecoPV = fTopoReconstructor->NPrimaryVertices();
1168  for ( int iP = 0; iP < NRecoPV; ++iP ) {
1169 
1170  const bool isBG = RtoMCPVId[iP].idsMI.size() != 0;
1171  const bool isGhost = !RtoMCPVId[iP].IsMatched();
1172 
1173  pvEff.IncReco(isGhost, isBG, "PV");
1174  pvEffMCReconstructable.IncReco(isGhost, isBG, "PV");
1175  }
1176 
1177  const int NMCPV = fPrimVertices.size();
1178  for ( int iV = 0; iV < NMCPV; ++iV ) {
1179  const KFMCVertex &pvMC = fPrimVertices[iV];
1180  if ( pvMC.IsReconstructable() )
1181  {
1182  const bool isReco = pvMC.IsReconstructed();
1183  const int nClones = MCtoRPVId[iV].ids.size() - 1;
1184 
1185 
1186  pvEff.Inc(isReco, nClones, "PV");
1187  if ( pvMC.IsTriggerPV() )
1188  {
1189  pvEff.Inc(isReco, nClones, "PVtrigger");
1190  hPVefficiency[0][0]->Fill( pvMC.NDaughterTracks(), isReco );
1191  hPVefficiency[0][1]->Fill( NMCPV, isReco );
1192  hPVefficiency[0][2]->Fill( vMCTracks.size(), isReco );
1193  hPVefficiency[0][3]->Fill( pvMC.NReconstructedDaughterTracks(), isReco );
1194  hPVefficiency[0][4]->Fill( NRecoPV, isReco );
1195  hPVefficiency[0][5]->Fill( nTracks, isReco );
1196  }
1197  else
1198  {
1199  pvEff.Inc(isReco, nClones, "PVpileup");
1200  hPVefficiency[1][0]->Fill( pvMC.NDaughterTracks(), isReco );
1201  hPVefficiency[1][1]->Fill( NMCPV, isReco );
1202  hPVefficiency[1][2]->Fill( vMCTracks.size(), isReco );
1203  hPVefficiency[1][3]->Fill( pvMC.NReconstructedDaughterTracks(), isReco );
1204  hPVefficiency[1][4]->Fill( NRecoPV, isReco );
1205  hPVefficiency[1][5]->Fill( nTracks, isReco );
1206  }
1207  }
1208  if ( pvMC.IsMCReconstructable() )
1209  {
1210  const bool isReco = pvMC.IsReconstructed();
1211  const int nClones = MCtoRPVId[iV].ids.size() - 1;
1212 
1213 
1214  pvEffMCReconstructable.Inc(isReco, nClones, "PV");
1215  if ( pvMC.IsTriggerPV() )
1216  {
1217  pvEffMCReconstructable.Inc(isReco, nClones, "PVtrigger");
1218  hPVefficiency[2][0]->Fill( pvMC.NDaughterTracks(), isReco );
1219  hPVefficiency[2][1]->Fill( NMCPV, isReco );
1220  hPVefficiency[2][2]->Fill( vMCTracks.size(), isReco );
1221  hPVefficiency[2][3]->Fill( pvMC.NReconstructedDaughterTracks(), isReco );
1222  hPVefficiency[2][4]->Fill( NRecoPV, isReco );
1223  hPVefficiency[2][5]->Fill( nTracks, isReco );
1224  }
1225  else
1226  {
1227  pvEffMCReconstructable.Inc(isReco, nClones, "PVpileup");
1228  hPVefficiency[3][0]->Fill( pvMC.NDaughterTracks(), isReco );
1229  hPVefficiency[3][1]->Fill( NMCPV, isReco );
1230  hPVefficiency[3][2]->Fill( vMCTracks.size(), isReco );
1231  hPVefficiency[3][3]->Fill( pvMC.NReconstructedDaughterTracks(), isReco );
1232  hPVefficiency[3][4]->Fill( NRecoPV, isReco );
1233  hPVefficiency[3][5]->Fill( nTracks, isReco );
1234  }
1235  }
1236  }
1237 
1238  fPVeff += pvEff;
1239 
1240  pvEff.CalcEff();
1241  fPVeff.CalcEff();
1242 
1243  fPVeffMCReconstructable += pvEffMCReconstructable;
1244  pvEffMCReconstructable.CalcEff();
1245  fPVeffMCReconstructable.CalcEff();
1246 
1247  if(fNEvents%fPrintEffFrequency == 0)
1248  {
1249  std::cout << " ---- KF PV finder --- " << std::endl;
1250  std::cout << "ACCUMULATED STAT : " << fNEvents << " EVENTS " << std::endl << std::endl;
1251  std::cout << "PV with at least 2 reconstructed tracks is reconstructable:" << std::endl;
1252  fPVeff.PrintEff();
1253  std::cout << std::endl;
1254  std::cout << "PV with at least 2 MC tracks with 15 MC points is reconstructable:" << std::endl;
1255  fPVeffMCReconstructable.PrintEff();
1256 
1257  std::cout<<std::endl;
1258  }
1259 }
1260 
1261 void KFTopoPerformance::FillParticleParameters(KFParticle& TempPart,
1262  int iParticle,
1263  int iP,
1264  int iPV,
1265  TH1F* histoParameters[nParametersSet][KFPartEfficiencies::nParticles][nHistoPartParam],
1266  TH2F* histoParameters2D[nParametersSet][KFPartEfficiencies::nParticles][nHistoPartParam2D],
1267  TH3F* histoParameters3D[1][KFPartEfficiencies::nParticles][nHistoPartParam3D],
1268  TH1F* histoFit[KFPartEfficiencies::nParticles][nFitQA],
1269  TH1F* histoFitDaughtersQA[KFPartEfficiencies::nParticles][nFitQA],
1270  TH1F* histoDSToParticleQA[KFPartEfficiencies::nParticles][nDSToParticleQA],
1271  vector<int>* multiplicities)
1272 {
1275  const std::map<int,bool>& decays = fTopoReconstructor->GetKFParticleFinder()->GetReconstructionList();
1276  if(!(decays.empty()) && (iParticle < fParteff.fFirstStableParticleIndex || iParticle > fParteff.fLastStableParticleIndex))
1277  if(decays.find(TempPart.GetPDG()) == decays.end()) return;
1278 
1279  float M, M_t, ErrM;
1280  float dL, ErrdL; // decay length
1281  float cT, ErrcT; // c*tau
1282  float P, ErrP;
1283  float Pt;
1284  float Rapidity;
1285  float Theta;
1286  float Phi;
1287  float X,Y,Z,R;
1288  float QtAlpha[2];
1289 
1290  TempPart.GetMass(M,ErrM);
1291  TempPart.GetMomentum(P,ErrP);
1292  Pt = TempPart.GetPt();
1293  Rapidity = TempPart.GetRapidity();
1294 
1295  KFParticle TempPartTopo = TempPart;
1296  TempPartTopo.SetProductionVertex(fTopoReconstructor->GetPrimVertex(0));
1297  TempPartTopo.GetDecayLength(dL,ErrdL);
1298  TempPartTopo.GetLifeTime(cT,ErrcT);
1299 
1300  float chi2 = TempPart.GetChi2();
1301  Int_t ndf = TempPart.GetNDF();
1302  float prob = TMath::Prob(chi2, ndf);//(TDHelper<float>::Chi2IProbability( ndf, chi2 ));
1303  Theta = TempPart.GetTheta();
1304  Phi = TempPart.GetPhi();
1305  X = TempPart.GetX();
1306  Y = TempPart.GetY();
1307  Z = TempPart.GetZ();
1308 #ifdef CBM
1309  if(Z>=1. && iParticle>=54 && iParticle<=64) return;
1310 #endif
1311  R = sqrt(X*X+Y*Y);
1312  M_t = sqrt(Pt*Pt+fParteff.GetMass(iParticle)*fParteff.GetMass(iParticle))-fParteff.GetMass(iParticle);
1313 
1314  KFParticleSIMD tempSIMDPart(TempPart);
1315  float_v l,dl;
1316  KFParticleSIMD pv(fTopoReconstructor->GetPrimVertex(iPV));
1317  tempSIMDPart.GetDistanceToVertexLine(pv, l, dl);
1318 #ifdef __ROOT__
1319  if( (l[0] > 0.2f || Pt < 0.f) && (abs( TempPart.GetPDG() ) == 4122 ||
1320  abs( TempPart.GetPDG() ) == 104122 ||
1321  abs( TempPart.GetPDG() ) == 204122 ||
1322  abs( TempPart.GetPDG() ) == 304122 ||
1323  abs( TempPart.GetPDG() ) == 404122 ||
1324  abs( TempPart.GetPDG() ) == 504122 ) ) return;
1325  if( (l[0] > 0.2f || Pt < 0.f) && (abs( TempPart.GetPDG() ) == 421 ||
1326  abs( TempPart.GetPDG() ) == 420 ||
1327  abs( TempPart.GetPDG() ) == 425 ||
1328  abs( TempPart.GetPDG() ) == 426 ||
1329  abs( TempPart.GetPDG() ) == 427 ||
1330  abs( TempPart.GetPDG() ) == 429) ) return;
1331  if( (l[0] > 0.4f || Pt < 0.f) && (abs( TempPart.GetPDG() ) == 411 ||
1332  abs( TempPart.GetPDG() ) == 100411 ||
1333  abs( TempPart.GetPDG() ) == 200411 ||
1334  abs( TempPart.GetPDG() ) == 300411) ) return;
1335  if( (l[0] > 0.2f || Pt < 0.f) && (abs( TempPart.GetPDG() ) == 431 ||
1336  abs( TempPart.GetPDG() ) == 100431 ||
1337  abs( TempPart.GetPDG() ) == 200431 ||
1338  abs( TempPart.GetPDG() ) == 300431 ||
1339  abs( TempPart.GetPDG() ) == 400431) ) return;
1340 
1341 // if(Pt < 2. && (abs( TempPart.GetPDG() ) == 443 ||
1342 // abs( TempPart.GetPDG() ) == 100443 ||
1343 // abs( TempPart.GetPDG() ) == 200443 ||
1344 // abs( TempPart.GetPDG() ) == 300443 ||
1345 // abs( TempPart.GetPDG() ) == 400443 ||
1346 // abs( TempPart.GetPDG() ) == 500443) ) return;
1347 
1348  if(Pt < 0.5f && (abs( TempPart.GetPDG() ) == 3000 ||
1349  abs( TempPart.GetPDG() ) == 3001) ) return;
1350 #endif
1351  float parameters[17] = {M, P, Pt, Rapidity, dL, cT, chi2/ndf, prob, Theta, Phi, X, Y, Z, R, l[0], l[0]/dl[0], M_t };
1352 
1353  //for all particle-candidates
1354  for(int iParam=0; iParam<17; iParam++)
1355  histoParameters[0][iParticle][iParam]->Fill(parameters[iParam]);
1356 
1357  if(multiplicities)
1358  multiplicities[0][iParticle]++;
1359 
1360  histoParameters2D[0][iParticle][0]->Fill(Rapidity,Pt,1);
1361  histoParameters2D[0][iParticle][3]->Fill(Rapidity,M_t,1);
1362 
1363  const bool drawZR = IsCollectZRHistogram(iParticle);
1364  if(histoParameters2D[0][iParticle][1] && drawZR)
1365  {
1366  histoParameters2D[0][iParticle][1]->Fill(Z,R,1);
1367  }
1368 
1369  if(TempPart.NDaughters() == 2 && IsCollectArmenteros(iParticle))
1370  {
1371  int index1 = TempPart.DaughterIds()[0];
1372  int index2 = TempPart.DaughterIds()[1];
1373  if(index1 >= int(fTopoReconstructor->GetParticles().size()) ||
1374  index2 >= int(fTopoReconstructor->GetParticles().size()) ||
1375  index1 < 0 || index2 < 0 )
1376  return;
1377 
1378  KFParticle posDaughter, negDaughter;
1379  if(int(fTopoReconstructor->GetParticles()[index1].Q()) > 0)
1380  {
1381  posDaughter = fTopoReconstructor->GetParticles()[index1];
1382  negDaughter = fTopoReconstructor->GetParticles()[index2];
1383  }
1384  else
1385  {
1386  negDaughter = fTopoReconstructor->GetParticles()[index1];
1387  posDaughter = fTopoReconstructor->GetParticles()[index2];
1388  }
1389  float vertex[3] = {TempPart.GetX(), TempPart.GetY(), TempPart.GetZ()};
1390  posDaughter.TransportToPoint(vertex);
1391  negDaughter.TransportToPoint(vertex);
1392  KFParticle::GetArmenterosPodolanski(posDaughter, negDaughter, QtAlpha );
1393 
1394  histoParameters2D[0][iParticle][2]->Fill(QtAlpha[1],QtAlpha[0],1);
1395  }
1396 
1397  //Fill 3D histograms for multi differential analysis
1398  if( histoParameters3D && IsCollect3DHistogram(iParticle))
1399  {
1400  histoParameters3D[0][iParticle][0]->Fill(Rapidity,Pt,M,1);
1401  histoParameters3D[0][iParticle][1]->Fill(Rapidity,M_t,M,1);
1402  if(fCentralityBin>=0)
1403  {
1404  histoParameters3D[0][iParticle][2]->Fill(fCentralityBin, Pt, M, fCentralityWeight);
1405  histoParameters3D[0][iParticle][3]->Fill(fCentralityBin, Rapidity, M, fCentralityWeight);
1406  histoParameters3D[0][iParticle][4]->Fill(fCentralityBin, M_t, M, fCentralityWeight);
1407  }
1408  histoParameters3D[0][iParticle][5]->Fill(cT, Pt, M, 1);
1409  }
1410 
1411  //Fill histograms for the side bands analysis
1412  if(histoDSToParticleQA && IsCollect3DHistogram(iParticle))
1413  {
1414  if(fabs(fParteff.GetMass(iParticle)-M) < 3.f*fParteff.GetMassSigma(iParticle))//SignalReco
1415  {
1416  for(int iParam=0; iParam<17; iParam++)
1417  histoParameters[4][iParticle][iParam]->Fill(parameters[iParam]);
1418 
1419  if(multiplicities)
1420  multiplicities[4][iParticle]++;
1421 
1422  histoParameters2D[4][iParticle][0]->Fill(Rapidity,Pt,1);
1423  histoParameters2D[4][iParticle][3]->Fill(Rapidity,M_t,1);
1424 
1425  if(drawZR)
1426  {
1427  if(histoParameters2D[4][iParticle][1])
1428  histoParameters2D[4][iParticle][1]->Fill(Z,R,1);
1429  if(histoParameters2D[4][iParticle][2])
1430  histoParameters2D[4][iParticle][2]->Fill(QtAlpha[1],QtAlpha[0],1);
1431  }
1432  }
1433 
1434  if( fabs(fParteff.GetMass(iParticle)-M) > 3.f*fParteff.GetMassSigma(iParticle) &&
1435  fabs(fParteff.GetMass(iParticle)-M) <= 6.f*fParteff.GetMassSigma(iParticle) )//BGReco
1436  {
1437  for(int iParam=0; iParam<17; iParam++)
1438  histoParameters[5][iParticle][iParam]->Fill(parameters[iParam]);
1439 
1440  if(multiplicities)
1441  multiplicities[5][iParticle]++;
1442 
1443  histoParameters2D[5][iParticle][0]->Fill(Rapidity,Pt,1);
1444  histoParameters2D[5][iParticle][3]->Fill(Rapidity,M_t,1);
1445 
1446  if(drawZR)
1447  {
1448  if(histoParameters2D[5][iParticle][1])
1449  histoParameters2D[5][iParticle][1]->Fill(Z,R,1);
1450  if(histoParameters2D[5][iParticle][2])
1451  histoParameters2D[5][iParticle][2]->Fill(QtAlpha[1],QtAlpha[0],1);
1452  }
1453  }
1454  }
1455 
1456  if(!fStoreMCHistograms) return;
1457 
1458  int iSet = 1;
1459  if(!RtoMCParticleId[iP].IsMatchedWithPdg()) //background
1460  {
1461  if(!RtoMCParticleId[iP].IsMatched()) iSet = 3; // for ghost particles - combinatorial background
1462  else iSet = 2; // for physical background
1463  }
1464 
1465  //Check if PV association is correct
1466  if(!histoDSToParticleQA && iSet == 1)
1467  {
1468  int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
1469  KFMCParticle &mcPart = vMCParticles[iMCPart];
1470  int iMCTrack = mcPart.GetMCTrackID();
1471  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
1472  int motherId = mcTrack.MotherId();
1473  bool isSecondaryParticle = motherId >= 0;
1474 
1475  if(iPV >=0)
1476  {
1477  if(isSecondaryParticle)
1478  iSet = 4;
1479  else
1480  {
1481  int iMCPV = -1;
1482  if(RtoMCPVId[iPV].IsMatchedWithPdg())
1483  iMCPV = RtoMCPVId[iPV].GetBestMatch();
1484 
1485  int iMCPVFromParticle = fMCTrackToMCPVMatch[iMCTrack];
1486  if(iMCPV != iMCPVFromParticle)
1487  iSet = 4;
1488  }
1489  }
1490  else
1491  {
1492  if(!isSecondaryParticle)
1493  iSet = 4;
1494  }
1495  }
1496 
1497  //for signal particles
1498  for(int iParam=0; iParam<17; iParam++)
1499  histoParameters[iSet][iParticle][iParam]->Fill(parameters[iParam]);
1500 
1501  if(multiplicities)
1502  multiplicities[iSet][iParticle]++;
1503 
1504  histoParameters2D[iSet][iParticle][0]->Fill(Rapidity,Pt,1);
1505  if(drawZR)
1506  {
1507  if(histoParameters2D[iSet][iParticle][1])
1508  histoParameters2D[iSet][iParticle][1]->Fill(Z,R,1);
1509  if(histoParameters2D[iSet][iParticle][2])
1510  histoParameters2D[iSet][iParticle][2]->Fill(QtAlpha[1],QtAlpha[0],1);
1511  }
1512  histoParameters2D[iSet][iParticle][3]->Fill(Rapidity,M_t,1);
1513 
1514  if(iSet != 1) return;
1515 
1516  int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
1517  KFMCParticle &mcPart = vMCParticles[iMCPart];
1518  // Fit quality of the mother particle
1519  if(histoFit)
1520  {
1521  int iMCTrack = mcPart.GetMCTrackID();
1522  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
1523  int mcDaughterId = -1;
1524  if(iParticle >= fParteff.fFirstStableParticleIndex && iParticle <= fParteff.fLastStableParticleIndex)
1525  mcDaughterId = iMCTrack;
1526  else if(mcTrack.PDG() == 22 && TempPart.NDaughters() == 1)
1527  mcDaughterId = iMCTrack;
1528  else if(iParticle >= fParteff.fFirstMissingMassParticleIndex && iParticle <= fParteff.fLastMissingMassParticleIndex)
1529  mcDaughterId = mcPart.GetDaughterIds()[1]; //the charged daughter
1530  else
1531  mcDaughterId = mcPart.GetDaughterIds()[0];
1532 
1533  KFMCTrack &mcDaughter = vMCTracks[mcDaughterId];
1534 
1535  float mcX = mcTrack.X();
1536  float mcY = mcTrack.Y();
1537  float mcZ = mcTrack.Z();
1538  if(histoDSToParticleQA || hPartParamPrimary == histoParameters)
1539  {
1540  mcX = mcDaughter.X();
1541  mcY = mcDaughter.Y();
1542  mcZ = mcDaughter.Z();
1543  }
1544  const float mcPx = mcTrack.Par(3);
1545  const float mcPy = mcTrack.Par(4);
1546  const float mcPz = mcTrack.Par(5);
1547 
1548  float decayVtx[3] = { mcTrack.X(), mcTrack.Y(), mcTrack.Z() };
1549  float recParam[8] = { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f };
1550  float errParam[8] = { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f };
1551 
1552  for(int iPar=0; iPar<3; iPar++)
1553  {
1554  recParam[iPar] = TempPart.GetParameter(iPar);
1555  Double_t error = TempPart.GetCovariance(iPar,iPar);
1556  if(error < 0.) { error = 1.e20;}
1557  errParam[iPar] = TMath::Sqrt(error);
1558  }
1559  TempPart.TransportToPoint(decayVtx);
1560  for(int iPar=3; iPar<7; iPar++)
1561  {
1562  recParam[iPar] = TempPart.GetParameter(iPar);
1563  Double_t error = TempPart.GetCovariance(iPar,iPar);
1564  if(error < 0.) { error = 1.e20;}
1565  errParam[iPar] = TMath::Sqrt(error);
1566  }
1567 
1568  int jParticlePDG = fParteff.GetParticleIndex(mcTrack.PDG());
1569  Double_t massMC = (jParticlePDG>=0) ? fParteff.partMass[jParticlePDG] :0.13957;
1570 
1571  Double_t Emc = sqrt(mcTrack.P()*mcTrack.P() + massMC*massMC);
1572  Double_t res[8] = {0},
1573  pull[8] = {0},
1574  mcParam[8] = { mcX, mcY, mcZ,
1575  mcPx, mcPy, mcPz, Emc, massMC };
1576  for(int iPar=0; iPar < 7; iPar++ )
1577  {
1578  res[iPar] = recParam[iPar] - mcParam[iPar];
1579  if(fabs(errParam[iPar]) > 1.e-20) pull[iPar] = res[iPar]/errParam[iPar];
1580  }
1581 
1582  res[7] = M - mcParam[7];
1583  if(fabs(ErrM) > 1.e-20) pull[7] = res[7]/ErrM;
1584 
1585  for(int iPar=0; iPar < 8; iPar++ )
1586  {
1587  histoFit[iParticle][iPar]->Fill(res[iPar]);
1588  histoFit[iParticle][iPar+8]->Fill(pull[iPar]);
1589  }
1590  }
1591  // Fit quality of daughters
1592  int daughterIndex[2] = {-1, -1};
1593 
1594  if(histoFitDaughtersQA)
1595  {
1596  for(int iD=0; iD<mcPart.NDaughters(); ++iD)
1597  {
1598  int mcDaughterId = mcPart.GetDaughterIds()[iD];
1599  // if(!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg()) continue;
1600  if(!MCtoRParticleId[mcDaughterId].IsMatched()) continue;
1601  KFMCTrack &mcTrack = vMCTracks[mcDaughterId];
1602  // int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
1603  int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatch();
1604  KFParticle Daughter = fTopoReconstructor->GetParticles()[recDaughterId];
1605  Daughter.GetMass(M,ErrM);
1606 
1607  const float mcX = mcTrack.X();
1608  const float mcY = mcTrack.Y();
1609  const float mcZ = mcTrack.Z();
1610  const float mcPx = mcTrack.Px();
1611  const float mcPy = mcTrack.Py();
1612  const float mcPz = mcTrack.Pz();
1613 
1614  float_v decayVtx[3] = {mcX, mcY, mcZ};
1615 
1616  KFParticleSIMD DaughterSIMD(Daughter);
1617  DaughterSIMD.TransportToPoint(decayVtx);
1618 
1619  int jParticlePDG = fParteff.GetParticleIndex(mcTrack.PDG());
1620  Double_t massMC = (jParticlePDG>=0) ? fParteff.partMass[jParticlePDG] :0.13957;
1621  Double_t Emc = sqrt(mcTrack.P()*mcTrack.P() + massMC*massMC);
1622 
1623  Double_t res[8] = {0},
1624  pull[8] = {0},
1625  mcParam[8] = { mcX, mcY, mcZ,
1626  mcPx, mcPy, mcPz, Emc, massMC };
1627  for(int iPar=0; iPar < 7; iPar++ )
1628  {
1629  Double_t error = DaughterSIMD.GetCovariance(iPar,iPar)[0];
1630  if(error < 0.) { error = 1.e20;}
1631  error = TMath::Sqrt(error);
1632  Double_t recoPar = DaughterSIMD.GetParameter(iPar)[0];
1633  res[iPar] = recoPar - mcParam[iPar];
1634  if(fabs(error) > 1.e-20) pull[iPar] = res[iPar]/error;
1635  }
1636  res[7] = M - mcParam[7];
1637  if(fabs(ErrM) > 1.e-20) pull[7] = res[7]/ErrM;
1638 
1639  for(int iPar=0; iPar < 8; iPar++ )
1640  {
1641  histoFitDaughtersQA[iParticle][iPar]->Fill(res[iPar]);
1642  histoFitDaughtersQA[iParticle][iPar+8]->Fill(pull[iPar]);
1643  }
1644 
1645  //fill Histos for GetDStoParticle
1646  if(iD == 0)
1647  daughterIndex[0] = recDaughterId;
1648  if(iD == 1 && daughterIndex[0] > -1 && histoDSToParticleQA)
1649  {
1650  daughterIndex[1] = recDaughterId;
1651  KFParticle d1 = fTopoReconstructor->GetParticles()[daughterIndex[0]];
1652  KFParticle d2 = fTopoReconstructor->GetParticles()[daughterIndex[1]];
1653 
1654  KFParticleSIMD daughters[2] = {d2, d1};
1655 
1656  float_v dS[2] = {0.f, 0.f};
1657  float_v dsdr[4][6];
1658  for(int i1=0; i1<4; i1++)
1659  for(int i2=0; i2<6; i2++)
1660  dsdr[i1][i2] = 0.f;
1661 
1662  daughters[0].GetDStoParticle(daughters[1], dS, dsdr);
1663  float_v pD[2][8], cD[2][36], corrPD[2][36], corrCD[2][36];
1664 
1665  for(int iDR=0; iDR<2; iDR++)
1666  {
1667  for(int iPD = 0; iPD<8; iPD++)
1668  {
1669  pD[iDR][iPD] = 0;
1670  corrPD[iDR][iPD] = 0;
1671  }
1672  for(int iCD=0; iCD<36; iCD++)
1673  {
1674  cD[iDR][iCD] = 0;
1675  corrCD[iDR][iCD] = 0;
1676  }
1677  }
1678 
1679  float_v F[4][36];
1680  {
1681  for(int i1=0; i1<4; i1++)
1682  for(int i2=0; i2<36; i2++)
1683  F[i1][i2] = 0;
1684  }
1685  daughters[0].Transport(dS[0], dsdr[0], pD[0], cD[0], dsdr[1], F[0], F[1]);
1686  daughters[1].Transport(dS[1], dsdr[3], pD[1], cD[1], dsdr[2], F[3], F[2]);
1687 
1688  daughters[0].MultQSQt( F[1], daughters[1].CovarianceMatrix(), corrCD[0], 6);
1689  daughters[0].MultQSQt( F[2], daughters[0].CovarianceMatrix(), corrCD[1], 6);
1690  for(int iDR=0; iDR<2; iDR++)
1691  for(int iC=0; iC<6; iC++)
1692  cD[iDR][iC] += corrCD[iDR][iC];
1693 
1694  for(int iDR=0; iDR<2; iDR++)
1695  {
1696  cD[iDR][1] = cD[iDR][2];
1697  cD[iDR][2] = cD[iDR][5];
1698  for(int iPar=0; iPar<3; iPar++)
1699  {
1700  res[iPar] = pD[iDR][iPar][0] - decayVtx[iPar][0];
1701 
1702  Double_t error = cD[iDR][iPar][0];
1703  if(error < 0.) { error = 1.e20;}
1704  error = sqrt(error);
1705 
1706  pull[iPar] = res[iPar] / error;
1707 
1708  histoDSToParticleQA[iParticle][iPar]->Fill(res[iPar]);
1709  histoDSToParticleQA[iParticle][iPar+3]->Fill(pull[iPar]);
1710  }
1711  }
1712 
1713  Double_t dXds = pD[0][0][0] - pD[1][0][0];
1714  Double_t dYds = pD[0][1][0] - pD[1][1][0];
1715  Double_t dZds = pD[0][2][0] - pD[1][2][0];
1716 
1717  Double_t dRds = sqrt(dXds*dXds + dYds*dYds + dZds*dZds);
1718  histoDSToParticleQA[iParticle][6]->Fill(dRds);
1719  }
1720  }
1721  }
1722 }
1723 
1724 void KFTopoPerformance::FillHistos()
1725 {
1727  vector<int> multiplicities[6];
1728  for(int iV=0; iV<6; iV++)
1729  multiplicities[iV].resize(KFPartEfficiencies::nParticles, 0);
1730 
1731  //fill histograms for found short-lived particles
1732  for(unsigned int iP=0; iP<fTopoReconstructor->GetParticles().size(); iP++)
1733  {
1734  int iParticle = fParteff.GetParticleIndex(fTopoReconstructor->GetParticles()[iP].GetPDG());
1735  if(iParticle < 0) continue;
1736  KFParticle TempPart = fTopoReconstructor->GetParticles()[iP];
1737 
1738  FillParticleParameters(TempPart,iParticle, iP, 0, hPartParam, hPartParam2D, hPartParam3D,
1739  hFitQA, hFitDaughtersQA, hDSToParticleQA, multiplicities);
1740  }
1741 
1742  if(fStoreMCHistograms)
1743  {
1744  for(int iSet=0; iSet<KFParticleFinder::GetNSecondarySets(); iSet++)
1745  {
1746  const std::vector<KFParticle>& SecondaryCandidates = fTopoReconstructor->GetKFParticleFinder()->GetSecondaryCandidates()[iSet];
1747  for(unsigned int iP=0; iP<SecondaryCandidates.size(); iP++)
1748  {
1749  KFParticle TempPart = SecondaryCandidates[iP];
1750  int iParticle = fParteff.GetParticleIndex(TempPart.GetPDG());
1751  if(iParticle < 0) continue;
1752 
1753  const int id = TempPart.Id();
1754  FillParticleParameters(TempPart, iParticle, id, 0, hPartParamSecondaryMass, hPartParam2DSecondaryMass, 0);
1755 
1756  TempPart = fTopoReconstructor->GetParticles()[id];
1757  FillParticleParameters(TempPart, iParticle, id, 0, hPartParamSecondary, hPartParam2DSecondary, 0);
1758  }
1759  }
1760 
1761  for(int iSet=0; iSet<KFParticleFinder::GetNPrimarySets(); iSet++)
1762  {
1763  for(int iPV=0; iPV<fTopoReconstructor->NPrimaryVertices(); iPV++)
1764  {
1765  const std::vector<KFParticle>& PrimaryCandidates = fTopoReconstructor->GetKFParticleFinder()->GetPrimaryCandidates()[iSet][iPV];
1766  for(unsigned int iP=0; iP<PrimaryCandidates.size(); iP++)
1767  {
1768  KFParticle TempPart = PrimaryCandidates[iP];
1769  int iParticle = fParteff.GetParticleIndex(TempPart.GetPDG());
1770  if(iParticle < 0) continue;
1771 
1772  const int id = TempPart.Id();
1773  FillParticleParameters(TempPart,iParticle, id, iPV, hPartParamPrimaryMass, hPartParam2DPrimaryMass, 0, hFitQAMassConstraint);
1774 
1775  TempPart = fTopoReconstructor->GetParticles()[id];
1776  FillParticleParameters(TempPart,iParticle, id, iPV, hPartParamPrimary, hPartParam2DPrimary, 0, hFitQANoConstraint);
1777  }
1778 
1779  const std::vector<KFParticle>& PrimaryCandidatesTopo = fTopoReconstructor->GetKFParticleFinder()->GetPrimaryTopoCandidates()[iSet][iPV];
1780  for(unsigned int iP=0; iP<PrimaryCandidatesTopo.size(); iP++)
1781  {
1782  KFParticle TempPart = PrimaryCandidatesTopo[iP];
1783  int iParticle = fParteff.GetParticleIndex(TempPart.GetPDG());
1784  if(iParticle < 0) continue;
1785 
1786  FillParticleParameters(TempPart,iParticle, TempPart.Id(), iPV, hPartParamPrimaryTopo, hPartParam2DPrimaryTopo, 0, hFitQATopoConstraint);
1787  }
1788 
1789  const std::vector<KFParticle>& PrimaryCandidatesTopoMass = fTopoReconstructor->GetKFParticleFinder()->GetPrimaryTopoMassCandidates()[iSet][iPV];
1790  for(unsigned int iP=0; iP<PrimaryCandidatesTopoMass.size(); iP++)
1791  {
1792  KFParticle TempPart = PrimaryCandidatesTopoMass[iP];
1793  int iParticle = fParteff.GetParticleIndex(TempPart.GetPDG());
1794  if(iParticle < 0) continue;
1795 
1796  FillParticleParameters(TempPart,iParticle, TempPart.Id(), iPV, hPartParamPrimaryTopoMass, hPartParam2DPrimaryTopoMass, 0, hFitQATopoMassConstraint);
1797  }
1798  }
1799  }
1800  }
1801  //fill histograms with ChiPrim for every particle
1802  for(unsigned int iP=0; iP<fTopoReconstructor->GetParticles().size(); iP++)
1803  {
1804  KFParticle TempPart = fTopoReconstructor->GetParticles()[iP];
1805  KFParticle vtx = fTopoReconstructor->GetPrimVertex(0);
1806 
1807  if(RtoMCParticleId[iP].IsMatched())
1808  {
1809  int iMCPV = vMCParticles[RtoMCParticleId[iP].GetBestMatch()].GetMotherId();
1810  if(iMCPV<0.)
1811  {
1812  iMCPV = -iMCPV - 1;
1813  if(MCtoRPVId[iMCPV].IsMatched())
1814  {
1815  vtx = fTopoReconstructor->GetPrimVertex(MCtoRPVId[iMCPV].GetBestMatch());
1816  }
1817  }
1818  }
1819 // else
1820 // KFParticle & vtx = fTopoReconstructor->GetPrimVertex(0);
1821 
1822 
1823  float chi2 = TempPart.GetDeviationFromVertex(vtx);
1824  int ndf = 2;
1825 
1826  hTrackParameters[KFPartEfficiencies::nParticles]->Fill(chi2);
1827  hTrackParameters[KFPartEfficiencies::nParticles+4]->Fill(TMath::Prob(chi2, ndf));
1828 
1829  if(!RtoMCParticleId[iP].IsMatched())
1830  {
1831  hTrackParameters[KFPartEfficiencies::nParticles+3]->Fill(chi2);
1832  hTrackParameters[KFPartEfficiencies::nParticles+7]->Fill(TMath::Prob(chi2, ndf));
1833  continue;
1834  }
1835 
1836  int iMCPart = RtoMCParticleId[iP].GetBestMatch();
1837  KFMCParticle &mcPart = vMCParticles[iMCPart];
1838  if(mcPart.GetMotherId() < 0)
1839  {
1840  hTrackParameters[KFPartEfficiencies::nParticles+1]->Fill(chi2 );
1841  hTrackParameters[KFPartEfficiencies::nParticles+5]->Fill(TMath::Prob(chi2, ndf));
1842  }
1843  else
1844  {
1845  hTrackParameters[KFPartEfficiencies::nParticles+2]->Fill(chi2 );
1846  hTrackParameters[KFPartEfficiencies::nParticles+6]->Fill(TMath::Prob(chi2, ndf));
1847  }
1848  int iParticle = fParteff.GetParticleIndex(fTopoReconstructor->GetParticles()[iP].GetPDG());
1849  if(iParticle > -1 && iParticle<KFPartEfficiencies::nParticles)
1850  hTrackParameters[iParticle]->Fill(chi2 );
1851 
1852  }
1853 
1854 
1855  //fill histograms of the primary vertex quality
1856  for(int iPV = 0; iPV<fTopoReconstructor->NPrimaryVertices(); iPV++)
1857  {
1858  KFParticle & vtx = fTopoReconstructor->GetPrimVertex(iPV);
1859  vector<int> &tracks = fTopoReconstructor->GetPVTrackIndexArray(iPV);
1860 
1861  Double_t probPV = TMath::Prob(vtx.Chi2(), vtx.NDF());//(TDHelper<float>::Chi2IProbability( ndf, chi2 ));
1862  vector<Double_t> dzPV;
1863  if(RtoMCPVId[iPV].IsMatched())
1864  {
1865  int iCurrMCPV = RtoMCPVId[iPV].GetBestMatch();
1866  for(int iPV2 = iPV+1; iPV2 < fTopoReconstructor->NPrimaryVertices(); iPV2++)
1867  {
1868  if(!RtoMCPVId[iPV2].IsMatched()) continue;
1869  int iCurrMCPV2 = RtoMCPVId[iPV2].GetBestMatch();
1870  if(iCurrMCPV != iCurrMCPV2) continue;
1871  KFParticle & vtx2 = fTopoReconstructor->GetPrimVertex(iPV2);
1872 
1873  dzPV.push_back(fabs(vtx.Z() - vtx2.Z()));
1874  }
1875  }
1876 
1877  hPVParam[ 0]->Fill(vtx.X());
1878  hPVParam[ 1]->Fill(vtx.Y());
1879  hPVParam[ 2]->Fill(vtx.Z());
1880  hPVParam[ 3]->Fill(sqrt(vtx.X()*vtx.X() + vtx.Y()*vtx.Y()));
1881  hPVParam[ 4]->Fill(tracks.size());
1882  hPVParam[ 5]->Fill(vtx.Chi2());
1883  hPVParam[ 6]->Fill(vtx.NDF());
1884  hPVParam[ 7]->Fill(vtx.Chi2()/vtx.NDF());
1885  hPVParam[ 8]->Fill(probPV);
1886  hPVParam[ 9]->Fill(fPVPurity[iPV]);
1887  hPVParam[10]->Fill(fPVTracksRate[0][iPV]);
1888  hPVParam[11]->Fill(fPVTracksRate[1][iPV]);
1889  hPVParam[12]->Fill(fPVTracksRate[2][iPV]);
1890  hPVParam[13]->Fill(fPVTracksRate[3][iPV]);
1891  for(unsigned int iZ=0; iZ<dzPV.size(); iZ++)
1892  hPVParam[14]->Fill(dzPV[iZ]);
1893 
1894  hPVParam2D[0]->Fill(vtx.X(),vtx.Y());
1895 
1896 
1897  if(!RtoMCPVId[iPV].IsMatchedWithPdg())
1898  {
1899  if(!RtoMCPVId[iPV].IsMatched())
1900  {
1901  hPVParamGhost[ 0]->Fill(vtx.X());
1902  hPVParamGhost[ 1]->Fill(vtx.Y());
1903  hPVParamGhost[ 2]->Fill(vtx.Z());
1904  hPVParamGhost[ 3]->Fill(sqrt(vtx.X()*vtx.X() + vtx.Y()*vtx.Y()));
1905  hPVParamGhost[ 4]->Fill(tracks.size());
1906  hPVParamGhost[ 5]->Fill(vtx.Chi2());
1907  hPVParamGhost[ 6]->Fill(vtx.NDF());
1908  hPVParamGhost[ 7]->Fill(vtx.Chi2()/vtx.NDF());
1909  hPVParamGhost[ 8]->Fill(probPV);
1910  hPVParamGhost[ 9]->Fill(fPVPurity[iPV]);
1911  hPVParamGhost[10]->Fill(fPVTracksRate[0][iPV]);
1912  hPVParamGhost[11]->Fill(fPVTracksRate[1][iPV]);
1913  hPVParamGhost[12]->Fill(fPVTracksRate[2][iPV]);
1914  hPVParamGhost[13]->Fill(fPVTracksRate[3][iPV]);
1915  for(unsigned int iZ=0; iZ<dzPV.size(); iZ++)
1916  hPVParamGhost[14]->Fill(dzPV[iZ]);
1917  }
1918  else
1919  {
1920  hPVParamBG[ 0]->Fill(vtx.X());
1921  hPVParamBG[ 1]->Fill(vtx.Y());
1922  hPVParamBG[ 2]->Fill(vtx.Z());
1923  hPVParamBG[ 3]->Fill(sqrt(vtx.X()*vtx.X() + vtx.Y()*vtx.Y()));
1924  hPVParamBG[ 4]->Fill(tracks.size());
1925  hPVParamBG[ 5]->Fill(vtx.Chi2());
1926  hPVParamBG[ 6]->Fill(vtx.NDF());
1927  hPVParamBG[ 7]->Fill(vtx.Chi2()/vtx.NDF());
1928  hPVParamBG[ 8]->Fill(probPV);
1929  hPVParamBG[ 9]->Fill(fPVPurity[iPV]);
1930  hPVParamBG[10]->Fill(fPVTracksRate[0][iPV]);
1931  hPVParamBG[11]->Fill(fPVTracksRate[1][iPV]);
1932  hPVParamBG[12]->Fill(fPVTracksRate[2][iPV]);
1933  hPVParamBG[13]->Fill(fPVTracksRate[3][iPV]);
1934  for(unsigned int iZ=0; iZ<dzPV.size(); iZ++)
1935  hPVParamBG[14]->Fill(dzPV[iZ]);
1936  }
1937  continue;
1938  }
1939 
1940  int iMCPV = RtoMCPVId[iPV].GetBestMatch();
1941  KFMCVertex &mcPV = fPrimVertices[iMCPV]; // primary vertex positions (currently only one vertex is implemented)
1942 
1943  int iPVType = 0;
1944  if(mcPV.IsTriggerPV())
1945  {
1946  hPVParamSignal[ 0]->Fill(vtx.X());
1947  hPVParamSignal[ 1]->Fill(vtx.Y());
1948  hPVParamSignal[ 2]->Fill(vtx.Z());
1949  hPVParamSignal[ 3]->Fill(sqrt(vtx.X()*vtx.X() + vtx.Y()*vtx.Y()));
1950  hPVParamSignal[ 4]->Fill(tracks.size());
1951  hPVParamSignal[ 5]->Fill(vtx.Chi2());
1952  hPVParamSignal[ 6]->Fill(vtx.NDF());
1953  hPVParamSignal[ 7]->Fill(vtx.Chi2()/vtx.NDF());
1954  hPVParamSignal[ 8]->Fill(probPV);
1955  hPVParamSignal[ 9]->Fill(fPVPurity[iPV]);
1956  hPVParamSignal[10]->Fill(fPVTracksRate[0][iPV]);
1957  hPVParamSignal[11]->Fill(fPVTracksRate[1][iPV]);
1958  hPVParamSignal[12]->Fill(fPVTracksRate[2][iPV]);
1959  hPVParamSignal[13]->Fill(fPVTracksRate[3][iPV]);
1960  for(unsigned int iZ=0; iZ<dzPV.size(); iZ++)
1961  hPVParamSignal[14]->Fill(dzPV[iZ]);
1962  }
1963  else
1964  {
1965  hPVParamPileup[ 0]->Fill(vtx.X());
1966  hPVParamPileup[ 1]->Fill(vtx.Y());
1967  hPVParamPileup[ 2]->Fill(vtx.Z());
1968  hPVParamPileup[ 3]->Fill(sqrt(vtx.X()*vtx.X() + vtx.Y()*vtx.Y()));
1969  hPVParamPileup[ 4]->Fill(tracks.size());
1970  hPVParamPileup[ 5]->Fill(vtx.Chi2());
1971  hPVParamPileup[ 6]->Fill(vtx.NDF());
1972  hPVParamPileup[ 7]->Fill(vtx.Chi2()/vtx.NDF());
1973  hPVParamPileup[ 8]->Fill(probPV);
1974  hPVParamPileup[ 9]->Fill(fPVPurity[iPV]);
1975  hPVParamPileup[10]->Fill(fPVTracksRate[0][iPV]);
1976  hPVParamPileup[11]->Fill(fPVTracksRate[1][iPV]);
1977  hPVParamPileup[12]->Fill(fPVTracksRate[2][iPV]);
1978  hPVParamPileup[13]->Fill(fPVTracksRate[3][iPV]);
1979  for(unsigned int iZ=0; iZ<dzPV.size(); iZ++)
1980  hPVParamPileup[14]->Fill(dzPV[iZ]);
1981  iPVType = 1;
1982  }
1983  //Find MC parameters of the primary vertex
1984  float mcPVx[3]={mcPV.X(), mcPV.Y(), mcPV.Z()};
1985 
1986  float errPV[3] = {vtx.CovarianceMatrix()[0], vtx.CovarianceMatrix()[2], vtx.CovarianceMatrix()[5]};
1987  for(int iErr=0; iErr<3; iErr++)
1988  if(fabs(errPV[iErr]) < 1.e-8f) errPV[iErr] = 1.e8;
1989 
1990  float dRPVr[3] = {vtx.X()-mcPVx[0],
1991  vtx.Y()-mcPVx[1],
1992  vtx.Z()-mcPVx[2]};
1993  float dRPVp[3] = {static_cast<float>(dRPVr[0]/sqrt(errPV[0])),
1994  static_cast<float>(dRPVr[1]/sqrt(errPV[1])),
1995  static_cast<float>(dRPVr[2]/sqrt(errPV[2]))};
1996 
1997  for(unsigned int iHPV=0; iHPV<3; ++iHPV)
1998  hPVFitQa[iPVType][iHPV]->Fill(dRPVr[iHPV]);
1999  for(unsigned int iHPV=3; iHPV<6; ++iHPV)
2000  hPVFitQa[iPVType][iHPV]->Fill(dRPVp[iHPV-3]);
2001 
2002  for(unsigned int iHPV=0; iHPV<3; ++iHPV)
2003  hPVFitQa2D[iPVType][1][iHPV]->Fill(fNCorrectPVTracks[iPV],dRPVr[iHPV]);
2004  for(unsigned int iHPV=3; iHPV<6; ++iHPV)
2005  hPVFitQa2D[iPVType][1][iHPV]->Fill(fNCorrectPVTracks[iPV],dRPVp[iHPV-3]);
2006 
2007  for(unsigned int iHPV=0; iHPV<3; ++iHPV)
2008  hPVFitQa2D[iPVType][0][iHPV]->Fill(mcPV.NReconstructedDaughterTracks(),dRPVr[iHPV]);
2009  for(unsigned int iHPV=3; iHPV<6; ++iHPV)
2010  hPVFitQa2D[iPVType][0][iHPV]->Fill(mcPV.NReconstructedDaughterTracks(),dRPVp[iHPV-3]);
2011 
2012  hPVFitQa[iPVType][6]->Fill( double(mcPV.NReconstructedDaughterTracks() - fNCorrectPVTracks[iPV])/double(mcPV.NReconstructedDaughterTracks()) );
2013  }
2014 
2015  //fill histograms with quality of input tracks from PV
2016  for(unsigned int iP=0; iP<fTopoReconstructor->GetParticles().size(); iP++)
2017  {
2018  KFParticle TempPart = fTopoReconstructor->GetParticles()[iP];
2019  int nDaughters = TempPart.NDaughters();
2020  if(nDaughters > 1) continue; //use only tracks, not short lived particles
2021 
2022  if(!RtoMCParticleId[iP].IsMatchedWithPdg()) continue; //ghost
2023 
2024  int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
2025  KFMCParticle &mcPart = vMCParticles[iMCPart];
2026 
2027  int iMCTrack = mcPart.GetMCTrackID();
2028  KFMCTrack &mcTrack = vMCTracks[iMCTrack];
2029 
2030  if( mcTrack.MotherId() > -1 ) continue; // select only PV tracks
2031 
2032  const float mcX = mcTrack.X();
2033  const float mcY = mcTrack.Y();
2034  const float mcZ = mcTrack.Z();
2035  const float mcPx = mcTrack.Px();
2036  const float mcPy = mcTrack.Py();
2037  const float mcPz = mcTrack.Pz();
2038 
2039  float decayVtx[3] = {mcX, mcY, mcZ};
2040  TempPart.TransportToPoint(decayVtx);
2041 
2042 
2043  Double_t res[6] = {0},
2044  pull[6] = {0},
2045  mcParam[6] = { mcX, mcY, mcZ,
2046  mcPx, mcPy, mcPz };
2047  for(int iPar=0; iPar < 6; iPar++ )
2048  {
2049  Double_t error = TempPart.GetCovariance(iPar,iPar);
2050  if(error < 0.) { error = 1.e20;}
2051  error = TMath::Sqrt(error);
2052  res[iPar] = TempPart.GetParameter(iPar) - mcParam[iPar];
2053  if(fabs(error) > 1.e-20) pull[iPar] = res[iPar]/error;
2054 
2055  hFitPVTracksQA[iPar]->Fill(res[iPar]);
2056  hFitPVTracksQA[iPar+6]->Fill(pull[iPar]);
2057  }
2058  }
2059 
2060  if(fStoreMCHistograms)
2061  {
2062  for(int iV=0; iV<6; iV++)
2063  for(int iP=0; iP < KFPartEfficiencies::nParticles; iP++)
2064  if(hPartParam[iV][iP][17])
2065  hPartParam[iV][iP][17]->Fill(multiplicities[iV][iP]);
2066  FillMCHistos();
2067  }
2068  else
2069  for(int iP=0; iP < KFPartEfficiencies::nParticles; iP++)
2070  if(hPartParam[0][iP][17])
2071  hPartParam[0][iP][17]->Fill(multiplicities[0][iP]);
2072 } // void KFTopoPerformance::FillHistos()
2073 
2074 void KFTopoPerformance::FillMCHistos()
2075 {
2077  for(unsigned int iMCTrack=0; iMCTrack<vMCTracks.size(); iMCTrack++)
2078  {
2079  int iPDG = fParteff.GetParticleIndex(vMCTracks[iMCTrack].PDG());
2080  if(iPDG < 0) continue;
2081 
2082  if(vMCTracks[iMCTrack].MotherId()>=0) continue;
2083  KFMCParticle &part = vMCParticles[iMCTrack];
2084 
2085  float M = fParteff.partMass[iPDG];
2086  float P = vMCTracks[iMCTrack].P();
2087  float Pt = vMCTracks[iMCTrack].Pt();
2088  float E = sqrt(M*M+P*P);
2089  float Rapidity = 0.5*log((E+vMCTracks[iMCTrack].Pz())/(E-vMCTracks[iMCTrack].Pz()));
2090  float M_t = sqrt(Pt*Pt+M*M)-M;
2091 
2092  float X;
2093  float Y;
2094  float Z;
2095  float R;
2096 
2097  if (part.NDaughters()>0)
2098  {
2099  X = vMCTracks[part.GetDaughterIds()[0]].X();
2100  Y = vMCTracks[part.GetDaughterIds()[0]].Y();
2101  Z = vMCTracks[part.GetDaughterIds()[0]].Z();
2102  }
2103  else
2104  {
2105  X = vMCTracks[iMCTrack].X();
2106  Y = vMCTracks[iMCTrack].Y();
2107  Z = vMCTracks[iMCTrack].Z();
2108  }
2109  R = sqrt(X*X+Y*Y);
2110 
2111 
2112  float parameters[17] = {M, P, Pt, Rapidity, 0, 0, 0, 0, 0, 0, X, Y, Z, R, 0, 0, M_t};
2113  //for all particle-candidates
2114  for(int iParam=0; iParam<17; iParam++)
2115  if(hPartParam[6][iPDG][iParam]) hPartParam[6][iPDG][iParam]->Fill(parameters[iParam]);
2116 
2117  if(hPartParam2D[6][iPDG][0]) hPartParam2D[6][iPDG][0]->Fill(Rapidity,Pt,1);
2118  if(hPartParam2D[6][iPDG][3]) hPartParam2D[6][iPDG][3]->Fill(Rapidity,M_t,1);
2119 
2120  if(IsCollectZRHistogram(iPDG))
2121  if(hPartParam2D[6][iPDG][1]) hPartParam2D[6][iPDG][1]->Fill(Z,R,1);
2122 
2123  if( part.IsReconstructable(2) && IsCollectArmenteros(iPDG))
2124  {
2125  int index1 = part.GetDaughterIds()[0];
2126  int index2 = part.GetDaughterIds()[1];
2127  KFMCTrack positive, negative;
2128  if(vMCTracks[index1].Par(6) > 0)
2129  {
2130  positive = vMCTracks[index1];
2131  negative = vMCTracks[index2];
2132  }
2133  else
2134  {
2135  negative = vMCTracks[index1];
2136  positive = vMCTracks[index2];
2137  }
2138 
2139  float alpha = 0., qt = 0.;
2140  float spx = positive.Px() + negative.Px();
2141  float spy = positive.Py() + negative.Py();
2142  float spz = positive.Pz() + negative.Pz();
2143  float sp = sqrt(spx*spx + spy*spy + spz*spz);
2144  float pn, pln, plp;
2145  pn = sqrt(negative.Px()*negative.Px() + negative.Py()*negative.Py() + negative.Pz()*negative.Pz());
2146  pln = (negative.Px()*spx+negative.Py()*spy+negative.Pz()*spz)/sp;
2147  plp = (positive.Px()*spx+positive.Py()*spy+positive.Pz()*spz)/sp;
2148  float ptm = (1.-((pln/pn)*(pln/pn)));
2149  qt= (ptm>=0.)? pn*sqrt(ptm) :0;
2150  alpha = (plp-pln)/(plp+pln);
2151 
2152  if(hPartParam2D[6][iPDG][2]) hPartParam2D[6][iPDG][2]->Fill(alpha,qt,1);
2153  }
2154  }
2155 }
2156 
2157 void KFTopoPerformance::AddV0Histos()
2158 {
2160  int iV0 = fParteff.nParticles - 1;
2161  int iK0 = fParteff.GetParticleIndex(310);
2162 
2163  for(int iH=0; iH<nFitQA; iH++)
2164  {
2165  hFitDaughtersQA[iV0][iH]->Add(hFitDaughtersQA[iK0][iH]);
2166  hFitQA[iV0][iH]->Add(hFitQA[iK0][iH]);
2167  }
2168 
2169  for(int iV=0; iV<4; iV++)
2170  for(int iH=0; iH<nHistoPartParam; iH++)
2171  hPartParam[iV][iV0][iH]->Add(hPartParam[iV][iK0][iH]);
2172 }
2173 
2174 void KFTopoPerformance::FillHistos(const KFPHistogram* histograms)
2175 {
2177  for(int iParticle=0; iParticle<KFPartEfficiencies::nParticles; iParticle++)
2178  {
2179  const int& nHistograms = histograms->GetHistogramSet(0).GetNHisto1D();
2180  for(int iHistogram=0; iHistogram<nHistograms; iHistogram++)
2181  {
2182  const KFPHistogram1D& histogram = histograms->GetHistogram(iParticle,iHistogram);
2183  for(int iBin=0; iBin<histogram.Size(); iBin++)
2184  hPartParam[0][iParticle][iHistogram]->SetBinContent( iBin, histogram.GetHistogram()[iBin] );
2185  }
2186  }
2187 }
2188 
2189 #endif //DO_TPCCATRACKER_EFF_PERFORMANCE