StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StKFVerticesCollection.cxx
1 // $Id: StKFVerticesCollection.cxx,v 2.5 2018/04/10 11:32:10 smirnovd Exp $
2 #include "StKFVerticesCollection.h"
3 #include "TArrayI.h"
4 #include "TArrayD.h"
5 #include "TRSymMatrix.h"
6 #include "TRVector.h"
7 #include "TPolyMarker.h"
8 #include "TList.h"
9 using namespace std;
10 Double_t StKFVerticesCollection::fgVxPenaltyFactor = 1000;
11 //________________________________________________________________________________
12 StKFVerticesCollection::StKFVerticesCollection(Int_t NoPeaks, Double_t *zOfPeaks, Double_t sigmaXY, Double_t sigmaZ) :
13  fVertices(NoPeaks,0) {
14  fVertices.SetOwner(kTRUE);
15  for (Int_t peak = 0; peak < NoPeaks; peak++) {
16  AddVertex(0.,0., zOfPeaks[peak], sigmaXY, sigmaZ);
17  }
18 }
19 //________________________________________________________________________________
20 void StKFVerticesCollection::AddVertex(Double_t x, Double_t y, Double_t z, Double_t sigmaXY, Double_t sigmaZ) {
21  StKFVertex *vtx = new StKFVertex(fVertices.GetEntriesFast() + 1);
22  vtx->Vertex().SetBeamConstraint(x, y, z, sigmaXY, sigmaXY, sigmaZ);
23  fVertices.AddLast(vtx);
24 }
25 //________________________________________________________________________________
26 void StKFVerticesCollection::operator +=(StKFVerticesCollection &col) {
27  Int_t N2 = col.NoVertices();
28  for (Int_t i = N2-1; i >= 0; i--) {
29  fVertices.AddLast(col.Remove(i));
30  }
31  col.Compress();
32 }
33 //________________________________________________________________________________
34 void StKFVerticesCollection::SetMc(Int_t NoMuMcVertex, Int_t NoMuMcTrack, const Float_t *time,
35  const Float_t *x,const Float_t *y,const Float_t *z,
36  const Int_t *NoDaughters,const Int_t *IdParTrk,const Int_t *gePid) {
37  Int_t Nvtx = NoVertices();
38  for (Int_t l = 0; l < Nvtx; l++) {
39  StKFVertex *V = Vertex(l);
40  if (! V) continue;
41  Int_t kv = V->IdTruth();
42  if (kv > 0 && kv <= NoMuMcVertex) {
43  if (time && x && y && z && NoDaughters && IdParTrk) {
44  Int_t kvp = IdParTrk[kv-1];
45  Int_t ge = 0;
46  if (kvp > 0 && kvp <= NoMuMcTrack) ge = gePid[kvp-1];
47  V->SetMc(time[kv-1],x[kv-1],y[kv-1],z[kv-1],NoDaughters[kv-1],ge);
48  }
49  }
50  }
51 }
52 //________________________________________________________________________________
53 ostream& operator<<(ostream& os, const StKFVerticesCollection& vc) {
54  Int_t Nvtx = vc.NoVertices();
55  for (Int_t l = 0; l < Nvtx; l++) {
56  const StKFVertex *V = vc.Vertex(l);
57  if (! V) continue;
58  os << Form("Vtx: %3i",l) << *V << endl;
59  }
60  return os;
61 }
62 //________________________________________________________________________________
63 Double_t StKFVerticesCollection::DoTrack2VertexAssociation(const TObjArray &particles) {
64  Double_t chi2Total = 0;
65  Int_t NVtx = NoVertices();
66  Int_t Ntracks = particles.GetEntriesFast();
67  TArrayI Idx(Ntracks);
68  TArrayD Chi2s(Ntracks);
69  for (Int_t l = 0; l < NVtx; l++) {
70  StKFVertex *vtx = Vertex(l);
71  if (! vtx) continue;
72  // if (! vtx->Vertex()) continue;
73  // vtx->PrintW();
74  KFParticle *particle = 0;
75  for (Int_t k = 0; k < Ntracks; k++) {
76  Chi2s[k] = 1e10;
77  particle = (KFParticle *) particles.UncheckedAt(k);
78  if (! particle) continue;
79  if (particle->GetID() > 100000) continue;
80  Double_t chi2il = particle->GetDeviationFromVertex(vtx->Vertex());
81  chi2il *= 2*chi2il;
82  Chi2s[k] = chi2il;
83  }
84  vtx->Clear();
85  TMath::Sort(Ntracks,Chi2s.GetArray(),Idx.GetArray(),0);
86  Double_t chi2Vx = 0;
87  for (Int_t j = 0; j < Ntracks; j++) {
88  Int_t k = Idx[j];
89  particle = (KFParticle *) particles.UncheckedAt(k);
90  if (! particle) continue;
91  if (Chi2s[k] > StAnneling::Chi2Cut()) break;
92  StKFTrack *track = new StKFTrack(k,particle,Chi2s[k]);
93  chi2Vx += track->Chi2()/2 + TMath::Log(track->Weight() + StAnneling::Weight());
94  vtx->AddTrack(track);
95  }
96  if (vtx->NoTracks() < 2) {
97  delete vtx; Vertex(l) = 0;
98  continue;
99  }
100  if (StKFVertex::Debug()) vtx->PrintW("DoTrack2VertexAssociation ");
101  chi2Total += chi2Vx;
102  // vtx->PrintW();
103  }
104  fVertices.Compress();
105  if (NoVertices()) UpdateWeights();
106  return chi2Total;
107 }
108 //________________________________________________________________________________
109 Double_t StKFVerticesCollection::UpdateStVertexTrackAssociation() {
110  Double_t chi2Total = 0;
111  Int_t NVtx = NoVertices();
112  for (Int_t l = 0; l < NVtx; l++) {
113  StKFVertex *vtx = Vertex(l);
114  if (! vtx) continue;
115  chi2Total += vtx->UpdateVertex2TrackChi2();
116  }
117  return chi2Total;
118 }
119 //________________________________________________________________________________
120 void StKFVerticesCollection::CleanDuplicatedVertices() {
121  Int_t NVtx = NoVertices();
122  // Check that vertices are the same
123  for (Int_t l = 1; l < NVtx; l++) {
124  StKFVertex *vtxl = Vertex(l);
125  if (! vtxl) continue;
126  Int_t NL = vtxl->NoTracks();
127  TRVector vL(3,vtxl->Vertex().Parameters());
128  TRSymMatrix CL(3,vtxl->Vertex().CovarianceMatrix());
129  for (Int_t m = 0; m < l; m++) {
130  StKFVertex *vtxm = Vertex(m);
131  if (! vtxm) continue;
132  // compliete overlap by an other vertex
133  Int_t NM = vtxm->NoTracks();
134  Int_t Nmatched = 0;
135  for (Int_t i = 0; i < NL; i++)
136  for (Int_t j = 0; j < NM; j++)
137  if (vtxl->Track(i)->OrigParticle() == vtxm->Track(j)->OrigParticle()) Nmatched++;
138  if (Nmatched == TMath::Min(NL,NM)) {
139  TRVector vM(3,vtxm->Vertex().Parameters());
140  TRSymMatrix CM(3,vtxm->Vertex().CovarianceMatrix());
141  vM -= vL;
142  CM += CL;
143  TRSymMatrix G(CM,TRArray::kInverted);
144  Double_t chi2 = G.Product(vM,TRArray::kATxSxA);
145  Double_t prob = TMath::Prob(chi2,3);
146  if (prob > 0.10) {
147  if ((NL > NM) || ((NL == NM) && (vtxl->Vertex().GetChi2() < vtxm->Vertex().GetChi2()))) {
148  if (StKFVertex::Debug()) {
149  vtxm->Print(Form("Cleaned Vertex prob %7.2f M %3i keep L %3i L(%3i/%7.2f) M (%3i/%7.2f)\t",
150  prob,m,l,NL,vtxl->Vertex().GetChi2(),NM,vtxm->Vertex().GetChi2()));
151  }
152  delete vtxm; Vertex(m) = 0;
153  continue;
154  } else {
155  if (StKFVertex::Debug()) {
156  vtxl->Print(Form("Cleaned Vertex prob %7.2f L %3i keep M %3i M(%3i/%7.2f) L (%3i/%7.2f)",
157  prob,l,m,NM,vtxm->Vertex().GetChi2(),NL,vtxl->Vertex().GetChi2()));
158  }
159  delete vtxl; Vertex(l) = 0;
160  break;
161  }
162  }
163  }
164  }
165  }
166  fVertices.Compress();
167 }
168 //________________________________________________________________________________
169 void StKFVerticesCollection::MergeDuplicatedVertices() {
170  // Check that vertices are the same
171  for (Int_t l = 1; l < NoVertices(); l++) {
172  StKFVertex *vtxl = Vertex(l);
173  if (! vtxl) continue;
174  if (! vtxl->NoTracks()) {delete vtxl; vtxl = Vertex(l) = 0; continue;}
175  TRVector vL(3,vtxl->Vertex().Parameters());
176  TRSymMatrix CL(3,vtxl->Vertex().CovarianceMatrix());
177  for (Int_t m = 0; m < l; m++) {
178  StKFVertex *vtxm = Vertex(m);
179  if (! vtxm) continue;
180  if (! vtxm->NoTracks()) {delete vtxm; vtxm = Vertex(m) = 0; continue;}
181  TRVector vM(3,vtxm->Vertex().Parameters());
182  vM -= vL;
183  if (vM.Mag() > 5.) continue;
184  TRSymMatrix CM(3,vtxm->Vertex().CovarianceMatrix());
185  CM += CL;
186  CM[0] += 0.0001; // 100 mkm tolerance
187  CM[2] += 0.0001;
188  CM[5] += 0.0001;
189  TRSymMatrix G(CM,TRArray::kInverted);
190  Double_t chi2 = G.Product(vM,TRArray::kATxSxA);
191  Double_t prob = TMath::Prob(chi2,3);
192  if (prob > 1e-4) {
193  *vtxm += *vtxl;
194  // if (vtxl->NoTracks()) vtxl->Clear("keep");
195  delete vtxl; vtxl = Vertex(l) = 0;
196  vtxm->Fit();
197  break;
198  }
199  }
200  }
201  fVertices.Compress();
202 }
203 //________________________________________________________________________________
204 void StKFVerticesCollection::UpdateWeights() {
205  Int_t NVtx = NoVertices();
206  for (Int_t l = 0; l < NVtx; l++) {
207  StKFVertex *vtxl = Vertex(l);
208  if (! vtxl) continue;
209  // recalculate weights
210  if (StKFVertex::Debug()) vtxl->PrintW("Weights to Update ");
211  for (Int_t i = 0; i < vtxl->NoTracks(); i++) {
212  Double_t Dominator = TMath::Exp(-StAnneling::Chi2Cut()/(2*StAnneling::Temperature())) + vtxl->Track(i)->Weight();
213  for (Int_t m = 0; m < NVtx; m++) {
214  if (l == m) continue;
215  StKFVertex *vtxm = Vertex(m);
216  if (! vtxm) continue;
217  for (Int_t j = 0; j < vtxm->NoTracks(); j++) {
218  if (vtxl->Track(i)->OrigParticle() == vtxm->Track(j)->OrigParticle()) {
219  Dominator += vtxm->Track(j)->Weight();
220  break;
221  }
222  }
223  }
224  vtxl->Track(i)->W() = vtxl->Track(i)->Weight()/Dominator;
225  vtxl->Track(i)->Particle() = *(vtxl->Track(i)->OrigParticle());
226  Double_t *CovXyz = vtxl->Track(i)->Particle().CovarianceMatrix();
227  for (Int_t j = 0; j < 36; j++) CovXyz[j] = CovXyz[j]/vtxl->Track(i)->W();
228  }
229  if (StKFVertex::Debug()) vtxl->PrintW("Updated Weights ");
230  }
231 }
232 //________________________________________________________________________________
233 void StKFVerticesCollection::UniqueTracks2VertexAssociation(){
234  // Make track associated with only vertex (by maximum weight to the vertex)
235  Int_t NVtx = NoVertices();
236  for (Int_t l = 0; l < NVtx; l++) {
237  StKFVertex *vtxl = Vertex(l);
238  if (! vtxl) continue;
239  // recalculate weights
240  for (Int_t i = vtxl->NoTracks()-1; i >= 0; i--) {
241  if (! vtxl->Track(i)) continue;
242  Double_t WMax = vtxl->Track(i)->Weight();
243  Int_t iMax = i;
244  Int_t lMax = l;
245  Int_t nPart = 1;
246  KFParticle *particleMax = vtxl->Track(i)->OrigParticle();
247  for (Int_t m = 0; m < NVtx; m++) {
248  if (l == m) continue;
249  StKFVertex *vtxm = Vertex(m);
250  if (! vtxm) continue;
251  for (Int_t j = 0; j < vtxm->NoTracks(); j++) {
252  if (particleMax == vtxm->Track(j)->OrigParticle()) {
253  nPart++;
254  if (vtxm->Track(j)->Weight() > WMax) {
255  WMax = vtxm->Track(j)->Weight();
256  iMax = j;
257  lMax = m;
258  particleMax = vtxm->Track(j)->OrigParticle();
259  break;
260  }
261  }
262  }
263  }
264  if (WMax < 0.01) { // Particle weight is too small => remove the particle from all vertices
265  for (Int_t m = 0; m < NVtx; m++) {
266  StKFVertex *vtxm = Vertex(m);
267  if (! vtxm) continue;
268  delete vtxm->Remove(particleMax);
269  vtxm->Compress();
270  }
271  if (! vtxl->NoTracks()) break;
272  continue;
273  }
274  if (nPart > 1) {
275  for (Int_t m = 0; m < NVtx; m++) {
276  StKFVertex *vtxm = Vertex(m);
277  if (! vtxm) continue;
278  if (m != lMax) {
279  if (particleMax->GetID()%100000) { // beam track is not in the game
280  delete vtxm->Remove(particleMax);
281  vtxm->Compress();
282  if (vtxm->NoTracks() == 0) {delete vtxm; Vertex(m) = 0;}
283  }
284  } else vtxm->Track(iMax)->W() = vtxm->Track(iMax)->Weight();
285  }
286  }
287  }
288  }
289  for (Int_t l = 0; l < NVtx; l++) {
290  StKFVertex *vtxl = Vertex(l);
291  if (! vtxl) continue;
292  if (vtxl->NoTracks() == 0) {delete vtxl; Vertex(l) = 0;}
293  }
294  fVertices.Compress();
295  NVtx = NoVertices();
296  // Set particle ID
297  for (Int_t l = 0; l < NVtx; l++) {
298  StKFVertex *vtxl = Vertex(l);
299  if (! vtxl) continue;
300  Int_t N = vtxl->NoTracks();
301  for (Int_t i = 0; i < N; i++) {
302  KFParticle *particle = vtxl->Track(i)->OrigParticle();;
303  Int_t ID = particle->GetID()%100000 + 100000*vtxl->ID();;
304  particle->SetID(ID);
305  }
306  }
307 }
308 //________________________________________________________________________________
309 Double_t StKFVerticesCollection::Fit(Int_t marker, TCanvas *c1, TH1 *Vtx) {
310  // Primary Vertex fit
311  Double_t chi2Total = 1e10;
312  fVertices.Compress();
313  Int_t NVtx = NoVertices();
314  if (! NVtx) return chi2Total;
315  for (Int_t l = 0; l < NVtx; l++) {
316  StKFVertex *vtx = Vertex(l);
317  if (! vtx) continue;
318  vtx->Vertex().SetBeamConstraintOff();
319  vtx->Fit();
320  if (vtx->Vertex().GetNDF() < 1 || vtx->NoTracks() < 2) {
321  delete vtx;
322  Vertex(l) = 0;
323  continue;
324  } else {
325  chi2Total += fgVxPenaltyFactor;
326  }
327  }
328  fVertices.Compress();
329  CleanDuplicatedVertices();
330  fVertices.Compress();
331  chi2Total = UpdateStVertexTrackAssociation();
332  UpdateWeights();
333  //
334  if (StKFVertex::Debug()) {
335  cout << "chi2Total = " << chi2Total
336  << " at Temperature " << StAnneling::Temperature()
337  << " and Log(Temperature) " << TMath::Log(StAnneling::Temperature())
338  << " no. vertices " << NoVertices()
339  << endl;
340  }
341  if (Vtx) {
342  Double_t ymax = Vtx->GetMaximum();
343  for (Int_t i = 0; i < NVtx; i++) {
344  StKFVertex *vtx = Vertex(i);
345  if (! vtx) continue;
346  Double_t X = vtx->Vertex().GetParameter(2);
347  Double_t Y = vtx->NoTracks();
348  if (Y > ymax) ymax = Y;
349  TPolyMarker * pm = new TPolyMarker(1, &X, &Y);
350  Vtx->GetListOfFunctions()->Add(pm);
351  pm->SetMarkerColor(TMath::Log(StAnneling::Temperature())+2);
352  Int_t m = 22;
353  if (marker) {
354  m = marker;
355  pm->SetMarkerColor(4);
356  if (StKFVertex::Debug()) vtx->PrintW();
357  }
358  pm->SetMarkerStyle(m);
359  pm->SetMarkerSize(2);
360  // chi2NDF->Fill(TMath::Log10(vtx->NoTracks()),vtx->Vertex().GetChi2()/vtx->Vertex().GetNDF());
361  }
362  if (c1) {
363  Vtx->SetMaximum(ymax);
364  Vtx->Draw("same");
365  c1->Update();
366  }
367  }
368  return chi2Total;
369 }
370 // $Log: StKFVerticesCollection.cxx,v $
371 // Revision 2.5 2018/04/10 11:32:10 smirnovd
372 // Minor corrections across multiple files
373 //
374 // - Remove ClassImp macro
375 // - Change white space
376 // - Correct windows newlines to unix
377 // - Remove unused debugging
378 // - Correct StTpcRTSHitMaker header guard
379 // - Remove unused preprocessor directives in StiCA
380 // - Minor changes in status and debug print out
381 // - Remove using std namespace from StiKalmanTrackFinder
382 // - Remove includes for unused headers
383 //
384 // Revision 2.4 2013/04/10 22:14:20 fisyak
385 // Roll back to version 04/04/2013
386 //
387 // Revision 2.2 2012/06/11 15:33:41 fisyak
388 // std namespace
389 //
390 // Revision 2.1 2012/05/07 14:56:14 fisyak
391 // Add StKFVertexMaker
392 //
393 // Revision 1.5 2012/03/29 23:35:47 fisyak
394 // Fix problem with multiple beam tracks
395 //
396 // Revision 1.4 2012/03/26 23:42:35 fisyak
397 // Add beam constrain
398 //
399 // Revision 1.3 2012/02/07 19:38:26 fisyak
400 // Repackage
401 //
Generated by   doxygen 1.8.5