StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPPVertexFinder.cxx
1 /************************************************************
2  *
3  * $Id: StPPVertexFinder.cxx,v 1.120 2017/10/03 21:25:15 genevb Exp $
4  *
5  * Author: Jan Balewski
6  ************************************************************
7  *
8  * Description: does not fear any pileup
9  ************************************************************
10  *
11  * $Log: StPPVertexFinder.cxx,v $
12  * Revision 1.120 2017/10/03 21:25:15 genevb
13  * Proper reporting/usage of useBTOFmatchOnly, plus LOGGER usage bugs fixed
14  *
15  * Revision 1.119 2017/08/04 21:14:55 genevb
16  * Remove memory leak of StEmcCollection (RT 3303)
17  *
18  * Revision 1.118 2017/05/30 18:27:19 smirnovd
19  * StPPvertexFinder: Removed overlooked reference to a debug histogram
20  *
21  * See commit 312bb0f6 "StPPVertexFinder: Do not fill debug histograms"
22  *
23  * Revision 1.117 2017/05/24 05:02:05 genevb
24  * Options for number of unqualified verts to store, and using only BTOF-matched tracks
25  *
26  *
27  ************************************************************/
28 
29 #include <St_base/StMessMgr.h>
30 
31 #include "TFile.h"
32 #include "TH1D.h"
33 #include "TH1F.h"
34 #include "TH2F.h"
35 #include "TObjArray.h"
36 
37 #include <tables/St_g2t_vertex_Table.h> // tmp for Dz(vertex)
38 
39 #include "StGenericVertexMaker/StiPPVertex/StPPVertexFinder.h"
40 #include "StGenericVertexMaker/StGenericVertexMaker.h"
41 #include "StGenericVertexMaker/Minuit/St_VertexCutsC.h"
42 #include "StEvent/StDcaGeometry.h"
43 #include "StEvent/StEventTypes.h"
44 #include "StEvent/StEnumerations.h"
45 #include "StEvent/StTrack.h"
46 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
47 #include "StMuDSTMaker/COMMON/StMuDst.h"
48 #include "StMuDSTMaker/COMMON/StMuTrack.h"
49 #include "StMuDSTMaker/COMMON/StMuEmcCollection.h"
50 #include "StMuDSTMaker/COMMON/StMuEmcUtil.h"
51 
52 #include <Sti/StiToolkit.h>
53 #include <Sti/StiKalmanTrack.h>
54 #include <Sti/StiKalmanTrackNode.h>
55 #include <Sti/StiTrackContainer.h>
56 #include "Sti/StiTrack.h"
57 
58 #include <St_db_Maker/St_db_Maker.h>
59 #include <StIOMaker/StIOMaker.h> // to save local histos
60 #include <StBFChain/StBFChain.h>
61 
62 #include "StGenericVertexMaker/StiPPVertex/BtofHitList.h"
63 #include "StGenericVertexMaker/StiPPVertex/CtbHitList.h"
64 #include "StGenericVertexMaker/StiPPVertex/BemcHitList.h"
65 #include "StGenericVertexMaker/StiPPVertex/EemcHitList.h"
66 
67 #include "StEvent/StEmcCollection.h"
68 #include "StEvent/StBTofCollection.h"
69 #include "StBTofUtil/StBTofGeometry.h"
70 
71 //==========================================================
72 //==========================================================
73 
74 StPPVertexFinder::StPPVertexFinder(VertexFit_t fitMode) :
75  StGenericVertexFinder(SeedFinder_t::PPVLikelihood, fitMode),
76  mTrackData(), mVertexData(),
77  mTotEve(0), eveID(0), nBadVertex(0),
78  mAlgoSwitches(kSwitchOneHighPT),
79  hA{}, hACorr(nullptr), hL(nullptr), hM(nullptr), hW(nullptr),
80  HList(),
81  ntrk{},
82  mMinTrkPt(0.2),
83  mMaxTrkDcaRxy(3.),
84  mMaxZradius(3.),
85  mMinMatchTr(5),
86  mMaxZrange(200.),
87  mMinZBtof(-3.),
88  mMaxZBtof(3.),
89  mMinAdcBemc(8),
90  mMinAdcEemc(5),
91  mMinFitPfrac(0.51),
92  mFitPossWeighting(false),
93  mDropPostCrossingTrack(true), // default PCT rejection on
94  mStoreUnqualifiedVertex(5),
95  mCut_oneTrackPT(10.),
96  mUseBTOFmatchOnly(false),
97  mToolkit(nullptr),
98  btofList(nullptr),
99  ctbList(nullptr),
100  bemcList(new BemcHitList()),
101  eemcList(new EemcHitList()),
102  mStMuDst(nullptr)
103 {
104  mDebugLevel = 1;
105  mUseCtb = true; // default CTB is in the data stream
106  mVertexOrderMethod = orderByRanking; // change ordering by ranking
107 
108  // special histogram for finding the vertex, not to be saved
109  int nb = 5000;
110  float zRange = 250; // (cm)
111 
112  hL = new TH1D("ppvL", "Vertex likelyhood; Z /cm", nb, -zRange, zRange);
113  // needed only for better errZ calculation
114  hM = new TH1D("ppvM", "cumulative track multiplicity; Z /cm", nb, -zRange, zRange);
115  hW = new TH1D("ppvW", "cumulative track weight; Z /cm", nb, -zRange, zRange);
116 }
117 
118 
119 //==========================================================
120 //==========================================================
121 void StPPVertexFinder::Init()
122 {
123  assert(mTotEve==0); // can't be called twice
124  LOG_INFO << Form("PPV-algo switches=0x%0x, following cuts have been activated:",mAlgoSwitches)<<endm;
125 
126  //get pointer to Sti toolkit
127  mToolkit = StiToolkit::instance();
128 
129  // BTOF and/or CTB hits can be requested after the finder is constructed but
130  // before the Init() is called. In this case we need to create the
131  // corresponding hit lists if they are not available
132  if (mUseBtof && !btofList) {
133  btofList = new BtofHitList();
134  }
135 
136  if (mUseCtb && !ctbList) {
137  ctbList = new CtbHitList();
138  }
139 
140  LOG_INFO << "Finished Init" << endm;
141 }
142 
143 //==========================================================
144 //==========================================================
145 void StPPVertexFinder::InitRun(int run_number, const St_db_Maker* db_maker)
146 {
147  StGenericVertexFinder::InitRun(run_number, db_maker);
148 
149  int dateY = db_maker->GetDateTime().GetYear();
150 
151  //.. set various params
152  // It is not clear why one would hard code cuts for any specific run or
153  // a period since they can be set in the database. Here we'll assume that for
154  // Runs 5 to 12 the PPV cuts are optimized and there is no need to access the
155  // values from the database.
156  if (run_number >= 6000000 && run_number < 13000000) {
157  // old defaults, pre-Run12
158  // (important if we want to reprocess old data with different cuts!)
159  LOG_INFO << "PPV InitRun() using old, hardwired cuts" << endm;
160  mMaxTrkDcaRxy = 3.0; // cm
161  mMinTrkPt = 0.20; // GeV/c //was 0.2 in 2005 prod
162  mMinFitPfrac = 0.7; // nFit /nPossible points on the track
163  mMaxZradius = 3.0; //+sigTrack, to match tracks to Zvertex
164  mMinMatchTr = 2; // required to accept vertex
165  } else {
166  St_VertexCutsC* vtxCuts = St_VertexCutsC::instance();
167  mMaxTrkDcaRxy = vtxCuts->RImpactMax();
168  mMinTrkPt = vtxCuts->MinTrackPt();
169  mMinFitPfrac = vtxCuts->MinFracOfPossFitPointsOnTrack();
170  mMaxZradius = vtxCuts->DcaZMax(); //+sigTrack, to match tracks to Zvertex
171  mMinMatchTr = vtxCuts->MinTrack(); // required to accept vertex
172  mFitPossWeighting = true;
173  }
174 
175  if(dateY<2006) {
176  mMinAdcBemc = 15; // BTOW used calibration of maxt Et @ ~27Gev
177  } else {
178  mMinAdcBemc = 8; // BTOW used calibration of maxt Et @ ~60Gev
179  }
180 
181  // Unfortunately, forced to remove const...
182  St_db_Maker* st_db_maker = const_cast<St_db_Maker*>(db_maker);
183 
184  if (mUseBtof) btofList->initRun(st_db_maker);
185  if (mUseCtb) ctbList->initRun();
186 
187  bemcList->initRun(st_db_maker);
188  eemcList->initRun(st_db_maker);
189 
190  LOG_INFO << "PPV::cuts "
191  << "\n MinNumberOfFitPointsOnTrack = unused"
192  << "\n MinFitPfrac=nFit/nPos = " << mMinFitPfrac
193  << "\n MaxTrkDcaRxy/cm= " << mMaxTrkDcaRxy
194  << "\n MinTrkPt GeV/c = " << mMinTrkPt
195  << "\n MinMatchTr of prim tracks = " << mMinMatchTr
196  << "\n MaxZrange (cm)for glob tracks = " << mMaxZrange
197  << "\n MaxZradius (cm) for prim tracks &Likelihood = " << mMaxZradius
198  << "\n Min/Max Z position for BTOF hit = " << mMinZBtof << " " << mMaxZBtof
199  << "\n MinAdcBemc for MIP = " << mMinAdcBemc
200  << "\n MinAdcEemc for MIP = " << mMinAdcEemc
201  << "\n bool useCtb = " << mUseCtb
202  << "\n bool useBtof = " << mUseBtof
203  << "\n bool useBTOFmatchOnly = " << mUseBTOFmatchOnly
204  << "\n bool nFit/nPoss weighting = " << mFitPossWeighting
205  << "\n bool DropPostCrossingTrack = " << mDropPostCrossingTrack
206  << "\n Store # of UnqualifiedVertex = " << mStoreUnqualifiedVertex
207  << "\n Store=" << (mAlgoSwitches & kSwitchOneHighPT)
208  << " oneTrack-vertex if track PT/GeV>" << mCut_oneTrackPT << endm;
209 }
210 
211 
212 //==========================================================
213 //==========================================================
214 void StPPVertexFinder::initHisto()
215 {
216  hA[0]=new TH1F("ppvStat","event types; 1=inp, 2=trg, 3=-, 4=w/trk, 5=anyMch, 6=Bmch 7=Emch 8=anyV, 9=mulV",10,0.5,10.5);
217  hA[1]=new TH1F("ch1","chi2/Dof, ppv pool tracks",100,0,10);
218  hA[2]=new TH1F("nP","No. of fit points, ppv pool tracks",30,-.5,59.5);
219  hA[3]=new TH1F("zV","reconstructed vertices ; Z (cm)",100,-200,200);
220  hA[4]=new TH1F("nV","No. of vertices per eve",20,-0.5,19.5);
221 
222  hA[5]=new TH1F("rxyDca","Rxy to beam @ DCA ; (cm)",40,-0.1,3.9);
223  hA[6]=new TH1F("nTpcM","No. tracks: tpcMatch /eve ",60,-.5,59.5);
224  hA[7]=new TH1F("nTpcV","No. tracks: tpcVeto /eve ",60,-.5,59.5);
225 
226  hA[8]=0; // (TH1F*) new TH2F ("xyE","Y vs. X of match tracks in EEMC; X (cm); Y(cm)",200,-250.,250,200,-250,250);
227 
228  hA[9]=new TH1F("zDca","Z DCA for all accepted tracks; Z (cm)",100,-200,200);
229 
230  hA[10]=new TH1F("zCtb","Z @CTB for all accepted tracks; Z (cm)",50,-250,250);
231  hA[11]=new TH1F("zBemc","Z @Bemc for all accepted tracks; Z (cm)",100,-400,400);
232  hA[12]=new TH1F("dzVerTr","zVerGeant - zDca of tracks used by any vertex ; (cm)",100,-5,5);
233  hA[13]=new TH1F("dzVerVer","zVerGeant - best reco vertex ; (cm)",100,-5,5);
234 
235  hA[14]=new TH1F("EzDca","Error of Z DCA for all accepted tracks; Z (cm)",100,-0.,4);
236  hA[15]=new TH1F("nTpcT","No. tracks: accepted Dca /eve ",201,-.5,200.5);
237  hA[16]=new TH1F("ptTr","pT, ppv pool tracks; pT (GeV/c) ",50,0.,10.);
238  hA[17]=new TH1F("vRL","PPV Vertex rank, 'funny' X-axis; X=Log10(rank-1e6)+offset", 150, -11,25);
239 
240  hACorr=new TH2F("BTOFvsBEMC","BTOF vs BEMC", 5,-2.5,2.5,5,-2.5,2.5);
241 
242  for (int i=0; i<mxH; i++) if(hA[i]) HList.Add(hA[i]);
243 
244  HList.Add(hACorr);
245 }
246 
247 
248 void StPPVertexFinder::findSeeds_TSpectrum()
249 {
250  std::vector<double> vertexZs = StGenericVertexFinder::FindSeeds_TSpectrum();
251 
252  // Loop over the seeds and associate tracks with it. Then for each seed
253  // create a vertex candidate of VertexData type
254  for (double vertexZ : vertexZs)
255  {
256  VertexData vertex( TVector3(0, 0, vertexZ) );
257  // Need to add a findMatchingTracks(vertex, mTrackData) method...
258  mVertexData.push_back( vertex );
259  }
260 }
261 
262 
273 void StPPVertexFinder::findSeeds_PPVLikelihood()
274 {
275  const float par_rankOffset = 1e6; // to separate class of vertices (approximately)
276 
277  int vertexID=0;
278 
279  while(1)
280  {
281  if ( !buildLikelihoodZ() ) break;
282 
283  VertexData vertex(++vertexID);
284 
285  if ( !findVertexZ(vertex) ) break;
286 
287  bool trigV = evalVertexZ(vertex); // vertex.print();
288 
289  //bump up rank of 2+ track all vertices
290  if (vertex.nAnyMatch >= mMinMatchTr) vertex.Lmax += par_rankOffset;
291 
292  if (!trigV) {
293  // Ignore this "bad" vertex
294  if (nBadVertex >= mStoreUnqualifiedVertex) continue;
295  // ... or keep it
296  nBadVertex++;
297  // ... and bump down rank of sub-prime vertices
298  vertex.Lmax -= par_rankOffset;
299  }
300 
301  mVertexData.push_back(vertex);
302  }
303 }
304 
305 
306 //==========================================================
307 //==========================================================
308 void StPPVertexFinder::Clear()
309 {
310  LOG_INFO << "StPPVertexFinder::Clear(): Finished event " << mTotEve
311  << ": Found " << mVertexData.size()-nBadVertex << " \"good\" and "
312  << nBadVertex << " \"bad\" vertices" << endm;
313 
314  StGenericVertexFinder::Clear();
315 
316  if (btofList) btofList->clear();
317  if (ctbList) ctbList->clear();
318  bemcList->clear();
319  eemcList->clear();
320 
321  mTrackData.clear();
322  mVertexData.clear();
323  eveID = -1;
324  nBadVertex = 0;
325 
326  ntrk.fill(0);
327 
328  // the clear below is not needed but cleans up stale result
329  hL->Reset();
330  hM->Reset();
331  hW->Reset();
332 }
333 
334 
335 //======================================================
336 //======================================================
337 void StPPVertexFinder::printInfo(ostream& os) const
338 {
339  LOG_INFO << "\n"
340  << Form("PPV:: # of input track = %d\n", ntrk[0])
341  << Form("PPV:: dropped due to flag/dummy = %d\n", ntrk[1])
342  << Form("PPV:: dropped due to pt = %d\n", ntrk[2])
343  << Form("PPV:: dropped due to PCT check = %d\n", ntrk[3])
344  << Form("PPV:: dropped due to DCA check = %d\n", ntrk[4])
345  << Form("PPV:: dropped due to NHit check = %d\n", ntrk[5])
346  << Form("PPV:: dropped due to TOF check = %d\n", ntrk[6])
347  << Form("PPV:: # of track after all cuts = %d", ntrk[7]) << endm;
348 
349  if(mUseBtof) btofList->print();
350  if(mUseCtb) ctbList->print();
351 
352  bemcList->print();
353  eemcList->print();
354 
355 
356  os << "StPPVertexFinder ver=1 - Fit Statistics:\n"
357  << "StPPVertexFinder::result " << mVertexData.size() << " vertices found" << std::endl;
358 
359  int nTpcM=0, nTpcV=0;
360  int k=0;
361  for (const TrackData &track : mTrackData) {
362  if( track.mTpc>0) nTpcM++;
363  else if ( track.mTpc<0) nTpcV++;
364  if(track.vertexID<=0) continue; // skip not used or pileup vertex
365  k++;
366  LOG_DEBUG
367  << Form("%d track@z0=%.2f +/- %.2f gPt=%.3f vertID=%d match: bin,Fired,Track:\n",
368  k,track.zDca,track.ezDca,track.gPt,track.vertexID);
369 
370  if (mUseBtof) { LOG_DEBUG
371  << Form(" Btof %3d,%d,%d",track.btofBin,btofList->getFired(track.btofBin),btofList->getTrack(track.btofBin));
372  }
373 
374  if (mUseCtb) { LOG_DEBUG
375  << Form(" CTB %3d,%d,%d",track.ctbBin,ctbList->getFired(track.ctbBin),ctbList->getTrack(track.ctbBin));
376  }
377 
378  LOG_DEBUG
379  << Form(" Bemc %3d,%d,%d",track.bemcBin,bemcList->getFired(track.bemcBin),bemcList->getTrack(track.bemcBin))
380  << Form(" Eemc %3d,%d,%d",track.eemcBin,eemcList->getFired(track.eemcBin),eemcList->getTrack(track.eemcBin))
381  << Form(" TPC %d",track.mTpc)
382  <<endm;
383  }
384 
385  LOG_INFO << Form("PPVend eveID=%d, list of found %d vertices from pool of %d tracks\n",
386  eveID, mVertexData.size(), mTrackData.size()) << endm;
387 
388  for (const VertexData &vertex : mVertexData)
389  vertex.print(os);
390 }
391 
392 
393 //==========================================================
394 //==========================================================
395 int StPPVertexFinder::fit(StEvent* event)
396 {
397  mTotEve++;
398  eveID=event->id();
399 
400  LOG_INFO << "***** START FIT\n"
401  << " @@@@@@ PPVertex::Fit START nEve=" << mTotEve
402  << " eveID=" << eveID << endm;
403 
404  hL->SetTitle("Vertex likelihood, eveID=" + TString(eveID) );
405 
406  // get BTOF info
407  if(mUseBtof) {
408  StBTofCollection *btofColl = (StBTofCollection*)event->btofCollection();
409  if(btofColl==0) {
410  LOG_WARN << "no btofCollection, continue THE SAME eve" << endm;
411  } else {
412  btofList->build(btofColl);
413  }
414  }
415 
416  // get CTB info, does not work for embeding
417  if(mUseCtb) {// CTB could be off since 2006
418  StTriggerData *trgD=event->triggerData ();
419  ctbList->buildFromData(trgD); // use real data
420  }
421 
422 
423  StEmcCollection* emcC =(StEmcCollection*)event->emcCollection();
424  if(emcC==0) {
425  LOG_WARN << "No StEmcCollection found, continue with this event" << endm;
426  } else {
427  StEmcDetector* btow = emcC->detector( kBarrelEmcTowerId);
428  if(btow==0) {
429  LOG_WARN << "No BEMC found in StEmcCollection, continue with this event" << endm;
430  } else {
431  bemcList->build(btow, mMinAdcBemc);
432  }
433 
434  StEmcDetector* etow = emcC->detector(kEndcapEmcTowerId);
435  if(etow==0) {
436  LOG_WARN << "No EEMC found in StEmcCollection, continue with this event" << endm;
437  } else {
438  eemcList->build(etow, mMinAdcEemc);
439  }
440  }
441 
442  //get the Sti track container...
443  assert(mToolkit); // internal error of Sti
444  StiTrackContainer* stiTracks = mToolkit->getTrackContainer();
445  if(stiTracks==0) {
446  LOG_WARN << "No Sti tracks found, skipping this event" << endm;
447  return 0 ;
448  }
449 
450  //select reasonable tracks and add them to my list
451  int kBtof=0,kCtb=0,kBemc=0, kEemc=0,kTpc=0;
452  int nTracksMatchingAnyFastDetector=0;
453 
454  for (const StiTrack* stiTrack : *stiTracks)
455  {
456  const StiKalmanTrack* stiKalmanTrack = static_cast<const StiKalmanTrack*>(stiTrack);
457 
458  TrackDataT<StiKalmanTrack> track(*stiKalmanTrack);
459 
460  ntrk[0]++;
461 
462  if (stiKalmanTrack->getFlag() <0) { ntrk[1]++; continue; }
463  if (stiKalmanTrack->getPt() < mMinTrkPt) { ntrk[2]++; continue; }
464  if (mDropPostCrossingTrack &&
465  isPostCrossingTrack(stiKalmanTrack)) { ntrk[3]++; continue; } // kill if it has hits in wrong z
466  if (!examinTrackDca(stiKalmanTrack, track)) { ntrk[4]++; continue; } // drop from DCA
467  if (!matchTrack2Membrane(track)) { ntrk[5]++; continue; } // kill if nFitP too small
468 
469 
470  // Match to various detectors
471  if (mUseBtof) matchTrack2BTOF(stiKalmanTrack, track); // matching track to btofGeometry
472  if (mUseBTOFmatchOnly && (track.mBtof <= 0)) { ntrk[6]++; continue; }
473 
474  if (mUseCtb) matchTrack2CTB(stiKalmanTrack, track);
475  matchTrack2BEMC(track);
476  matchTrack2EEMC(track);
477 
478  ntrk[7]++;
479 
480  // ...all test done on this track
481  mTrackData.push_back(track);
482 
483 
484  if (track.mBtof > 0) kBtof++;
485  if (track.mCtb > 0) kCtb++;
486  if (track.mBemc > 0) kBemc++;
487  if (track.mEemc > 0) kEemc++;
488  if (track.mTpc > 0) kTpc++;
489 
490  if (track.mBtof>0 || track.mCtb>0 || track.mBemc>0 || track.mEemc>0 || track.mTpc>0)
491  nTracksMatchingAnyFastDetector++;
492  }
493 
494  LOG_INFO << Form("PPV::TpcList size=%d nMatched=%d\n\n",mTrackData.size(),kTpc)
495  << "PPV::fit() nEve=" << mTotEve << " , "
496  << nTracksMatchingAnyFastDetector << " traks with good DCA, matching: BTOF="
497  << kBtof << " CTB=" << kCtb << " BEMC=" << kBemc << " EEMC=" << kEemc << endm;
498 
499  if (nTracksMatchingAnyFastDetector >= mMinMatchTr || mStoreUnqualifiedVertex > 0)
500  {
501  seed_fit_export();
502  } else {
503  LOG_INFO << "StPPVertexFinder::fit() nEve=" << mTotEve << " Quit, to few matched tracks" << endm;
504  }
505 
506  return size();
507 }
508 
509 
510 int StPPVertexFinder::fit(const StMuDst& muDst)
511 {
512  mTotEve++;
513 
514  mStMuDst = &muDst;
515 
516  // Similar to fit() we need to populate bemcList
517  StMuEmcCollection *muEmcCollection = muDst.muEmcCollection();
518 
519  StMuEmcUtil muEmcUtil;
520 
521  StEmcCollection* emcC = muEmcUtil.getEmc(muEmcCollection);
522 
523  StEmcDetector* btow = emcC->detector(kBarrelEmcTowerId);
524  bemcList->build(btow, mMinAdcBemc);
525 
526  StEmcDetector* etow = emcC->detector(kEndcapEmcTowerId);
527  eemcList->build(etow, mMinAdcEemc);
528 
529  delete emcC; emcC=0;
530 
531  // Access btof data from ... branch
532  //TClonesArray* muBTofHits = muDst.btofArray(muBTofHit);
533  //btofList->build(*muBTofHits);
534 
535  // Access array of all StDcaGeometry objects (i.e. tracks)
536  TObjArray* globalTracks = muDst.globalTracks();
537  TClonesArray* covGlobTracks = muDst.covGlobTrack();
538 
539 
540  for (const TObject* obj : *globalTracks)
541  {
542  ntrk[0]++;
543 
544  const StMuTrack& stMuTrack = static_cast<const StMuTrack&>(*obj);
545 
546  if (stMuTrack.pt() < mMinTrkPt) { ntrk[2]++; continue; }
547 
548  // Supposedly equivalent to isPostCrossingTrack()
549  if ( (stMuTrack.flagExtension() & kPostXTrack) != 0 ) { ntrk[3]++; continue; }
550 
551  // Supposedly equivalent to DCA check with examinTrackDca()
552  if (stMuTrack.index2Cov() < 0) { ntrk[4]++; continue; }
553 
554  StDcaGeometry* dca = static_cast<StDcaGeometry*>(covGlobTracks->At(stMuTrack.index2Cov()));
555 
556  if ( std::fabs(dca->z()) > mMaxZrange ) { ntrk[4]++; continue; }
557  if ( std::fabs(dca->impact()) > mMaxTrkDcaRxy) { ntrk[4]++; continue; }
558 
559  // Condition similar to one in matchTrack2Membrane
560  double fracFit2PossHits = static_cast<double>(stMuTrack.nHitsFit(kTpcId)) / stMuTrack.nHitsPoss(kTpcId);
561  if (fracFit2PossHits < mMinFitPfrac) { ntrk[5]++; continue; } // kill if nFitP too small
562 
563  // Test TOF match if required
564  if (mUseBTOFmatchOnly && (stMuTrack.tofHit() == 0)) { ntrk[6]++; continue; }
565 
566  ntrk[7]++;
567 
568  TrackDataT<StMuTrack> track(stMuTrack, dca);
569 
570  // Modify track weights
571  matchTrack2BEMC(track);
572  matchTrack2EEMC(track);
573  matchTrack2Membrane(track);
574 
575  mTrackData.push_back(track);
576  }
577 
578  seed_fit_export();
579 
580  return size();
581 }
582 
583 
584 void StPPVertexFinder::seed_fit_export()
585 {
586  // Select a method to find vertex candidates/seeds. The methods work using the
587  // `mTrackData` and `mDCAs` containers as input whereas the reconstructed
588  // vertices are put in the private container `mVertexData`
589  switch (mSeedFinderType)
590  {
591  case SeedFinder_t::TSpectrum:
592  findSeeds_TSpectrum();
593  break;
594 
595  case SeedFinder_t::PPVLikelihood:
596  default:
597  findSeeds_PPVLikelihood();
598  break;
599  }
600 
601  // Refit vertex position for all cases (currently NoBeamline, Beamline1D, and
602  // Beamline3D) except when the BeamlineNoFit option is specified. This is
603  // done to keep backward compatible behavior when by default the vertex was
604  // placed on the beamline
605  if (mVertexFitMode == VertexFit_t::BeamlineNoFit)
606  {
607  for (VertexData &vertex : mVertexData) {
608  const double& z = vertex.r.Z();
609  vertex.r.SetXYZ( beamX(z), beamY(z), z);
610  }
611  }
612  else
613  {
614  size_t n_seeds = mVertexData.size();
615 
616  auto cannot_fit = [this] (VertexData &vertex) { return fitTracksToVertex(vertex) != 0; };
617  mVertexData.erase( std::remove_if(mVertexData.begin(), mVertexData.end(), cannot_fit), mVertexData.end() );
618  // Update the "bad" vertex counter as some vertices could have been
619  // removed in the previous step
620  nBadVertex -= n_seeds - mVertexData.size();
621  }
622 
623  exportVertices();
624 
625  if (mDebugLevel) printInfo();
626 }
627 
628 
629 //==========================================================
630 //==========================================================
631 bool StPPVertexFinder::buildLikelihoodZ()
632 {
633  hL->Reset();
634  hM->Reset();
635  hW->Reset();
636 
637  float dzMax2 = mMaxZradius*mMaxZradius;
638 
639  int nt=mTrackData.size();
640  LOG_DEBUG<< Form("PPV::buildLikelihood() pool of nTracks=%d",nt)<<endm;
641  if(nt<=0) return false;
642 
643  int nTracksMatchingAnyFastDetector = 0;
644 
645  double *La = hL->GetArray(); // PPV main likelihood histogram
646  double *Ma = hM->GetArray(); // track multiplicity histogram
647  double *Wa = hW->GetArray(); // track weight histogram
648 
649  // Loop over pre-selected tracks only
650  for (const TrackData &track : mTrackData)
651  {
652  // Skip tracks already assigned to a vertex
653  if (track.vertexID != 0) continue;
654 
655  if (track.anyMatch) nTracksMatchingAnyFastDetector++;
656 
657  float z0 = track.zDca; // z coordinate at DCA
658  float ez = track.ezDca; // error on z coordinate at DCA
659  float ez2 = ez*ez;
660  int j1 = hL->FindBin(z0-mMaxZradius-.1);
661  int j2 = hL->FindBin(z0+mMaxZradius+.1);
662  float base = dzMax2/2/ez2;
663  float totW = track.weight;
664 
665  for (int j=j1; j<=j2; j++)
666  {
667  float z = hL->GetBinCenter(j);
668  float dz = z-z0;
669  float xx = base - dz*dz/2/ez2;
670  if (xx <= 0) continue;
671  La[j] += xx*totW;
672  Ma[j] += 1.;
673  Wa[j] += totW;
674  }
675  }
676 
677  LOG_DEBUG << Form("PPV::buildLikelihood() %d tracks w/ matched @ Lmax=%f", nTracksMatchingAnyFastDetector, hL->GetMaximum()) << endm;
678 
679  return (nTracksMatchingAnyFastDetector >= mMinMatchTr) || (mStoreUnqualifiedVertex > 0);
680 }
681 
682 //==========================================================
683 //==========================================================
684 bool StPPVertexFinder::findVertexZ(VertexData &vertex)
685 {
686  if(hL->GetMaximum()<=0) return false; // no more tracks left
687 
688  int iMax = hL->GetMaximumBin();
689  float z0 = hL->GetBinCenter(iMax);
690  float Lmax = hL->GetBinContent(iMax);
691  float accM = hM->GetBinContent(iMax);
692  float accW = hW->GetBinContent(iMax);
693  assert(accM>0);
694  float avrW = accW/accM;
695 
696  // search for sigma of the vertex
697  float Llow = 0.9* Lmax;
698  if ((Lmax-Llow) < 8*avrW ) Llow = Lmax - 8*avrW; // to be at least 4 sigma
699 
700  double *L = hL->GetArray(); // likelyhood
701 
702  int iLow = -1, iHigh = -1;
703  for(int i=iMax; i<=hL->GetNbinsX(); i++) {
704  if(L[i] > Llow) continue;
705  iHigh = i;
706  break;
707  }
708  for(int i=iMax; i>=1; i--) {
709  if(L[i] > Llow) continue;
710  iLow = i;
711  break;
712  }
713 
714  float zLow = hL->GetBinCenter(iLow);
715  float zHigh = hL->GetBinCenter(iHigh);
716 
717  float kSig = std::sqrt(2*(Lmax-Llow)/avrW);
718  float sigZ = (zHigh-zLow)/2/kSig;
719 
720  LOG_DEBUG << Form("PPV:: iLow/iMax/iHigh=%d/%d/%d\n",iLow,iMax,iHigh)
721  << Form(" accM=%f accW=%f avrW=%f\n",accM,accW,avrW)
722  << Form(" Z low/max/high=%f %f %f, kSig=%f, sig=%f\n",zLow,z0,zHigh,kSig,sigZ)
723  << Form(" found PPVertex(ID=%d,neve=%d) z0 =%.2f +/- %.2f\n",vertex.id,mTotEve,z0,sigZ)<<endm;
724 
725  if (sigZ < 0.1) sigZ = 0.1; // tmp, make it not smaller than the bin size
726 
727  // For approximate seed position we use (x,y)=(0,0) because the tracks are
728  // extrapolated to (0,0) anyway. The x and y coordinates can be updated later
729  // in a proper fit.
730  vertex.r = TVector3(0, 0, z0);
731  vertex.er = TVector3(0.1, 0.1, sigZ);
732  vertex.Lmax = Lmax;
733 
734  return true;
735 }
736 
737 //==========================================================
738 //==========================================================
739 bool StPPVertexFinder::evalVertexZ(VertexData &vertex) // and tag used tracks
740 {
741  // returns true if vertex is accepted accepted
742  int nHiPt=0;
743 
744  for (TrackData &track : mTrackData)
745  {
746  // Skip tracks already matched to a vertex
747  if (track.vertexID != 0) continue;
748 
749  // Do not match tracks to vertex if they are too far in Z
750  // (i.e. |delta_z| > (mMaxZradius + sigma_z))
751  if ( !track.matchVertex(vertex, mMaxZradius) ) continue;
752 
753  // Otherwise, this track belongs to this vertex
754  track.vertexID = vertex.id;
755  vertex.gPtSum += track.gPt;
756  vertex.nUsedTrack++;
757 
758  if( track.gPt>mCut_oneTrackPT && ( track.mBemc>0|| track.mEemc>0) ) nHiPt++;
759 
760  if( track.mTpc>0) vertex.nTpc++;
761  else if ( track.mTpc<0) vertex.nTpcV++;
762 
763  if( track.mBtof>0) vertex.nBtof++;
764  else if ( track.mBtof<0) vertex.nBtofV++;
765 
766  if( track.mCtb>0) vertex.nCtb++;
767  else if ( track.mCtb<0) vertex.nCtbV++;
768 
769  if( track.mBemc>0) vertex.nBemc++;
770  else if ( track.mBemc<0) vertex.nBemcV++;
771 
772  if( track.mEemc>0) vertex.nEemc++;
773  else if ( track.mEemc<0) vertex.nEemcV++;
774 
775  if( track.anyMatch) vertex.nAnyMatch++;
776  else if (track.anyVeto) vertex.nAnyVeto++;
777  }
778 
779  vertex.isTriggered = (vertex.nAnyMatch >= mMinMatchTr) || ( (mAlgoSwitches & kSwitchOneHighPT) && nHiPt>0 );
780  if (!vertex.isTriggered) { // discrad vertex
781  //no match tracks in this vertex, tag vertex ID in tracks differently
782  //vertex.print(cout);
783  LOG_DEBUG << "StPPVertexFinder::evalVertex Vid="<<vertex.id<<" rejected"<<endm;
784  for (TrackData &track : mTrackData) {
785  if(track.vertexID!=vertex.id) continue;
786  track.vertexID=-vertex.id;
787  }
788  return false;
789  }
790 
791  LOG_INFO << "StPPVertexFinder::evalVertex Vid=" << vertex.id
792  << " accepted, nAnyMatch=" << vertex.nAnyMatch
793  << " nAnyVeto=" << vertex.nAnyVeto << endm;
794 
795  return true;
796 }
797 
798 
807 void StPPVertexFinder::createTrackDcas(const VertexData &vertex)
808 {
809  // Consider muDst case
810  if (mStMuDst)
811  {
812  // Just clean the pointers owned by something else
813  mDCAs.clear();
814 
815  for (const TrackData & track : mTrackData) {
816  if ( std::fabs(track.vertexID) != vertex.id) continue;
817  mDCAs.push_back(track.dca);
818  }
819 
820  return;
821  }
822 
823  // Fill member array of pointers to StDcaGeometry objects for selected tracks
824  // in mTrackData corresponding to this vertex. These will be used in static
825  // minimization function
826  while (!mDCAs.empty()) delete mDCAs.back(), mDCAs.pop_back();
827 
828 
829  for (const TrackData & track : mTrackData)
830  {
831  if ( std::fabs(track.vertexID) != vertex.id) continue;
832  if (!track.mother) continue;
833 
834  // This code is adopted from StiStEventFiller::fillDca()
835  StiKalmanTrack tmpTrack = *track.getMother<StiKalmanTrack>();
836  StiKalmanTrackNode *tNode = tmpTrack.extrapolateToBeam();
837 
838  if (!tNode) continue;
839 
840  const StiNodePars &pars = tNode->fitPars();
841  const StiNodeErrs &errs = tNode->fitErrs();
842  float alfa = tNode->getAlpha();
843  float setp[7] = {(float)pars.y(), (float)pars.z(), (float)pars.phi(),
844  (float)pars.ptin(), (float)pars.tanl(), (float)pars.curv(), (float)pars.hz()};
845  setp[2] += alfa;
846  float sete[15];
847 
848  for (int i=1, li=1, jj=0; i<kNPars; li += ++i) {
849  for (int j=1;j<=i;j++) {
850  sete[jj++] = errs.G()[li+j];
851  }
852  }
853 
854  StDcaGeometry* dca = new StDcaGeometry();
855  dca->set(setp, sete);
856  mDCAs.push_back(dca);
857  }
858 }
859 
860 
870 int StPPVertexFinder::fitTracksToVertex(VertexData &vertex)
871 {
872  createTrackDcas(vertex);
873 
874  bool fitRequiresBeamline = star_vertex::requiresBeamline(mVertexFitMode);
875 
876  bool prerequisites = mDCAs.size() > 1 || (mDCAs.size() > 0 && fitRequiresBeamline);
877 
878  if ( !prerequisites )
879  {
880  LOG_WARN << "StPPVertexFinder::fitTracksToVertex: At least two tracks required "
881  << "OR one with beam line. This vertex (id = " << vertex.id
882  << ") coordinates will not be updated" << endm;
883  return 5;
884  }
885 
886  // Recalculate vertex seed coordinates to be used as initial point in the fit
887  StThreeVectorD vertexSeed = CalcVertexSeed(mDCAs);
888 
889  // For fits with beamline force the seed to be on the beamline
890  if ( fitRequiresBeamline )
891  {
892  vertexSeed.setX( beamX(vertexSeed.z()) );
893  vertexSeed.setY( beamY(vertexSeed.z()) );
894  }
895 
896  // Make sure the global pointer points to valid object so Minuit uses correct data
898 
899  int minuitStatus;
900 
901  mMinuit->mnexcm("clear", 0, 0, minuitStatus);
902 
903  static double step[3] = {0.01, 0.01, 0.01};
904 
905  double x_lo = vertexSeed.x() - mMaxTrkDcaRxy;
906  double y_lo = vertexSeed.y() - mMaxTrkDcaRxy;
907  double z_lo = vertexSeed.z() - mMaxZradius;
908 
909  double x_hi = vertexSeed.x() + mMaxTrkDcaRxy;
910  double y_hi = vertexSeed.y() + mMaxTrkDcaRxy;
911  double z_hi = vertexSeed.z() + mMaxZradius;
912 
913  if (mVertexFitMode == VertexFit_t::Beamline1D)
914  {
915  mMinuit->mnparm(0, "z", vertexSeed.z(), step[2], z_lo, z_hi, minuitStatus);
916  } else {
917  mMinuit->mnparm(0, "x", vertexSeed.x(), step[0], x_lo, x_hi, minuitStatus);
918  mMinuit->mnparm(1, "y", vertexSeed.y(), step[1], y_lo, y_hi, minuitStatus);
919  mMinuit->mnparm(2, "z", vertexSeed.z(), step[2], z_lo, z_hi, minuitStatus);
920  }
921 
922  mMinuit->mnexcm("minimize", 0, 0, minuitStatus);
923 
924  // Check fit result
925  if (minuitStatus) {
926  LOG_WARN << "StPPVertexFinder::fitTracksToVertex: Fit did not converge. "
927  << "Check TMinuit::mnexcm() status flag: " << minuitStatus << ". "
928  << "This vertex (id = " << vertex.id << ") coordinates will not be updated" << endm;
929 
930  // The fit has failed but let's keep the vertex anyway. For cases with
931  // beam line we put the vertex on the beam line
932  if ( fitRequiresBeamline ) {
933  const double& z = vertex.r.Z();
934  vertex.r.SetXYZ( beamX(z), beamY(z), z);
935  vertex.er.SetXYZ( mBeamline.err_x0, mBeamline.err_y0, vertex.er.Z());
936  }
937  // Return 0 (=success) in order to keep the vertex
938  return 0;
939  }
940 
941  double chisquare, fedm, errdef;
942  int npari, nparx;
943 
944  mMinuit->mnstat(chisquare, fedm, errdef, npari, nparx, minuitStatus);
945  mMinuit->mnhess();
946 
947  // The default dimension 9=3*3 should work for 1D case (npar = 1) as well
948  double emat[9];
949  /* 0 1 2
950  3 4 5
951  6 7 8 */
952  mMinuit->mnemat(emat, 3);
953 
954  if (mVertexFitMode == VertexFit_t::Beamline1D)
955  {
956  const double& z = mMinuit->fU[0];
957  vertex.r.SetXYZ( beamX(z), beamY(z), z);
958  vertex.er.SetXYZ( mBeamline.err_x0, mBeamline.err_y0, std::sqrt(emat[0]) );
959  } else {
960  vertex.r.SetXYZ(mMinuit->fU[0], mMinuit->fU[1], mMinuit->fU[2]);
961  vertex.er.SetXYZ( std::sqrt(emat[0]), std::sqrt(emat[4]), std::sqrt(emat[8]) );
962  }
963 
964  return 0;
965 }
966 
967 
972 void StPPVertexFinder::exportVertices()
973 {
974  for (VertexData &vertex : mVertexData)
975  {
976  StThreeVectorD r(vertex.r.x(), vertex.r.y(), vertex.r.z());
977 
978  float cov[6]{};
979 
980  cov[0] = vertex.er.x() * vertex.er.x();
981  cov[2] = vertex.er.y() * vertex.er.y();
982  cov[5] = vertex.er.z() * vertex.er.z(); // [5] is correct,JB
983 
984  StPrimaryVertex primV;
985  primV.setPosition(r);
986  primV.setCovariantMatrix(cov);
987  primV.setVertexFinderId(mUseCtb ? ppvVertexFinder : ppvNoCtbVertexFinder);
988  primV.setNumTracksUsedInFinder(vertex.nUsedTrack);
989  primV.setNumMatchesWithBTOF(vertex.nBtof);
990  primV.setNumMatchesWithCTB(vertex.nCtb);
991  primV.setNumMatchesWithBEMC(vertex.nBemc);
992  primV.setNumMatchesWithEEMC(vertex.nEemc);
993  primV.setNumTracksCrossingCentralMembrane(vertex.nTpc);
994  primV.setSumOfTrackPt(vertex.gPtSum);
995  primV.setRanking(vertex.Lmax);
996  primV.setFlag(1); //??? is it a right value?
997 
998  if (mStMuDst)
999  {
1000  for (TrackData &track : mTrackData)
1001  {
1002  StThreeVectorF v_position(vertex.r.x(), vertex.r.y(), vertex.r.z());
1003  StThreeVectorF dist = v_position - track.dca->origin();
1004 
1005  // Calculate total error as fully correlated between DCA and vertex
1006  float total_err_perp = std::sqrt( vertex.er.Perp2() + track.dca->errMatrix()[0] ); // fully uncorrelated
1007  float total_err_z = std::sqrt( vertex.er.z()*vertex.er.z() + track.dca->errMatrix()[2] );
1008 
1009  bool is_daughter = (track.vertexID == vertex.id ||
1010  (std::fabs(dist.perp())/total_err_perp < 3 && std::fabs(dist.z())/total_err_z < 3) );
1011 
1012  if ( !is_daughter ) continue;
1013 
1014  track.vertexID = vertex.id;
1015 
1016  StMuTrack* stMuTrack = const_cast<StMuTrack*>( track.getMother<StMuTrack>() );
1017  stMuTrack->setType(primary);
1018 
1019  // Create StTrack from StMuTrack so idTruth can be calculated for this vertex
1020  StTrack* primTrack = StMuDst::createStTrack(stMuTrack);
1021  primV.addDaughter(primTrack);
1022  }
1023 
1024  primV.setIdTruth();
1025  vertex.mIdTruth = primV.idTruth();
1026 
1027  // The daughter StTrack-s are removed at this time because we don't save them
1028  while ( !primV.daughters().empty() )
1029  delete primV.daughters().back(), primV.daughters().pop_back();
1030  }
1031 
1032  //..... add vertex to the list
1033  addVertex(primV);
1034  }
1035 }
1036 
1037 //-------------------------------------------------
1038 //-------------------------------------------------
1039 void StPPVertexFinder::Finish()
1040 {
1041  LOG_INFO << "StPPVertexFinder::Finish() done, seen eve=" << mTotEve << endm;
1042 }
1043 
1044 //-------------------------------------------------
1045 //-------------------------------------------------
1046 void StPPVertexFinder::saveHisto(TString fname)
1047 {
1048  TString outName = fname + ".hist.root";
1049  TFile f(outName, "recreate");
1050  assert(f.IsOpen());
1051  printf("%d histos are written to '%s' ...\n", HList.GetEntries(), outName.Data());
1052  HList.Write();
1053  f.Close();
1054 }
1055 
1056 //==========================================================
1057 //==========================================================
1058 void StPPVertexFinder::dumpKalmanNodes(const StiKalmanTrack*track)
1059 {
1060  //.................... print all nodes ...........
1062  int in=0,nh=0,nTpc=0;
1063  float zL=999, zH=-999;
1064  for (it=track->begin();it!=track->end();it++,in++) {
1065  StiKalmanTrackNode& ktn = (*it);
1066  if(!ktn.isValid()) continue;
1067  if(ktn.getHit() && ktn.getChi2() >1000) continue;
1068  const StiDetector * det=ktn.getDetector();
1069  assert(!(ktn.x()) || det);
1070  float rxy=ktn.getX();
1071  bool actv= !det || det->isActive(ktn.getY(), ktn.getZ());
1072  if(!det || (rxy>58 && rxy < 190)){
1073  float z=ktn.z_g();
1074  if(zL>z) zL=z;
1075  if(zH<z) zH=z;
1076  if(actv) {
1077  nTpc++;
1078  if(ktn.getHit()) nh++;
1079  }
1080  }
1081  }
1082  int nn=in;
1083  TString tagPlp=" "; if((nTpc-nh)>10) tagPlp=" plp";
1084  TString tagMemb=" "; if(zH*zL<0) tagMemb=" memb";
1085 
1086  LOG_INFO
1087  <<"\n#e dumpKalmanNodes nNodes="<<nn<<" actv: nTPC="<<nTpc<<" nHit="<<nh
1088  <<" zL="<<zL<<" zH="<<zH <<tagPlp<<tagMemb
1089  <<endm;
1090 
1091  // ........................print both ends ....................
1092  LOG_INFO << "#e |P|="<<track->getP()<<" pT="<<track->getPt()<<" eta="<<track->getPseudoRapidity()<<" nFitP="<<track->getFitPointCount() << endm;
1093  StiKalmanTrackNode* inNode=track->getInnerMostNode();
1094  LOG_INFO << "#e @InnerMostNode x:"<< inNode->x_g()<<" y:"<< inNode->y_g()<<" z:"<< inNode->z_g()<<" Eta="<<inNode->getEta()<<" |P|="<<inNode->getP() << endm;
1095  StiKalmanTrackNode* ouNode=track->getOuterMostNode();
1096  LOG_INFO << "#e @OuterMostNode g x:"<< ouNode->x_g()<<" y:"<< ouNode->y_g()<<" z:"<< ouNode->z_g()<<" Eta="<<ouNode->getEta()<<" |P|="<<ouNode->getP() << endm;
1097 
1098  in=0;
1099  for (it=track->begin();it!=track->end();it++,in++) {
1100  // if(in>=2 && in<nn-5) continue; // print only ends of the track
1101  StiKalmanTrackNode& ktn = (*it);
1102  if(!ktn.isValid()) continue;
1103  if(ktn.getHit() && ktn.getChi2() >1000) continue;
1104  float sy = std::sqrt(ktn.getCyy());
1105  float sz = std::sqrt(ktn.getCzz());
1106  const StiDetector * det=ktn.getDetector();
1107  assert(!(ktn.x()) || det);
1108 
1109  LOG_INFO << "#e in="<<in<<" |P|="<<ktn.getP()<<" Local: x="<<ktn.getX()<<" y="<<ktn.getY()<<" +/- "<<sy<<" z="<<ktn.getZ()<<" +/- "<<sz;
1110 
1111  if(ktn.getHit()) {LOG_INFO <<" hit=1";}
1112  else {LOG_INFO <<" hit=0";}
1113 
1114  if(det==0) {LOG_INFO <<" noDet ";}
1115  else {LOG_INFO <<" detActv="<<(!det || det->isActive(ktn.getY(), ktn.getZ()));}
1116  LOG_INFO << endm;
1117  }
1118 }
1119 
1120 //==========================================================
1121 //==========================================================
1122 bool StPPVertexFinder::examinTrackDca(const StiKalmanTrack* stiTrack, TrackData &track)
1123 {
1124  // .......... test DCA to beam .............
1125  StiKalmanTrackNode * bmNode = stiTrack->getInnerMostNode();
1126  if (!bmNode) return 0;
1127  if (!bmNode->isDca()) return 0;
1128 
1129  float rxy = std::sqrt(bmNode->x_g()*bmNode->x_g() + bmNode->y_g()*bmNode->y_g());
1130 
1131  if(rxy>mMaxTrkDcaRxy) return false;
1132  if( std::fabs(bmNode->z_g())> mMaxZrange ) return false ;
1133 
1134  track.zDca = bmNode->getZ();
1135  track.ezDca = std::sqrt(bmNode->getCzz());
1136  track.rxyDca = rxy;
1137  track.gPt = bmNode->getPt();
1138 
1139  return true;
1140 }
1141 
1142 
1143 //==========================================================
1144 //==========================================================
1145 void StPPVertexFinder::matchTrack2BTOF(const StiKalmanTrack* stiTrack, TrackData &track)
1146 {
1147  StBTofGeometry* geom = btofList->Geometry();
1148 
1149  StiKalmanTrackNode* ouNode=stiTrack->getOuterMostNode();
1150 
1151  StThreeVectorD posTOF;
1152  // helix extrapolation:
1153  StThreeVectorD ou(ouNode->getX(),ouNode->getY(),ouNode->getZ());
1154  ou.rotateZ(ouNode->getAlpha());
1155  StPhysicalHelixD phys_helix(std::fabs(ouNode->getCurvature()),
1156  ouNode->getDipAngle(),
1157  ouNode->getPhase(),
1158  ou,
1159  ouNode->getHelicity());
1160  IntVec idVec;
1161  DoubleVec pathVec;
1162  PointVec crossVec;
1163 
1164  IntVec iBinVec;
1165  if(geom->HelixCrossCellIds(phys_helix,idVec,pathVec,crossVec)) {
1166  for(size_t i=0;i<idVec.size();i++) {
1167  int itray, imodule, icell;
1168  geom->DecodeCellId(idVec[i], icell, imodule, itray);
1169 
1170  Double_t local[3], global[3];
1171  for(int j=0;j<3;j++) local[j] = 0;
1172  global[0] = crossVec[i].x();
1173  global[1] = crossVec[i].y();
1174  global[2] = crossVec[i].z();
1175  StBTofGeomSensor *sensor = geom->GetGeomSensor(imodule,itray);
1176  if(!sensor) {
1177  LOG_WARN << "no sensitive module in this projection??? - weird" << endm;
1178  continue;
1179  }
1180  sensor->Master2Local(&global[0],&local[0]);
1181  if(local[2]<mMinZBtof||local[2]>mMaxZBtof) continue;
1182  int iBin = btofList->cell2bin(itray, imodule, icell);
1183  iBinVec.push_back(iBin);
1184  btofList->addBtofTrack(itray, imodule, icell);
1185  LOG_DEBUG << " !!! Push to the list ...tray/module/cell " << itray << "/" << imodule << "/" << icell << endm;
1186  }
1187  }
1188 
1189  bool btofMatch=btofList->isMatched(iBinVec);
1190  bool btofVeto =btofList->isVetoed(iBinVec);
1191  float btofW =btofList->getWeight(iBinVec);
1192  btofList->addBtofMatch(iBinVec); // update the nMatch statistics
1193 
1194  LOG_DEBUG << " ** BTOF ** match/veto/weight = " << btofMatch << " " << btofVeto << " " << btofW << endm;
1195 
1196  track.updateAnyMatch(btofMatch,btofVeto,track.mBtof);
1197  track.weight*=btofW;
1198  track.btofBin= iBinVec.size() ? iBinVec[0] : -1;
1199 }
1200 
1201 //==========================================================
1202 //==========================================================
1203 void
1204 StPPVertexFinder::matchTrack2CTB(const StiKalmanTrack* stiTrack,TrackData &track){
1205  const double Rctb=213.6; // (cm) radius of the CTB
1206 
1207  StiKalmanTrackNode* ouNode=stiTrack->getOuterMostNode();
1208 
1209  StThreeVectorD posCTB;
1210  float path=-1;
1211  //alternative helix extrapolation:
1212  if(1){
1213  StiKalmanTrackNode * inNode = ouNode;
1214  StThreeVectorD in(inNode->getX(),inNode->getY(),inNode->getZ());
1215  in.rotateZ(inNode->getAlpha());
1216  StPhysicalHelixD phys_helix(std::fabs(inNode->getCurvature()),
1217  inNode->getDipAngle(),
1218  inNode->getPhase(),
1219  in,
1220  inNode->getHelicity());
1221  pairD d2;
1222  d2 = phys_helix.pathLength(Rctb);
1223  path=d2.second;
1224  if(d2.first>=0 || d2.second<=0) {
1225  LOG_DEBUG <<Form("WARN MatchTrk , unexpected solution for track crossing CTB\n")<<
1226  Form(" d2.firts=%f, second=%f, try first", d2.first, d2.second)<<endm;
1227  path=d2.first;
1228  }
1229  posCTB = phys_helix.at(path);
1230  }
1231 
1232  // official Sti node extrapolation
1233  if(0){
1234  StiKalmanTrack track2=*stiTrack;
1235  StiKalmanTrackNode* ctbNode=track2.extrapolateToRadius(Rctb);
1236 
1237  if(ctbNode==0) {
1238  LOG_INFO <<"#e @ctbNode NULL"<<endm;
1239  LOG_INFO <<"#e @track dump"<< *stiTrack << endm;
1240  LOG_INFO <<"#e @OuterMostNode dump"<< *ouNode <<endm;
1241  return;
1242  }
1243 
1244  posCTB=StThreeVectorD( ctbNode->x_g(),ctbNode->y_g(),ctbNode->z_g());
1245  }
1246 
1247  float phi=atan2(posCTB.y(),posCTB.x());
1248  if(phi<0) phi+=2*M_PI;// now phi is [0,2Pi] as for CTB slats
1249  float eta=posCTB.pseudoRapidity();
1250 
1251  int iBin=ctbList->addTrack(eta,phi);
1252 
1253  bool ctbMatch=ctbList->isMatched(iBin);
1254  bool ctbVeto =ctbList->isVetoed(iBin);
1255  float ctbW =ctbList->getWeight(iBin);
1256 
1257  track.updateAnyMatch(ctbMatch,ctbVeto,track.mCtb);
1258  track.weight*=ctbW;
1259  track.ctbBin=iBin;
1260 }
1261 
1262 //==========================================================
1263 //==========================================================
1264 void StPPVertexFinder::matchTrack2BEMC(TrackDataT<StiKalmanTrack> &track)
1265 {
1266  const StiKalmanTrack* stiTrack = track.getMother();
1267 
1268  StiKalmanTrackNode* ouNode=stiTrack->getOuterMostNode();
1269 
1270  //alternative helix extrapolation:
1271  StThreeVectorD ou(ouNode->getX(),ouNode->getY(),ouNode->getZ());
1272  ou.rotateZ(ouNode->getAlpha());
1273  StPhysicalHelixD phys_helix(std::fabs(ouNode->getCurvature()),
1274  ouNode->getDipAngle(),
1275  ouNode->getPhase(),
1276  ou,
1277  ouNode->getHelicity());
1278 
1279  matchTrack2BEMC(phys_helix, track);
1280 }
1281 
1282 
1283 void StPPVertexFinder::matchTrack2BEMC(TrackDataT<StMuTrack> &track)
1284 {
1285  const StMuTrack& muTrack = *track.getMother();
1286  matchTrack2BEMC(muTrack.outerHelix(), track);
1287 }
1288 
1289 
1290 void StPPVertexFinder::matchTrack2BEMC(const StPhysicalHelixD& phys_helix, TrackData &track)
1291 {
1292  const double Rxy = 242.; // middle of tower in Rxy
1293 
1294  pairD d2;
1295  d2 = phys_helix.pathLength(Rxy);
1296  float path = d2.second;
1297 
1298  if(d2.first>=0 || d2.second<=0) {
1299  LOG_DEBUG <<Form("WARN MatchTrk , unexpected solution for track crossing BEMC Cyl\n")<<
1300  Form(" d2.firts=%f, second=%f, try first\n", d2.first, d2.second)<<endm;
1301  path=d2.first;
1302  }
1303 
1304  StThreeVectorD posCyl = phys_helix.at(path);
1305 
1306  float phi = atan2(posCyl.y(), posCyl.x());
1307  if (phi < 0) phi += 2*M_PI; // now phi is [0,2Pi] as for Cyl slats
1308  float eta = posCyl.pseudoRapidity();
1309 
1310  int iBin = bemcList->addTrack(eta, phi);
1311  bool bemcMatch = bemcList->isMatched(iBin);
1312  bool bemcVeto = bemcList->isVetoed(iBin);
1313  float bemcW = bemcList->getWeight(iBin);
1314 
1315  track.updateAnyMatch(bemcMatch, bemcVeto, track.mBemc);
1316  track.weight *= bemcW;
1317  track.bemcBin = iBin;
1318 }
1319 
1320 
1321 //==========================================================
1322 //==========================================================
1323 void StPPVertexFinder::matchTrack2EEMC(TrackDataT<StiKalmanTrack> &track)
1324 {
1325  const StiKalmanTrack* stiTrack = track.getMother();
1326 
1327  const double minEta=0.7 ;// tmp cut
1328 
1329  StiKalmanTrackNode* ouNode=stiTrack->getOuterMostNode();
1330  StiKalmanTrackNode* inNode=stiTrack->getInnerMostNode();
1331 
1332  //direction of extrapolation must be toward West (Z+ axis)
1333  if(inNode->getZ()> ouNode->getZ()) return;
1334 
1335  // droop too steep tracks
1336  if(stiTrack->getPseudoRapidity()<minEta) return;
1337 
1338  StThreeVectorD ou(ouNode->getX(),ouNode->getY(),ouNode->getZ());
1339  ou.rotateZ(ouNode->getAlpha());
1340  StPhysicalHelixD phys_helix(std::fabs(ouNode->getCurvature()),
1341  ouNode->getDipAngle(),ouNode->getPhase(),
1342  ou,ouNode->getHelicity());
1343 
1344  // path length at intersection with plane
1345  // double pathLength(const StThreeVectorD& r,
1346  // const StThreeVectorD& n) const;
1347 
1348  matchTrack2EEMC(phys_helix, track);
1349 }
1350 
1351 
1352 void StPPVertexFinder::matchTrack2EEMC(TrackDataT<StMuTrack> &track)
1353 {
1354  const StMuTrack& muTrack = *track.getMother();
1355 
1356  const double minEta = 0.7;
1357 
1358  //direction of extrapolation must be toward West (Z+ axis)
1359  if (muTrack.firstPoint().z() > muTrack.lastPoint().z()) return;
1360 
1361  // drop too steep tracks
1362  if(muTrack.eta() < minEta) return;
1363 
1364  matchTrack2EEMC(muTrack.outerHelix(), track);
1365 }
1366 
1367 
1368 void StPPVertexFinder::matchTrack2EEMC(const StPhysicalHelixD& phys_helix, TrackData &track)
1369 {
1370  const double eemc_z_position = 288.; // middle of tower in Z
1371  const double maxPath=200 ;// tmp, cut too long extrapolation
1372 
1373  StThreeVectorD rSmd=StThreeVectorD(0,0,eemc_z_position);
1374  StThreeVectorD n=StThreeVectorD(0,0,1);
1375 
1376  double path = phys_helix.pathLength(rSmd,n);
1377  if(path>maxPath) return; // too long extrapolation
1378 
1379  StThreeVectorD r = phys_helix.at(path);
1380  double periodL=phys_helix. period();
1381 
1382  if(periodL<2*path) {
1383  LOG_DEBUG <<Form(" Warn, long path fac=%.1f ",path/periodL)<<
1384  Form(" punchEEMC1 x,y,z=%.1f, %.1f, %.1f path=%.1f period=%.1f\n",r.x(),r.y(),r.z(),path,periodL)<<endm;
1385  }
1386 
1387  double phi=atan2(r.y(),r.x());
1388  if(phi<0) phi+=2*M_PI;// now phi is [0,2Pi] as for Cyl slats
1389  double eta=r.pseudoRapidity();
1390 
1391  int iBin=eemcList->addTrack(eta,phi);
1392  bool eemcMatch=eemcList->isMatched(iBin);
1393  bool eemcVeto =eemcList->isVetoed(iBin);
1394  float eemcW =eemcList->getWeight(iBin);
1395 
1396  track.updateAnyMatch(eemcMatch,eemcVeto,track.mEemc);
1397  track.weight*=eemcW;
1398  track.eemcBin=iBin;
1399 
1400 }
1401 
1402 
1403 //==========================================================
1404 //==========================================================
1405 bool StPPVertexFinder::matchTrack2Membrane(TrackDataT<StiKalmanTrack> &track)
1406 {
1407  const StiKalmanTrack* stiTrack = track.getMother();
1408 
1409  const double RxyMin=59, RxyMax=199, zMax=200;
1410  const double zMembraneDepth=1; // (cm) ignore signe change for nodes so close to membrane
1411 
1412  //generate bitt pattern for TPC nodes with hits
1413  std::vector<int> hitPatt;
1414  int nPos=0,nFit=0;
1415  int in=0;
1416  double lastRxy=9999;
1417  double lastZ=9999;
1418 
1419  int jz0=0;
1421  for (it=stiTrack->begin();it!=stiTrack->end();it++) {
1422  StiKalmanTrackNode* ktnp=& (*it);
1423  if(!ktnp->isValid()) continue;
1424  //if(ktnp->getHit() && ktnp->getChi2() >1000) continue; // ---> those track need to be counted as npossiblehit, commented out
1425  double rxy = std::sqrt(ktnp->x_g()*ktnp->x_g() + ktnp->y_g()*ktnp->y_g()); //ktn.getX();
1426  double z=ktnp->z_g(); //ktn.z_g();
1427  if(rxy<RxyMin) continue;
1428  if(rxy>RxyMax) continue;
1429  if(std::fabs(z)>zMax) continue;
1430  // .........node is within TPC fiducial volume
1431  if(lastRxy<=rxy){
1432  LOG_WARN << "StPPVertexFinder::matchTrack2Membrane() \n the "<<in<<" node of the kalmanTrack below is out of order and is ignorred in (some) of vertex finder analysis"<<"\n Z="<<z<<" rXY="<<rxy<<" last_rxy="<<lastRxy<<endm;
1433  continue;
1434  }
1435  lastRxy=rxy;
1436  if(in==0) lastZ=z;
1437  in++;
1438  if(fabsf(z)>zMembraneDepth) { //ignore hits too close to z=0
1439  if(lastZ*z<0) { // track just crossed Z=0 plane
1440  if(jz0>0) {
1441  LOG_WARN << "StPPVertexFinder::matchTrack2Membrane() \n the "<<in<<" node of the kalmanTrack crosses Z=0 for the 2nd time, this track has a strange list of nodes - continue"<<endm;
1442  }
1443  //assert(jz0==0); // only one crosss point is expected
1444  jz0 = hitPatt.size();
1445  }
1446  lastZ=z;
1447  }
1448  const StiDetector * det=ktnp->getDetector();
1449  assert(!(ktnp->x()) || det);
1450  bool active = !det || det->isActive(ktnp->getY(), ktnp->getZ());
1451  int hit = ktnp->getHit() ? 1 : 0;
1452  if (active) {
1453  hitPatt.push_back(hit);
1454  nPos++;
1455  if(hit && ktnp->getChi2() <=1000 ) nFit++;
1456  }
1457  }
1458 
1459  if (nFit < mMinFitPfrac * nPos) return false; // too short fragment of a track
1460 
1461  if( mFitPossWeighting)
1462  track.weight *= 1.*nFit/nPos; // introduced in 2012 for pp510 to differentiate between global track quality, together with lowering the overall threshold from 0.7 to 0.51, Jan
1463 
1464  track.scanNodes(hitPatt, jz0); // if central membrane is crossed, scale weight inside
1465 
1466 
1467  return true;
1468 }
1469 
1470 
1471 void StPPVertexFinder::matchTrack2Membrane(TrackDataT<StMuTrack> &track)
1472 {
1473  const StMuTrack& muTrack = *track.getMother();
1474 
1475  // Code from matchTrack2Membrane
1476  if (mFitPossWeighting) { // introduced in 2012 for pp510 to differentiate between global track quality, together with lowering the overall threshold from 0.7 to 0.51
1477  double fracFit2PossHits = static_cast<double>(muTrack.nHitsFit(kTpcId)) / muTrack.nHitsPoss(kTpcId);
1478  track.weight *= fracFit2PossHits;
1479  }
1480 
1481  const StThreeVectorF& firstPoint = muTrack.firstPoint();
1482  const StThreeVectorF& lastPoint = muTrack.lastPoint();
1483 
1484  // Require the track to be within TPC volume (approximately)
1485  const float RxyMin = 59, RxyMax = 199, zMax = 200;
1486 
1487  bool isTrackInside = firstPoint.perp() < RxyMin || lastPoint.perp() < RxyMin ||
1488  firstPoint.perp() > RxyMax || lastPoint.perp() > RxyMax ||
1489  std::fabs(firstPoint.z()) > zMax ||
1490  std::fabs(lastPoint.z()) > zMax;
1491 
1492  // Require start and end to be on different sides of the z=0 plane
1493  bool crossMembrane = firstPoint.z() * lastPoint.z() < 0;
1494 
1495  if ( !isTrackInside || !crossMembrane) return;
1496 
1497  // Find crossing point of the track with z=0 membrane then identify jz0
1498  double t = firstPoint.z() / (firstPoint.z() - lastPoint.z());
1499 
1500  double intersect_x = firstPoint.x() + t * (lastPoint.x() - firstPoint.x());
1501  double intersect_y = firstPoint.y() + t * (lastPoint.y() - firstPoint.y());
1502  //intersect_z = firstPoint.z() + t * (lastPoint.z() - firstPoint.z()); // == 0 (=z)
1503 
1504  double intersect_r = std::sqrt(intersect_x*intersect_x + intersect_y*intersect_y);
1505 
1506  // TPC padrow radii
1507  const std::array<double, 45> padrow_radii
1508  {
1509  60.0, 64.8, 69.6, 74.4, 79.2, 84.0, 88.8, 93.6, 98.8, 104.0,
1510  109.2, 114.4, 119.6, 127.195, 129.195, 131.195, 133.195, 135.195, 137.195, 139.195,
1511  141.195, 143.195, 145.195, 147.195, 149.195, 151.195, 153.195, 155.195, 157.195, 159.195,
1512  161.195, 163.195, 165.195, 167.195, 169.195, 171.195, 173.195, 175.195, 177.195, 179.195,
1513  181.195, 183.195, 185.195, 187.195, 189.195
1514  };
1515 
1516  std::vector<int> hitPatt;
1517  int jz0 = 0;
1518 
1519  for (auto padrow_r : padrow_radii)
1520  {
1521  int curr_padrow_index = hitPatt.size();
1522 
1523  if (intersect_r > padrow_r)
1524  jz0 = curr_padrow_index;
1525 
1526  int hit = muTrack.topologyMap().hasHitInRow(kTpcId, curr_padrow_index+1) ? 1 : 0;
1527 
1528  hitPatt.push_back(hit);
1529  }
1530 
1531  track.scanNodes(hitPatt, jz0); // if central membrane is crossed, scale weight inside
1532 }
1533 
1534 
1538 bool StPPVertexFinder::isPostCrossingTrack(const StiKalmanTrack* stiTrack)
1539 {
1540  const float RxyMin=59, RxyMax=199, zMax=200;
1541  const float zMembraneDepth=1.0;
1542  const int nWrongZHitCut=2;
1543  int nWrongZHit=0;
1545  for (it=stiTrack->begin();it!=stiTrack->end();it++)
1546  {
1547  StiKalmanTrackNode* ktnp=& (*it);
1548  if(!ktnp->isValid() || ktnp->getChi2()>1000 ) continue;
1549 
1550  StiHit* stihit=ktnp->getHit();
1551 
1552  if (!stihit) continue;
1553 
1554  StHit* sthit=(StHit*)stihit->stHit();
1555 
1556  if (!sthit) continue;
1557  if (sthit->detector() != kTpcId) continue;
1558 
1559  StTpcHit* hit=(StTpcHit*) sthit;
1560  float r=hit->position().perp();
1561  if (r < RxyMin) continue;
1562  if (r > RxyMax) continue;
1563 
1564  float z=hit->position().z();
1565  if (std::fabs(z) > zMax) continue;
1566 
1567  if ((z < -zMembraneDepth && hit->sector() <= 12) ||
1568  (z > zMembraneDepth && hit->sector() > 12))
1569  {
1570  nWrongZHit++;
1571  if (nWrongZHit >= nWrongZHitCut) {return true;}
1572  }
1573  }
1574  return false;
1575 }
Bool_t HelixCrossCellIds(const StHelixD &helix, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec) const
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
Definition of Kalman Track.
double getP() const
Calculates and returns the momentum of the track at this node.
bool isTriggered
Indicates whether the vertex potentially belongs to triggered event.
Definition: VertexData.h:12
Definition: StHit.h:125
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
Abstract definition of a Track.
Definition: StiTrack.h:59
SeedFinder_t mSeedFinderType
The type of vertex seed finder to use in derived concrete implementation.
double beamX(double z) const
Definition: StiHit.h:51
StiKalmanTrackNode * getOuterMostNode(int qua=0) const
Accessor method returns the outer most node associated with the track.
double beamY(double z) const
double getP() const
Calculates and returns the momentum of the track at the inner most node.
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
int getFitPointCount(int detectorId=0) const
Returns the number of hits associated and used in the fit of this track.
std::vector< double > FindSeeds_TSpectrum()
double getPt() const
Calculates and returns the transverse momentum of the track at this node.
StiKalmanTrackNode * getInnerMostNode(int qua=0) const
Accessor method returns the inner most node associated with the track.
StThreeVectorD CalcVertexSeed(const StDcaList &trackDcas)
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
Definition of Kalman Track.
const StThreeVectorF & firstPoint() const
Returns positions of first measured point.
Definition: StMuTrack.h:261
double getPseudoRapidity() const
Return the pseudorapidity of the track.
StiKalmanTrackNode * extrapolateToBeam()
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
VertexFit_t mVertexFitMode
The type of vertex fit to use in derived concrete implementation.
static StTrack * createStTrack(const StMuTrack *)
creates a StTrack from an StMuTrack and return pointer to it
Definition: StMuDst.cxx:936
const StThreeVectorF & lastPoint() const
Returns positions of last measured point.
Definition: StMuTrack.h:262
Abstract interface for a STI toolkit.
const StiKTNBidirectionalIterator & end() const
double getCurvature() const
Calculates and returns the tangent of the track pitch angle at this node.
static StGenericVertexFinder * sSelf
By default point to invalid object.
const StMeasuredPoint * stHit() const
Definition: StiHit.h:104
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
StTrackTopologyMap topologyMap() const
Returns topology map.
Definition: StMuTrack.h:254
StiKTNBidirectionalIterator begin() const
double getPt() const
Calculates and returns the transverse momentum of the track at the inner most node.