StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StKFVertexMaker.cxx
1 // $Id: StKFVertexMaker.cxx,v 2.8 2018/04/10 11:32:09 smirnovd Exp $
2 #include "RVersion.h"
3 #if ROOT_VERSION_CODE < 331013
4 #include "TCL.h"
5 #else
6 #include "TCernLib.h"
7 #endif
8 #include <map>
9 using std::map;
10 #include "TMath.h"
11 #include "TH1.h"
12 #include "TCanvas.h"
13 #include "StDcaGeometry.h"
14 #include "KFParticle.h"
15 #include "KFVertex.h"
16 #include "MTrack.h"
17 #include "VVertex.h"
18 #include "TH1K.h"
19 #include "StAnneling.h"
20 #include "StKFEvent.h"
21 #include "StKFTrack.h"
22 #include "StKFVertex.h"
23 #include "StKFVerticesCollection.h"
24 #include "StVertexP.h"
25 #include "StVertexT.h"
26 #include "TDirectory.h"
27 #include "StEventTypes.h"
28 #include "Stypes.h"
29 #include "SystemOfUnits.h"
30 #include "StKFVertexMaker.h"
31 #include "StDetectorDbMaker/St_vertexSeedC.h"
32 #include "Sti/StiHit.h"
33 #include "Sti/StiKalmanTrack.h"
34 #include "Sti/StiKalmanTrackNode.h"
35 #include "StiStEventFiller.h"
36 #include "TRMatrix.h"
37 #include "TRSymMatrix.h"
38 #include "Sti/StiToolkit.h"
39 #include "TArrayI.h"
40 
41 StKFVerticesCollection *StKFVertexMaker::fcVertices = 0;
42 //________________________________________________________________________________
43 StKFVertexMaker::~StKFVertexMaker() {
44  SafeDelete(fVtxM);
45  for (Int_t pass = 0; pass < fNPasses; pass++) {
46  SafeDelete(fVtxKs[pass]);
47  SafeDelete(fVtxKs[pass]);
48  }
49  SafeDelete(fSpectrum);
50  SafeDelete(func);
51  SafeDelete(fminBrent);
52  delete [] fVerticesPass; fVerticesPass = 0;
53  SafeDelete(fParticles);
54 }
55 //________________________________________________________________________________
56 void StKFVertexMaker::Clear(Option_t *option) {
57  for (Int_t pass = 0; pass < fNPasses; pass++) {
58  fVtxKs[pass]->Reset();
59  fVtxKs[pass]->SetMaximum();
60  fVtxs[pass]->Reset();
61  fVtxs[pass]->SetMaximum();
62  }
63  fVtx = fVtxs[0]; // << switch between types Vtx = fVtxKs[0];
64  fVtxM->Reset();
65  fcVertices = 0;
66  fParticles->Clear("C");
67 }
68 //________________________________________________________________________________
69 StKFVertexMaker::StKFVertexMaker(const char *name) : StMaker(name),
70  fNzBins(2500), fNPasses(2), fSpectrum(0), fzWindow(2),
71  fVtxM(0), fVerticesPass(0), fTempLog(2), fminBrent(0), func(0),
72  mBeamLine(kFALSE), fc1(0)
73 {
74  Int_t npeaks = 100;
75  Double_t zmin = -250;
76  Double_t zmax = 250;
77  // StKFVertex::_debug = 1;
78  for (Int_t pass = 0; pass < fNPasses; pass++) {
79  fVtxs[pass] = new TH1F(Form("Vtx%1i",pass),Form("z-dca distribution for pass = %1i",pass),fNzBins,zmin,zmax);
80  fVtxs[pass]->SetDirectory(0);
81  if (pass) fVtxs[pass]->SetLineColor(5);
82  fVtxs[pass]->SetDefaultSumw2();
83  fVtxs[pass]->SetStats(0);
84  fVtxKs[pass] = new TH1K(Form("VtxK%1i",pass),Form("z-dca distribution for pass = %1i",pass),fNzBins,zmin,zmax);
85  fVtxKs[pass]->SetDirectory(0);
86  fVtxKs[pass]->SetStats(0);
87  fVtxKs[pass]->SetLineColor(2);
88  }
89  fVtxM = new TH1F("VtxM","MuDst reconstructed multiplicities versus Z",fNzBins,zmin,zmax);
90  fVtxM->SetDirectory(0);
91  fSpectrum = new TSpectrum(2*npeaks);
92  func = new ROOT::Math::Functor1D(&StKFVertexMaker::AnnelingFcn);
93  fminBrent = new ROOT::Math::GSLMinimizer1D();
94  fVerticesPass = new StKFVerticesCollection *[fNPasses+1];
95  memset (fVerticesPass, 0, (fNPasses+1)*sizeof(StKFVerticesCollection *));
96  fParticles = new TObjArray();
97  fParticles->SetOwner(kTRUE);
98  mVertexOrderMethod = orderByRanking; // change ordering by ranking
99 }
100 //_____________________________________________________________________________
101 Int_t StKFVertexMaker::Init(){
102  mBeamLine = IAttr("beamLine");
103  return StMaker::Init();
104 }
105 //_____________________________________________________________________________
106 Int_t StKFVertexMaker::InitRun(Int_t runumber){
107  return StMaker::InitRun(runumber);
108 }
109 //________________________________________________________________________________
110 Int_t StKFVertexMaker::Make() {
111  StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
112  if (! pEvent) {
113  LOG_WARN << "StKFVertexMaker::fit: no StEvent " << endm;
114  return kStOK; // if no event, we're done
115  }
116  Double_t bField = 0;
117  if (pEvent->runInfo()) bField = pEvent->runInfo()->magneticField();
118  KFParticle::SetField(bField);
119  if (mBeamLine) {
120  St_vertexSeedC* vSeed = St_vertexSeedC::instance();
121  Double_t x0 = vSeed->x0() ; Double_t err_x0 = vSeed->err_x0();
122  Double_t y0 = vSeed->y0() ; Double_t err_y0 = vSeed->err_y0();
123  Double_t dxdz = vSeed->dxdz(); Double_t err_dxdz = vSeed->err_dxdz();
124  Double_t dydz = vSeed->dydz(); Double_t err_dydz = vSeed->err_dydz();
125  Double_t weight = vSeed->weight();
126  if (err_x0 < 0.010) err_x0 = 0.010;
127  if (err_y0 < 0.010) err_y0 = 0.010;
128  static Bool_t firstTime = kTRUE;
129  if (firstTime) {
130  firstTime = kFALSE;
131  LOG_INFO << "BeamLine Constraint: weight = " << weight << endm;
132  LOG_INFO << "x(z) = (" << x0 << " +- " << err_x0 << ") + (" << dxdz << " +- " << err_dxdz << ") * z" << endm;
133  LOG_INFO << "y(z) = (" << y0 << " +- " << err_y0 << ") + (" << dydz << " +- " << err_dydz << ") * z" << endm;
134  }
135  static Double_t pZ = 1000;
136  static MTrack track;
137  Double_t xyzP[6] = { x0, y0, 0.,
138  pZ*dxdz, pZ*dydz, pZ};
139  Double_t CovXyzP[21] = {
140  err_x0*err_x0,
141  0 ,err_y0*err_y0,
142  0 ,0 , 0,
143  0 ,0 , 0, (err_dxdz*pZ)*(err_dxdz*pZ),
144  0 ,0 , 0, 0, (err_dydz*pZ)*(err_dydz*pZ)
145  };
146  track.SetParameters(xyzP);
147  track.SetCovarianceMatrix(CovXyzP);
148  track.SetNDF(1);
149  track.SetID(0);
150  track.SetCharge(1);
151  KFParticle *beam = new KFParticle(track, 321);
152  fParticles->AddAt(beam, 0);
153  }
154  StSPtrVecTrackNode& trackNode = pEvent->trackNodes();
155  UInt_t nTracks = trackNode.size();
156  StTrackNode *node=0;
157  Int_t NGoodGlobals = 0;
158  map<Int_t,StTrackNode*> TrackNodeMap;
159  for (UInt_t i=0; i < nTracks; i++) {
160  node = trackNode[i];
161  if (!node) continue;
162  StGlobalTrack *gTrack = static_cast<StGlobalTrack *>(node->track(global));
163  if (! gTrack) continue;
164  const StDcaGeometry* dca = gTrack->dcaGeometry();
165  if (! dca) continue;
166  if (gTrack->flag() < 0) continue; // Bad fit
167  if (gTrack->flag() > 700) continue; // FTPC
168  if (gTrack->flag()%100 == 11) continue; // Short track pointing to EEMC
169  if ((gTrack->isWestTpcOnly() || gTrack->isEastTpcOnly()) && gTrack->isPostXTrack()) continue; // wrong TPC side track
170  Int_t kg = gTrack->key();
171  TrackNodeMap[kg] = node;
172  KFParticle *particle = AddTrackAt(dca,kg);
173  if (Debug() > 1) {
174  if (Debug() > 2) {LOG_INFO << Form("particle: %4i/%4i ",NGoodGlobals,kg) << *particle << endm;}
175  LOG_INFO << "Add to map[" << kg << "] node = " << TrackNodeMap[kg] << endm;
176  }
177  NGoodGlobals++;
178  }
179  if (NGoodGlobals < 2) return 0;
180  Fit();
181  if (! Vertices()) return 0;
182  //
183  // In case there are no tracks left we better quit
184  //
185  StSPtrVecTrackDetectorInfo& detInfoVec = pEvent->trackDetectorInfo();
186  Int_t Nvtx = Vertices()->NoVertices();
187  for (Int_t l = 0; l < Nvtx; l++) {
188  const StKFVertex *V = Vertices()->Vertex(l);
189  if (! V) continue;
190  if (Debug() > 2) V->PrintW();
191  // Store vertex
192  StPrimaryVertex *primV = new StPrimaryVertex;
193  StThreeVectorF XVertex(&V->Vertex().X());
194  primV->setPosition(XVertex);
195  primV->setChiSquared(V->Vertex().Chi2()/V->Vertex().GetNDF());
196  primV->setProbChiSquared(TMath::Prob(V->Vertex().GetChi2(),V->Vertex().GetNDF()));
197  Float_t cov[6];
198  TCL::ucopy(&V->Vertex().Covariance(0),cov,6);
199  primV->setCovariantMatrix(cov);
200  primV->setVertexFinderId(KFVertexFinder);
201  primV->setFlag(1); // was not set earlier by this vertex finder ?? Jan
202  primV->setRanking(333);
203  primV->setNumTracksUsedInFinder(V->NoTracks());
204  Bool_t beam = kFALSE;
205  Double_t Pars[6];
206  TCL::ucopy(&V->Vertex().X(), Pars, 6);
207  Double_t Cov[21];
208  TCL::ucopy(&V->Vertex().Covariance(0), Cov, 21);
209  StiHit *Vertex = StiToolkit::instance()->getHitFactory()->getInstance();
210  Vertex->setGlobal(0, 0, V->Vertex().X(), V->Vertex().Y(), V->Vertex().Z(), 0);
211  Vertex->setError(cov);
212  Int_t NoTracks = V->NoTracks();
213  TArrayI indexT(NoTracks); Int_t *indexes = indexT.GetArray();
214  TArrayI IdT(NoTracks); Int_t *Ids = IdT.GetArray();
215  for (Int_t itk = 0; itk < NoTracks; itk++) {
216  Ids[itk] = 999999;
217  const StKFTrack* track = V->Track(itk);
218  if (! track) continue;
219  const KFParticle &P = track->Particle();
220  Int_t kg = P.GetID()%100000;
221  Ids[itk] = kg;
222  }
223  TMath::Sort(NoTracks,Ids,indexes,0);
224  for (Int_t i = 0; i < NoTracks; i++) {
225  Int_t itk = indexes[i];
226  const StKFTrack* track = V->Track(itk);
227  if (! track) continue;
228  const KFParticle &P = track->Particle();
229  Int_t kg = P.GetID()%100000;
230  if (kg == 0) {
231  assert(!beam);
232  beam = kTRUE;
233  continue;
234  }
235  if (Debug() > 2) {
236  const KFParticle *PO = track->OrigParticle();
237  const KFParticle *PS[2] = {PO, &P};
238  for (Int_t m = 0; m < 2; m++) {
239  if (! m) cout << "Original";
240  else cout << "Fitted ";
241  static const Char_t *names[6] = {"x","y","z","px","py","pz"};
242  for (Int_t j = 0; j < 6; j++) {
243  cout << Form(" %2s: %8.3f +/- %8.3f",names[j],
244  PS[m]->GetParameter(j),
245  PS[m]->GetCovariance(j,j) > 0 ? TMath::Sqrt(PS[m]->GetCovariance(j,j)) : -13);
246  }
247  cout << endl;
248  }
249  }
250  node = TrackNodeMap[kg];
251  if (! node) {
252  LOG_INFO << "Missing node in map[" << kg << "] node = " << TrackNodeMap[kg] << endm;
253  assert(node);
254  }
255  StiKalmanTrack* kTrack = (*StiStEventFiller::Node2TrackMap())[node];
256  assert(kTrack);
257  StGlobalTrack *gTrack = static_cast<StGlobalTrack *>(node->track(global));
258  assert(gTrack);
259 #ifdef ADD_NEW_NODE
260  // Replace dca node by a primary vertex
261  StiKalmanTrackNode *tNode = kTrack->getInnerMostNode();
262  if (! tNode->isDca()) continue;
263 #if 1
264  // tNode->print("XYZEPIJK");
265 #endif
266  tNode->rotate(-tNode->getAlpha());
267  // StiHit localVertex = *Vertex;
268  // localVertex.rotate(tNode->getAlpha());
269  // tNode->setHit(&localVertex);
270  tNode->setHit(Vertex);
271  tNode->setDetector(0);
272  Double_t Phi, dPhi;
273  P.GetPhi(Phi,dPhi);
274  Double_t pT, dpT;
275  P.GetPt(pT,dpT);
276  StiNodePars &FP = tNode->fitPars();
277  FP.x() = P.GetX();
278  FP.y() = P.GetY(); // local Y-coordinate of this track (reference plane)
279  FP.z() = P.GetZ(); // local Z-coordinate of this track (reference plane)
280  FP.eta() = Phi - tNode->getAlpha(); // (signed curvature)*(local Xc of helix axis - X current point on track)
281  FP.ptin() = - P.GetQ()/pT; // signed invert pt [sign = sign(-qB)]
282  FP.tanl() = P.GetPz()/pT; // tangent of the track momentum dip angle
283  // FP.curv() = FP.hz()/FP.ptin(); // signed curvature [sign = sign(-qB)]
284  FP.ready();
285  // FP.hz() = 0; // Z component magnetic field in units Pt(Gev) = Hz * RCurv(cm)
286  Double_t pzpT3 = - P.GetPz()/(pT*pT*pT);
287  Double_t f[6*6] = {
288  /* x, y, z, pX, pY, pZ */
289  /* x */ 1, 0, 0, 0, 0, 0,
290  /* y */ 0, 1, 0, 0, 0, 0,
291  /* z */ 0, 0, 1, 0, 0, 0,
292  /* eta */ 0, 0, 0, P.GetPy()/(pT*pT), -P.GetPx()/(pT*pT), 0,
293  /* q/pT */ 0, 0, 0, -FP.ptin()*P.GetPx()/(pT*pT), -FP.ptin()*P.GetPy()/(pT*pT), 0,
294  /* tanL */ 0, 0, 0, pzpT3*P.GetPx(), pzpT3*P.GetPy(), 1./pT};
295  TRMatrix F(6,6,f); if (Debug()) {LOG_INFO << "F\t" << F << endm;}
296  TRSymMatrix CovP(6,&((KFParticle *)&P)->Covariance(0)); if (Debug()) {LOG_INFO << "CovP\t" << CovP << endm;}
297  TRSymMatrix Covi(F,TRArray::kAxSxAT,CovP); if (Debug()) {LOG_INFO << "Covi\t" << Covi << endm;}
298  StiNodeErrs &FE = tNode->fitErrs();
299  TCL::ucopy(Covi.GetArray(), FE.A, 21);
300  tNode->setReady();
301 #if 0
302  // tNode->print("XYZEPIJK");
303  StiKalmanTrackNode *test = kTrack->getInnerMostHitNode();
304  assert(test == tNode);
305  // Int_t status = kTrack->refit(); // refit with primary vertex
306  Int_t status = kTrack->fit(kInsideOut);
307  if (status) continue; // failed to refit
308 #endif
309  kTrack->setPrimary(l+1);
310 #else /* ! ADD_NEW_NODE */
311  StiKalmanTrackNode *extended = (StiKalmanTrackNode*) kTrack->extendToVertex(Vertex);
312  if (extended) {
313  kTrack->add(extended,kOutsideIn);
314  if (extended && !extended->isValid()) extended=0;
315  if (extended && extended->getChi2()>1000) extended=0;
316  }
317  kTrack->reduce();
318  if (! extended) continue;
319  //? kTrack->add(extended,kOutsideIn);
320  kTrack->setPrimary(l+1);
321  extended->setUntouched();
322  Int_t ifail = kTrack->refit();
323  ifail |= (kTrack->getInnerMostHitNode(3)!=extended);
324  kTrack->reduce();
325  // something is wrong. It is not a primary
326  if (ifail) {
327  kTrack->removeLastNode();
328  kTrack->setPrimary(0);
329  continue;
330  }
331 #endif /* ADD_NEW_NODE */
332  //________________________________________________________________________________
333  StTrackDetectorInfo* detInfo = new StTrackDetectorInfo;
334  StiStEventFiller::instance()->fillDetectorInfo(detInfo,kTrack,kFALSE); //3d argument used to increase/not increase the refCount. MCBS oct 04.
335  // StiStEventFiller::instance()->fillPulls(kTrack,1);
336  StPrimaryTrack* pTrack = new StPrimaryTrack;
337  node->addTrack(pTrack); // StTrackNode::addTrack() calls track->setNode(this);
338  pTrack->setKey( gTrack->key());
339  pTrack->setFlagExtension( gTrack->flagExtension());
340  StiStEventFiller::instance()->fillTrack(pTrack,kTrack, detInfo);
341  // set up relationships between objects
342  detInfoVec.push_back(detInfo);
343  primV->addDaughter(pTrack);
344  //________________________________________________________________________________
345  }
346  if (beam ) primV->setBeamConstrained();
347  //..... add vertex to the list
348  if (primV->numberOfDaughters() < 1) {
349  delete primV;
350  } else {
351  primV->setTrackNumbers();
352  calculateRank(primV);
353  pEvent->addPrimaryVertex(primV,orderByRanking);
354  }
355  }
356  return kStOK;
357 }
358 //________________________________________________________________________________
359 void StKFVertexMaker::calculateRank(StPrimaryVertex *primV) {
360  // Calculation of veretx ranks to select 'best' (i.e. triggered) vertex
361  // Simpilfied version (w/o weighting)
362  Float_t rank = primV->probChiSquared();
363  static Float_t Wveto = 1;
364  static Float_t Wmatch = 4;
365  if (primV->isBeamConstrained()) rank += Wmatch;
366  rank -= Wveto*primV->numPostXTracks();
367  rank += Wmatch*primV->numTracksWithPromptHit();
368  rank += Wmatch*primV->numTracksCrossingCentralMembrane();
369  rank += Wmatch*primV->numMatchesWithCTB()
370  - Wveto*primV->numNotMatchesWithCTB();
371  rank += Wmatch*primV->numMatchesWithBTOF()
372  - Wveto*primV->numNotMatchesWithBTOF();
373  rank += Wmatch*(primV->numMatchesWithBEMC() + primV->numMatchesWithEEMC());
374  // It is unclear why the next line was uncommented but apparently it is not
375  // used. So, either this line or the semicolumn at the end of the preceeding
376  // line should be removed
377  //- Wveto*(primV->numNotMatchesWithBEMC() + primV->numNotMatchesWithEEMC());
378  if (primV->numTracksTpcWestOnly() > 0 && primV->numTracksTpcEastOnly() > 0)
379  rank += Wmatch*TMath::Min(primV->numTracksTpcWestOnly(),primV->numTracksTpcEastOnly());
380  primV->setRanking(rank);
381  if (Debug()) primV->Print();
382 }
383 //________________________________________________________________________________
384 KFParticle *StKFVertexMaker::AddTrackAt(const StDcaGeometry *dca, Int_t kg) {
385  fParticles->AddAtAndExpand (0,kg);
386  if (! dca) return 0;
387  Double_t xyzp[6], CovXyzp[21];
388  dca->GetXYZ(xyzp,CovXyzp);
389  static MTrack track;
390  track.SetParameters(xyzp);
391  track.SetCovarianceMatrix(CovXyzp);
392  track.SetNDF(1);
393  // track.SetChi2(GlobalTracks_mChiSqXY[k]);
394  track.SetID(kg);
395  Int_t q = 1;
396  Int_t pdg = 211;
397  if (dca->charge() < 0) {
398  q = -1;
399  pdg = -211;
400  }
401  track.SetCharge(q);
402  KFParticle *particle = new KFParticle(track, pdg);
403  particle->SetID(kg);
404  fParticles->AddAt(particle,kg);
405  return particle;
406 }
407 //________________________________________________________________________________
408 void StKFVertexMaker::Fit() {
409  if (Debug() != 2) StKFVertex::SetDebug(Debug());
410  fcVertices = 0;
411  for (Int_t i = 0; i < fNPasses+1; i++) {
412  SafeDelete(fVerticesPass[i]);
413  }
414  Int_t NGoodGlobals = Particles().GetLast();
415 
416  Double_t TempLog = fTempLog; // default Temperature Log
417  for (Int_t pass = 0; pass < fNPasses; pass++) {
418  Int_t nAccepted = 0;
419  Double_t dZ = fVtxs[pass]->GetBinWidth(1);
420  for (Int_t k = 0; k < NGoodGlobals; k++) {
421  KFParticle *particle = (KFParticle *) Particles()[k];
422  if (! particle) continue;
423  Double_t pT;
424  Double_t dpT;
425  particle->GetPt(pT,dpT);
426  Double_t offset = 0.5*particle->GetPz()/pT;
427  Double_t SigmaZ = TMath::Sqrt(particle->Covariance(2,2) + offset*offset);
428  SigmaZ += dZ;
429  Double_t Z = particle->GetZ();
430  fVtxKs[pass]->Fill(Z);
431  Int_t bin1 = fVtxs[pass]->FindBin(Z - 5*SigmaZ);
432  if (bin1 < 1) bin1 = 1;
433  Int_t bin2 = fVtxs[pass]->FindBin(Z + 5*SigmaZ);
434  if (bin2 > fNzBins) bin2 = fNzBins;
435  Double_t z = fVtxs[pass]->GetBinCenter(bin1);
436  for (Int_t bin = bin1; bin <= bin2; bin++, z += dZ) {
437  fVtxs[pass]->Fill(z,(TMath::Erfc((z - Z - fzWindow)/SigmaZ) - TMath::Erfc((z - Z + fzWindow)/SigmaZ))/2.);
438  }
439  nAccepted++;
440  }
441  Double_t F = fVtxKs[pass]->GetEntries();
442  if (F < 1) continue;
443  fVtxKs[pass]->SetNormFactor(F/dZ);
444  fVtx = fVtxs[0]; // << switch between types Vtx = fVtxKs[0];
445  TString opt("new");
446  if (! Canvas()) opt = "goff";
447  Int_t nfound = fSpectrum->Search(fVtx,3,opt,TMath::Min(0.1,5./NGoodGlobals));
448  if (! nfound) continue;
449  if (Canvas()) {
450  Canvas()->cd();
451  fVtxs[0]->Draw(); fVtxKs[0]->Draw("same");
452  fVtxM->Draw("same");
453  if (pass) fVtx->Draw("same");
454  Canvas()->Update();
455  }
456  if (StKFVertex::Debug() > 1) {
457  LOG_INFO << "Found " << nfound
458  << " candidate peaks to fit with " << NGoodGlobals
459  << " good globals from with " << nAccepted << " accepted" << endm;
460  }
461  Double_t *zOfPeaks = new Double_t[nfound];
462  Int_t npeaks = 0;
463 #if ROOT_VERSION_CODE > 336641 /* ROOT_VERSION(5,35,1) */
464  Double_t *xpeaks = fSpectrum->GetPositionX();
465 #else
466  Float_t *xpeaks = fSpectrum->GetPositionX();
467 #endif
468  for (Int_t p = 0; p < nfound; p++) {
469 #if ROOT_VERSION_CODE > 336641 /* ROOT_VERSION(5,35,1) */
470  Double_t xp = xpeaks[p];
471 #else
472  Float_t xp = xpeaks[p];
473 #endif
474  Int_t bin = fVtx->GetXaxis()->FindBin(xp);
475  Double_t yp = fVtx->GetBinContent(bin);
476  Double_t ep = fVtx->GetBinError(bin);
477  if (yp-1.25*ep < 0) continue;
478  zOfPeaks[npeaks] = xp;
479  npeaks++;
480  }
481  if (StKFVertex::Debug() > 1) {
482  LOG_INFO << "Found " << npeaks << " useful peaks to fit" << endm;
483  }
484  if (! npeaks) {delete [] zOfPeaks; break; }
485  if (fVerticesPass[pass]) {delete fVerticesPass[pass]; fVerticesPass[pass] = 0;}
486  fVerticesPass[pass] = new StKFVerticesCollection(npeaks, zOfPeaks);
487  delete [] zOfPeaks;
488  fcVertices = fVerticesPass[pass];
489  fcVertices->DoTrack2VertexAssociation(Particles());
490  if (! fcVertices->NoVertices()) continue;
491  if (AnnelingFcn(TMath::Exp(-TempLog)) <= 0) continue;
492  if (! fcVertices->NoVertices()) continue;
493  fcVertices->UniqueTracks2VertexAssociation(); // Make track associated with only vertex
494  // fcVertices->PrintV(NoMuMcVertex,NoMuMcTrack,StMuMcVertex_time,
495  // StMuMcVertex_xyzV_mX1,StMuMcVertex_xyzV_mX2,StMuMcVertex_xyzV_mX3,
496  // StMuMcVertex_NoDaughters,StMuMcVertex_IdParTrk,StMuMcTrack_gePid);
497  }
498  if (! fVerticesPass[0]) return;
499  if (fNPasses > 1 && Canvas()) {
500  Canvas()->cd();
501  fVtxs[1]->Draw("same");
502  Canvas()->Update();
503  }
504  Int_t N1 = fVerticesPass[0]->NoVertices();
505  if (! N1) return;
506  if (fVerticesPass[1]) {
507  *fVerticesPass[0] += *fVerticesPass[1];
508  }
509  fcVertices = fVerticesPass[0];
510  fcVertices->MergeDuplicatedVertices();
511  if (! fcVertices->NoVertices()) return;
512  // Double_t Temperature = TMath::Exp(TempLog);
513  TempLog = 5;
514  Double_t Temperature = TMath::Exp(TempLog);
515 #if 1
516  // secondary vertices
517  Int_t pass = fNPasses;
518  if (fVerticesPass[pass]) {delete fVerticesPass[pass]; fVerticesPass[pass] = 0;}
519  fVerticesPass[pass] = new StKFVerticesCollection();
520  fcVertices = fVerticesPass[pass];
521  StAnneling::SetTemperature(Temperature);
522  for (Int_t k = 1; k < NGoodGlobals; k++) {
523  KFParticle *particleK = (KFParticle *) Particles()[k];
524  if (! particleK) continue;
525  if (particleK->GetID() > 100000) continue;
526  StKFVertex *vtx = 0;
527  for (Int_t l = k+1; l < NGoodGlobals; l++) {
528  KFParticle *particleL = (KFParticle *) Particles()[l];
529  if (! particleL) continue;
530  if (particleL->GetID() > 100000) continue;
531  Double_t dist = particleK->GetDistanceFromParticle(*particleL);
532  if (dist > 5.0) continue;
533  if (! vtx) {
534  vtx = new StKFVertex(fcVertices->NoVertices() + 1);
535  vtx->AddTrack(new StKFTrack(k,particleK));
536  }
537  vtx->AddTrack(new StKFTrack(l,particleL));
538  }
539  if (! vtx) continue;
540  vtx->Fit();
541  Int_t N = vtx->NoTracks();
542  if (! N) {delete vtx; vtx = 0; continue;}
543  Double_t X = vtx->Vertex().X();
544  Double_t Y = vtx->Vertex().Y();
545  Double_t R = TMath::Sqrt(X*X + Y*Y);
546  if (R > 200 ) {delete vtx; vtx = 0; continue;}
547  Double_t prob = TMath::Prob(vtx->Vertex().GetChi2(),vtx->Vertex().GetNDF());
548  if (N > 2 || prob > 1.e-3) {// Allow V2 to share tracks
549  for (Int_t i = 0; i < N; i++) {
550  KFParticle *particle = vtx->Track(i)->OrigParticle();;
551  Int_t ID = particle->GetID()%100000 + 100000*vtx->ID();;
552  particle->SetID(ID);
553  }
554  }
555  fcVertices->AddVertex(vtx);
556  }
557  if (StKFVertex::Debug() > 1) {
558  LOG_INFO << "Candidate for secondary vertices: " << fcVertices->NoVertices() << endm;
559  }
560  if ( fcVertices->NoVertices() ) {
561  // fcVertices->PrintV(NoMuMcVertex,NoMuMcTrack,StMuMcVertex_time,
562  // StMuMcVertex_xyzV_mX1,StMuMcVertex_xyzV_mX2,StMuMcVertex_xyzV_mX3,
563  // StMuMcVertex_NoDaughters,StMuMcVertex_IdParTrk,StMuMcTrack_gePid);
564  *fVerticesPass[0] += *fVerticesPass[fNPasses];
565  }
566  // end of loop for secondary vertices
567 #endif
568  fcVertices = fVerticesPass[0];
569  fcVertices->Compress();
570  if (! fcVertices->NoVertices()) return;
571  fcVertices->MergeDuplicatedVertices();
572  fminBrent->SetFunction(*func,TMath::Exp(-0.5*(TempLog)),TMath::Exp(-TempLog),1);
573  if (! fminBrent->Minimize(10,0.1,0.1)) {
574  LOG_WARN << "Temperature fit has failed" << endm;
575  Temperature = 1;
576  } else {
577  Temperature = 1./fminBrent->XMinimum();
578  }
579  StAnneling::SetTemperature(Temperature);
580  fcVertices->UniqueTracks2VertexAssociation(); // Make track associated with only vertex
581  fcVertices->Fit(29,Canvas(),fVtx);
582  if (Canvas()) Canvas()->Update();
583 }
584 //________________________________________________________________________________
585 Double_t StKFVertexMaker::AnnelingFcn(Double_t TInv) {
586  if (! fcVertices) return 0;
587  Double_t Temperature = 1./TInv;
588  StAnneling::SetTemperature(Temperature);
589  Double_t Chi2 = fcVertices->Fit();
590  if (StKFVertex::Debug())
591  LOG_INFO << "StKFVertexMaker::AnnelingFcn\tTemperature = " << Temperature << " Chi2 = " << Chi2 << endm;
592  return Chi2;
593 }
594 //________________________________________________________________________________
595 // $Log: StKFVertexMaker.cxx,v $
596 // Revision 2.8 2018/04/10 11:32:09 smirnovd
597 // Minor corrections across multiple files
598 //
599 // - Remove ClassImp macro
600 // - Change white space
601 // - Correct windows newlines to unix
602 // - Remove unused debugging
603 // - Correct StTpcRTSHitMaker header guard
604 // - Remove unused preprocessor directives in StiCA
605 // - Minor changes in status and debug print out
606 // - Remove using std namespace from StiKalmanTrackFinder
607 // - Remove includes for unused headers
608 //
609 // Revision 2.7 2016/06/08 23:32:46 smirnovd
610 // Integration of StiCA
611 //
612 // This is a squashed commit with all changes combined. To see individual
613 // modifications check out the ds-StiCA_2016 branch in star-sti repository.
614 // Alternatively, one can explore the StiCA_2016 branch in the STAR's CVS
615 // repository.
616 //
617 // Revision 2.6 2013/04/10 22:14:20 fisyak
618 // Roll back to version 04/04/2013
619 //
620 // Revision 2.4 2013/01/28 21:51:17 fisyak
621 // Correct ranking
622 //
623 // Revision 2.3 2013/01/17 15:57:25 fisyak
624 // Add handles for debugging
625 //
626 // Revision 2.2 2012/09/16 21:38:42 fisyak
627 // use of Tpc West Only and East Only tracks, clean up
628 //
629 // Revision 2.1 2012/05/07 14:56:14 fisyak
630 // Add StKFVertexMaker
631 //
632 // Revision 1.5 2012/04/13 14:42:58 fisyak
633 // Freeze
634 //
635 // Revision 1.4 2012/03/29 23:35:47 fisyak
636 // Fix problem with multiple beam tracks
637 //
638 // Revision 1.3 2012/03/26 23:42:35 fisyak
639 // Add beam constrain
640 //
641 // Revision 1.2 2012/02/20 22:38:34 fisyak
642 // Freeze before go for ranking
643 //
644 // Revision 1.1 2012/02/18 23:20:52 fisyak
645 // Rename StKFVertexFitter => StKFVertexMaker
646 //
647 // Revision 1.3 2012/02/07 19:38:26 fisyak
648 // Repackage
649 //
StKFVertexMaker::Clear
void Clear(Option_t *option="")
User defined functions.
Definition: StKFVertexMaker.cxx:56
StiKalmanTrack::getInnerMostHitNode
StiKalmanTrackNode * getInnerMostHitNode(int qua=0) const
Accessor method returns the inner most hit node associated with the track.
Definition: StiKalmanTrack.cxx:1229
StiKalmanTrack
Definition of Kalman Track.
Definition: StiKalmanTrack.h:82
StiHit::setGlobal
void setGlobal(const StiDetector *detector, const StMeasuredPoint *stHit, Float_t x, Float_t y, Float_t z, Float_t energy)
Definition: StiHit.cxx:133
StiKalmanTrack::extendToVertex
StiTrackNode * extendToVertex(StiHit *vertex)
Definition: StiKalmanTrack.cxx:1352
StiHit
Definition: StiHit.h:51
StiKalmanTrackNode::rotate
int rotate(double alpha)
Definition: StiKalmanTrackNode.cxx:1840
StKFVertexMaker::Make
virtual Int_t Make()
Definition: StKFVertexMaker.cxx:110
Factory::getInstance
virtual Abstract * getInstance()=0
Get a pointer to instance of objects served by this factory.
StiKalmanTrack::getInnerMostNode
StiKalmanTrackNode * getInnerMostNode(int qua=0) const
Accessor method returns the inner most node associated with the track.
Definition: StiKalmanTrack.h:594
StiKalmanTrack.h
Definition of Kalman Track.
MTrack
Definition: MTrack.h:4
kStOK
Definition: Stypes.h:40
beam
Definition: beam.h:43
StiToolkit.h
Abstract interface for a STI toolkit.
StiKalmanTrack::add
virtual void add(StiTrackNode *node, int direction, StiTrackNode *near=0)
Definition: StiKalmanTrack.cxx:1518
St_vertexSeedC
Definition: St_vertexSeedC.h:7
KFParticle::SetField
static void SetField(Double_t Bz)
Definition: KFParticle.h:335
StiStEventFiller::fillDetectorInfo
void fillDetectorInfo(StTrackDetectorInfo *detInfo, StiKalmanTrack *kTrack, bool refCountIncr)
Definition: StiStEventFiller.cxx:929
Generated by   doxygen 1.8.5