13 #include "StEmcUtil/geometry/StEmcGeom.h"
14 #include "StEvent/StCtbTriggerDetector.h"
15 #include "StEvent/StDcaGeometry.h"
16 #include "StEvent/StEmcCollection.h"
17 #include "StEvent/StEmcDetector.h"
18 #include "StEvent/StEmcModule.h"
19 #include "StEvent/StEmcRawHit.h"
20 #include "StEvent/StEvent.h"
21 #include "StEvent/StGlobalTrack.h"
22 #include "StEvent/StTrackDetectorInfo.h"
23 #include "StEvent/StTrackNode.h"
24 #include "StEvent/StTriggerDetectorCollection.h"
25 #include "StGenericVertexMaker/Minuit/StMinuitVertexFinder.h"
26 #include "StGenericVertexMaker/Minuit/St_VertexCutsC.h"
27 #include "StGenericVertexMaker/StCtbMatcher.h"
28 #include "St_base/StMessMgr.h"
31 std::vector<StPhysicalHelixD> StMinuitVertexFinder::mHelices;
32 std::vector<UShort_t> StMinuitVertexFinder::mHelixFlags;
33 std::vector<Double_t > StMinuitVertexFinder::mZImpact;
34 Bool_t StMinuitVertexFinder::requireCTB;
35 Int_t StMinuitVertexFinder::nCTBHits;
38 void StMinuitVertexFinder::Clear(){
39 StGenericVertexFinder::Clear();
47 mExternalSeedPresent = kTRUE;
52 StMinuitVertexFinder::StMinuitVertexFinder(VertexFit_t fitMode) :
55 mExternalSeedPresent = kFALSE;
59 mUseOldBEMCRank = kFALSE;
60 mLowerSplitVtxRank = kFALSE;
67 StMinuitVertexFinder::~StMinuitVertexFinder()
69 LOG_WARN <<
"Skipping delete Minuit in StMinuitVertexFinder::~StMinuitVertexFinder()" << endm;
75 void StMinuitVertexFinder::InitRun(
int run_number,
const St_db_Maker* db_maker)
77 StGenericVertexFinder::InitRun(run_number, db_maker);
80 mMinNumberOfFitPointsOnTrack = cuts->MinNumberOfFitPointsOnTrack();
81 mDcaZMax = cuts->DcaZMax();
82 mMinTrack = (mMinTrack<0 ? cuts->MinTrack() : mMinTrack);
83 mRImpactMax = cuts->RImpactMax();
91 LOG_INFO <<
"Set cuts: MinNumberOfFitPointsOnTrack = " << mMinNumberOfFitPointsOnTrack
92 <<
" DcaZMax = " << mDcaZMax
93 <<
" MinTrack = " << mMinTrack
94 <<
" RImpactMax = " << mRImpactMax
95 <<
" ZMin = " << mZMin
96 <<
" ZMax = " << mZMax
103 StMinuitVertexFinder::setFlagBase(){
111 Int_t StMinuitVertexFinder::findSeeds() {
114 int zIdxMax = (int) (mZMax - mZMin);
115 Int_t zImpactArr[zIdxMax];
116 for (Int_t i=0; i < zIdxMax; i++)
119 Int_t nTrk = mZImpact.size();
120 for (Int_t iTrk=0; iTrk < nTrk; iTrk++) {
121 if ((mZImpact[iTrk] > mZMin) &&
122 (mZImpact[iTrk] < mZMax))
123 zImpactArr[int(mZImpact[iTrk]-mZMin)]++;
130 for (Int_t iBin=0; iBin < zIdxMax - nBinZ; iBin++) {
132 for (Int_t iBin2=0; iBin2 < nBinZ; iBin2++) {
133 nTrkBin += zImpactArr[iBin + iBin2];
135 if (nTrkBin > nOldBin)
137 else if (nTrkBin < nOldBin) {
139 if (mNSeed < maxSeed) {
140 Float_t seed_z = mZMin + iBin + (Float_t)nBinZ / 2 - 1;
143 for (Int_t iTrk = 0; iTrk < nTrk; iTrk ++ ) {
144 if ( fabs(mZImpact[iTrk] - seed_z ) < mDcaZMax ) {
145 meanZ += mZImpact[iTrk];
149 if (nTrkZ >= mMinTrack) {
151 LOG_INFO <<
"Seed " << mNSeed <<
", z " << seed_z <<
" nTrk " << nTrkZ <<
" meanZ/nTrkZ " << meanZ/nTrkZ << endm;
153 seed_z = meanZ/nTrkZ;
154 mSeedZ[mNSeed] = seed_z;
159 LOG_WARN <<
"Too many vertex seeds, limit=" << maxSeed << endm;
164 if (mDebugLevel > 1) {
165 LOG_INFO <<
"iBin " << iBin <<
" nTrkBin " << nTrkBin <<
" nOldBin " << nOldBin <<
", slope " << slope <<
" mNSeed " << mNSeed << endm;
170 LOG_INFO <<
"Found " << mNSeed <<
" seeds" << endm;
174 void StMinuitVertexFinder::fillBemcHits(
StEvent *event){
175 static Int_t nMod = 120;
176 static Int_t nEta = 20;
177 static Int_t nSub = 2;
178 static Float_t mEmcThres = 0.15;
179 for (Int_t m=0; m < nMod; m++)
180 for (Int_t e=0; e < nEta; e++)
181 for (Int_t s=0; s < nSub; s++)
185 if (event->emcCollection() &&
event->emcCollection()->detector(kBarrelEmcTowerId)) {
186 StEmcDetector* bemcDet =
event->emcCollection()->detector(kBarrelEmcTowerId);
187 for (Int_t iMod=0; iMod < nMod; iMod++) {
188 if (!bemcDet->module(iMod))
191 StSPtrVecEmcRawHit& hits = mod->hits();
192 for (StEmcRawHitIterator hitIter = hits.begin(); hitIter != hits.end(); hitIter++) {
194 if (hit->energy() > mEmcThres) {
195 mBemcHit[hit->module()-1][hit->eta()-1][hit->sub()-1]=1;
202 LOG_INFO <<
"Found " << n_emc_hit <<
" emc hits" << endm;
207 StMinuitVertexFinder::matchTrack2BEMC(
const StTrack *
track){
208 static const Double_t rBemc = 242;
209 static StEmcGeom *bemcGeom = StEmcGeom::getEmcGeom(
"bemc");
213 if (track->outerGeometry()==0) {
215 LOG_INFO <<
"No outer track geom" << endm;
222 if (!helix.
valid()) {
224 LOG_INFO <<
"Invalid helix" << endm;
231 if (d2.first > 99999999 && d2.second > 99999999)
233 Double_t path=d2.second;
234 if (d2.first >= 0 && d2.first < d2.second)
236 else if(d2.first>=0 || d2.second<=0) {
237 LOG_WARN <<
"Unexpected solution for track crossing Cyl:"
239 << track->geometry()->momentum().mag()
240 <<
", d2.first=" << d2.first
241 <<
", second=" << d2.second <<
", trying first" << endm;
246 Float_t phi=posCyl.phi();
247 Float_t eta=posCyl.pseudoRapidity();
250 Int_t mod, ieta, sub;
251 if (bemcGeom->
getBin(phi, eta, mod, ieta, sub) == 0) {
254 if (sub > 0 && mBemcHit[mod-1][ieta-1][sub-1]) {
262 Int_t StMinuitVertexFinder::checkCrossMembrane(
const StTrack *track) {
263 Int_t n_pnt_neg=0, n_pnt_pos=0;
265 StPtrVecHit hits = track->detectorInfo()->hits(kTpcId);
266 for (StHitIterator hitIter = hits.begin(); hitIter != hits.end(); hitIter++) {
267 if ((*hitIter)->position().z() < 0)
269 if ((*hitIter)->position().z() > 0)
272 return (n_pnt_pos > 5 && n_pnt_neg > 5);
275 void StMinuitVertexFinder::calculateRanks() {
297 Int_t nBemcMatchTot = 0;
298 Int_t nVtxTrackTot = 0;
299 for (Int_t iVertex=0; iVertex < size(); iVertex++) {
301 nVtxTrackTot += primV->numTracksUsedInFinder();
302 nBemcMatchTot += primV->numMatchesWithBEMC();
307 for (Int_t iVertex=0; iVertex < size(); iVertex++) {
310 Float_t avg_dip_expected = -0.0033*primV->position().z();
311 Float_t n_bemc_expected = 0;
314 n_bemc_expected = (1-0.25*(1-(float)primV->numTracksUsedInFinder()/nVtxTrackTot))*nBemcMatchTot;
316 n_bemc_expected = (float)primV->numTracksUsedInFinder()/nVtxTrackTot*nBemcMatchTot;
319 Float_t n_cross_expected = fabs(primV->position().z())*0.0020*primV->numTracksUsedInFinder();
322 LOG_INFO <<
"vertex z " << primV->position().z() <<
" dip expected " << avg_dip_expected <<
" bemc " << n_bemc_expected <<
" cross " << n_cross_expected << endm;
325 Float_t rank_avg_dip = 1 - fabs(primV->meanDip() - avg_dip_expected)*sqrt((
float)primV->numTracksUsedInFinder())/0.67;
327 Float_t rank_bemc = 0;
328 if (n_bemc_expected >= 1) {
330 Float_t sigma = 0.5*sqrt(n_bemc_expected);
333 sigma = ( sigma < 0.75 ? 0.75 : sigma);
335 rank_bemc = (primV->numMatchesWithBEMC() - n_bemc_expected)/sigma;
340 Float_t rank_cross = 0;
341 if ( n_cross_expected >= 1 ) {
342 Float_t sigma=1.1*sqrt(n_cross_expected);
343 rank_cross = (primV->numTracksCrossingCentralMembrane() - n_cross_expected)/sigma;
347 rank_avg_dip = ( rank_avg_dip < -5 ? -5 : rank_avg_dip );
348 rank_bemc = ( rank_bemc < -5 ? -5 : (rank_bemc > 1 ? 1 : rank_bemc) );
349 rank_cross = ( rank_cross < -5 ? -5 : (rank_cross > 1 ? 1 : rank_cross) );
352 LOG_INFO <<
"rankings: " << rank_avg_dip <<
" " << rank_bemc <<
" " << rank_cross << endm;
356 if (mLowerSplitVtxRank && mCTBSum > 6000.0*primV->numTracksUsedInFinder()/80. + 2000)
357 primV->setRanking(rank_cross+rank_bemc+rank_avg_dip-3);
359 primV->setRanking(rank_cross+rank_bemc+rank_avg_dip);
361 if (primV->ranking() > mBestRank) {
362 mBestRank = primV->ranking();
369 int StMinuitVertexFinder::fit(
StEvent* event)
375 std::vector<ctbHit> ctbHits;
380 ctbDet = &(trigCol->ctb());
382 for (UInt_t slat = 0; slat < ctbDet->numberOfSlats(); slat++) {
383 for (UInt_t tray = 0; tray < ctbDet->numberOfTrays(); tray++) {
385 curHit.adc = ctbDet->mips(tray,slat,0);
387 mCTBSum += curHit.adc;
388 ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
389 ctbHits.push_back(curHit);
406 Int_t n_ctb_match_tot = 0;
407 Int_t n_bemc_match_tot = 0;
408 Int_t n_cross_tot = 0;
410 for (
const StTrackNode* stTrack : event->trackNodes())
413 if (!accept(g))
continue;
415 if (! gDCA)
continue;
416 if (TMath::Abs(gDCA->impact()) > mRImpactMax)
continue;
417 mDCAs.push_back(gDCA);
420 mHelices.push_back(g->geometry()->helix());
421 mHelixFlags.push_back(1);
422 Double_t z_lin = gDCA->z();
423 mZImpact.push_back(z_lin);
425 Bool_t shouldHitCTB = kFALSE;
426 Double_t etaInCTBFrame = -999;
427 bool ctb_match = EtaAndPhiToOrriginAtCTB(g,&ctbHits,shouldHitCTB,etaInCTBFrame);
429 mHelixFlags[mHelixFlags.size()-1] |= kFlagCTBMatch;
433 if (matchTrack2BEMC(g)) {
434 mHelixFlags[mHelixFlags.size()-1] |= kFlagBEMCMatch;
438 if (checkCrossMembrane(g)) {
439 mHelixFlags[mHelixFlags.size()-1] |= kFlagCrossMembrane;
445 LOG_INFO <<
"Found " << n_ctb_match_tot <<
" ctb matches, " << n_bemc_match_tot <<
" bemc matches, " << n_cross_tot <<
" tracks crossing central membrane" << endm;
449 if (mHelices.empty()) {
450 LOG_WARN <<
"StMinuitVertexFinder::fit: no tracks to fit." << endm;
455 LOG_INFO <<
"StMinuitVertexFinder::fit size of helix vector: " << mHelices.size() << endm;
458 if (mRequireCTB) requireCTB = kTRUE;
461 mMinuit->mnexcm(
"CLEar", 0, 0, mStatusMin);
473 static Double_t step[3] = {0.03, 0.03, 0.03};
477 if (!mExternalSeedPresent) {
484 Float_t old_vtx_z = -999;
485 Double_t seed_z = -999;
486 Double_t chisquare = 0;
488 for (Int_t iSeed = 0; iSeed < mNSeed; iSeed++) {
491 mMinuit->mnexcm(
"CLEar", 0, 0, mStatusMin);
493 seed_z= mSeedZ[iSeed];
495 if (mExternalSeedPresent)
496 seed_z = mExternalSeed.z();
498 LOG_INFO <<
"Vertex seed = " << seed_z << endm;
502 mMinuit->mnparm(0,
"z", seed_z, step[2], 0, 0, mStatusMin);
505 mMinuit->mnparm(0,
"x", 0, step[0], 0, 0, mStatusMin);
506 mMinuit->mnparm(1,
"y", 0, step[1], 0, 0, mStatusMin);
507 mMinuit->mnparm(2,
"z", seed_z, step[2], 0, 0, mStatusMin);
514 Int_t n_helix = mHelices.size();
518 for (Int_t i=0; i < n_helix; i++) {
519 if (fabs(mZImpact[i]-seed_z) < mDcaZMax) {
520 mHelixFlags[i] |= kFlagDcaz;
524 mHelixFlags[i] &= ~kFlagDcaz;
528 LOG_INFO << n_trk_vtx <<
" tracks within dcaZ cut (iter " << iter <<
" )" << endm;
530 if (n_trk_vtx < mMinTrack) {
532 LOG_INFO <<
"Less than mMinTrack (=" << mMinTrack <<
") tracks, skipping vtx" << endm;
536 mMinuit->mnexcm(
"MINImize", 0, 0, mStatusMin);
541 LOG_WARN <<
"StMinuitVertexFinder::fit: error in Minuit::mnexcm(), check status flag. ( iter=" << iter <<
")" << endm;
545 Double_t fedm, errdef;
548 mMinuit->mnstat(chisquare, fedm, errdef, npari, nparx, mStatusMin);
550 if (mStatusMin != 3) {
551 LOG_INFO <<
"Warning: Minuit Status: " << mStatusMin <<
", func val " << chisquare<< endm;
556 Double_t new_z, zerr;
558 mMinuit->GetParameter(0, new_z, zerr);
561 mMinuit->GetParameter(2, new_z, zerr);
564 if (fabs(new_z - seed_z) > 1)
568 for (Int_t i=0; i < n_helix; i++) {
569 if ( fabs(mZImpact[i] - new_z) < mDcaZMax ) {
573 if ( 10 * abs(n_trk - n_trk_vtx) >= n_trk_vtx )
578 }
while (!done && iter < 5 && n_trk_vtx >= mMinTrack);
580 if (n_trk_vtx < mMinTrack)
584 LOG_WARN <<
"Vertex unstable: no convergence after " << iter <<
" iterations. Skipping vertex " << endm;
588 if (!mExternalSeedPresent && fabs(seed_z-mSeedZ[iSeed]) > mDcaZMax) {
589 LOG_WARN <<
"Vertex walks during fits: seed was " << mSeedZ[iSeed] <<
", vertex at " << seed_z << endm;
592 if (fabs(seed_z - old_vtx_z) < mDcaZMax) {
594 LOG_INFO <<
"Vertices too close (<mDcaZMax). Skipping" << endm;
602 memset(cov,0,
sizeof(cov));
606 mMinuit->GetParameter(0, val, verr);
607 XVertex.setZ(val); cov[5]=verr*verr;
611 XVertex.setX(
beamX(val)); cov[0]=0.1;
612 XVertex.setY(
beamY(val)); cov[2]=0.1;
615 XVertex =
StThreeVectorD(mMinuit->fU[0],mMinuit->fU[1],mMinuit->fU[2]);
620 mMinuit->mnemat(emat,3);
629 primV.setPosition(XVertex);
630 primV.setChiSquared(chisquare);
631 primV.setCovariantMatrix(cov);
632 primV.setVertexFinderId(minuitVertexFinder);
634 primV.setRanking(333);
635 primV.setNumTracksUsedInFinder(n_trk_vtx);
637 Int_t n_ctb_match = 0;
638 Int_t n_bemc_match = 0;
642 Double_t mean_dip = 0;
644 for (UInt_t i=0; i < mHelixFlags.size(); i++) {
645 if (!(mHelixFlags[i] & kFlagDcaz))
648 if (mHelixFlags[i] & kFlagCTBMatch)
650 if (mHelixFlags[i] & kFlagBEMCMatch)
652 if (mHelixFlags[i] & kFlagCrossMembrane)
654 mean_dip += mHelices[i].dipAngle();
655 sum_pt += mHelices[i].momentum(0).perp();
658 mean_dip /= n_trk_vtx;
661 LOG_INFO <<
"check n_trk_vtx " << n_trk_vtx <<
", found "
662 << n_ctb_match <<
" ctb matches, "
663 << n_bemc_match <<
" bemc matches, "
664 << n_cross <<
" tracks crossing central membrane\n"
665 <<
"mean dip " << mean_dip << endm;
667 primV.setNumMatchesWithCTB(n_ctb_match);
668 primV.setNumMatchesWithBEMC(n_bemc_match);
669 primV.setNumTracksCrossingCentralMembrane(n_cross);
670 primV.setMeanDip(mean_dip);
671 primV.setSumOfTrackPt(sum_pt);
676 old_vtx_z = XVertex.z();
683 mExternalSeedPresent = kFALSE;
692 double StMinuitVertexFinder::CalcChi2DCAs(
const StThreeVectorD &vtx) {
696 if (fabs(vtx.x())> 10)
return 1e6;
697 if (fabs(vtx.y())> 10)
return 1e6;
698 if (fabs(vtx.z())>300)
return 1e6;
699 for (UInt_t i=0; i<
mDCAs.size(); i++) {
701 if ( !(mHelixFlags[i] & kFlagDcaz) || (requireCTB && !(mHelixFlags[i] & kFlagCTBMatch)) )
705 if (! gDCA)
continue;
711 Double_t chi2 = gDCA->thelix().
Dca(&(vtx.x()),&err2);
714 static double scale = 100;
715 f += scale*(1. - TMath::Exp(-chi2/scale));
717 if((mHelixFlags[i] & kFlagCTBMatch) && e<3.0) nCTBHits++;
726 bool StMinuitVertexFinder::accept(
StTrack* track)
const
729 track->flag() >= 0 &&
730 track->fitTraits().numberOfFitPoints() >= mMinNumberOfFitPointsOnTrack &&
731 !track->topologyMap().trackFtpc() &&
732 finite(track->length()) &&
733 track->geometry()->helix().
valid());
740 mMinuit->SetPrintLevel(level);
744 void StMinuitVertexFinder::printInfo(ostream& os)
const
746 os <<
"StMinuitVertexFinder - Statistics:" << endl;
747 os <<
"Number of vertices found ......." << size() << endl;
748 os <<
"Rank of best vertex ............" << mBestRank << endl;
750 os <<
"Properties of best vertex:" << endl;
751 os <<
"position ..................... " << mBestVtx->position() << endl;
752 os <<
"position errors .............. " << mBestVtx->positionError()<< endl;
753 os <<
"# of used tracks ............. " << mBestVtx->numTracksUsedInFinder() << endl;
754 os <<
"Chisquare .................... " << mBestVtx->chiSquared() << endl;
756 os <<
"min # of fit points for tracks . " << mMinNumberOfFitPointsOnTrack << endl;
760 Int_t StMinuitVertexFinder::NCtbMatches() {
bool valid(double world=1.e+5) const
checks for valid parametrization
double beamX(double z) const
double beamY(double z) const
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
void setPrintLevel(Int_t=0)
Use mMinuit print level.
VertexFit_t mVertexFitMode
The type of vertex fit to use in derived concrete implementation.
static StGenericVertexFinder * sSelf
By default point to invalid object.
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
double Dca(const double point[3], double *dcaErr=0) const
DCA to given space point (with error matrix)