StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StKFVertex.cxx
1 // $Id: StKFVertex.cxx,v 2.5 2018/04/10 11:32:09 smirnovd Exp $
2 #include "StKFVertex.h"
3 #include "StKFTrack.h"
4 #include "TArrayC.h"
5 
6 #define __DEBUG__
7 #if defined(__DEBUG__)
8 #define PrPP(A,B) if (Debug()) {cout << "StKFVertex::" << (#A) << "\t" << (#B) << " = \t" << (B) << endl;}
9 #else
10 #define PrPP(A,B)
11 #endif
12 using namespace std;
13 Int_t StKFVertex::_debug = 0;
14 const Char_t *StKFVertex::GeNames[52] = {
15  // 1 2 3 4 5 6 7 8 9 10
16  "",
17  "gamma" ,"e+" ,"e-" ,"nu" ,"mu+" ,"mu-" ,"pi0" ,"pi+" ,"pi-" ,"K0L",
18  "K+" ,"K-" ,"N" ,"P" ,"Pbar" ,"K0S" ,"eta" ,"Lambda","Sigma+" ,"Sigma0",
19  "S-" ,"Xi0" ,"Xi-" ,"Omega","Nbar" ,"LamBar","SBar-","SBar0" ,"SBar+" ,"XiBar0",
20  "XiBar+" ,"OmBar","tau+","tau-" ,"D+" ,"D-" ,"D0" ,"Dbar0" ,"Ds+" ,"Ds-" ,
21  "LambC+" ,"W+" ,"W-" ,"Z0" ,"H2" ,"H3" ,"alpha","geanti","He3" ,"Cerenk",
22  "??????"};
23 //________________________________________________________________________________
24 ostream& operator<<(ostream& os, const StKFVertex& v) {
25  Int_t iWest = v.MultW();
26  Int_t iEast = v.MultE();
27  os << Form(" with %4i tracks Q:%3i W/E = %3i/%3i",v.NoTracks(),v.Charge(),iWest,iEast);
28  for (Int_t i = 0; i < 3; i++) {
29  os << Form("%9.3f +/- %5.3f",v.Vertex().GetParameter(i),TMath::Sqrt(v.Vertex().GetCovariance(i,i)));
30  }
31  Double_t prob = TMath::Prob(v.Vertex().GetChi2(),v.Vertex().GetNDF());
32  os << Form("\tchi2/NDF = %8.2f/%4i",v.Vertex().GetChi2(),v.Vertex().GetNDF())
33  << Form(" prob = %6.4f",prob);
34  Double_t M, dM;
35  if (! v.Vertex().GetMass(M,dM)) {
36  os << " M = " << Form("%7.3f +/- %5.3f",M,dM);
37  }
38  Int_t kv = v.IdTruth();
39  if (kv > 0) {
40  os << Form(" Mc/QA/t:%4i/%3i/%6.0f xyz: %8.3f%8.3f%8.3f m:%4i %6s",kv, v.QaTruth(),
41  v.TimeMc(), v.XyzMc().X(), v.XyzMc().Y(), v.XyzMc().Z(),
42  v.NoDaughtersMc(),StKFVertex::StKFVertex::GeNames[v.gePidMc()]);
43  }
44  return os;
45 }
46 //________________________________________________________________________________
47 void StKFVertex::Fit() {
48  Compress();
49  Int_t N = NoTracks();
50  if (! N) return;
51  KFParticle **particles = new KFParticle*[N];
52  for (Int_t i = 0; i < N; i++) {
53  // Track(i)->Reset();
54  particles[i] = &(Track(i)->Particle());
55  }
56  TArrayC Flag(N);
57  Vertex().ConstructPrimaryVertex((const KFParticle **) particles, N,
58  (Bool_t*) Flag.GetArray(),TMath::Sqrt(StAnneling::Chi2Cut()/2));
59  //Check Covariance Matrix
60  // Double_t prob = TMath::Prob(Vertex().GetChi2(),Vertex().GetNDF());
61  TRSymMatrix CL(3,Vertex().CovarianceMatrix());
62  if (CL[0] <= 0 || CL[2] <= 0 || CL[5] <= 0) {
63  for (Int_t i = N-1; i >= 0; i--) delete Remove(i);
64  } else {
65  for (Int_t i = N-1; i >= 0; i--) if (! Flag[i]) delete Remove(i);
66  }
67  delete [] particles;
68  Compress();
69  N = NoTracks();
70  // Assign MC and RC
71  struct vertexPing {
72  Int_t Id;
73  Int_t nPings;
74  };
75  static vertexPing candidates[20];
76  memset(candidates,0,sizeof(candidates));
77  Int_t NC = 0;
78  for (Int_t i = 0; i < N; i++) {
79  const StKFTrack *pTrack = Track(i);
80  if (! pTrack) continue;
81  Int_t IdVx = pTrack->Particle().IdParentVx();
82  if (IdVx <= 0) continue;
83  Int_t J = -1;
84  for (Int_t j = 0; j < NC; j++) if (candidates[j].Id == IdVx) {J = j; break;}
85  if (J < 0) {J = NC; if (NC < 18) NC++;}
86  candidates[J].Id = IdVx;
87  candidates[J].nPings++;
88  }
89  Int_t dominant = -1;
90  Int_t J = -1;
91  for (Int_t j = 0; j < NC; j++) if (candidates[j].nPings > dominant) {dominant = candidates[j].nPings; J = j;}
92  if (J > -1) {
93  Int_t Id = candidates[J].Id;
94  Int_t QA = (100*dominant)/N;
95  SetIdTruth(Id,QA);
96  }
97  for (Int_t i = 0; i < N; i++) {
98  Track(i)->Particle().SetProductionVertex(Vertex());
99  }
100 }
101 //________________________________________________________________________________
102 void StKFVertex::AddTrack(const StKFTrack *track) {
103  Int_t k2 = track->K()%100000;
104  Int_t N1 = NoTracks();
105  for (Int_t j = 0; j < N1; j++) {// protect from multiple copies of beam track
106  StKFTrack* t1 = Track(j);
107  Int_t k1 = t1->K()%100000;
108  if (k1 == k2) {
109  PrintW("AddTrack");
110  assert(0);
111  }
112  }
113  fKFTracks.AddLast((TObject *)track);
114 }
115 //________________________________________________________________________________
116 Double_t StKFVertex::UpdateVertex2TrackChi2() {
117  Int_t Ntracks = NoTracks();
118  Double_t chi2Vx = 0;
119  PrPP(UpdateVertex2TrackChi2,fVertex);
120  if (_debug) PrintW("old Weghts ");
121  for (Int_t k = Ntracks - 1; k >= 0; k--) {
122  KFVertex vTmp = fVertex;
123  // PrPP(UpdateVertex2TrackChi2,vTmp);
124  StKFTrack &track = *Track(k);
125  // PrPP(UpdateVertex2TrackChi2,track.Particle());
126  vTmp -= track.Particle();
127  // PrPP(UpdateVertex2TrackChi2,vTmp);
128  KFParticle *particle = track.OrigParticle();
129  if (! particle) continue;
130  // PrPP(UpdateVertex2TrackChi2,*particle);
131  Double_t chi2il = particle->GetDeviationFromVertex(vTmp);
132  chi2il *= 2*chi2il;
133  if (chi2il > StAnneling::Chi2Cut()) {
134  Remove(k);
135  continue;
136  }
137  track.SetChi2(chi2il);
138  chi2Vx += track.Chi2()/2 + TMath::Log(track.Weight() + StAnneling::Weight());
139  }
140  Compress();
141  if (_debug) PrintW("new Weights ");
142  return chi2Vx;
143 }
144 //________________________________________________________________________________
145 StKFTrack* StKFVertex::Remove(KFParticle *particle) {
146  Int_t N = NoTracks();
147  for (Int_t k = 0; k < N; k++) if (particle == Track(k)->OrigParticle()) return Remove(k);
148  return 0;
149 }
150 //________________________________________________________________________________
151 Int_t StKFVertex::MultWE(Int_t k) const {
152  Int_t N = NoTracks();
153  Int_t iWE = 0;
154  for (Int_t i = 0; i < N; i++) {
155  const StKFTrack* t = Track(i);
156  if (t) {
157  Int_t id = (t->OrigParticle()->GetID()/100000)%10;
158  if (id == k) iWE++;
159  }
160  }
161  return iWE;
162 }
163 //________________________________________________________________________________
164 Int_t StKFVertex::Q() const {
165  Int_t iQ = 0;
166  Int_t N = NoTracks();
167  for (Int_t i = 0; i < N; i++) {
168  const StKFTrack* t = Track(i);
169  if (t) {iQ += t->OrigParticle()->GetQ();}
170  }
171  return iQ;
172 }
173 //________________________________________________________________________________
174 void StKFVertex::operator +=(StKFVertex &vtx) {
175  if (_debug) {
176  PrintW("Before Merge 1");
177  vtx.PrintW("Before Merge 2");
178  }
179  Int_t N2 = vtx.NoTracks();
180  for (Int_t i = N2-1; i >= 0; i--) {
181  StKFTrack* t2 = (StKFTrack* ) vtx.Remove(i);
182  Int_t k2 = t2->K()%100000;
183  Int_t N1 = NoTracks();
184  for (Int_t j = 0; j < N1; j++) {// protect from multiple copies of beam track
185  StKFTrack* t1 = Track(j);
186  Int_t k1 = t1->K()%100000;
187  if (k1 == k2) {SafeDelete(t2); break;}
188  }
189  if (t2) fKFTracks.AddLast(t2);
190  }
191  vtx.Compress();
192  if (_debug) {
193  PrintW("After Merge 1");
194  vtx.PrintW("After Merge 2");
195  }
196 }
197 //________________________________________________________________________________
198 void StKFVertex::PrintW(Option_t *option) const {
199  Int_t N = NoTracks();
200  cout << Form("Vertex %5i with %5i tracks\t",fID,N);
201  Print(option);
202  cout << Form(" i k Weight W Z chi2") << endl;
203  for (Int_t i = 0; i < N; i++) {
204  const StKFTrack* t = Track(i);
205  cout << Form("%6i",i) << *t << endl;
206  }
207 }
208 //________________________________________________________________________________
209 void StKFVertex::SetMc(Float_t time, Float_t x, Float_t y, Float_t z, Int_t NoDaughters, Int_t gePid) {
210  fTimeMc = 1e9*time;
211  fXyzMc = TVector3(x,y,z);
212  fNoDaughtersMc = NoDaughters;
213  fgePidMc = StKFTrack::CorrectGePid(gePid);
214 }
215 #undef PrPP
216 // $Log: StKFVertex.cxx,v $
217 // Revision 2.5 2018/04/10 11:32:09 smirnovd
218 // Minor corrections across multiple files
219 //
220 // - Remove ClassImp macro
221 // - Change white space
222 // - Correct windows newlines to unix
223 // - Remove unused debugging
224 // - Correct StTpcRTSHitMaker header guard
225 // - Remove unused preprocessor directives in StiCA
226 // - Minor changes in status and debug print out
227 // - Remove using std namespace from StiKalmanTrackFinder
228 // - Remove includes for unused headers
229 //
230 // Revision 2.4 2013/04/10 22:14:20 fisyak
231 // Roll back to version 04/04/2013
232 //
233 // Revision 2.2 2012/06/11 15:33:41 fisyak
234 // std namespace
235 //
236 // Revision 2.1 2012/05/07 14:56:14 fisyak
237 // Add StKFVertexMaker
238 //
239 // Revision 1.3 2012/03/29 23:35:47 fisyak
240 // Fix problem with multiple beam tracks
241 //
242 // Revision 1.2 2012/02/07 19:38:26 fisyak
243 // Repackage
244 //
Track
C++ STL includes.
Definition: AgUStep.h:47
Generated by   doxygen 1.8.5