00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 #include <assert.h>
00137 #include "StMinuitVertexFinder.h"
00138 #include "StEventTypes.h"
00139 #include "StEnumerations.h"
00140 #include "TMinuit.h"
00141 #if 0
00142 #include "TH1K.h"
00143 #include "TSpectrum.h"
00144 #endif
00145 #include "StGlobals.hh"
00146 #include "SystemOfUnits.h"
00147 #include "StCtbMatcher.h"
00148 #include "StMessMgr.h"
00149 #include <math.h>
00150 #include "StEmcUtil/geometry/StEmcGeom.h"
00151 #include "StDcaGeometry.h"
00152 #include "St_VertexCutsC.h"
00153 #include "StMaker.h"
00154 vector<StDcaGeometry*> StMinuitVertexFinder::mDCAs;
00155 vector<StPhysicalHelixD> StMinuitVertexFinder::mHelices;
00156 vector<UShort_t> StMinuitVertexFinder::mHelixFlags;
00157 vector<Double_t > StMinuitVertexFinder::mSigma;
00158 vector<Double_t > StMinuitVertexFinder::mZImpact;
00159 Double_t StMinuitVertexFinder::mWidthScale = 1;
00160 Double_t StMinuitVertexFinder::mX0;
00161 Double_t StMinuitVertexFinder::mY0;
00162 Double_t StMinuitVertexFinder::mdxdz;
00163 Double_t StMinuitVertexFinder::mdydz;
00164 Bool_t StMinuitVertexFinder::requireCTB;
00165 Int_t StMinuitVertexFinder::nCTBHits;
00166
00167
00168 void StMinuitVertexFinder::Clear(){
00169 StGenericVertexFinder::Clear();
00170 mStatusMin = 0;
00171 mNSeed = 0;
00172 }
00173
00174 void
00175 StMinuitVertexFinder::setExternalSeed(const StThreeVectorD& s)
00176 {
00177 mExternalSeedPresent = kTRUE;
00178 mExternalSeed = s;
00179 }
00180
00181
00182 StMinuitVertexFinder::StMinuitVertexFinder() {
00183 LOG_INFO << "StMinuitVertexFinder::StMinuitVertexFinder is in use." << endm;
00184 mBeamHelix =0;
00185
00186 mMinuit = new TMinuit(3);
00187 mMinuit->SetFCN(&StMinuitVertexFinder::fcn);
00188 mMinuit->SetPrintLevel(-1);
00189 mMinuit->SetMaxIterations(1000);
00190 mExternalSeedPresent = kFALSE;
00191 mVertexConstrain = kFALSE;
00192 mRequireCTB = kFALSE;
00193 requireCTB = kFALSE;
00194 mUseITTF = kFALSE;
00195 mUseOldBEMCRank = kFALSE;
00196 mLowerSplitVtxRank = kFALSE;
00197 mVertexOrderMethod = orderByRanking;
00198 mMinTrack = -1;
00199 mCTBSum = 0;
00200 }
00201
00202
00203 StMinuitVertexFinder::~StMinuitVertexFinder()
00204 {
00205 delete mBeamHelix; mBeamHelix=0;
00206 LOG_WARN << "Skipping delete Minuit in StMinuitVertexFinder::~StMinuitVertexFinder()" << endm;
00207
00208 mDCAs.clear();
00209 mHelices.clear();
00210 mHelixFlags.clear();
00211 mZImpact.clear();
00212 mSigma.clear();
00213 }
00214
00215 void StMinuitVertexFinder::InitRun(Int_t runumber) {
00216 St_VertexCutsC *cuts = St_VertexCutsC::instance();
00217 mMinNumberOfFitPointsOnTrack = cuts->MinNumberOfFitPointsOnTrack();
00218 mDcaZMax = cuts->DcaZMax();
00219 mMinTrack = (mMinTrack<0 ? cuts->MinTrack() : mMinTrack);
00220 mRImpactMax = cuts->RImpactMax();
00221 LOG_INFO << "Set cuts: MinNumberOfFitPointsOnTrack = " << mMinNumberOfFitPointsOnTrack
00222 << " DcaZMax = " << mDcaZMax
00223 << " MinTrack = " << mMinTrack
00224 << " RImpactMax = " << mRImpactMax
00225 << endm;
00226 }
00227
00228
00229
00230 void
00231 StMinuitVertexFinder::setFlagBase(){
00232 if(mUseITTF){
00233 mFlagBase = 8000;
00234 } else {
00235 mFlagBase = 1000;
00236 }
00237 }
00238
00239 Int_t StMinuitVertexFinder::findSeeds() {
00240 mNSeed = 0;
00241
00242 Int_t zImpactArr[400];
00243 for (Int_t i=0; i < 400; i++)
00244 zImpactArr[i]=0;
00245
00246 Int_t nTrk = mZImpact.size();
00247 for (Int_t iTrk=0; iTrk < nTrk; iTrk++) {
00248 if (fabs(mZImpact[iTrk]) < 200)
00249 zImpactArr[int(mZImpact[iTrk]+200)]++;
00250 }
00251
00252
00253 Int_t nOldBin = 0;
00254 Int_t slope = 0;
00255 Int_t nBinZ = 3;
00256 for (Int_t iBin=0; iBin < 400 - nBinZ; iBin++) {
00257 Int_t nTrkBin = 0;
00258 for (Int_t iBin2=0; iBin2 < nBinZ; iBin2++) {
00259 nTrkBin += zImpactArr[iBin + iBin2];
00260 }
00261 if (nTrkBin > nOldBin)
00262 slope = 1;
00263 else if (nTrkBin < nOldBin) {
00264 if (slope == 1) {
00265 if (mNSeed < maxSeed) {
00266 Float_t seed_z = -200 + iBin + (Float_t)nBinZ / 2 - 1;
00267 Double_t meanZ = 0;
00268 Int_t nTrkZ = 0;
00269 for (Int_t iTrk = 0; iTrk < nTrk; iTrk ++ ) {
00270 if ( fabs(mZImpact[iTrk] - seed_z ) < mDcaZMax ) {
00271 meanZ += mZImpact[iTrk];
00272 nTrkZ++;
00273 }
00274 }
00275 if (nTrkZ >= mMinTrack) {
00276 if (mDebugLevel) {
00277 LOG_INFO << "Seed " << mNSeed << ", z " << seed_z << " nTrk " << nTrkZ << " meanZ/nTrkZ " << meanZ/nTrkZ << endm;
00278 }
00279 seed_z = meanZ/nTrkZ;
00280 mSeedZ[mNSeed] = seed_z;
00281 mNSeed ++;
00282 }
00283 }
00284 else {
00285 LOG_WARN << "Too many vertex seeds, limit=" << maxSeed << endm;
00286 }
00287 }
00288 slope = -1;
00289 }
00290 if (mDebugLevel > 1) {
00291 LOG_INFO << "iBin " << iBin << " nTrkBin " << nTrkBin << " nOldBin " << nOldBin << ", slope " << slope << " mNSeed " << mNSeed << endm;
00292 }
00293 nOldBin = nTrkBin;
00294 }
00295
00296 LOG_INFO << "Found " << mNSeed << " seeds" << endm;
00297 return mNSeed;
00298 }
00299
00300 void StMinuitVertexFinder::fillBemcHits(StEvent *event){
00301 static Int_t nMod = 120;
00302 static Int_t nEta = 20;
00303 static Int_t nSub = 2;
00304 static Float_t mEmcThres = 0.15;
00305 for (Int_t m=0; m < nMod; m++)
00306 for (Int_t e=0; e < nEta; e++)
00307 for (Int_t s=0; s < nSub; s++)
00308 mBemcHit[m][e][s]=0;
00309
00310 Int_t n_emc_hit=0;
00311 if (event->emcCollection() && event->emcCollection()->detector(kBarrelEmcTowerId)) {
00312 StEmcDetector* bemcDet = event->emcCollection()->detector(kBarrelEmcTowerId);
00313 for (Int_t iMod=0; iMod < nMod; iMod++) {
00314 if (!bemcDet->module(iMod))
00315 continue;
00316 StEmcModule *mod = bemcDet->module(iMod);
00317 StSPtrVecEmcRawHit& hits = mod->hits();
00318 for (StEmcRawHitIterator hitIter = hits.begin(); hitIter != hits.end(); hitIter++) {
00319 StEmcRawHit *hit = *hitIter;
00320 if (hit->energy() > mEmcThres) {
00321 mBemcHit[hit->module()-1][hit->eta()-1][hit->sub()-1]=1;
00322 n_emc_hit++;
00323 }
00324 }
00325 }
00326 }
00327 if (mDebugLevel) {
00328 LOG_INFO << "Found " << n_emc_hit << " emc hits" << endm;
00329 }
00330 }
00331
00332 Int_t
00333 StMinuitVertexFinder::matchTrack2BEMC(const StTrack *track){
00334 static const Double_t rBemc = 242;
00335 static StEmcGeom *bemcGeom = StEmcGeom::getEmcGeom("bemc");
00336
00337
00338
00339 if (track->outerGeometry()==0) {
00340 if (mDebugLevel) {
00341 LOG_INFO << "No outer track geom" << endm;
00342 }
00343 return 0;
00344 }
00345
00346 StPhysicalHelixD helix = track->outerGeometry()->helix();
00347
00348 if (!helix.valid()) {
00349 if (mDebugLevel) {
00350 LOG_INFO << "Invalid helix" << endm;
00351 }
00352 return 0;
00353 }
00354
00355 pairD d2;
00356 d2 = helix.pathLength(rBemc);
00357 if (d2.first > 99999999 && d2.second > 99999999)
00358 return 0;
00359 Double_t path=d2.second;
00360 if (d2.first >= 0 && d2.first < d2.second)
00361 path = d2.first;
00362 else if(d2.first>=0 || d2.second<=0) {
00363 LOG_WARN << "Unexpected solution for track crossing Cyl:"
00364 << " track mom = "
00365 << track->geometry()->momentum().mag()
00366 << ", d2.first=" << d2.first
00367 << ", second=" << d2.second <<", trying first" << endm;
00368 path=d2.first;
00369 }
00370 StThreeVectorD posCyl = helix.at(path);
00371
00372 Float_t phi=posCyl.phi();
00373 Float_t eta=posCyl.pseudoRapidity();
00374
00375 if (fabs(eta) < 1) {
00376 Int_t mod, ieta, sub;
00377 if (bemcGeom->getBin(phi, eta, mod, ieta, sub) == 0) {
00378
00379
00380 if (sub > 0 && mBemcHit[mod-1][ieta-1][sub-1]) {
00381 return 1;
00382 }
00383 }
00384 }
00385 return 0;
00386 }
00387
00388 Int_t StMinuitVertexFinder::checkCrossMembrane(const StTrack *track) {
00389 Int_t n_pnt_neg=0, n_pnt_pos=0;
00390
00391 StPtrVecHit hits = track->detectorInfo()->hits(kTpcId);
00392 for (StHitIterator hitIter = hits.begin(); hitIter != hits.end(); hitIter++) {
00393 if ((*hitIter)->position().z() < 0)
00394 n_pnt_neg++;
00395 if ((*hitIter)->position().z() > 0)
00396 n_pnt_pos++;
00397 }
00398 return (n_pnt_pos > 5 && n_pnt_neg > 5);
00399 }
00400
00401 void StMinuitVertexFinder::calculateRanks() {
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 Int_t nBemcMatchTot = 0;
00424 Int_t nVtxTrackTot = 0;
00425 for (Int_t iVertex=0; iVertex < size(); iVertex++) {
00426 StPrimaryVertex *primV = getVertex(iVertex);
00427 nVtxTrackTot += primV->numTracksUsedInFinder();
00428 nBemcMatchTot += primV->numMatchesWithBEMC();
00429 }
00430
00431 mBestRank = -999;
00432 mBestVtx = 0;
00433 for (Int_t iVertex=0; iVertex < size(); iVertex++) {
00434 StPrimaryVertex *primV = getVertex(iVertex);
00435
00436 Float_t avg_dip_expected = -0.0033*primV->position().z();
00437 Float_t n_bemc_expected = 0;
00438 if (nVtxTrackTot) {
00439 if (mUseOldBEMCRank)
00440 n_bemc_expected = (1-0.25*(1-(float)primV->numTracksUsedInFinder()/nVtxTrackTot))*nBemcMatchTot;
00441 else
00442 n_bemc_expected = (float)primV->numTracksUsedInFinder()/nVtxTrackTot*nBemcMatchTot;
00443 }
00444
00445 Float_t n_cross_expected = fabs(primV->position().z())*0.0020*primV->numTracksUsedInFinder();
00446
00447 if (mDebugLevel) {
00448 LOG_INFO << "vertex z " << primV->position().z() << " dip expected " << avg_dip_expected << " bemc " << n_bemc_expected << " cross " << n_cross_expected << endm;
00449 }
00450 Float_t rank_avg_dip = 1 - fabs(primV->meanDip() - avg_dip_expected)*sqrt((float)primV->numTracksUsedInFinder())/0.67;
00451 if (rank_avg_dip < -5)
00452 rank_avg_dip = -5;
00453
00454 Float_t rank_bemc = 0;
00455 if (n_bemc_expected >= 1) {
00456
00457 Float_t sigma = 0.5*sqrt(n_bemc_expected);
00458 if ( sigma < 0.75 ) {
00459
00460 sigma = 0.75;
00461 }
00462 rank_bemc = (primV->numMatchesWithBEMC() - n_bemc_expected)/sigma;
00463 if (mUseOldBEMCRank)
00464 rank_bemc += 0.5;
00465 }
00466 if (rank_bemc < -5)
00467 rank_bemc = -5;
00468 if (rank_bemc > 1)
00469 rank_bemc = 1;
00470
00471 Float_t rank_cross = 0;
00472 if ( n_cross_expected >= 1 ) {
00473 Float_t sigma=1.1*sqrt(n_cross_expected);
00474 rank_cross = (primV->numTracksCrossingCentralMembrane() - n_cross_expected)/sigma;
00475 }
00476 if (rank_cross < -5)
00477 rank_cross = -5;
00478 if (rank_cross > 1)
00479 rank_cross = 1;
00480
00481 if (mDebugLevel) {
00482 LOG_INFO << "rankings: " << rank_avg_dip << " " << rank_bemc << " " << rank_cross << endm;
00483 }
00484
00485
00486 if (mLowerSplitVtxRank &&
00487 (mCTBSum > ((6000.0*(float)primV->numTracksUsedInFinder()/80)+2000)))
00488 primV->setRanking(rank_cross+rank_bemc+rank_avg_dip-3);
00489 else
00490 primV->setRanking(rank_cross+rank_bemc+rank_avg_dip);
00491 if (primV->ranking() > mBestRank) {
00492 mBestRank = primV->ranking();
00493 mBestVtx = primV;
00494 }
00495 }
00496 }
00497
00498 int
00499 StMinuitVertexFinder::fit(StEvent* event)
00500 {
00501 Double_t arglist[4];
00502
00503 setFlagBase();
00504
00505
00506 StCtbTriggerDetector* ctbDet = 0;
00507 vector<ctbHit> ctbHits;
00508
00509 StTriggerDetectorCollection* trigCol = event->triggerDetectorCollection();
00510 mCTBSum = 0;
00511 if(trigCol){
00512 ctbDet = &(trigCol->ctb());
00513
00514 for (UInt_t slat = 0; slat < ctbDet->numberOfSlats(); slat++) {
00515 for (UInt_t tray = 0; tray < ctbDet->numberOfTrays(); tray++) {
00516 ctbHit curHit;
00517 curHit.adc = ctbDet->mips(tray,slat,0);
00518 if(curHit.adc > 0){
00519 mCTBSum += curHit.adc;
00520 ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
00521 ctbHits.push_back(curHit);
00522 }
00523 }
00524 }
00525 }
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 mDCAs.clear();
00537 mHelices.clear();
00538 mHelixFlags.clear();
00539 mSigma.clear();
00540 mZImpact.clear();
00541
00542 fillBemcHits(event);
00543
00544 Bool_t ctb_match;
00545
00546 Int_t n_ctb_match_tot = 0;
00547 Int_t n_bemc_match_tot = 0;
00548 Int_t n_cross_tot = 0;
00549
00550 StSPtrVecTrackNode& nodes = event->trackNodes();
00551 UInt_t Nnodes = nodes.size();
00552 for (UInt_t k = 0; k < Nnodes; k++) {
00553 StGlobalTrack* g = ( StGlobalTrack*) nodes[k]->track(global);
00554 if (accept(g)) {
00555 mWidthScale = 0.1;
00556 StDcaGeometry* gDCA = g->dcaGeometry();
00557 if (! gDCA) continue;
00558 if (TMath::Abs(gDCA->impact()) > mRImpactMax) continue;
00559 mDCAs.push_back(gDCA);
00560
00561
00562 mHelices.push_back(g->geometry()->helix());
00563 mHelixFlags.push_back(1);
00564 Double_t z_lin = gDCA->z();
00565 mZImpact.push_back(z_lin);
00566 mSigma.push_back(-1);
00567
00568 Bool_t shouldHitCTB = kFALSE;
00569 Double_t etaInCTBFrame = -999;
00570 ctb_match = EtaAndPhiToOrriginAtCTB(g,&ctbHits,shouldHitCTB,etaInCTBFrame);
00571 if (ctb_match) {
00572 mHelixFlags[mHelixFlags.size()-1] |= kFlagCTBMatch;
00573 n_ctb_match_tot++;
00574 }
00575
00576 if (matchTrack2BEMC(g)) {
00577 mHelixFlags[mHelixFlags.size()-1] |= kFlagBEMCMatch;
00578 n_bemc_match_tot++;
00579 }
00580
00581 if (checkCrossMembrane(g)) {
00582 mHelixFlags[mHelixFlags.size()-1] |= kFlagCrossMembrane;
00583 n_cross_tot++;
00584 }
00585 }
00586 }
00587 if (mDebugLevel) {
00588 LOG_INFO << "Found " << n_ctb_match_tot << " ctb matches, " << n_bemc_match_tot << " bemc matches, " << n_cross_tot << " tracks crossing central membrane" << endm;
00589 }
00590
00591
00592
00593 if (mHelices.empty()) {
00594 LOG_WARN << "StMinuitVertexFinder::fit: no tracks to fit." << endm;
00595 mStatusMin = -1;
00596 return 0;
00597 }
00598 LOG_INFO << "StMinuitVertexFinder::fit size of helix vector: " << mHelices.size() << endm;
00599
00600
00601 if (mRequireCTB) requireCTB = kTRUE;
00602
00603
00604
00605
00606 mMinuit->mnexcm("CLEar", 0, 0, mStatusMin);
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 static Double_t step[3] = {0.03, 0.03, 0.03};
00618
00619
00620
00621
00622
00623 if (!mExternalSeedPresent) {
00624 if (!mVertexConstrain){
00625 arglist[0] = 3;
00626 }
00627 else {
00628 arglist[0]=1;
00629 }
00630 findSeeds();
00631 }
00632 else {
00633 mNSeed = 1;
00634 }
00635
00636 Float_t old_vtx_z = -999;
00637 Double_t seed_z = -999;
00638 Double_t chisquare = 0;
00639 for (Int_t iSeed = 0; iSeed < mNSeed; iSeed++) {
00640
00641
00642
00643 mMinuit->mnexcm("CLEar", 0, 0, mStatusMin);
00644
00645 seed_z= mSeedZ[iSeed];
00646
00647 if (mExternalSeedPresent)
00648 seed_z = mExternalSeed.z();
00649 if (mDebugLevel) {
00650 LOG_INFO << "Vertex seed = " << seed_z << endm;
00651 }
00652
00653 if (!mVertexConstrain){
00654 mMinuit->mnparm(0, "x", 0, step[0], 0, 0, mStatusMin);
00655 mMinuit->mnparm(1, "y", 0, step[1], 0, 0, mStatusMin);
00656 mMinuit->mnparm(2, "z", seed_z, step[2], 0, 0, mStatusMin);
00657 }
00658 else {
00659 mMinuit->mnparm(0, "z", seed_z, step[2], 0, 0, mStatusMin);
00660 }
00661
00662 if (!mVertexConstrain){
00663 arglist[0] = 3;
00664 }
00665 else {
00666 arglist[0] = 1;
00667 }
00668
00669 Int_t done = 0;
00670 Int_t iter = 0;
00671
00672 Int_t n_trk_vtx = 0;
00673 Int_t n_helix = mHelices.size();
00674 do {
00675
00676
00677 n_trk_vtx = 0;
00678 for (Int_t i=0; i < n_helix; i++) {
00679 if (fabs(mZImpact[i]-seed_z) < mDcaZMax) {
00680 mHelixFlags[i] |= kFlagDcaz;
00681 n_trk_vtx++;
00682 }
00683 else
00684 mHelixFlags[i] &= ~kFlagDcaz;
00685 }
00686
00687 if (mDebugLevel) {
00688 LOG_INFO << n_trk_vtx << " tracks within dcaZ cut (iter " << iter <<" )" << endm;
00689 }
00690 if (n_trk_vtx < mMinTrack) {
00691 if (mDebugLevel) {
00692 LOG_INFO << "Less than mMinTrack (=" << mMinTrack << ") tracks, skipping vtx" << endm;
00693 }
00694 continue;
00695 }
00696 mMinuit->mnexcm("MINImize", 0, 0, mStatusMin);
00697 done = 1;
00698
00699
00700
00701
00702
00703 if (mStatusMin) {
00704 LOG_WARN << "StMinuitVertexFinder::fit: error in Minuit::mnexcm(), check status flag. ( iter=" << iter << ")" << endm;
00705 done = 0;
00706 }
00707
00708 Double_t fedm, errdef;
00709 Int_t npari, nparx;
00710
00711 mMinuit->mnstat(chisquare, fedm, errdef, npari, nparx, mStatusMin);
00712
00713 if (mStatusMin != 3) {
00714 LOG_INFO << "Warning: Minuit Status: " << mStatusMin << ", func val " << chisquare<< endm;
00715 done = 0;
00716 }
00717 mMinuit->mnhess();
00718
00719 Double_t new_z, zerr;
00720 if (!mVertexConstrain) {
00721 mMinuit->GetParameter(2, new_z, zerr);
00722 }
00723 else {
00724 mMinuit->GetParameter(0, new_z, zerr);
00725 }
00726
00727 if (fabs(new_z - seed_z) > 1)
00728 done = 0;
00729
00730 Int_t n_trk = 0;
00731 for (Int_t i=0; i < n_helix; i++) {
00732 if ( fabs(mZImpact[i] - new_z) < mDcaZMax ) {
00733 n_trk++;
00734 }
00735 }
00736 if ( 10 * abs(n_trk - n_trk_vtx) >= n_trk_vtx )
00737 done = 0;
00738
00739 iter++;
00740 seed_z = new_z;
00741 } while (!done && iter < 5 && n_trk_vtx >= mMinTrack);
00742
00743 if (n_trk_vtx < mMinTrack)
00744 continue;
00745
00746 if (!done) {
00747 LOG_WARN << "Vertex unstable: no convergence after " << iter << " iterations. Skipping vertex " << endm;
00748 continue;
00749 }
00750
00751 if (!mExternalSeedPresent && fabs(seed_z-mSeedZ[iSeed]) > mDcaZMax) {
00752 LOG_WARN << "Vertex walks during fits: seed was " << mSeedZ[iSeed] << ", vertex at " << seed_z << endm;
00753 }
00754
00755 if (fabs(seed_z - old_vtx_z) < mDcaZMax) {
00756 if (mDebugLevel) {
00757 LOG_INFO << "Vertices too close (<mDcaZMax). Skipping" << endm;
00758 }
00759 continue;
00760 }
00761
00762
00763 StThreeVectorD XVertex;
00764 Float_t cov[6];
00765 memset(cov,0,sizeof(cov));
00766
00767 Double_t val, verr;
00768 if (!mVertexConstrain) {
00769 XVertex = StThreeVectorD(mMinuit->fU[0],mMinuit->fU[1],mMinuit->fU[2]);
00770 Double_t emat[9];
00771
00772
00773
00774 mMinuit->mnemat(emat,3);
00775 cov[0] = emat[0];
00776 cov[1] = emat[3];
00777 cov[2] = emat[4];
00778 cov[3] = emat[6];
00779 cov[4] = emat[7];
00780 cov[5] = emat[8];
00781 }
00782 else {
00783 mMinuit->GetParameter(0, val, verr);
00784 XVertex.setZ(val); cov[5]=verr*verr;
00785
00786
00787
00788 XVertex.setX(beamX(val)); cov[0]=0.1;
00789 XVertex.setY(beamY(val)); cov[2]=0.1;
00790 }
00791 StPrimaryVertex primV;
00792 primV.setPosition(XVertex);
00793 primV.setChiSquared(chisquare);
00794 primV.setCovariantMatrix(cov);
00795 primV.setVertexFinderId(minuitVertexFinder);
00796 primV.setFlag(1);
00797 primV.setRanking(333);
00798 primV.setNumTracksUsedInFinder(n_trk_vtx);
00799
00800 Int_t n_ctb_match = 0;
00801 Int_t n_bemc_match = 0;
00802 Int_t n_cross = 0;
00803 n_trk_vtx = 0;
00804
00805 Double_t mean_dip = 0;
00806 Double_t sum_pt = 0;
00807 for (UInt_t i=0; i < mHelixFlags.size(); i++) {
00808 if (!(mHelixFlags[i] & kFlagDcaz))
00809 continue;
00810 n_trk_vtx++;
00811 if (mHelixFlags[i] & kFlagCTBMatch)
00812 n_ctb_match++;
00813 if (mHelixFlags[i] & kFlagBEMCMatch)
00814 n_bemc_match++;
00815 if (mHelixFlags[i] & kFlagCrossMembrane)
00816 n_cross++;
00817 mean_dip += mHelices[i].dipAngle();
00818 sum_pt += mHelices[i].momentum(0).perp();
00819 }
00820
00821 mean_dip /= n_trk_vtx;
00822
00823 if (mDebugLevel) {
00824 LOG_INFO << "check n_trk_vtx " << n_trk_vtx << ", found " << n_ctb_match << " ctb matches, " << n_bemc_match << " bemc matches, " << n_cross << " tracks crossing central membrane" << endm;
00825 LOG_INFO << "mean dip " << mean_dip << endm;
00826 }
00827 primV.setNumMatchesWithCTB(n_ctb_match);
00828 primV.setNumMatchesWithBEMC(n_bemc_match);
00829 primV.setNumTracksCrossingCentralMembrane(n_cross);
00830 primV.setMeanDip(mean_dip);
00831 primV.setSumOfTrackPt(sum_pt);
00832
00833
00834 addVertex(&primV);
00835
00836 old_vtx_z = XVertex.z();
00837 }
00838
00839 calculateRanks();
00840
00841
00842
00843 mExternalSeedPresent = kFALSE;
00844
00845 requireCTB = kFALSE;
00846
00847 return 1;
00848 }
00849
00850 Double_t StMinuitVertexFinder::Chi2atVertex(StThreeVectorD &vtx) {
00851 static Int_t nCall=0; nCall++;
00852 Double_t f = 0;
00853 Double_t e;
00854 nCTBHits = 0;
00855 if (fabs(vtx.x())> 10) return 1e6;
00856 if (fabs(vtx.y())> 10) return 1e6;
00857 if (fabs(vtx.z())>300) return 1e6;
00858 for (UInt_t i=0; i<mDCAs.size(); i++) {
00859 if ((mHelixFlags[i] & kFlagDcaz) && (!requireCTB||(mHelixFlags[i] & kFlagCTBMatch))) {
00860 const StDcaGeometry* gDCA = mDCAs[i];
00861 if (! gDCA) continue;
00862 const StPhysicalHelixD helix = gDCA->helix();
00863 e = helix.distance(vtx, kFALSE);
00864
00865
00866 static Int_t nCall=0;
00867 nCall++;
00868 Double_t err2;
00869 Double_t chi2 = gDCA->thelix().Dca(&(vtx.x()),&err2);
00870 chi2*=chi2/err2;
00871
00872 Double_t scale = 1./(mWidthScale*mWidthScale);
00873 f += scale*(1. - TMath::Exp(-chi2/scale));
00874
00875 if((mHelixFlags[i] & kFlagCTBMatch) && e<3.0) nCTBHits++;
00876 }
00877 }
00878 return f;
00879 }
00880
00881
00882 void StMinuitVertexFinder::fcn1D(int& npar, double* gin, double& f, double* par, Int_t iflag)
00883 {
00884 Double_t z = par[0];
00885 Double_t x = beamX(z);
00886 Double_t y = beamY(z);
00887 StThreeVectorD vtx(x,y,z);
00888 f = Chi2atVertex(vtx);
00889 }
00890 void StMinuitVertexFinder::fcn(int& npar, double* gin, double& f, double* par, Int_t iflag)
00891 {
00892 StThreeVectorD vtx(par);
00893 f = Chi2atVertex(vtx);
00894 }
00895
00896
00897 bool
00898 StMinuitVertexFinder::accept(StTrack* track) const
00899 {
00900
00901
00902
00903
00904
00905 return (track &&
00906 track->flag() >= 0 &&
00907 track->fitTraits().numberOfFitPoints() >= mMinNumberOfFitPointsOnTrack &&
00908 !track->topologyMap().trackFtpc() &&
00909 finite(track->length()) &&
00910 track->geometry()->helix().valid());
00911 }
00912
00913
00915 void
00916 StMinuitVertexFinder::setPrintLevel(Int_t level)
00917 {
00918 mMinuit->SetPrintLevel(level);
00919 }
00920
00921 void
00922 StMinuitVertexFinder::printInfo(ostream& os) const
00923 {
00924
00925 os << "StMinuitVertexFinder - Statistics:" << endl;
00926 os << "Number of vertices found ......." << size() << endl;
00927 os << "Rank of best vertex ............" << mBestRank << endl;
00928 if(mBestVtx) {
00929 os << "Properties of best vertex:" << endl;
00930 os << "position ..................... " << mBestVtx->position() << endl;
00931 os << "position errors .............. " << mBestVtx->positionError()<< endl;
00932 os << "# of used tracks ............. " << mBestVtx->numTracksUsedInFinder() << endl;
00933 os << "Chisquare .................... " << mBestVtx->chiSquared() << endl;
00934 }
00935 os << "min # of fit points for tracks . " << mMinNumberOfFitPointsOnTrack << endl;
00936 os << "final potential width scale .... " << mWidthScale << endl;
00937 }
00938
00939 void StMinuitVertexFinder::UseVertexConstraint(Double_t x0, Double_t y0, Double_t dxdz, Double_t dydz, Double_t weight) {
00940 mVertexConstrain = kTRUE;
00941 mX0 = x0;
00942 mY0 = y0;
00943 mdxdz = dxdz;
00944 mdydz = dydz;
00945 mWeight = weight;
00946 LOG_INFO << "StMinuitVertexFinder::Using Constrained Vertex" << endm;
00947 LOG_INFO << "x origin = " << mX0 << endm;
00948 LOG_INFO << "y origin = " << mY0 << endm;
00949 LOG_INFO << "slope dxdz = " << mdxdz << endm;
00950 LOG_INFO << "slope dydz = " << mdydz << endm;
00951 LOG_INFO << "weight in fit = " << weight << endm;
00952 StThreeVectorD origin(mX0,mY0,0.0);
00953 Double_t pt = 88889999;
00954 Double_t nxy=::sqrt(mdxdz*mdxdz + mdydz*mdydz);
00955 if(nxy<1.e-5){
00956 LOG_WARN << "StMinuitVertexFinder:: Beam line must be tilted!" << endm;
00957 nxy=mdxdz=1.e-5;
00958 }
00959 Double_t p0=pt/nxy;
00960 Double_t px = p0*mdxdz;
00961 Double_t py = p0*mdydz;
00962 Double_t pz = p0;
00963 StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
00964 delete mBeamHelix;
00965 mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
00966
00967
00968 mMinuit = new TMinuit(1);
00969 mMinuit->SetFCN(&StMinuitVertexFinder::fcn1D);
00970 mMinuit->SetPrintLevel(1);
00971 mMinuit->SetMaxIterations(1000);
00972 mExternalSeedPresent = kFALSE;
00973
00974
00975 }
00976
00977
00978 Double_t StMinuitVertexFinder::beamX(Double_t z) {
00979 Float_t x = mX0 + mdxdz*z;
00980 return x;
00981 }
00982
00983 Double_t StMinuitVertexFinder::beamY(Double_t z) {
00984 Float_t y = mY0 + mdydz*z;
00985 return y;
00986 }
00987
00988 Int_t StMinuitVertexFinder::NCtbMatches() {
00989 return nCTBHits;
00990 }
00991 Int_t StMinuitVertexFinder::NCtbSlats() {
00992 return -777;
00993 }
00994
00995
00996
00997
00998
