Project:/VertexAlg/src/VertexFit.cc

CDF KITS
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
Architecture:	[ i386 ]
Version:	[ 4.10.4 ] [ 4.10.5 ] [ 4.8.4 ] [ 4.8.4l3s ] [ 4.8.5 ] [ 4.9.0 ] [ 4.9.1 ] [ 4.9.1.hpt3 ] [ 4.9.1hpt3 ] [ 4.9.1top1 ] [ 5.0.0 ] [ 5.1.0 ] [ 5.1.0beamonly ] [ 5.1.1 ] [ 5.2.0 ] [ 5.3.0 ] [ 5.3.1 ] [ 5.3.1dsp ] [ 5.3.3 ] [ 5.3.3_nt ] [ 5.3.4 ] [ 6.1.1 ] [ 6.1.1b ] [ 6.1.2 ] [ 6.1.3 ] [ 6.1.4 ] [ 6.1.4int3 ] [ 6.1.4mc ] [ 6.1.4mc_a ] [ 6.1.6 ] [ development ]
001  /************************************************************************

002  * VertexFit class implementation file                                   *

003  *                                                                       *

004  *                                                                       *

005  * Author: Craig Blocker, CDF Group                                      *

006  *         Phone 781-736-2920                                            *

007  *                                                                       *

008  * Description: Routines to implement VertexFit class                    *

009  *                                                                       *

010  * See VertexFit.hh in the VertexAlg/ area and                           *

011  * the file VertexFit.txt in the /doc area of VertexAlg for              *

012  * documentation.                                                        *

013  *                                                                       *

014  * Revision History:                                                     *

015  *    May 27, 1999   Craig Blocker   Initial creation                    * 

016  *    June 22, 2001  Craig Blocker   Fix a few bugs in checking on       *

017  *                                     vertices (such as pointing to     *

018  *                                     primary vertex wasn't allowed)    *

019  *                                                                       *

020  *    Sept 07, 2002 Joao Guimaraes   Added accessors for ijk errors and  *

021  *                                   track-id of tracks that produce     *

022  *                                   some of those errors                *

023  *                                   These methods are:                  *

024  *                                      getIJKErr(int&,int&,int&) const  *

025  *                                      getIJKErr0() const               *

026  *                                      getIJKErr1() const               *

027  *                                      getIJKErr2() const               *

028  *                                      getErrTrackId() const            *

029  *                                                                       *

030  *    Feb 10, 2002 Joao Guimaraes    Added two methods to allow for a    *

031  *                                   beam constrained fit with CTVMFT.   *

032  *                                   This only makes sense when fitting  *

033  *                                   a primary vertex and should not be  *

034  *                                   used for fitting secondary vertices *

035  *    Oct 15 2003 Sal Rappoccio (rappocc@fnal.gov)                       *

036  *                  Added accessor to add tracks to ctvmft via a track   *

037  *                  id and the track parameters and covariance           *

038  *                                                                       *

039  *  December 2 2003 Sal Rappoccio (rappocc@fnal.gov)                     *

040  *                  Added an interface to use extrapolated track errors  *

041  *                  in a "transparent" way. The new interface is used as:*

042  *                        VertexFit fit;                                 *

043  *                        fit.init();                                    *

044  *                        fit.setPrimaryVertex ( p1 ); // Note the order *

045  *                                                     // matters        *

046  *                        fit.setReferencePoint( p2 );                   *

047  *                        fit.addTrack( trk1, M_PION, 1 );               *

048  *                        fit.addTrack( trk2, M_PION, 1 );               *

049  *                        fit.fit();                                     *

050  *                        Hep3Vector & pos = fit.getVertex( 1 );         *

051  *                  "pos" will now hold the new vertex in CDF coordinates*

052  *                  constructed from tracks with errors referenced at the*

053  *                  new reference point "p2".                            *

054  *                  This track extrapolation code is OFF by default. It  *

055  *                  is turned ON by:   fit.setReferencePoint( p2 );      * 

056  *                  AddTrack was modified to use this feature.           *

057  *  Dec 2 2003 Joao Guimaraes                                            *

058  *                  Updated VertexFit::beamlineConstraint to use new     *

059  *                  reference point (see entry above).                   *

060  *                                                                       *

061  *************************************************************************/
062 
063 
064 #ifndef VERTEXFIT_CC
065 #define VERTEXFIT_CC
066 
067 //-----------------------

068 // This Class's Header --

069 //-----------------------

070 #include "VertexAlg/VertexFit.hh"
071 
072 extern "C" void ctvmft_(int&,int&,int&);
073 extern "C" float prob_(float&,int&);
074 extern "C" bool mcalc_(int&,int*,float&,float&,double*);
075 extern "C" void dcalc_(int&,int&,float*,float&,float&,float*);
076 
077 // C++ Headers --

078 //---------------

079 #if defined(DEFECT_OLD_IOSTREAM_HEADERS)
080 #include <iostream>
081 #else // Standard C++

082 #include <iostream>
083 #endif
084 
085 #if defined(DEFECT_NO_STDLIB_NAMESPACES)
086 // Do nothing

087 #else // Standard C++

088 using std::cerr ;
089 using std::cout ;
090 using std::endl ;
091 #endif
092 
093 #include <algorithm>
094 #include <math.h>
095 
096 //-------------------------------

097 // Collaborating Class Headers --

098 //-------------------------------

099 #include "TrackingObjects/Tracks/CdfTrack.hh"
100 #include "CLHEP/Matrix/Matrix.h"
101 #include "TrackingHL/Utility/TrackExtrapolator.hh"
102 
103 //for floating point exception handling Joe and Farrukh

104 #include <csignal>
105 #include <csetjmp>
106 
107 /******************************************************************************

108  * Default constructor                                                        *

109  *                                                                            *

110  *                                                                            *

111  *****************************************************************************/
112 
113 //the following will have scope throughout.

114 
115  jmp_buf env;
116  extern "C" void VertexFitSetStatus(int i) { 
117    std::cout << "Warning, you are handling a severe error in VertexFit" << std::endl;
118    longjmp(env,-66);
119  }
120 
121 VertexFit::VertexFit() :
122   _currentAllocatedVertexNumber(0),   // facilitates CandNode recursion

123   _referencePoint( 0,0,0 ), // Set Reference Point to (0,0,0) initially

124   _primaryVertex( 0,0,0 ), // Set primary vertex to (0,0,0) initially

125   _cdfPrimaryVertex( 0,0,0 ) // Set p.v. in CDF coords to (0,0,0) initially

126 {
127 // Set name and email of VertexFit expert.

128    _expert="Craig Blocker (blocker@fnal.gov)";
129 // First get pointers to various FORTAN common blocks.

130    _ctvmq_com = (CTVMQ*) ctvmq_address_();
131    _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
132    _fiddle_com = (FIDDLE*) fiddle_address_();
133    _trkprm_com = (TRKPRM*) trkprm_address_();
134 // Initialize various arrays

135    init();
136 // Don't bomb program on error.

137    _fiddle.excuse=1;
138 // Don't extrapolate track errors by default

139    _extrapolateTrackErrors = false;
140 }
141 
142 /******************************************************************************

143  * Destructor                                                                 *

144  *        Use default, i.e. do nothing                                        *

145  *****************************************************************************/
146 
147 VertexFit::~VertexFit(){ }
148 
149 
150 
151 // Particle masses

152 const float VertexFit::EL_MASS = 0.000510998902;
153 const float VertexFit::MU_MASS = 0.105658357;
154 const float VertexFit::PI_MASS = 0.13957018;
155 const float VertexFit::K_MASS = 0.493677;
156 const float VertexFit::PROTON_MASS = 0.938272;
157 
158 /*=======================================================================

159                          VertexFit::init

160   =======================================================================*/
161 
162 void VertexFit::init()
163 {
164   double bmag = 1.4116; // Magnetic field in Tesla

165    init(bmag);
166 }
167 
168 void VertexFit::init(double bmag)
169 {
170 // Now initialize CTVMQ common.  Run and trigger numbers are

171 // dummies - they are not used.

172    _ctvmq.runnum=1;
173    _ctvmq.trgnum=100;
174 // Eventually, we have to get the magnetic field from the right place.

175 // Origignal with      bmag = 14.116:

176 // _ctvmq.pscale=0.000149896*bmag;

177 // New bfield in Tesla bmag = 1.4116 [T]:

178    _ctvmq.pscale=0.00149896*bmag;
179 // Set the default maximum chi-square per dof.

180    _fiddle.chisqmax = 225.;
181 // We also need to get the primary vertex from the right place,

182 // but for now we put in (0,0,0).

183    setPrimaryVertex(0.,0.,0.);
184    float xverr[3][3];
185    for(int j = 0; j<3; j++) {
186       for(int k = 0; k<3; k++) {
187          xverr[j][k] = 0.;
188       }
189    }
190    xverr[0][0]=.005;
191    xverr[1][1]=.005;
192    xverr[2][2]=1.;
193    setPrimaryVertexError(xverr);
194 // Zero number of tracks, vertices and mass constraints.

195    _ctvmq.ntrack=0;
196    _ctvmq.nvertx=0;
197    _ctvmq.nmassc=0;
198 // Zero track list and arrays containing the vertex and

199 // mass constraint configuration.

200    for(int j=0; j<_maxtrk; ++j) {
201       _ctvmq.list[j]=0;
202       for(int jv=0; jv<_maxvtx; ++jv) {
203          _ctvmq.trkvtx[jv][j]=false;
204       }
205       for(int jmc=0; jmc<_maxmcn; ++jmc) {
206          _ctvmq.trkmcn[jmc][j]=false;
207       }
208    }
209 // Initialize the conversion and vertex pointing arrays.

210    for(int jv=0; jv<_maxvtx; ++jv) {
211       _ctvmq.cvtx[jv]=0;
212       _ctvmq.vtxpnt[0][jv]=-1;
213       _ctvmq.vtxpnt[1][jv]=0;
214    }
215    _ctvmq.drmax=2.;
216    _ctvmq.dzmax=20.;
217    _ctvmq.rvmax=70.;
218    _ctvmq.trnmax=0.5;
219    _ctvmq.dsmin=-2.;   
220 // Set stat to -999  and chisqq to -1 to symbolize that no fit has yet been done.

221    stat=-999;
222    _ctvmq.chisqr[0]=-1.;
223 
224    _primaryVertex = Hep3Vector(0,0,0);
225    _cdfPrimaryVertex = Hep3Vector(0,0,0);
226    _referencePoint = Hep3Vector(0,0,0);
227 }
228 
229 /*=======================================================================

230                          VertexFit::addTrack

231   =======================================================================*/
232 
233 bool VertexFit::addTrack(const CdfTrack* trk, float mass,
234                          vertexNumber jv)
235 {
236 // Check that this vertex number is within the allowed range.

237    if(jv<VERTEX_1 || jv>_maxvtx) return false;
238    _ctvmq.nvertx = jv>_ctvmq.nvertx ? jv : _ctvmq.nvertx;
239 
240 // Check that we have not exceeded the maximum number of tracks.

241    if(_ctvmq.ntrack>=_maxtrk) return false;
242 
243 
244 // Sal Rappoccio: Add track extrapolation, if the user set it

245    HepVector w = trk->par();
246    HepSymMatrix m = trk->getHelixFit().getErrorMatrix();
247    if ( _extrapolateTrackErrors ) {
248      // This will move the reference point

249      moveReferencePoint( w, m );
250    }
251 
252 // Add this track.

253    _ctvmq.list[_ctvmq.ntrack]=trk->id().value();
254    _ctvmq.tkbank[_ctvmq.ntrack][0]='Q';
255    _ctvmq.tkbank[_ctvmq.ntrack][1]='T';
256    _ctvmq.tkbank[_ctvmq.ntrack][2]='R';
257    _ctvmq.tkbank[_ctvmq.ntrack][3]='K';
258    _ctvmq.tmass[_ctvmq.ntrack]=mass;
259    _ctvmq.trkvtx[jv-1][_ctvmq.ntrack]=true;
260 
261 // Put this track's helix parameters and error matrix into

262 // a fortran common block so that they can  be accessed

263 // by gettrk.  This is a dummy for now.

264    _trkprm.trhelix[_ctvmq.ntrack][0]=w[0];
265    _trkprm.trhelix[_ctvmq.ntrack][1]=w[1];
266    _trkprm.trhelix[_ctvmq.ntrack][2]=w[2];
267    _trkprm.trhelix[_ctvmq.ntrack][3]=w[3];
268    _trkprm.trhelix[_ctvmq.ntrack][4]=w[4];
269 
270    for(int j=0; j<5; ++j) {
271       for(int k=0; k<5; ++k) {
272          _trkprm.trem[_ctvmq.ntrack][j][k]=m[j][k];
273       }
274    }
275    _ctvmq.ntrack++;
276 
277    return true;
278 }
279 
280 // Sal Rappoccio: 11/04/03: Added this method to add track using

281 // their track id and the explicit helix parameters.

282 // This is useful when for adding tracks with helical parameters relative 

283 // to a point other than (0,0)

284 
285 /*=======================================================================

286                          VertexFit::addTrack

287   =======================================================================*/
288 
289 bool VertexFit::addTrack( const HepVector & v,
290                           const HepSymMatrix & cov,
291                           int trackId,
292                           float mass,
293                           vertexNumber jv)
294 {
295 // Check that this vertex number is within the allowed range.

296    if(jv<VERTEX_1 || jv>_maxvtx) return false;
297    _ctvmq.nvertx = jv>_ctvmq.nvertx ? jv : _ctvmq.nvertx;
298 
299 // Sal Rappoccio: Add track extrapolation, if the user set it

300 
301    HepVector w = v;
302    HepSymMatrix m = cov;
303    if ( _extrapolateTrackErrors ) {
304      // This will move the reference point

305      moveReferencePoint( w, m );
306    }
307 
308 // Check that we have not exceeded the maximum number of tracks.

309    if(_ctvmq.ntrack>=_maxtrk) return false;
310 
311 // Add this track.

312    _ctvmq.list[_ctvmq.ntrack]=trackId;
313    _ctvmq.tkbank[_ctvmq.ntrack][0]='Q';
314    _ctvmq.tkbank[_ctvmq.ntrack][1]='T';
315    _ctvmq.tkbank[_ctvmq.ntrack][2]='R';
316    _ctvmq.tkbank[_ctvmq.ntrack][3]='K';
317    _ctvmq.tmass[_ctvmq.ntrack]=mass;
318    _ctvmq.trkvtx[jv-1][_ctvmq.ntrack]=true;
319 
320 // Put this track's helix parameters and error matrix into

321 // a fortran common block so that they can  be accessed

322 // by gettrk.  This is a dummy for now.

323    _trkprm.trhelix[_ctvmq.ntrack][0]=w[0];
324    _trkprm.trhelix[_ctvmq.ntrack][1]=w[1];
325    _trkprm.trhelix[_ctvmq.ntrack][2]=w[2];
326    _trkprm.trhelix[_ctvmq.ntrack][3]=w[3];
327    _trkprm.trhelix[_ctvmq.ntrack][4]=w[4];
328 
329    for(int j=0; j<5; ++j) {
330       for(int k=0; k<5; ++k) {
331          _trkprm.trem[_ctvmq.ntrack][j][k]=m[j][k];
332 
333       }
334    }
335    _ctvmq.ntrack++;
336 
337    return true;
338 }
339 
340 /*=======================================================================

341                          VertexFit::vertexPoint_2d

342   =======================================================================*/
343   
344 bool VertexFit::vertexPoint_2d(vertexNumber jv1, vertexNumber jv2)
345 {
346 // Check that these vertex numbers are within allowed range and

347 // that the vertices are unique.

348    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
349    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
350    if(jv1 <= jv2) return false;
351 
352 // Setup vertex pointing.

353    _ctvmq.vtxpnt[0][jv1-1]=jv2;
354    _ctvmq.vtxpnt[1][jv1-1]=1;           // 2d pointing.

355 
356    return true;
357 }
358 
359 /*=======================================================================

360                          VertexFit::vertexPoint_3d

361   =======================================================================*/
362 
363 bool VertexFit::vertexPoint_3d(vertexNumber jv1, vertexNumber jv2)
364 {
365 // Check that these vertex numbers are within allowed range and

366 // that the vertices are distinct.

367    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
368    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
369    if(jv1 <= jv2) return false;
370 // Setup vertex pointing.

371    _ctvmq.vtxpnt[0][jv1-1]=jv2;
372    _ctvmq.vtxpnt[1][jv1-1]=2;           // 3d pointing.

373 
374    return true;
375 }
376 
377 /*=======================================================================

378                          VertexFit::vertexPoint_1track

379   =======================================================================*/
380 
381 bool VertexFit::vertexPoint_1track(vertexNumber jv1, vertexNumber jv2)
382 {
383 // Check that these vertex numbers are within allowed range and are distinct.

384    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
385    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
386    if(jv1 <= jv2) return false;
387 
388 // Setup vertex pointing.

389    _ctvmq.vtxpnt[0][jv1-1]=jv2;
390    _ctvmq.vtxpnt[1][jv1-1]=3;           // Point to 1 track vertex.

391 
392    return true;
393 }
394 /*=======================================================================

395                          VertexFit::vertexPoint_0track

396   =======================================================================*/
397 
398 bool VertexFit::vertexPoint_0track(vertexNumber jv1, vertexNumber jv2)
399 {
400 
401   // jv2 is the zero track vertex. jv1 is the multi track vertex which points to jv2

402 
403   // Note: You must call this routine at least twice in order for ctvmft_zerotrackvtx to work

404   // ie there must be at least 2 vertices pointed at a zero track vertex. ctvmft

405   // for this, but the error message may not make it to your log file (look at the local

406   // variables in the stack frame especially IJKERR(2). The significance of this

407   // error code is documented at the top of whatever routine chucked you out (ctvmf00 in

408   // this case)

409 
410   // see ctvmft.f source file for discussion. See especially comments at the top

411   // of subroutines: ctvmft and ctvmfa

412 
413 // Check that these vertex numbers are within allowed range and are distinct.

414    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
415    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
416    if(jv1 <= jv2) return false;
417 
418 // Setup vertex pointing.

419    _ctvmq.vtxpnt[0][jv1-1]=jv2;
420    _ctvmq.vtxpnt[1][jv1-1]=4;           // Point to 0 track vertex.

421 
422    return true;
423 }//vertexPoint_0track

424 
425 /*=======================================================================

426                          VertexFit::conversion_2d

427   =======================================================================*/
428 
429 bool VertexFit::conversion_2d(vertexNumber jv)
430 {
431    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
432       return false;
433    }
434    _ctvmq.cvtx[jv-1]=1;
435    return true;
436 }
437 
438 /*=======================================================================

439                          VertexFit::conversion_3d

440   =======================================================================*/
441 
442 bool VertexFit::conversion_3d(vertexNumber jv)
443 {
444    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
445       return false;
446    }
447    _ctvmq.cvtx[jv-1]=2;
448    return true;
449 }
450 
451 /*=======================================================================

452                          VertexFit::massConstrain

453   =======================================================================*/
454 
455 bool VertexFit::massConstrain(const CdfTrack* trk1,
456                               const CdfTrack* trk2, float mass)
457 {
458    int ntrk=2;
459    const CdfTrack* tracks[2];
460    tracks[0]=trk1;
461    tracks[1]=trk2;
462    return massConstrain(ntrk,tracks,mass);
463 }
464 
465 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
466                                 const CdfTrack* trk3, float mass)
467 {
468    int ntrk=3;
469    const CdfTrack* tracks[3];
470    tracks[0]=trk1;
471    tracks[1]=trk2;
472    tracks[2]=trk3;
473    return massConstrain(ntrk,tracks,mass);
474 }
475 
476 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
477                     const CdfTrack* trk3, const CdfTrack* trk4,float mass)
478 {
479    int ntrk=4;
480    const CdfTrack* tracks[4];
481    tracks[0]=trk1;
482    tracks[1]=trk2;
483    tracks[2]=trk3;
484    tracks[3]=trk4;
485    return massConstrain(ntrk,tracks,mass);
486 }
487 
488 bool VertexFit::massConstrain(int ntrk, const CdfTrack* tracks[], float mass)
489 {
490 // Check that we have not exceeded the allowed number of mass constraints.

491    if(_ctvmq.nmassc>=_maxmcn) return false;
492 
493 // Set constraint mass

494    _ctvmq.cmass[_ctvmq.nmassc]=mass;
495      
496 // For each track in contraint, set trkmcn true.  Since the number in

497 // tracks[] is the track number, we have to find each track in the

498 // list of tracks.

499    for(int jt=0; jt<ntrk; ++jt) {
500       bool found=false;
501       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
502          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
503             _ctvmq.trkmcn[_ctvmq.nmassc][kt]=true;
504             found=true;
505          }
506       }
507       if(!found) return false;
508    }
509 // Increment number of mass constraints.   

510    _ctvmq.nmassc++;
511 
512    return true;
513 }
514 
515 
516 /*=======================================================================

517                      VertexFit::beamlineConstraint

518   =======================================================================*/
519 //  Turn ON this option only when fitting the primary with no additional 

520 // constraints. The routine is protected and will not function

521 // if those options are selected as well as the beamline constraint.

522 // In case this option is chosen the fit will calculate a result also with

523 // just one track in input.

524 // In order to enable this option the user must do the following:

525 // Fit the primary vertex as VERTEX_1 (without any other vertices).

526 // Enable beamlineConstraint by setting the pointing constraint variable

527 // vtxpnt[0][0] to -100 (no pointing constraints when <0)

528 // Note that index 0,0 corresponds to index 1,1 in fortran

529 // Provide the beamline parameters:

530 //  --> assume xv = xb + xzbslope*z, 

531 //             yv = yb + yzbslope*z, where (xv, yv) is the transverse

532 //      beam position at z, (xb, yb) is the transverse beam position at z = 0 

533 //      and (xzbslope, yzbslope) are the slopes of the beamline

534 //  --> set primary vertex at z=0 (xb, yb, 0)

535 //  --> set the diagonal elements of the primary vertex error matrix to the

536 //      sigma**2 of beam spread in x, y and z:

537 //      xverr[1][1] =  (~30 - 50 um)**2

538 //      xverr[2][2] =  (~30 - 50 um)**2

539 //      xverr[3][3] =  (~30 cm) **2

540 //      All other elements should be 0

541 // For help email: Joao Guimaraes  guima@huhepl.harvard.edu

542 //                 Franco Bedeschi bed@fnal.gov 

543 // See more information in the ctvmft.f

544 bool VertexFit::beamlineConstraint(float xb, float yb, HepSymMatrix berr,
545                                    float xzbslope, float yzbslope)
546 {
547   // Set beam position at z=0

548   setPrimaryVertex(xb,yb,0);
549   if (_extrapolateTrackErrors) {
550     float newXb = xb - _referencePoint.x() + 
551       _referencePoint.z() * xzbslope;
552     float newYb = yb - _referencePoint.y() + 
553       _referencePoint.z() * yzbslope;
554     setPrimaryVertex( newXb, newYb, 0);
555   }
556 
557   bool success = setPrimaryVertexError(berr);
558   //Set the beamline slope values

559   _ctvmq.xzslope= xzbslope;
560   _ctvmq.yzslope= yzbslope;
561   // Turn ON beamline constraint

562   _ctvmq.vtxpnt[0][0]=-100; 
563 
564   return success;
565 }
566 
567 bool VertexFit::beamlineConstraint(Hep3Vector pv, HepSymMatrix berr,
568                                    float xzbslope, float yzbslope)
569 {
570   // Check if input beam position coordinates are at z=0

571   if (pv.z() != 0) return false;
572 
573   return beamlineConstraint(pv.x(),pv.y(),berr,xzbslope,yzbslope);
574 }
575 
576 /*=======================================================================

577                          VertexFit::setPrimaryVertex

578   =======================================================================*/
579 
580 void VertexFit::setPrimaryVertex(float xv, float yv, float zv)
581 {
582 // Set x,y,z position of the primary vertex.

583    _ctvmq.xyzpv0[0]=xv;
584    _ctvmq.xyzpv0[1]=yv;
585    _ctvmq.xyzpv0[2]=zv;
586 
587    _primaryVertex = Hep3Vector( xv, yv, zv );
588    return;
589 }
590 
591 void VertexFit::setPrimaryVertex(Hep3Vector pv)
592 {
593 // Set x,y,z position of the primary vertex.

594    _ctvmq.xyzpv0[0]=pv.x();
595    _ctvmq.xyzpv0[1]=pv.y();
596    _ctvmq.xyzpv0[2]=pv.z();
597 
598    _primaryVertex = pv;
599    return;
600 }
601 
602 /*=======================================================================

603                          VertexFit::setPrimaryVertexError

604   =======================================================================*/
605 
606 void VertexFit::setPrimaryVertexError(const float xverr[3][3])
607 {
608 // Set the error matrix for the primary vertex.

609    for(int j=0; j<3; ++j) {
610       for(int k=0; k<3; ++k) {
611          _ctvmq.exyzpv[j][k]=xverr[j][k];
612       }
613    }
614    return;
615 }
616 
617 bool VertexFit::setPrimaryVertexError(const HepSymMatrix &xverr)
618 {
619 // Set the error matrix for the primary vertex using a HepSymMatrix.

620 // First check that the matrix is the correct size.

621    if(xverr.num_row() != 3) return false;
622    for(int j=0; j<3; j++) {
623      for(int k=0; k<3; k++) {
624        _ctvmq.exyzpv[j][k]=xverr[j][k];
625      }
626    }
627    return true;
628 }
629 
630 /*=======================================================================

631                          VertexFit::fit

632   =======================================================================*/
633 
634 bool VertexFit::fit()
635 {
636 // Check that the diagonal elements of all the track error matrices

637 //are positive.

638   bool mstat = true;
639   for(int trk=0; trk<_ctvmq.ntrack; ++trk) {
640     for(int j=0; j<5; ++j) {
641 // Check diagonal elements of error matrix.

642       if(_trkprm.trem[trk][j][j] < 0.) {
643 // Tell user that the covariance matrix could not be inverted.

644 // Set the error codes and fail this fit.

645         mstat = false;
646         _ctvmq.ijkerr[0] = 3;
647         _ctvmq.ijkerr[1] = 2;
648         _ctvmq.ijkerr[2] = trk + 1;
649       }
650     }
651 // Check that the curvature of the track is reasonable, that is,

652 // that the Pt is above ~10Mev/c.  If not, set the error codes

653 // and fail this fit.

654     if(fabs(_trkprm.trhelix[trk][1]) > 0.01) {
655       mstat = false;
656       _ctvmq.ijkerr[0] = 3;
657       _ctvmq.ijkerr[1] = 5;
658       _ctvmq.ijkerr[2] = trk + 1;
659     }