StRoot: StGenericVertexMaker/StppLMVVertexFinder.cxx Source File

00001 /************************************************************
00002  *
00003  * $Id: StppLMVVertexFinder.cxx,v 1.24 2010/01/26 21:01:49 fisyak Exp $
00004  *
00005  * Author: Jan Balewski
00006  ************************************************************
00007  *
00008  * Description: 
00009  *
00010  ************************************************************/
00011 
00012 #include <StEventTypes.h>
00013 #include <StEnumerations.h>
00014 #include <StGlobals.hh>
00015 #include <SystemOfUnits.h>
00016 #include <StMessMgr.h>
00017 #include <cmath>
00018 #include "math_constants.h"
00019 
00020 #include "tables/St_g2t_vertex_Table.h" // tmp for Dz(vertex)
00021 #include "StMaker.h"
00022 
00023 #include "StppLMVVertexFinder.h"
00024 #include "StGenericVertexMaker.h"
00025 
00026 //==========================================================
00027 //==========================================================
00028 StppLMVVertexFinder::StppLMVVertexFinder() {
00029   LOG_INFO << "StppLMVVertexFinder::StppLMVVertexFinder is in use." << endm;
00030   mBeamHelix           = 0;
00031   mVertexConstrain     = false;
00032   mTotEve              = 0;
00033   mdxdz=mdydz=mX0=mY0  = 0;
00034   mMode                = 1;
00035 }
00036 
00037 //==========================================================
00038 //==========================================================
00039 void 
00040 StppLMVVertexFinder::Init() {
00041 
00042   //jan default cuts
00043   if (mMode==0){
00044     LOG_INFO << "The  ppLMV4 cuts have been activated" << endm; 
00045     mMaxTrkDcaRxy    = 3.9;
00046     mMinTrkPt        = 0.2;
00047     mMinNumberOfFitPointsOnTrack = 15; // a typo (=10) was here before , JB
00048     mMaxZrange       = 180;            // for tracks;  typo (=250) was here before , JB
00049     mDVtxMax         = 4.0;            // max sigma multipl between tracks and current vertex, used for tracks rejection
00050     mMinMatchTr      = 1;              // minimal # of tracks matched to CTB 
00051     mBLequivNtr      = 100;            // equivalent # of tracks for BeamLine
00052     mMatchCtbMax_eta = mCtbEtaSeg/2.+0.02;
00053     mMatchCtbMax_phi = mCtbPhiSeg/2.+C_PI*1./180.;
00054   } else {
00055     LOG_INFO << "The  ppLMV5 cuts have been activated" << endm; 
00056     mMaxTrkDcaRxy    = 2.0;            
00057     mMinTrkPt        = 0.2;            
00058     mMinNumberOfFitPointsOnTrack = 15; 
00059     mMaxZrange       = 150;            
00060     mDVtxMax         = 4.0;            
00061     mMinMatchTr      = 2;              
00062     mBLequivNtr      = 20;             
00063     mMatchCtbMax_eta = mCtbEtaSeg/2.+0.02;
00064     mMatchCtbMax_phi = mCtbPhiSeg/2.+C_PI*0./180.;
00065   }
00066 }
00067 
00068 //==========================================================
00069 //==========================================================
00070 void 
00071 StppLMVVertexFinder::Clear(){  //  Reset vertex
00072   StGenericVertexFinder::Clear();
00073   mCtbHits.clear();
00074   mPrimCand.clear();
00075   LOG_INFO << "StppLMVVertexFinder::Clear() ..." << endm;
00076 }
00077 
00078 //==========================================================
00079 //==========================================================
00080 void 
00081 StppLMVVertexFinder::addFakeVerex(float z){
00082   StPrimaryVertex primV;
00083   Float_t cov[6] = {1,0.0,2,0.0,0.0,3  }; 
00084   StThreeVectorD XVertex(0.1,0.2,z);
00085   primV.setPosition(XVertex);
00086   primV.setCovariantMatrix(cov); 
00087   primV.setVertexFinderId(pplmvVertexFinder);
00088   primV.setFlag(1); 
00089   primV.setRanking(444);
00090   //..... add vertex to the list
00091   addVertex(&primV);
00092 }
00093 
00094 //==========================================================
00095 //==========================================================
00096 StppLMVVertexFinder::~StppLMVVertexFinder() {
00097    delete mBeamHelix;mBeamHelix=0;
00098 }
00099 
00100 //==========================================================
00101 //==========================================================
00102 int 
00103 StppLMVVertexFinder::fit(StEvent* event) {
00104   LOG_INFO << "StppLMVVertexFinder::fit() START ..." << endm;
00105  //
00106 
00107   n1=0,n2=0,n3=0, n4=0,n5=0,n6=0;
00108   //  mStatus = -888;
00109 
00110   mTotEve++;
00111   eveID=event->id();
00112 
00113   mBfield = event->runInfo()->magneticField(); 
00114   LOG_INFO << Form("ppLMV5:: mBfield[Tesla]=%e ",mBfield)<<endm;
00115   
00116 
00117   gMessMgr->Message("","I") << "ppLMV5::cuts"
00118                             <<"\n CtbThres_ch (real)="<<mCtbThres_ch
00119                             <<"\n CtbThres_mev (M-C)="<<mCtbThres_mev
00120                             <<"\n MinNumberOfFitPointsOnTrack="<<mMinNumberOfFitPointsOnTrack
00121                             <<"\n MaxTrkDcaRxy/cm="<<mMaxTrkDcaRxy
00122                             <<"\n MinTrkPt GeV/c ="<<mMinTrkPt
00123                             <<"\n MatchCtbMax_eta="<<mMatchCtbMax_eta
00124                             <<"\n MatchCtbMax_phi/deg="<<mMatchCtbMax_phi/3.1416*180
00125                             <<"\n BeamLequivNtr="<<mBLequivNtr
00126                             <<"\n DVtxMax (dz/sig)="<<mDVtxMax
00127                             <<"\n MinMatchTr ="<< mMinMatchTr
00128                             <<"\n MaxZrange (cm) ="<< mMaxZrange
00129                             <<"\n VertexConstrain="<<mVertexConstrain
00130     //                      <<"\n  ="<<
00131                             <<endm; 
00132    
00133   // get CTB info, does not  work for embeding 
00134   if(event->runId()<100000){
00135     St_DataSet *gds=StMaker::GetChain()->GetDataSet("geant");
00136     collectCTBhitsMC(gds);  // use M-C
00137   }  else {
00138     StTriggerData *trgD=event->triggerData ();
00139     collectCTBhitsData(trgD); // use real data
00140   }
00141   // printCtb(); //4ppv
00142 
00143   if(mCtbHits.size()<=0){
00144     LOG_WARN << "StppLMVVertexFinder::fit() no valid CTB hits found, quit" << endm;
00145     return 0;
00146   }   
00147    
00148   //  Loop all global tracks (TPC) and store the
00149   //  refering helices and their estimated DCA
00150   //  resolution in  vectors.
00151     
00152   StSPtrVecTrackNode& nodes = event->trackNodes();
00153   for (unsigned int k=0; k<nodes.size(); k++) {
00154     StTrack* g = nodes[k]->track(global);
00155     float sigma=0;
00156     if( !matchTrack2CTB(g,sigma)) continue;;    
00157     JHelix x;
00158     x.helix=g->geometry()->helix();
00159     x.sigma=sigma;
00160     mPrimCand.push_back(x);
00161     // printf("added tr %d w/ sigma=%f\n",mPrimCand.size(),sigma);
00162          
00163   } // end of track selection
00164 
00165   //  printf(", now n1=%d n2=%d n3=%d n4=%d  n6=%d\n",n1,n2,n3,n4,n6); //tmp
00166   LOG_DEBUG << ", now n1=" << n1 << " n2=" << n2 << " n3=" << n3 << " n4=" << n4 << "n6=" << n6 << endm;
00167 
00168 
00169   //
00170   //  In case there are no tracks left we better quit
00171   //
00172   if (mPrimCand.size() < mMinMatchTr){
00173     LOG_WARN << "StppLMVVertexFinder::fit() only "<< mPrimCand.size() << "tracks to fit - too few, quit ppLMV" << endm;
00174     return 0;
00175   }
00176   LOG_INFO << "StppLMVVertexFinder::fit() size of helix vector: " << mPrimCand.size() << endm;
00177  
00178   // use beam line constraint or not
00179   if(mVertexConstrain ) {
00180 
00181     float sigMin=100000; // pick large weight for this track 
00182     for( uint j=0;j<mPrimCand.size();j++) {
00183       float sig=mPrimCand[j].sigma;
00184       if(sigMin>sig) sigMin=sig;
00185     }
00186     float sigma=sigMin/::sqrt(mBLequivNtr); //<== assigne relative weight
00187     JHelix x;
00188     x.helix=*mBeamHelix;
00189     x.sigma=sigma;
00190     mPrimCand.push_back(x);
00191     //printf("WARN ppLMV: nominal beam line added with sigma=%f, now nTrack=%d \n",sigma, mPrimCand.size());
00192     LOG_WARN << "ppLMV: nominal beam line added with sigma=" << sigma 
00193                         << ", now nTrack=" << mPrimCand.size() << endm;
00194 
00195   } // end of beamLine
00196 
00197 
00198   //printf("PrimCand  before ppLMV for eveID=%d tracks at first point\n",eveID);
00199   LOG_DEBUG << "PrimCand  before ppLMV for eveID=" << eveID << "tracks at first point" << endm;
00200 
00201   //tmp - lits all matched tracks
00202   for( uint j=0;j<mPrimCand.size();j++) {
00203     StThreeVectorD p=mPrimCand[j].helix.momentum(mBfield*tesla);
00204     printf("j=%d  sig=%f pT=%f eta=%f phi/deg=%f\n",j,mPrimCand[j].sigma,p.perp(),p.pseudoRapidity(), p.phi()/C_PI*180);    
00205   } 
00206 
00207   //  ----------  D O   F I N D    V E R T E X
00208   int ret=ppLMV5();
00209   if(ret==false) {
00210     LOG_DEBUG << "ppLMV did not found vertex"<< endm;
00211     return 0;
00212   }
00213   
00214 
00215   LOG_INFO << "Prim tr used by ppLMV for eveID=" << eveID 
00216                    << ", tracks at vertex=" << mPrimCand.size()<<endm;
00217 
00218   assert(size());
00219   for( uint j=0;j<mPrimCand.size();j++) {
00220     double spath = mPrimCand[j].helix.pathLength( getVertex(0)->position().x(),getVertex(0)->position().y());
00221     StThreeVectorD p=mPrimCand[j].helix.momentumAt(spath,mBfield*tesla);
00222     // printf("j=%d  sig=%f pT=%f eta=%f phi/deg=%f cur=%f p=%f charg=%d spath=%f\n",j,mPrimCand[j].sigma,p.perp(),p.pseudoRapidity(), p.phi()/C_PI*180,mPrimCand[j].helix.curvature(), p.mag(),mPrimCand[j].helix.charge(mBfield*tesla),spath);
00223   }
00224          
00225   gMessMgr->Message("","I") << "Prim ppLMV Vertex at " << getVertex(0)->position()<<endm;
00226   return size();
00227 } 
00228 
00229 
00230 //======================================================
00231 //======================================================
00232 void
00233 StppLMVVertexFinder::printInfo(ostream& os) const
00234 {
00235     os << "StppLMVVertexFinder - Fit Statistics:" << endl;
00236     os << "found prim vertices ................ " <<size() << endl;
00237     os << "no more info printed at the moment, Jan"<<endl;
00238  }
00239 
00240 
00241 //======================================================
00242 //======================================================
00243 void 
00244 StppLMVVertexFinder::UseVertexConstraint(double x0, double y0, 
00245                                          double dxdz, double dydz, double weight) {
00246   mVertexConstrain = true;
00247   mX0 = x0;
00248   mY0 = y0;
00249   mdxdz = dxdz;
00250   mdydz = dydz;
00251   mWeight = weight;
00252   LOG_INFO << "StppLMVVertexFinder::Using Constrained Vertex" << endm;
00253   LOG_INFO << "x origin = " << mX0 << endm;
00254   LOG_INFO << "y origin = " << mY0 << endm;
00255   LOG_INFO << "slope dxdz = " << mdxdz << endm;
00256   LOG_INFO << "slope dydz = " << mdydz << endm;
00257   LOG_INFO << "NOT used (JB) weight in fit = " << weight <<  endm;
00258   StThreeVectorD origin(mX0,mY0,0.0);
00259   double pt  = 88889999;   
00260   double nxy=::sqrt(mdxdz*mdxdz +  mdydz*mdydz);
00261   if(nxy<1.e-5){ // beam line _MUST_ be tilted
00262     LOG_WARN << "StppLMVVertexFinder:: Beam line must be tilted!" << endm;
00263     nxy=mdxdz=1.e-5; 
00264   }
00265   double p0=pt/nxy;  
00266   double px   = p0*mdxdz;
00267   double py   = p0*mdydz;
00268   double pz   = p0; // approximation: nx,ny<<0
00269   StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
00270   delete mBeamHelix;
00271   mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
00272 }
00273 
00274 
00275 
00276 
00277 //==========================================================
00278 //==========================================================
00279 bool 
00280 StppLMVVertexFinder::matchTrack2CTB (StTrack* track, float & sigma) {
00281   /* upgrade:
00282      - used dE/dX for strag
00283      - use track length  
00284      - make clear()
00285      - check for total energy per slat for  geant
00286   */
00287 
00288   sigma=0;
00289   const double Rctb=213.6; // (cm) radius of the CTB 
00290   uint nPoss=0, nFitP=0, nSvtP=0;
00291 
00292   if (!track) return false; // it should never happen
00293   if(!finite(track->geometry()->helix().curvature())){
00294     LOG_WARN << "StppLMVVertexFinder::matchTrack2CTB: helix.curvature not finite, fix tracker, ppLMV aborting" << endm;
00295     mCtbHits.clear();
00296     mPrimCand.clear();    
00297     return false;
00298   }
00299 
00300   if( track->flag()<0) return false;
00301   if( track->topologyMap().trackFtpc() )return false;
00302 
00303   n1++; 
00304 
00305   StPhysicalHelixD TrkHlxIn=track->geometry()->helix();
00306 
00307   //           check Rxy_min condition  close to beam    
00308   double spath = TrkHlxIn.pathLength(mX0, mY0 );
00309   StThreeVectorD posDCA = TrkHlxIn.at(spath);
00310   //  cout<<" DCA Position: "<<posDCA<<endl;
00311   double x_m = posDCA.x(), y_m = posDCA.y();
00312   double dmin = ::sqrt(x_m*x_m + y_m*y_m);
00313   if( dmin > mMaxTrkDcaRxy ) return false;
00314   if (fabs(posDCA.z()) >mMaxZrange) return false;
00315   n2++;
00316   
00317   nFitP = track->detectorInfo()->numberOfPoints(kTpcId);
00318   nSvtP = track->detectorInfo()->numberOfPoints(kSvtId);
00319   nPoss=track->numberOfPossiblePoints(kTpcId);
00320 
00321   if(  nFitP <= mMinNumberOfFitPointsOnTrack) return false;
00322   n3++;
00323 
00324   if(track->geometry()->momentum().perp() <mMinTrkPt ) return false;
00325   n4++;
00326 
00327   //cout<<"\n\n DCA to beam at xy=00: "<<posDCA<<endl;
00328        
00329   //Find momentum direction at vertex point
00330   StThreeVectorD pmom = TrkHlxIn.momentumAt(spath,mBfield*tesla );
00331   double beta = pmom.mag()/::sqrt(pmom.mag()*pmom.mag()+0.139*0.139); //Assume pion 
00332 
00333   // old formula from TPT chain
00334   // float strag=0.0136/beta/pmom.mag()*fabs(spath); 
00335 
00336   // ppLMV face-lift : attenuated weights 
00337   float pmomM=pmom.mag();
00338   float spathL=fabs(spath); // reduce advantage of SVT matched tracks
00339 
00340   if(mMode ==1){
00341     // This is mode 1 (ppLMV5)
00342     if (pmomM >4 ) pmomM=4; //inhibit domination of high pT tracks
00343     if( spathL<40) spathL=40;
00344   } else {
00345     // This is mode 0 (ppLMV4)
00346     // ignore SVT contribution , ~aproximately
00347     if( spathL<60) spathL=60;
00348   } 
00349 
00350   float strag=0.0136/beta/pmomM*spathL; 
00351   if(fabs(mBfield)<0.01) strag=0.0136*spathL; // no field case, pT makes no sense
00352 
00353   //  printf("stragling=%f %f p=%f %f nFp=%d nPp=%d\n",strag,beta,pmom.mag(),spath, track->fitTraits().numberOfFitPoints(),nPoss );
00354   
00355 
00356   if ( !track->outerGeometry() )  return false;
00357   StPhysicalHelixD TrkHlxOut=track->outerGeometry()->helix();
00358   pairD  d2 = TrkHlxOut.pathLength(Rctb);
00359   //  printf(" path 1=%f, 2=%f, period=%f, R=%f\n",d2.first ,d2.second,TrkHlx.period(),1./TrkHlxOut.curvature());
00360   
00361   // assert(d2.first<0); // propagate backwards
00362   //assert(d2.second>0); // propagate forwards
00363   
00364   if(d2.first>=0 || d2.second<=0) {
00365     n5++;
00366     gMessMgr->Message("","W") << "ppLMV::matchTrack2CTB ()"  
00367       " unexpected solution for track crossing CTB, TrkHlxOut="<< TrkHlxOut<<endm;
00368     return false;
00369   }
00370   
00371   
00372   StThreeVectorD posCTB = TrkHlxOut.at(d2.second);
00373   //  double xmagn = ::sqrt( posCTB.x()*posCTB.x() + posCTB.y()*posCTB.y() );//tmp, out
00374   // printf(" punch2 x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f\n",posCTB.x(),posCTB.y(),posCTB.z(),xmagn);
00375   
00376   float phi=atan2(posCTB.y(),posCTB.x());
00377   if(phi<0) phi+=2*C_PI;// now phi is [0,2Pi] as for CTB slats
00378   
00379   // printf("posCTB.z()=%f posDCA.z()=%f\n",posCTB.z(),posDCA.z());
00380 
00381   uint ih;
00382   for(ih=0;ih<mCtbHits.size();ih++) {// loop over CTB hits
00383  
00384    // match to CTB slats in phi
00385     float del_phi=phi-mCtbHits[ih].phi;
00386     if(del_phi>C_PI) del_phi-=2*C_PI;
00387     if(del_phi<-C_PI) del_phi+=2*C_PI;
00388     //printf("phiRad trk=%f  CTB=%f del=%f\n",phi,mCtbHits[ih].phi,del_phi);
00389     //printf("match ih=%d del_phi=%f/deg\n",ih, del_phi/C_PI*180);
00390     if(fabs(del_phi) >mMatchCtbMax_phi) continue;
00391     
00392     // match to CTB slats in eta
00393     float eta=posCTB.pseudoRapidity();
00394     float del_eta=eta-mCtbHits[ih].eta;
00395     //printf("  match ih=%d del_eta=%f\n",ih, del_eta);
00396     if(fabs(del_eta) >mMatchCtbMax_eta) continue;
00397     
00398     // printf("  CTB match OK:  del_eta=%.2f, del_phi/deg=%.1f \n", del_eta,del_phi/C_PI*180);
00399     sigma=strag;
00400     n6++;    
00401     //  printf("add tr %d w/ sigma=%f p/GeV=%f spath/cm=%f nFitP=%d nPoss=%d nSVT=%d\n",mPrimCand.size(),sigma,pmom.mag(),spath,nFitP, nPoss,nSvtP);
00402     
00403     return true;
00404   }
00405   
00406   // match not found
00407   return false;
00408 }
00409 
00410 //==========================================================
00411 //==========================================================
00412 bool  
00413 StppLMVVertexFinder::ppLMV5() { 
00414   //  ----------  D O   F I N D    V E R T E X
00415   if(mMode == 0) LOG_WARN<<"ppLMV4 cuts have been activated"<<endm; 
00416 
00417   int totTr=mPrimCand.size();
00418   uint minTr=mMinMatchTr;
00419   if(mVertexConstrain) minTr++;
00420   //printf("passed %d tracks match to CTB,  BeamLine=%d\n",totTr,mVertexConstrain );
00421   LOG_DEBUG << "passed " << totTr << " tracks match to CTB,  BeamLine=" << mVertexConstrain << endm;
00422    
00423   double xo=0.0,yo=0.0;
00424   xo=mX0;
00425   yo=mY0; 
00426  
00427   //Do the actual vertex fitting, continue until good
00428   double A11=0.0,A12=0.0,A13=0.0,A21=0.0,A22=0.0,A23=0.0;
00429   double A31=0.0,A32=0.0,A33=0.0; // Matrix Elements
00430   double C11=0.0,C12=0.0,C13=0.0,C21=0.0,C22=0.0,C23=0.0;
00431   double C31=0.0,C32=0.0,C33=0.0; // C = A^-1
00432   int done = 0;
00433   int vertIter=0;
00434   double chi2=0;
00435   StThreeVectorD XVertex(999.,888.,777.);
00436   while( done != 1 ){
00437     vertIter++;
00438     // Check that there at least are 2 tracks
00439     if( mPrimCand.size() < minTr ){
00440       LOG_WARN  << "ppLMV5: below  "<<minTr<<" track remains. No vertex found." << endm;
00441       return false;
00442     }
00443   
00444     // Begin by doing a fit
00445     A11=0.0,A12=0.0,A13=0.0,A21=0.0,A22=0.0,A23=0.0;
00446     A31=0.0,A32=0.0,A33=0.0; // Matrix Elements
00447     C11=0.0,C12=0.0,C13=0.0,C21=0.0,C22=0.0,C23=0.0;
00448     C31=0.0,C32=0.0,C33=0.0; // C = A^-1
00449     double b1=0.0,b2=0.0,b3=0.0;
00450 
00451     // Compute matrix A and vector b
00452     for(unsigned int itr=0; itr < mPrimCand.size(); itr++){ 
00453 
00454       double spath = mPrimCand[itr].helix.pathLength(xo,yo);
00455       StThreeVectorD XClosest =  mPrimCand[itr].helix.at(spath);
00456       StThreeVectorD XMomAtClosest =  mPrimCand[itr].helix.momentumAt(spath,mBfield*tesla );
00457       double xp   = XClosest.x(); double yp= XClosest.y(); double zp= XClosest.z();  
00458       // printf("itr=%d  DCA x=%f y=%f z=%f  sig=%f\n",itr,xp,yp,zp,mPrimCand[itr].sigma);
00459       
00460       double xhat = XMomAtClosest.x()/XMomAtClosest.mag();
00461       double yhat = XMomAtClosest.y()/XMomAtClosest.mag();
00462       double zhat = XMomAtClosest.z()/XMomAtClosest.mag();
00463       float sig=mPrimCand[itr].sigma;
00464       A11=A11+(yhat*yhat+zhat*zhat)/sig;
00465       A12=A12-(xhat*yhat)/sig;
00466       A13=A13-(xhat*zhat)/sig;
00467       A22=A22+(xhat*xhat+zhat*zhat)/sig;
00468       A23=A23-(yhat*zhat)/sig;
00469       A33=A33+(xhat*xhat+yhat*yhat)/sig;
00470       b1=b1 + ( (yhat*yhat+zhat*zhat)*xp - xhat*yhat*yp - xhat*zhat*zp )/sig;
00471       b2=b2 + ( (xhat*xhat+zhat*zhat)*yp - xhat*yhat*xp - yhat*zhat*zp )/sig;
00472       b3=b3 + ( (xhat*xhat+yhat*yhat)*zp - xhat*zhat*xp - yhat*zhat*yp )/sig;
00473     }
00474     A21 = A12; A31=A13; A32=A23;
00475 
00476     // Invert A
00477     double detA =   A11*A22*A33 + A12*A23*A31 + A13*A21*A32;
00478     detA = detA   - A31*A22*A13 - A32*A23*A11 - A33*A21*A12;
00479     //    cout<<"Determinant= "<<detA<<endl;
00480     //    cout<<"A11,A12,A13: "<<A11<<" "<<A12<<" "<<A13<<endl;
00481     //    cout<<"A21,A22,A23: "<<A21<<" "<<A22<<" "<<A23<<endl;
00482     //    cout<<"A31,A32,A33: "<<A31<<" "<<A32<<" "<<A33<<endl;
00483     //    cout<<"b1,b2,b3 "<<b1<<" "<<b2<<" "<<b3<<endl;
00484     C11=(A22*A33-A23*A32)/detA; C12=(A13*A32-A12*A33)/detA; C13=(A12*A23-A13*A22)/detA;
00485     C21=C12;                    C22=(A11*A33-A13*A31)/detA; C23=(A13*A21-A11*A23)/detA;
00486     C31=C13;                    C32=C23;                    C33=(A11*A22-A12*A21)/detA;
00487 
00488     // Find Vertex Position
00489     double Xv = C11*b1 + C12*b2 + C13*b3;
00490     double Yv = C21*b1 + C22*b2 + C23*b3;
00491     double Zv = C31*b1 + C32*b2 + C33*b3;
00492     XVertex.setX(Xv); XVertex.setY(Yv); XVertex.setZ(Zv);
00493     LOG_DEBUG <<vertIter<<"  Vertex Position   : "<<XVertex.x()<<" "<<XVertex.y()<<" "<<XVertex.z()<<endm;
00494     LOG_DEBUG <<"Error in Position : "<<::sqrt(C11)<<" "<<::sqrt(C22)<<" "<<::sqrt(C33)<<endm;
00495     
00496 
00497     // Check if the fit is any good
00498     // Loop over tracks again to get Chi2 and check each track's deviation
00499     
00500     double dmax=0.0;
00501     vector<JHelix>::iterator itehlx, end, ikeep;
00502     end=mPrimCand.end();
00503     if(mVertexConstrain )end--; // to preserve beamL
00504 
00505     for( itehlx=mPrimCand.begin(); itehlx <end;  itehlx++) {
00506       StPhysicalHelixD hlx = (*itehlx).helix;
00507       double sig = (*itehlx).sigma;
00508       double spath = hlx.pathLength(XVertex.x(),XVertex.y()); 
00509       StThreeVectorD XHel = hlx.at(spath);
00510       double d=(XHel.x()-XVertex.x())*(XHel.x()-XVertex.x());
00511       d = d+(XHel.y()-XVertex.y())*(XHel.y()-XVertex.y());
00512       d = d+(XHel.z()-XVertex.z())*(XHel.z()-XVertex.z());
00513       d = ::sqrt(d);
00514       chi2 = chi2 + (d*d)/(sig*sig);
00515       double drel = d;  // do not use sig during track rejection ??? JB
00516       //printf(" DCA x=%f y=%f z=%f d=%f drel=%f dmax=%f sig=%f\n",XHel.x(),XHel.y(),XHel.z(),d,drel,dmax,sig);
00517       if( drel > dmax ){// save  track that deviates the most from vertex to be discarded later
00518         dmax = drel;
00519         ikeep = itehlx;
00520       }
00521     } 
00522 
00523     if( dmax > mDVtxMax ){
00524       LOG_INFO << "Removing a track! dmax=" << dmax << endm;
00525       mPrimCand.erase(ikeep);
00526       done=0;
00527     }
00528     else{
00529       done=1;
00530     }
00531   } // End While Loop
00532   
00533   //  double  chi2pdof = chi2/(mPrimCand.size()-1);
00534 
00535   StPrimaryVertex primV;
00536   Float_t cov[6] = {C11,0.0,C22,0.0,0.0,C33  };  //  cxx,?,cyy,?,?,czz
00537   
00538   primV.setPosition(XVertex);
00539   primV.setCovariantMatrix(cov); 
00540   primV.setVertexFinderId(pplmvVertexFinder);
00541   primV.setNumTracksUsedInFinder(mPrimCand.size());
00542   primV.setNumMatchesWithCTB(mPrimCand.size());// the same, only matched tracks are used
00543 
00544   /*  I do not understand if this elaborate scheme is applicable to multiple prim vertices for pp
00545       http://www.star.bnl.gov/STAR/html/all_l/html/dst_vertex_flags.html
00546       provide '1' as was
00547   */
00548   primV.setFlag(1); 
00549   primV.setRanking(555);
00550 
00551   //..... add vertex to the list
00552   addVertex(&primV);
00553 
00554 #if 0
00555   /*
00556     The fake vertex will never take precedence over the real one in 
00557     the muDst analysi since its rank is lower (if default ranking is used).
00558     But may pick a random match if luminosity (i.e. pileup is high).
00559     Unless the true vertex is not found. Then we would have a fake 
00560     vertex w/o prim tracks.
00561    Jan B.
00562   */
00563   if(eveID%2) addFakeVerex(XVertex.z()+20);
00564 #endif
00565 
00566 #if 0  
00567   mFitResult=XVertex;
00568   mFitError.setX(sqrt(C11));   mFitError.setY(sqrt(C22));   mFitError.setZ(sqrt(C33));
00569 #endif
00570 
00571 
00572   LOG_DEBUG <<"ppLMV5: Primary Vertex found!\nVert position: "<<XVertex<<", accepted tracks "<<mPrimCand.size()<<" of "<<totTr<<" eveID="<<eveID<<endm;
00573   printf("##V %6d %d %f %f %f    %d %d %d\n",eveID,mTotEve,XVertex.x(),XVertex.y(),XVertex.z(),mCtbHits.size(),n1,NCtbMatches());
00574 
00575 
00576   // get geant vertex
00577   St_DataSet *gds=StMaker::GetChain()->GetDataSet("geant");
00578   if(gds) {
00579     St_g2t_vertex  *g2t_ver=( St_g2t_vertex *)gds->Find("g2t_vertex");
00580     if(g2t_ver) {
00581       // --------------  A C C E S S    G E A N T   V E R T E X  (for histo)
00582       g2t_vertex_st *GVER= g2t_ver->GetTable();
00583       //  printf("GVER add=%d\n",GVER);
00584       // hPiFi[10]->Fill(GVER->ge_x[2]-rZver);
00585       // hPiFi[11]->Fill(GVER->ge_x[0]-rXver);
00586       // hPiFi[12]->Fill(GVER->ge_x[1]-rYver);
00587       printf("Z Geant-found=%.2f, dx=%.2f, dy=%.2f nCtbSl=%d n1=%d eveID=%d\n",GVER->ge_x[2]-XVertex.z(),GVER->ge_x[0]-XVertex.x(),GVER->ge_x[1]-XVertex.y(),mCtbHits.size(),n1,eveID);
00588       printf("##G %6d %d %f %f    %d %d %d\n",eveID,mTotEve,GVER->ge_x[2],XVertex.z(),mCtbHits.size(),n1,NCtbMatches());
00589       
00590     }
00591   }
00592   
00593   LOG_DEBUG << "StppLMVVertexFinder::ppLMV5() ends" << endm;
00594   return true;
00595 }
00596 
00597 
00598 //==========================================================
00599 //==========================================================
00600 int  StppLMVVertexFinder::NCtbMatches() { 
00601   int nTr=mPrimCand.size();
00602   if(mVertexConstrain ) nTr--; // drop beam line
00603   return nTr;
00604 }
00605 
00606 
00607 /*
00608  * $Log: StppLMVVertexFinder.cxx,v $
00609  * Revision 1.24  2010/01/26 21:01:49  fisyak
00610  * Clean up, switch from bit mask to attributes
00611  *
00612  * Revision 1.23  2006/05/05 18:35:39  balewski
00613  * block the fake second prim vertex
00614  *
00615  * Revision 1.22  2006/05/04 20:01:30  jeromel
00616  * Switched to logger
00617  *
00618  * Revision 1.21  2006/05/02 13:49:21  balewski
00619  * replace gufld() with  mBfield = event->runInfo()->magneticField();
00620  *
00621  * Revision 1.20  2005/07/19 21:56:58  perev
00622  * MultiVertex
00623  *
00624  * Revision 1.19  2005/06/21 02:16:36  balewski
00625  * multiple prim vertices are stored in StEvent
00626  *
00627  * Revision 1.18  2005/03/11 22:23:53  balewski
00628  * towards PPV
00629  *
00630  * Revision 1.17  2005/03/09 19:24:18  balewski
00631  * preparation for PPV vertex finder
00632  *
00633  * Revision 1.16  2004/12/16 17:01:36  balewski
00634  * 2 cuts in ppLMV4 were slightly off what was in the TPT version
00635  *
00636  * Revision 1.15  2004/12/13 20:39:58  fisyak
00637  * Add initaition of StGenericVertexFinder variables, replace mDumMaker by StMaker::GetChain() method
00638  *
00639  * Revision 1.14  2004/09/13 15:41:31  balewski
00640  * fix bug in ppLMV4/5 switch
00641  *
00642  * Revision 1.13  2004/09/03 14:24:15  jeromel
00643  * Fixed inverted comment
00644  *
00645  * Revision 1.12  2004/09/03 00:09:08  jeromel
00646  * Modified code to Implement Init() and SetMode() and allow passing a switch
00647  * to chose the vertex finder from within the same code implementation. Was
00648  * needed for ppLMV (one implementation, two algorithm)
00649  *
00650  * Revision 1.11  2004/09/01 18:45:01  balewski
00651  * ppLMV5/4 switch added
00652  *
00653  * Revision 1.10  2004/08/25 15:20:30  balewski
00654  * Z-vertex range increased to +/- 150 cm
00655  *
00656  * Revision 1.9  2004/08/18 20:08:04  balewski
00657  * outerGeometry()->helix() used by ppLMV for extrapolation to CTB
00658  *
00659  * Revision 1.8  2004/08/06 21:00:01  balewski
00660  * now should be fine for 2004 pp200 data,
00661  * high pT & small path cat-off included
00662  *  still no nFit/nPoss cut
00663  *
00664  * Revision 1.7  2004/08/06 04:49:14  balewski
00665  * that would be it
00666  *
00667  * Revision 1.6  2004/08/05 22:08:04  balewski
00668  * toward working point
00669  *
00670  * Revision 1.5  2004/08/04 21:57:56  balewski
00671  * toward smarter ppLMV5
00672  *
00673  * Revision 1.4  2004/07/29 00:44:30  balewski
00674  * ppLMV() returns false on error
00675  *
00676  * Revision 1.3  2004/07/24 02:57:40  balewski
00677  * clean up of ppLMV, CTB-util separated
00678  *
00679  * Revision 1.2  2004/07/23 01:01:33  jeromel
00680  * See .h + changed printf() / cout to gMessMgr + removed a ew asserts
00681  *
00682  * Revision 1.1  2004/07/21 01:53:18  balewski
00683  * first
00684  *
00685  **************************************************************************/
00686