CDF II
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  *************************************************************************/
036 
037 #ifndef VERTEXFIT_CC
038 #define VERTEXFIT_CC
039 
040 //-----------------------
041 // This Class's Header --
042 //-----------------------
043 #include "VertexAlg/VertexFit.hh"
044 
045 extern "C" void ctvmft_(int&,int&,int&);
046 extern "C" float prob_(float&,int&);
047 extern "C" bool mcalc_(int&,int*,float&,float&,double*);
048 extern "C" void dcalc_(int&,int&,float*,float&,float&,float*);
049 
050 // C++ Headers --
051 //---------------
052 #if defined(DEFECT_OLD_IOSTREAM_HEADERS)
053 #include <iostream>
054 #else // Standard C++
055 #include <iostream>
056 #endif
057 
058 #if defined(DEFECT_NO_STDLIB_NAMESPACES)
059 // Do nothing
060 #else // Standard C++
061 using std::cerr ;
062 using std::cout ;
063 using std::endl ;
064 #endif
065 
066 #include <algorithm>
067 #include <math.h>
068 
069 //-------------------------------
070 // Collaborating Class Headers --
071 //-------------------------------
072 #include "TrackingObjects/Tracks/CdfTrack.hh"
073 #include "CLHEP/Matrix/Matrix.h"
074 
075 //for floating point exception handling Joe and Farrukh
076 #include <csignal>
077 #include <csetjmp>
078 
079 /******************************************************************************
080  * Default constructor                                                        *
081  *                                                                            *
082  *                                                                            *
083  *****************************************************************************/
084 
085 //the following will have scope throughout.
086 
087  jmp_buf env;
088  extern "C" void VertexFitSetStatus(int i) { 
089    std::cout << "Warning, you are handling a severe error in VertexFit" << std::endl;
090    longjmp(env,-66);
091  }
092 
093 VertexFit::VertexFit() :
094   _currentAllocatedVertexNumber(0)   // facilitates CandNode recursion
095 {
096 // Set name and email of VertexFit expert.
097    _expert="Craig Blocker (blocker@fnal.gov)";
098 // First get pointers to various FORTAN common blocks.
099    _ctvmq_com = (CTVMQ*) ctvmq_address_();
100    _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
101    _fiddle_com = (FIDDLE*) fiddle_address_();
102    _trkprm_com = (TRKPRM*) trkprm_address_();
103 // Initialize various arrays
104    init();
105 // Don't bomb program on error.
106    _fiddle.excuse=1;
107 }
108 
109 /******************************************************************************
110  * Destructor                                                                 *
111  *        Use default, i.e. do nothing                                        *
112  *****************************************************************************/
113 
114 VertexFit::~VertexFit(){ }
115 
116 
117 
118 // Particle masses
119 const float VertexFit::EL_MASS = 0.000510998902;
120 const float VertexFit::MU_MASS = 0.105658357;
121 const float VertexFit::PI_MASS = 0.13957018;
122 const float VertexFit::K_MASS = 0.493677;
123 const float VertexFit::PROTON_MASS = 0.938272;
124 
125 /*=======================================================================
126                          VertexFit::init
127   =======================================================================*/
128 
129 void VertexFit::init()
130 {
131   double bmag = 1.4116; // Magnetic field in Tesla
132    init(bmag);
133 }
134 
135 void VertexFit::init(double bmag)
136 {
137 // Now initialize CTVMQ common.  Run and trigger numbers are
138 // dummies - they are not used.
139    _ctvmq.runnum=1;
140    _ctvmq.trgnum=100;
141 // Eventually, we have to get the magnetic field from the right place.
142 // Origignal with      bmag = 14.116:
143 // _ctvmq.pscale=0.000149896*bmag;
144 // New bfield in Tesla bmag = 1.4116 [T]:
145    _ctvmq.pscale=0.00149896*bmag;
146 // Set the default maximum chi-square per dof.
147    _fiddle.chisqmax = 225.;
148 // We also need to get the primary vertex from the right place,
149 // but for now we put in (0,0,0).
150    setPrimaryVertex(0.,0.,0.);
151    float xverr[3][3];
152    for(int j = 0; j<3; j++) {
153       for(int k = 0; k<3; k++) {
154          xverr[j][k] = 0.;
155       }
156    }
157    xverr[0][0]=.005;
158    xverr[1][1]=.005;
159    xverr[2][2]=1.;
160    setPrimaryVertexError(xverr);
161 // Zero number of tracks, vertices and mass constraints.
162    _ctvmq.ntrack=0;
163    _ctvmq.nvertx=0;
164    _ctvmq.nmassc=0;
165 // Zero track list and arrays containing the vertex and
166 // mass constraint configuration.
167    for(int j=0; j<_maxtrk; ++j) {
168       _ctvmq.list[j]=0;
169       for(int jv=0; jv<_maxvtx; ++jv) {
170          _ctvmq.trkvtx[jv][j]=false;
171       }
172       for(int jmc=0; jmc<_maxmcn; ++jmc) {
173          _ctvmq.trkmcn[jmc][j]=false;
174       }
175    }
176 // Initialize the conversion and vertex pointing arrays.
177    for(int jv=0; jv<_maxvtx; ++jv) {
178       _ctvmq.cvtx[jv]=0;
179       _ctvmq.vtxpnt[0][jv]=-1;
180       _ctvmq.vtxpnt[1][jv]=0;
181    }
182    _ctvmq.drmax=2.;
183    _ctvmq.dzmax=20.;
184    _ctvmq.rvmax=70.;
185    _ctvmq.trnmax=0.5;
186    _ctvmq.dsmin=-2.;   
187 // Set stat to -999  and chisqq to -1 to symbolize that no fit has yet been done.
188    stat=-999;
189    _ctvmq.chisqr[0]=-1.;
190 }
191 
192 /*=======================================================================
193                          VertexFit::addTrack
194   =======================================================================*/
195 
196 bool VertexFit::addTrack(const CdfTrack* trk, float mass,
197                          vertexNumber jv)
198 {
199 // Check that this vertex number is within the allowed range.
200    if(jv<VERTEX_1 || jv>_maxvtx) return false;
201    _ctvmq.nvertx = jv>_ctvmq.nvertx ? jv : _ctvmq.nvertx;
202 
203 // Check that we have not exceeded the maximum number of tracks.
204    if(_ctvmq.ntrack>=_maxtrk) return false;
205 
206 // Add this track.
207    _ctvmq.list[_ctvmq.ntrack]=trk->id().value();
208    _ctvmq.tkbank[_ctvmq.ntrack][0]='Q';
209    _ctvmq.tkbank[_ctvmq.ntrack][1]='T';
210    _ctvmq.tkbank[_ctvmq.ntrack][2]='R';
211    _ctvmq.tkbank[_ctvmq.ntrack][3]='K';
212    _ctvmq.tmass[_ctvmq.ntrack]=mass;
213    _ctvmq.trkvtx[jv-1][_ctvmq.ntrack]=true;
214 
215 // Put this track's helix parameters and error matrix into
216 // a fortran common block so that they can  be accessed
217 // by gettrk.  This is a dummy for now.
218    _trkprm.trhelix[_ctvmq.ntrack][0]=trk->lambda();
219    _trkprm.trhelix[_ctvmq.ntrack][1]=trk->curvature();
220    _trkprm.trhelix[_ctvmq.ntrack][2]=trk->z0();
221    _trkprm.trhelix[_ctvmq.ntrack][3]=trk->d0();
222    _trkprm.trhelix[_ctvmq.ntrack][4]=trk->phi0();
223 
224    for(int j=0; j<5; ++j) {
225       for(int k=0; k<5; ++k) {
226          _trkprm.trem[_ctvmq.ntrack][j][k]=
227                          trk->getHelixFit().getErrorMatrix()[j][k];
228       }
229    }
230    _ctvmq.ntrack++;
231 
232    return true;
233 }
234 
235 /*=======================================================================
236                          VertexFit::vertexPoint_2d
237   =======================================================================*/
238   
239 bool VertexFit::vertexPoint_2d(vertexNumber jv1, vertexNumber jv2)
240 {
241 // Check that these vertex numbers are within allowed range and
242 // that the vertices are unique.
243    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
244    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
245    if(jv1 <= jv2) return false;
246 
247 // Setup vertex pointing.
248    _ctvmq.vtxpnt[0][jv1-1]=jv2;
249    _ctvmq.vtxpnt[1][jv1-1]=1;           // 2d pointing.
250 
251    return true;
252 }
253 
254 /*=======================================================================
255                          VertexFit::vertexPoint_3d
256   =======================================================================*/
257 
258 bool VertexFit::vertexPoint_3d(vertexNumber jv1, vertexNumber jv2)
259 {
260 // Check that these vertex numbers are within allowed range and
261 // that the vertices are distinct.
262    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
263    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
264    if(jv1 <= jv2) return false;
265 // Setup vertex pointing.
266    _ctvmq.vtxpnt[0][jv1-1]=jv2;
267    _ctvmq.vtxpnt[1][jv1-1]=2;           // 3d pointing.
268 
269    return true;
270 }
271 
272 /*=======================================================================
273                          VertexFit::vertexPoint_1track
274   =======================================================================*/
275 
276 bool VertexFit::vertexPoint_1track(vertexNumber jv1, vertexNumber jv2)
277 {
278 // Check that these vertex numbers are within allowed range and are distinct.
279    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
280    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
281    if(jv1 <= jv2) return false;
282 
283 // Setup vertex pointing.
284    _ctvmq.vtxpnt[0][jv1-1]=jv2;
285    _ctvmq.vtxpnt[1][jv1-1]=3;           // Point to 1 track vertex.
286 
287    return true;
288 }
289 /*=======================================================================
290                          VertexFit::vertexPoint_0track
291   =======================================================================*/
292 
293 bool VertexFit::vertexPoint_0track(vertexNumber jv1, vertexNumber jv2)
294 {
295 
296   // jv2 is the zero track vertex. jv1 is the multi track vertex which points to jv2
297 
298   // Note: You must call this routine at least twice in order for ctvmft_zerotrackvtx to work
299   // ie there must be at least 2 vertices pointed at a zero track vertex. ctvmft
300   // for this, but the error message may not make it to your log file (look at the local
301   // variables in the stack frame especially IJKERR(2). The significance of this
302   // error code is documented at the top of whatever routine chucked you out (ctvmf00 in
303   // this case)
304 
305   // see ctvmft.f source file for discussion. See especially comments at the top
306   // of subroutines: ctvmft and ctvmfa
307 
308 // Check that these vertex numbers are within allowed range and are distinct.
309    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
310    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
311    if(jv1 <= jv2) return false;
312 
313 // Setup vertex pointing.
314    _ctvmq.vtxpnt[0][jv1-1]=jv2;
315    _ctvmq.vtxpnt[1][jv1-1]=4;           // Point to 0 track vertex.
316 
317    return true;
318 }//vertexPoint_0track
319 
320 /*=======================================================================
321                          VertexFit::conversion_2d
322   =======================================================================*/
323 
324 bool VertexFit::conversion_2d(vertexNumber jv)
325 {
326    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
327       return false;
328    }
329    _ctvmq.cvtx[jv-1]=1;
330    return true;
331 }
332 
333 /*=======================================================================
334                          VertexFit::conversion_3d
335   =======================================================================*/
336 
337 bool VertexFit::conversion_3d(vertexNumber jv)
338 {
339    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
340       return false;
341    }
342    _ctvmq.cvtx[jv-1]=2;
343    return true;
344 }
345 
346 /*=======================================================================
347                          VertexFit::massConstrain
348   =======================================================================*/
349 
350 bool VertexFit::massConstrain(const CdfTrack* trk1,
351                               const CdfTrack* trk2, float mass)
352 {
353    int ntrk=2;
354    const CdfTrack* tracks[2];
355    tracks[0]=trk1;
356    tracks[1]=trk2;
357    return massConstrain(ntrk,tracks,mass);
358 }
359 
360 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
361                                 const CdfTrack* trk3, float mass)
362 {
363    int ntrk=3;
364    const CdfTrack* tracks[3];
365    tracks[0]=trk1;
366    tracks[1]=trk2;
367    tracks[2]=trk3;
368    return massConstrain(ntrk,tracks,mass);
369 }
370 
371 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
372                     const CdfTrack* trk3, const CdfTrack* trk4,float mass)
373 {
374    int ntrk=4;
375    const CdfTrack* tracks[4];
376    tracks[0]=trk1;
377    tracks[1]=trk2;
378    tracks[2]=trk3;
379    tracks[3]=trk4;
380    return massConstrain(ntrk,tracks,mass);
381 }
382 
383 bool VertexFit::massConstrain(int ntrk, const CdfTrack* tracks[], float mass)
384 {
385 // Check that we have not exceeded the allowed number of mass constraints.
386    if(_ctvmq.nmassc>=_maxmcn) return false;
387 
388 // Set constraint mass
389    _ctvmq.cmass[_ctvmq.nmassc]=mass;
390      
391 // For each track in contraint, set trkmcn true.  Since the number in
392 // tracks[] is the track number, we have to find each track in the
393 // list of tracks.
394    for(int jt=0; jt<ntrk; ++jt) {
395       bool found=false;
396       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
397          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
398             _ctvmq.trkmcn[_ctvmq.nmassc][kt]=true;
399             found=true;
400          }
401       }
402       if(!found) return false;
403    }
404 // Increment number of mass constraints.   
405    _ctvmq.nmassc++;
406 
407    return true;
408 }
409 
410 
411 /*=======================================================================
412                      VertexFit::beamlineConstraint
413   =======================================================================*/
414 //  Turn ON this option only when fitting the primary with no additional 
415 // constraints. The routine is protected and will not function
416 // if those options are selected as well as the beamline constraint.
417 // In case this option is chosen the fit will calculate a result also with
418 // just one track in input.
419 // In order to enable this option the user must do the following:
420 // Fit the primary vertex as VERTEX_1 (without any other vertices).
421 // Enable beamlineConstraint by setting the pointing constraint variable
422 // vtxpnt[0][0] to -100 (no pointing constraints when <0)
423 // Note that index 0,0 corresponds to index 1,1 in fortran
424 // Provide the beamline parameters:
425 //  --> assume xv = xb + xzbslope*z, 
426 //             yv = yb + yzbslope*z, where (xv, yv) is the transverse
427 //      beam position at z, (xb, yb) is the transverse beam position at z = 0 
428 //      and (xzbslope, yzbslope) are the slopes of the beamline
429 //  --> set primary vertex at z=0 (xb, yb, 0)
430 //  --> set the diagonal elements of the primary vertex error matrix to the
431 //      sigma**2 of beam spread in x, y and z:
432 //      xverr[1][1] =  (~30 - 50 um)**2
433 //      xverr[2][2] =  (~30 - 50 um)**2
434 //      xverr[3][3] =  (~30 cm) **2
435 //      All other elements should be 0
436 // For help email: Joao Guimaraes  guima@huhepl.harvard.edu
437 //                 Franco Bedeschi bed@fnal.gov 
438 // See more information in the ctvmft.f
439 bool VertexFit::beamlineConstraint(float xb, float yb, HepSymMatrix berr,
440                                    float xzbslope, float yzbslope)
441 {
442   // Set beam position at z=0
443   setPrimaryVertex(xb,yb,0);
444   bool success = setPrimaryVertexError(berr);
445   //Set the beamline slope values
446   _ctvmq.xzslope= xzbslope;
447   _ctvmq.yzslope= yzbslope;
448   // Turn ON beamline constraint
449   _ctvmq.vtxpnt[0][0]=-100; 
450 
451   return success;
452 }
453 
454 bool VertexFit::beamlineConstraint(Hep3Vector pv, HepSymMatrix berr,
455                                    float xzbslope, float yzbslope)
456 {
457   // Check if input beam position coordinates are at z=0
458   if (pv.z() != 0) return false;
459 
460   return beamlineConstraint(pv.x(),pv.y(),berr,xzbslope,yzbslope);
461 }
462 
463 /*=======================================================================
464                          VertexFit::setPrimaryVertex
465   =======================================================================*/
466 
467 void VertexFit::setPrimaryVertex(float xv, float yv, float zv)
468 {
469 // Set x,y,z position of the primary vertex.
470    _ctvmq.xyzpv0[0]=xv;
471    _ctvmq.xyzpv0[1]=yv;
472    _ctvmq.xyzpv0[2]=zv;
473 
474    return;
475 }
476 
477 void VertexFit::setPrimaryVertex(Hep3Vector pv)
478 {
479 // Set x,y,z position of the primary vertex.
480    _ctvmq.xyzpv0[0]=pv.x();
481    _ctvmq.xyzpv0[1]=pv.y();
482    _ctvmq.xyzpv0[2]=pv.z();
483 
484    return;
485 }
486 
487 /*=======================================================================
488                          VertexFit::setPrimaryVertexError
489   =======================================================================*/
490 
491 void VertexFit::setPrimaryVertexError(const float xverr[3][3])
492 {
493 // Set the error matrix for the primary vertex.
494    for(int j=0; j<3; ++j) {
495       for(int k=0; k<3; ++k) {
496          _ctvmq.exyzpv[j][k]=xverr[j][k];
497       }
498    }
499    return;
500 }
501 
502 bool VertexFit::setPrimaryVertexError(const HepSymMatrix &xverr)
503 {
504 // Set the error matrix for the primary vertex using a HepSymMatrix.
505 // First check that the matrix is the correct size.
506    if(xverr.num_row() != 3) return false;
507    for(int j=0; j<3; j++) {
508      for(int k=0; k<3; k++) {
509        _ctvmq.exyzpv[j][k]=xverr[j][k];
510      }
511    }
512    return true;
513 }
514 
515 /*=======================================================================
516                          VertexFit::fit
517   =======================================================================*/
518 
519 bool VertexFit::fit()
520 {
521 // Check that the diagonal elements of all the track error matrices
522 //are positive.
523   bool mstat = true;
524   for(int trk=0; trk<_ctvmq.ntrack; ++trk) {
525     for(int j=0; j<5; ++j) {
526 // Check diagonal elements of error matrix.
527       if(_trkprm.trem[trk][j][j] < 0.) {
528         mstat = false;
529 // Record offending track
530         _ctvmq.ijkerr[2] = trk + 1;
531       }
532     }
533   }
534   if(!mstat) {
535 // Tell user that the covariance matrix could not be inverted.
536     _ctvmq.ijkerr[0] = 3;
537     _ctvmq.ijkerr[1] = 2;
538     stat = 1;
539     return false;
540   }
541 
542 // First copy information into CTVMFT common blocks
543    *_ctvmq_com=_ctvmq;
544    *_ctvmfr_com=_ctvmfr;
545    *_fiddle_com=_fiddle;
546    *_trkprm_com=_trkprm;
547 // Do the vertex fit.
548    int print=0;
549    int level=0;
550    ctvmft_(print,level,stat);
551 // Now copy information from CTVMFT common blocks to local storage
552    _ctvmq=*_ctvmq_com;
553    _ctvmfr=*_ctvmfr_com;
554    _fiddle=*_fiddle_com;
555    _trkprm=*_trkprm_com;
556    
557    return stat==0;
558 }
559 
560 /*=======================================================================
561                          VertexFit::print
562   =======================================================================*/
563 
564 void VertexFit::print() const
565 {
566    print(std::cout);
567 }
568 
569 void VertexFit::print(std::ostream& os) const
570 {
571    os << "****************************** VertexFit "
572       << "******************************" << std::endl;
573    os << "Number of tracks: " << _ctvmq.ntrack << std::endl;
574    os << "   Tracks: ";
575    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
576       if(jt != 0) os << ",  ";
577       CdfTrack::Id trkId(_ctvmq.list[jt]);
578       os << trkId;
579    }
580    os << std::endl;
581    os << "Number of vertices: " << _ctvmq.nvertx << std::endl;
582    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
583       os << "   Vertex " << jv+1 << " tracks: ";
584       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
585          if(_ctvmq.trkvtx[jv][jt]) {
586             os << " " << _ctvmq.list[jt];
587          }
588       }
589       os << std::endl;
590    }
591    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
592       if(_ctvmq.vtxpnt[0][jv]==0) {
593          os << "   Vertex " << jv+1 << " points to the primary vertex ";
594       }
595       else if(_ctvmq.vtxpnt[0][jv]>0) {
596          os << "   Vertex " << jv+1 << " points to vertex "
597               << _ctvmq.vtxpnt[0][jv];
598       }
599       if(_ctvmq.vtxpnt[1][jv]==1) {
600          os << " in 2 dimensions" << std::endl;
601       }
602       else if(_ctvmq.vtxpnt[1][jv]==2) {
603          os << " in 3 dimensions" << std::endl;
604       }
605       else if(_ctvmq.vtxpnt[1][jv]==3) {
606          os << ", a single track vertex" << std::endl;
607       }
608       if(_ctvmq.cvtx[jv]>0) {
609          os << "   Vertex " << jv+1 << " is a conversion" << std::endl;
610       }
611    }
612    os << "Number of mass constraints: " << _ctvmq.nmassc << std::endl;
613    for(int jmc=0; jmc<_ctvmq.nmassc; ++jmc) {
614       os << "   Tracks ";
615       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
616          if(_ctvmq.trkmcn[jmc][jt]) {
617             os << " " << _ctvmq.list[jt];
618          }
619       }
620       os << " constrained to mass " << _ctvmq.cmass[jmc]
621          << " Gev/c^2" << std::endl;
622    }
623    if(stat==-999) {
624       os << "No fit has been done." << std::endl;
625    }
626    else {
627       os << "***** Results of Fit *****" << std::endl;
628       printErr(os);
629       os << "   Status = " << stat << std::endl;
630       os << "   Chi-square = " << _ctvmq.chisqr[0]
631          << " for " << _ctvmq.ndof << " degrees of freedom." << std::endl;
632       os << "   => probability = " << prob() << std::endl;
633       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
634          os << "Vertex " << jv+1 << " position: "
635             << _ctvmq.xyzvrt[jv+1][0] << " "
636             << _ctvmq.xyzvrt[jv+1][1] << " "
637             << _ctvmq.xyzvrt[jv+1][2] << std::endl;
638       }
639       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
640          os << "Track " << _ctvmq.list[jt] << " P4: " 
641             << _ctvmq.trkp4[0][jt] << " "
642             << _ctvmq.trkp4[1][jt] << " "
643             << _ctvmq.trkp4[2][jt] << " "
644             << _ctvmq.trkp4[3][jt] << " " << std::endl;
645       }
646    }
647    os << "****************************************"
648       << "**************************" << std::endl;
649 
650    return;
651 }
652 
653 /*=======================================================================
654                          VertexFit::printErr
655   =======================================================================*/
656 
657 void VertexFit::printErr() const
658 {
659    printErr(std::cout);
660 }
661 
662 void VertexFit::printErr(std::ostream& os) const
663 {
664    os << "VertexFit: IJKERR = " << _ctvmq.ijkerr[0] << ", "
665                         << _ctvmq.ijkerr[1] << ", "
666                         << _ctvmq.ijkerr[2] << std::endl;
667    if(status()==0 && _ctvmq.ijkerr[0]==0) return;
668    if(_ctvmq.ijkerr[0] == -1) {
669       os << "   Problem with GETTRK:  track requested is not in list."
670          << std::endl
671          << "   This should not happen - Contact VertexFit expert "
672          << _expert << "." <<std::endl;
673    }
674    else if(_ctvmq.ijkerr[0]==1) {
675       os << "   Problem in CTVM00:" << std::endl;
676       if(_ctvmq.ijkerr[1]==1) {
677          os << "      Number of tracks is " << _ctvmq.ntrack
678             << "." << std::endl;
679          if(_ctvmq.ntrack < 2) {
680             os << ", which is too few (must be at least 2)." << std::endl;
681          }
682          else if(_ctvmq.ntrack > _maxtrk) {
683             os << ", which is too many (maximum is " << _maxtrk
684                << ")." << std::endl;
685          }
686          else {
687             os << "      Problem with number of tracks"
688                << " for unknown reasons." << std::endl;
689          }
690       }
691       else if(_ctvmq.ijkerr[1]==2) {
692          os << "      Number of vertices is " << _ctvmq.nvertx
693             << "." << std::endl;
694          if(_ctvmq.nvertx < 1) {
695             os << ", which is too few (must be at least 1)." << std::endl;
696          }
697          else if(_ctvmq.nvertx > _maxvtx) {
698             os << ", which is too many (maximum is " << _maxvtx
699                << ")." << std::endl;
700          }
701          else {
702             os << std::endl << "      Problem with number of vertices"
703                << " for unknown reasons." << std::endl;
704          }
705       }
706       else if(_ctvmq.ijkerr[1]==3) {
707          os << "      Number of mass constraints is " << _ctvmq.nmassc
708             << "." << std::endl;
709          if(_ctvmq.nmassc < 0) {
710             os << ", which is negative." << std::endl;
711          }
712          else if(_ctvmq.nmassc > _maxmcn) {
713             os << ", which is too many (maximum is " << _maxmcn
714                << ")." << std::endl;
715          }
716          else {
717             os << std::endl << "      Problem with number of mass"
718                << " constraints for unknown reasons." << std::endl;
719          }
720       }
721       else if(_ctvmq.ijkerr[1]==11) {
722          os << "      Vertex " << _ctvmq.ijkerr[2]
723             << " has less than one track." << std::endl;
724       }
725       else if(_ctvmq.ijkerr[1]==12) {
726          os << "      Vertex " << _ctvmq.ijkerr[2]
727             << " is a conversion vertex with a number of tracks"
728             << " different than two." << std::endl;
729       }
730       else if(_ctvmq.ijkerr[1]==13) {
731          os << "      Vertex " << _ctvmq.ijkerr[2]
732             << " is a one track vertex that has no multi-track"
733             << " descendents." << std::endl;
734       }
735       else if(_ctvmq.ijkerr[1]==14) {
736          os << "      Vertex " << _ctvmq.ijkerr[2]
737             << " does not point at a vertex with a lower number."
738             << std::endl;
739       }
740       else if(_ctvmq.ijkerr[1]==15) {
741          os << "      Vertex " << _ctvmq.ijkerr[2]
742             << " has a parent vertex that is a conversion." << std::endl;
743       }
744       else if(_ctvmq.ijkerr[1]==16) {
745          os << "      Vertex " << _ctvmq.ijkerr[2]
746             << " does 1 track pointing to a vertex with"
747             << " more than 1 track." << std::endl;
748       }
749       else if(_ctvmq.ijkerr[1]==17) {
750          os << "      Vertex " << _ctvmq.ijkerr[2]
751             << " does 0 track pointing to a vertex with"
752             << " more than 0 track (?)." << std::endl;
753       }
754       else if(_ctvmq.ijkerr[1]==19) {
755          os << "      Primary vertex error matrix is singular." << std::endl;
756       }
757       else if(_ctvmq.ijkerr[1]==21) {
758          os << "      Track with Id " << _ctvmq.ijkerr[2]
759             << "is not in any vertex." << std::endl;
760       }
761       else if(_ctvmq.ijkerr[1]==22) {
762          os << "      Track with Id " << _ctvmq.ijkerr[2]
763             << "is in multiple vertices." << std::endl;
764       }
765       else if(_ctvmq.ijkerr[1]==23) {
766          os << "      Track with Id " << _ctvmq.ijkerr[2]
767             << "occurs more than once." << std::endl;
768       }
769       else if(_ctvmq.ijkerr[1]==31) {
770          os << "      A mass constraint has less than 2 tracks." << std::endl;
771       }
772       else if(_ctvmq.ijkerr[1]==32) {
773          os << "      The sum masses of the tracks in a mass constraint"
774             << " exceeds the constraint mass." << std::endl;
775       }
776       else if(_ctvmq.ijkerr[1]==33) {
777          os << "      Beamline constraint. Beam covariance not set properly."
778             << " Negative diagonal elements." << std::endl;
779       }
780       else if(_ctvmq.ijkerr[1]==34) {
781          os << "      Beamline constraint. Beam covariance not set properly."
782             << " Off-diagonal elements not zero." << std::endl;
783       }
784       else if(_ctvmq.ijkerr[1]==36) {
785          os << "      Beamline constraint. Number of vertices = " 
786             << _ctvmq.nvertx << " Should be 1." << std::endl;
787       }
788    }
789    else if(_ctvmq.ijkerr[0] == 2) {
790       if(_ctvmq.ijkerr[1] == 20) {
791          os << "   Problem in CTVM00: " << std::endl;
792          os << "      Track has negative Id = "
793             << _ctvmq.list[_ctvmq.ijkerr[2]-1] << "." << std::endl;
794       }
795       else {
796          os << "   Problem in CTVMFA with vertex "
797             << _ctvmq.ijkerr[2] << ": " << std::endl;
798          os << "      Failure in vertex first approximation." << std::endl;
799          if(_ctvmq.ijkerr[1] == 1) {
800             os << "      Tracks are concentric circles." << std::endl;
801          }
802          if(_ctvmq.ijkerr[1] == 2) {
803             os << "      Conversion vertex has widely separated"
804                << " exterior circles at midpoint." << std::endl;
805          }
806          if(_ctvmq.ijkerr[1] == 3) {
807             os << "      Conversion vertex has widely separated"
808                << " interior circles at midpoint." << std::endl;
809          }
810          if(_ctvmq.ijkerr[1] == 4) {
811             os << "      Vertex has widely separated"
812                << " exterior circles at approximate vertex." << std::endl;
813          }
814          if(_ctvmq.ijkerr[1] == 5) {
815             os << "      Vertex has widely separated"
816                << " interior circles at approximate vertex." << std::endl;
817          }
818          if(_ctvmq.ijkerr[1] == 6) {
819             os << "      Rv is too large at the chosen"
820                << " intersection point." << std::endl;
821          }
822          if(_ctvmq.ijkerr[1] == 7) {
823             os << "      Delta z is too large at the chosen"
824                << " intersection point." << std::endl;
825          }
826          if(_ctvmq.ijkerr[1] == 8) {
827             os << "      A track's turning to the chosen vertex"
828                << " is too large." << std::endl;
829          }
830          if(_ctvmq.ijkerr[1] == 9) {
831             os << "      There is no solution with an adequately"
832                << " positive arc length." << std::endl;
833          }
834          if(_ctvmq.ijkerr[1] == 21) {
835             os << "      zero-track vertexing: either/both vertex "
836                << " momenta are too small (<0.01 MeV)." << std::endl;
837          }
838          if(_ctvmq.ijkerr[1] == 22) {
839             os << "      zero-track vertexing: Two lines (tracks) are "
840                << " parallel/antiparallel." << std::endl;
841          }
842 
843       }
844    }
845    else if(_ctvmq.ijkerr[0] == 3) {
846       os << "   Problem in CTVM01 with track with Id = "
847          << _ctvmq.list[_ctvmq.ijkerr[2]-1] << ": " << std::endl;
848       if(_ctvmq.ijkerr[1] == 1) {
849          os << "      GETTRK cannot find Id in list." << std::endl;
850       }
851       if(_ctvmq.ijkerr[1] == 2) {
852          os << "      Covariance matrix could not be inverted."  << endl 
853             << "      Offending track number (in order addded) is "
854             << _ctvmq.ijkerr[2] << "." << std::endl;
855       }
856       if(_ctvmq.ijkerr[1] == 3) {
857          os << "      Track turns through too large an angle"
858             << " to the vertex." << std::endl;
859       }
860       if(_ctvmq.ijkerr[1] == 4) {
861          os << "      Track moves too far backward to vertex." << std::endl;
862       }
863    }
864    else if(status() == 9) {
865       os << "   General fit problem: " << std::endl;
866       if(_ctvmq.ijkerr[1] == 1) {
867          os << "      Singular solution matrix." << std::endl;
868       }
869       if(_ctvmq.ijkerr[1] == 2 || _ctvmq.ijkerr[1] == 3) {
870          os << "      Too many iterations ( "
871             << _ctvmq.ijkerr[2] << "(." << std::endl;
872       }
873       if(_ctvmq.ijkerr[1] == 4) {
874          os << "      Convergence failure." << std::endl;
875       }
876       if(_ctvmq.ijkerr[1] == 5) {
877          os << "      Bad convergence." << std::endl;
878       }
879       if(_ctvmq.ijkerr[1] == 9) {
880          os << "      Ill-formed  covariance matrix." << std::endl;
881       }
882    }
883    else {
884       os << "   The error codes above are not recognized." << std::endl
885          << "   Contact VertexFit expert " << _expert << "." << std::endl;
886    }
887    return;
888 }
889 
890 
891 /*=======================================================================
892                          VertexFit::getIJKErr
893   =======================================================================*/
894 
895 void VertexFit::getIJKErr(int& err0, int& err1, int& err2) const
896 {
897   err0 = _ctvmq.ijkerr[0];
898   err1 = _ctvmq.ijkerr[1];
899   err2 = _ctvmq.ijkerr[2];
900   return;
901 }
902 
903 int VertexFit::getIJKErr0() const
904 {
905   return _ctvmq.ijkerr[0];
906 }
907 int VertexFit::getIJKErr1() const
908 {
909   return _ctvmq.ijkerr[1];
910 }
911 int VertexFit::getIJKErr2() const
912 {
913   return _ctvmq.ijkerr[2];
914 }
915 
916 /*=======================================================================
917                        VertexFit::getErrTrackId()
918   =======================================================================*/
919 
920 int VertexFit::getErrTrackId() const
921  {
922    if (status() == 0) return 0;
923    int trkId = 0;
924    // Problems with track in CTVM01 or track has negative id in CTVM00
925    // See PrintErr() for a more detailed list of error codes.
926    if ((_ctvmq.ijkerr[0] == 2 && _ctvmq.ijkerr[1] == 20) ||
927        _ctvmq.ijkerr[0] == 3) {
928      trkId = _ctvmq.list[_ctvmq.ijkerr[2]-1];
929    }
930    
931    return trkId;
932  }
933 
934 /*=======================================================================
935                          VertexFit::expert
936   =======================================================================*/
937 
938 string VertexFit::expert() const
939 {
940    return _expert;
941 }
942 
943 /*=======================================================================
944                          VertexFit::status
945   =======================================================================*/
946 
947 int VertexFit::status() const
948 {
949 
950    return stat;
951 }
952 
953 /*=======================================================================
954                          VertexFit::chisq
955   =======================================================================*/
956 
957 float VertexFit::chisq() const
958 {
959 // Chi-square of fit
960       return _ctvmq.chisqr[0];
961 }
962 
963 /*=======================================================================
964                          VertexFit::ndof
965   =======================================================================*/
966 
967 int VertexFit::ndof() const
968 {
969 // Number of degrees of freedom of fit.
970    if(_ctvmq.chisqr[0] >= 0) {
971       return _ctvmq.ndof; }
972    else {
973       return 0;}
974 }
975 
976 /*=======================================================================
977                          VertexFit::prob
978   =======================================================================*/
979 
980 float VertexFit::prob() const
981 {
982 // Probability of chi-square of fit
983    if(_ctvmq.chisqr[0]>=0.) {
984       float chisq=_ctvmq.chisqr[0];
985       int nd=_ctvmq.ndof;
986       return prob_(chisq,nd);
987    }
988    else {
989       return -999.;
990    }
991 }
992 
993 /*=======================================================================
994                          VertexFit::chisq(CdfTrack*)
995   =======================================================================*/
996 
997 float VertexFit::chisq(const CdfTrack* trk) const
998 {
999 // This method returns the chisquare contribution for one track 
1000 // If fit not successful or not done yet, return -1.
1001    if(_ctvmq.chisqr[0] < 0) return -1.;
1002 // Look for this track in the list of tracks.
1003    for(int jt = 0; jt < _ctvmq.ntrack; ++jt) {
1004      if(trk->id().value() == _ctvmq.list[jt]) {
1005 // Found the track, so return its chisquare contribution.
1006        return _ctvmq.chit[jt];
1007      }
1008    }
1009 // If track is not in list, return -1.
1010    return -1.;
1011 }
1012 
1013 /*=======================================================================
1014                          VertexFit::chisq_rphi()
1015   =======================================================================*/
1016 
1017 float VertexFit::chisq_rphi() const
1018 {
1019 // This method returns the chisquare contribution in the r-phi plane.
1020    int index[3] = {0,1,3};
1021 // If fit not successful or not done yet, return -1.
1022    if(_ctvmq.chisqr[0] < 0) return -1.;
1023 // Loop over the tracks in the event.
1024    float chisq = 0.;
1025    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1026 // Double loop over the relevant parameter indices.
1027      for(int k1=0; k1<3; ++k1) {
1028        for(int k2=0; k2<3; ++k2) {
1029 // Add contribution to chisquare.
1030          chisq += _ctvmq.pardif[jt][index[k1]] *
1031                   _ctvmq.g[jt][index[k1]][index[k2]] *
1032                   _ctvmq.pardif[jt][index[k2]];
1033        }
1034      }
1035    }
1036 // Return the chisquare.
1037    return chisq;
1038 }
1039 
1040 
1041 /*=======================================================================
1042                          VertexFit::chisq_z()
1043   =======================================================================*/
1044 
1045 float VertexFit::chisq_z() const
1046 {
1047 // This method returns the chisquare contribution in the z direction.
1048    int index[2] = {2,4};
1049 // If fit not successful or not done yet, return -1.
1050    if(_ctvmq.chisqr[0] < 0) return -1.;
1051 // Loop over the tracks in the event.
1052    float chisq = 0.;
1053    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1054 // Double loop over the relevant parameter indices.
1055      for(int k1=0; k1<2; ++k1) {
1056        for(int k2=0; k2<2; ++k2) {
1057 // Add contribution to chisquare.
1058          chisq += _ctvmq.pardif[jt][index[k1]] *
1059                   _ctvmq.g[jt][index[k1]][index[k2]] *
1060                   _ctvmq.pardif[jt][index[k2]];
1061        }
1062      }
1063    }
1064 // Return the chisquare.
1065    return chisq;
1066 }
1067 
1068 
1069 /*=======================================================================
1070                          VertexFit::chisq_rphiz()
1071   =======================================================================*/
1072 
1073 float VertexFit::chisq_rphiz() const
1074 {
1075 // This method returns the chisquare contribution of the cross
1076 // terms in the r-phi and z directions.
1077    int index1[2] = {2,4};
1078    int index2[3] = {0,1,3};
1079 // If fit not successful or not done yet, return -1.
1080    if(_ctvmq.chisqr[0] < 0) return -1.;
1081 // Loop over the tracks in the event.
1082    float chisq = 0.;
1083    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1084 // Double loop over the relevant parameter indices.
1085      for(int k1=0; k1<2; ++k1) {
1086        for(int k2=0; k2<3; ++k2) {
1087 // Add contribution to chisquare.
1088          chisq += _ctvmq.pardif[jt][index1[k1]] *
1089                   _ctvmq.g[jt][index1[k1]][index2[k2]] *
1090                   _ctvmq.pardif[jt][index2[k2]];
1091        }
1092      }
1093    }
1094 // Return the chisquare.
1095    return 2.*chisq;
1096 }
1097 
1098 /*=======================================================================
1099                          VertexFit::chisq_rphi(CdfTrack*)
1100   =======================================================================*/
1101 
1102 float VertexFit::chisq_rphi(const CdfTrack* trk) const
1103 {
1104 // This method returns the chisquare contribution in the r-phi plane.
1105    int index[3] = {0,1,3};
1106 // If fit not successful or not done yet, return -1.
1107    if(_ctvmq.chisqr[0] < 0) return -1.;
1108 // Loop over the tracks in the event, looking for the one we want
1109    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1110      if(trk->id().value() == _ctvmq.list[jt]) {
1111 // Found the track, so calculate its chisquare contribution.
1112        float chisq = 0.;
1113 // Double loop over the relevant parameter indices.
1114        for(int k1=0; k1<3; ++k1) {
1115          for(int k2=0; k2<3; ++k2) {
1116 // Add contribution to chisquare.
1117            chisq += _ctvmq.pardif[jt][index[k1]] *
1118                     _ctvmq.g[jt][index[k1]][index[k2]] *
1119                     _ctvmq.pardif[jt][index[k2]];
1120          }
1121        }
1122        return chisq;
1123      }
1124    }
1125 // Track not found, return -1.
1126    return -1.;
1127 }
1128 
1129 
1130 /*=======================================================================
1131                          VertexFit::chisq_z(CdfTrack*)
1132   =======================================================================*/
1133 
1134 float VertexFit::chisq_z(const CdfTrack* trk) const
1135 {
1136 // This method returns the chisquare contribution in the z direction.
1137    int index[2] = {2,4};
1138 // If fit not successful or not done yet, return -1.
1139    if(_ctvmq.chisqr[0] < 0) return -1.;
1140 // Loop over the tracks in the event, looking for the one we want.
1141    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1142      if(trk->id().value() == _ctvmq.list[jt]) {
1143 // Found the track, so calculate its chisquare contribution.
1144        float chisq = 0.;
1145 // Double loop over the relevant parameter indices.
1146        for(int k1=0; k1<2; ++k1) {
1147          for(int k2=0; k2<2; ++k2) {
1148 // Add contribution to chisquare.
1149            chisq += _ctvmq.pardif[jt][index[k1]] *
1150                     _ctvmq.g[jt][index[k1]][index[k2]] *
1151                     _ctvmq.pardif[jt][index[k2]];
1152          }
1153        }
1154        return chisq;
1155      }
1156    }
1157 // Track not found - return -1.
1158    return -1.;
1159 }
1160 
1161 
1162 /*=======================================================================
1163                          VertexFit::chisq_rphiz(CdfTrack*)
1164   =======================================================================*/
1165 
1166 float VertexFit::chisq_rphiz(const CdfTrack* trk) const
1167 {
1168 // This method returns the chisquare contribution of the cross
1169 // terms in the r-phi and z directions.
1170    int index1[2] = {2,4};
1171    int index2[3] = {0,1,3};
1172 // If fit not successful or not done yet, return -1.
1173    if(_ctvmq.chisqr[0] < 0) return -1.;
1174 // Loop over the tracks in the event.
1175    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1176      if(trk->id().value() == _ctvmq.list[jt]) {
1177 // Found the track, so calculate its chisquare contribution.
1178        float chisq = 0.;
1179 // Double loop over the relevant parameter indices.
1180        for(int k1=0; k1<2; ++k1) {
1181          for(int k2=0; k2<3; ++k2) {
1182 // Add contribution to chisquare.
1183            chisq += _ctvmq.pardif[jt][index1[k1]] *
1184                     _ctvmq.g[jt][index1[k1]][index2[k2]] *
1185                     _ctvmq.pardif[jt][index2[k2]];
1186          }
1187        }
1188        return 2.*chisq;
1189      }
1190    }
1191 // Track not found so return -1.
1192    return -1.;
1193 }
1194 
1195 /*=======================================================================
1196                          VertexFit::getTrackP4
1197   =======================================================================*/
1198 
1199 HepLorentzVector VertexFit::getTrackP4(const CdfTrack* trk) const
1200 {
1201    if(stat != 0) return HepLorentzVector(0,0,0,0);
1202 // return four momentum of fit track 
1203    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1204 //   Find which track matches this Id.
1205       if(trk->id().value()==_ctvmq.list[jt]) {
1206          HepLorentzVector p4((double)_ctvmq.trkp4[0][jt],
1207                              (double)_ctvmq.trkp4[1][jt],
1208                              (double)_ctvmq.trkp4[2][jt],
1209                              (double)_ctvmq.trkp4[3][jt]);
1210          return p4;
1211       }
1212    }
1213    return HepLorentzVector(0,0,0,0);
1214 }
1215 
1216 /*=======================================================================
1217                          VertexFit::getHelix
1218   =======================================================================*/
1219 
1220 Helix VertexFit::getHelix(const CdfTrack* trk) const
1221 {
1222    if(stat != 0) return Helix(0,0,0,0,0);
1223 // return helix parameters of fit track 
1224    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1225 //   Find which track matches this Id.
1226       if(trk->id().value()==_ctvmq.list[jt]) {
1227          Helix trhelix((double)_ctvmq.par[jt][2],
1228                        (double)_ctvmq.par[jt][0],
1229                        (double)_ctvmq.par[jt][4],
1230                        (double)_ctvmq.par[jt][3],
1231                        (double)_ctvmq.par[jt][1]);
1232          return trhelix;
1233       }
1234    }
1235    return Helix(0,0,0,0,0);
1236 }
1237 
1238 /*=======================================================================
1239                          VertexFit::getMass
1240   =======================================================================*/
1241 
1242 float VertexFit::getMass(const CdfTrack* trk1,
1243                          const CdfTrack* trk2, float& dmass) const
1244 {
1245 
1246 //
1247 // Vertex Fitting is prone to floating point errors.  Here, we
1248 // turn those off during the duration of this routine: Joe and Farrukh
1249 //
1250 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1251   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1252   sigaction(SIGFPE, &myaction, &oldaction);
1253   if (setjmp(env)!=0) {
1254     sigaction(SIGFPE, &oldaction,0);
1255     return -999;
1256   }
1257 #endif
1258 
1259    dmass=-999.;
1260    if(stat!=0) return -999.;
1261 // Get the fit invariant mass of tracks trk1 and trk2.
1262 // dmass is the error on the mass.
1263    int ntrk=2;
1264    const CdfTrack* trks[2];
1265    trks[0]=trk1;
1266    trks[1]=trk2;
1267 
1268   double result = getMass(ntrk,trks,dmass);
1269    //return something for FPE Protection-> Joe and Farrukh
1270 
1271 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1272   sigaction(SIGFPE, &oldaction,0);
1273 #endif
1274 
1275   return result;
1276 }
1277 
1278 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
1279                                 const CdfTrack* trk3,float& dmass) const
1280 {
1281 
1282 //
1283 // Vertex Fitting is prone to floating point errors.  Here, we
1284 // turn those off during the duration of this routine: Joe and Farrukh
1285 //
1286 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1287   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1288   sigaction(SIGFPE, &myaction, &oldaction);
1289   if (setjmp(env)!=0) {
1290     sigaction(SIGFPE, &oldaction,0);
1291     return -999;
1292   }
1293 #endif
1294 
1295 
1296 
1297    dmass=-999.;
1298    if(stat!=0) return -999.;
1299 // Get the fit invariant mass of tracks trk1, trk2, and trk3.
1300 // dmass is the error on the mass.
1301    int ntrk=3;
1302    const CdfTrack* trks[3];
1303    trks[0]=trk1;
1304    trks[1]=trk2;
1305    trks[2]=trk3;
1306 
1307 
1308    double result = getMass(ntrk,trks,dmass);
1309 
1310 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1311   sigaction(SIGFPE, &oldaction,0);
1312 #endif
1313 
1314   return result;
1315 }
1316 
1317 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
1318      const CdfTrack* trk3, const CdfTrack* trk4, float& dmass) const
1319 {
1320    dmass=-999.;
1321    if(stat!=0) return -999.;
1322 // Get the fit invariant mass of tracks trk1, trk2, trk3, and trk4.
1323 // dmass is the error on the mass.
1324    int ntrk=4;
1325    const CdfTrack* trks[4];
1326    trks[0]=trk1;
1327    trks[1]=trk2;
1328    trks[2]=trk3;
1329    trks[3]=trk4;
1330 //
1331 // Vertex Fitting is prone to floating point errors.  Here, we
1332 // turn those off during the duration of this routine: Joe and Farrukh
1333 //
1334 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1335   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1336   sigaction(SIGFPE, &myaction, &oldaction);
1337   if (setjmp(env)!=0) {
1338     sigaction(SIGFPE, &oldaction,0);
1339     return -999;
1340   }
1341 #endif
1342 
1343   double result = getMass(ntrk,trks,dmass);
1344 
1345    //return something for FPE Protection-> Joe and Farrukh
1346 
1347 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1348   sigaction(SIGFPE, &oldaction,0);
1349 #endif
1350 
1351   return result;
1352 }
1353 
1354 float VertexFit::getMass(int ntrk, const CdfTrack* tracks[],
1355                                                 float& dmass) const
1356 {
1357 
1358 // Vertex Fitting is prone to floating point errors.  Here, we
1359 // turn those off during the duration of this routine: Joe and Farrukh
1360 //
1361 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1362   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1363   sigaction(SIGFPE, &myaction, &oldaction);
1364   if (setjmp(env)!=0) {
1365     sigaction(SIGFPE, &oldaction,0);
1366     return -999;
1367   }
1368 #endif
1369 
1370    dmass=-999.;
1371    if(stat!=0) return -999.;
1372 // Get fit invariant mass of ntrk tracks listed in array tracks.
1373 // dmass is the error on the mass.
1374    dmass=-999.;
1375    if(ntrk <= 0) return 0;
1376    int jtrks[_maxtrk];
1377    for(int jt=0; jt<ntrk; ++jt) {
1378       bool found=false;
1379       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
1380          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
1381             found=true;
1382             jtrks[jt]=kt+1;
1383          }
1384       }
1385       if(!found) return 0;
1386    }
1387 // Copy information into CTVMFT common blocks
1388    *_ctvmq_com=_ctvmq;
1389    *_ctvmfr_com=_ctvmfr;
1390    int ntr=ntrk;
1391    float mass;
1392    double p4[4];
1393    mcalc_(ntr, jtrks, mass, dmass, p4);
1394 
1395    //return something for FPE-> Joe and Farrukh
1396 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1397   sigaction(SIGFPE, &oldaction,0);
1398 #endif
1399 
1400    return mass;
1401 }
1402 
1403 /*=======================================================================
1404                          VertexFit::getDecayLength
1405   =======================================================================*/
1406 
1407 float VertexFit::getDecayLength(vertexNumber nv, vertexNumber mv, 
1408                                 const Hep3Vector& dir, float& dlerr) const
1409 {
1410    dlerr=-999.;
1411    if(stat!=0) return -999.;
1412 // Get the signed decay length from vertex nv to vertex mv
1413 // along the x-y direction defined by the 3 vector dir.
1414 // dlerr is the error on the decay length including the
1415 // full error matrix.
1416 // Check that vertices are within range.
1417    if(nv<0 || nv>=_ctvmq.nvertx) return -999.;
1418    if(mv<1 || mv>_ctvmq.nvertx) return -999.;
1419    if(nv>=mv) return -999.;
1420    float dir_t=sqrt(dir.x()*dir.x()+dir.y()*dir.y());
1421    if(dir_t <= 0.) return -999.;
1422    Hep3Vector dv=getVertex(mv)-getVertex(nv);
1423    float dl=(dv.x()*dir.x()+dv.y()*dir.y())/dir_t;
1424 // Set up the column matrix of derivatives
1425    HepMatrix A(4,1);
1426    A(1,1)=dir.x()/dir_t;
1427    A(2,1)=dir.y()/dir_t;
1428    A(3,1)=-dir.x()/dir_t;
1429    A(4,1)=-dir.y()/dir_t;
1430 // Check if first vertex (nv) is primary vertex.  If it is,
1431 // check if it was used in the primary vertex.  If not,
1432 // all of the corresponding error matrix elements are those
1433 // supplied for the primary vertex.
1434    int nvf=0;
1435    if(nv==0) {
1436       nvf=-1;
1437       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
1438          if(_ctvmq.vtxpnt[jv][0]==0) nvf=0;
1439       }
1440    }
1441 // Get the relevant submatrix of the full error matrix.
1442    HepMatrix V(4,4,0);
1443    if(nvf < 0) {
1444       V(1,1)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+1);
1445       V(1,2)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2);
1446       V(2,1)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+1);
1447       V(2,2)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+2);
1448       V(3,3)=_ctvmq.exyzpv[0][0];
1449       V(3,4)=_ctvmq.exyzpv[0][1];
1450       V(4,3)=_ctvmq.exyzpv[1][0];
1451       V(4,4)=_ctvmq.exyzpv[1][1];
1452    }
1453    else {
1454 //    Get the indices into the error matrix vmat
1455       int index[4]={_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2,0,0};
1456       if(nv==0) {
1457          index[2]=1;
1458          index[3]=2;
1459       }
1460       else {
1461         index[2] = _ctvmq.voff[nv-1]+1;
1462         index[3] = _ctvmq.voff[nv-1]+2;
1463       }
1464       for(int j=0; j<4; ++j) {
1465          for(int k=0; k<4; ++k) {
1466             V[j][k]=getErrorMatrix(index[j],index[k]);
1467          }
1468       }
1469    }
1470 // Calculate square of dlerr
1471    dlerr=(A.T()*V*A)(1,1);
1472    if(dlerr >= 0.) {
1473       dlerr=sqrt(dlerr);
1474    }
1475    else {
1476       dlerr=-sqrt(-dlerr);
1477    }
1478 
1479    return dl;
1480 }
1481 
1482 /*=======================================================================
1483                          VertexFit::get_dr
1484   =======================================================================*/
1485 
1486 float VertexFit::get_dr(vertexNumber nv, vertexNumber mv,
1487                                            float& drerr) const
1488 {
1489    drerr=-999.;
1490    if(stat!=0) return -999.;
1491 // Get the transvese distance between vertices nv and mv
1492 // and return it as the function value.  drerr is the
1493 // uncertainty on the transverse distance, calculated from the
1494 // full error matrix including correlations.
1495    float dxyz[3];
1496    float dr;
1497    float dz;
1498    float dl[3];
1499 
1500    int mvert=mv;
1501    int nvert=nv;
1502 // Copy information into CTVMFT common blocks
1503    *_ctvmq_com=_ctvmq;
1504    *_ctvmfr_com=_ctvmfr;
1505 // Do calculation
1506    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
1507    drerr=dl[0];
1508    return dr;
1509 }
1510 
1511 /*=======================================================================
1512                          VertexFit::get_dz
1513   =======================================================================*/
1514 
1515 float VertexFit::get_dz(vertexNumber nv, vertexNumber mv,
1516                                              float& dzerr) const
1517 {
1518    dzerr=-999.;
1519    if(stat!=0) return -999.;
1520 // Get the longitudinal distance between vertices nv and mv
1521 // and return it as the function value.  dzerr is the
1522 // uncertainty on the longitudinal distance, calculated from the
1523 // full error matrix including correlations.
1524    float dxyz[3];
1525    float dr;
1526    float dz;
1527    float dl[3];
1528 
1529    int mvert=mv;
1530    int nvert=nv;
1531 // Copy information into CTVMFT common blocks
1532    *_ctvmq_com=_ctvmq;
1533    *_ctvmfr_com=_ctvmfr;
1534 // Do calculation
1535    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
1536    dzerr=dl[1];
1537    return dz;
1538 }
1539 
1540 /*=======================================================================
1541                          VertexFit::getVertex
1542   =======================================================================*/
1543 
1544 Hep3Vector VertexFit::getVertex(vertexNumber nv) const
1545 {
1546    if(stat!=0) return -999.;
1547 // Return x,y,z position of vertex nv.
1548    if(nv<0 || nv>_ctvmq.nvertx) return Hep3Vector(0,0,0);
1549 // Check if first vertex (nv) is primary vertex.  If it is,
1550 // check if it was used in the primary vertex.  If not,
1551 // all of the corresponding error matrix elements are those
1552 // supplied for the primary vertex.
1553    int nvf=0;
1554    if(nv==0) {
1555       nvf=-1;
1556       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
1557          if(_ctvmq.vtxpnt[jv][0]==0) nvf=0;
1558       }
1559    }
1560    Hep3Vector vertex;
1561 // If primary vertex was not part of fit, take vertex
1562 // location as supplied.
1563    if(nvf < 0) {
1564       vertex.setX(_ctvmq.xyzpv0[0]);
1565       vertex.setY(_ctvmq.xyzpv0[1]);
1566       vertex.setZ(_ctvmq.xyzpv0[2]);
1567    }
1568    else{
1569       vertex.setX(_ctvmq.xyzvrt[nv][0]);
1570       vertex.setY(_ctvmq.xyzvrt[nv][1]);
1571       vertex.setZ(_ctvmq.xyzvrt[nv][2]);
1572    }
1573    return vertex;
1574 }
1575 
1576 /*=======================================================================
1577                          VertexFit::getErrorMatrix
1578   =======================================================================*/
1579 
1580 double VertexFit::getErrorMatrix(int j, int k) const
1581 {
1582    if(stat!=0) return -999.;
1583 // Return the j,k element of the full error matrix VMAT.
1584 // Since the CTVMFT documentation assumes the indices
1585 // start from 1 (ala Fortran), we will also assume this
1586 // and convert the C++ indices.  Note also the that order of
1587 // Fortran and C++ indices is different.  We assume that
1588 // j and k are in the Fortran order.
1589    if(j<1 || k<1 || j>_maxdim+1 || k>_maxdim) {
1590       return -999.;
1591    }
1592    return _ctvmfr.vmat[k-1][j-1];
1593 }
1594 
1595 HepSymMatrix VertexFit::getErrorMatrix(VertexFit::vertexNumber nv) const 
1596 {  
1597   HepSymMatrix cov(3,0);
1598   if(stat!=0) return cov;
1599   // Return x,y,z position of vertex nv.
1600   if(nv<VERTEX_1 || nv>_ctvmq.nvertx) return cov;
1601   // get offset for vertex nv
1602   int voff = _ctvmq.voff[nv-1];
1603   // fill matrix
1604   for(int i = 0 ; i < 3 ; ++i)
1605     for(int j = i ; j < 3 ; ++j)
1606       cov[i][j] = _ctvmfr.vmat[voff+i][voff+j];
1607   return cov;
1608 }
1609 
1610 HepSymMatrix VertexFit::getErrorMatrix(const CdfTrack* trk) const 
1611 {
1612   HepSymMatrix cov(3,0);
1613   if (stat != 0) return cov;
1614 
1615   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
1616 
1617     // Find which track matches this Id
1618     if (trk->id().value() == _ctvmq.list[nt]) {
1619 
1620       // Position of track nt get offset for track nt
1621       int toff = _ctvmq.toff[nt];
1622 
1623       // Fill matrix -- Crv,Phi,Ctg
1624       for(int i = 0 ; i < 3 ; ++i)
1625         for(int j = i ; j < 3 ; ++j)
1626           cov[i][j] = _ctvmfr.vmat[toff+i][toff+j];
1627     }
1628   }
1629   return cov;
1630 }
1631 
1632 float VertexFit::getPtError(const CdfTrack* trk) const
1633 {  
1634   if (stat != 0) return 0;
1635   
1636   int   toff;
1637   float pt,curv,curvErr,ptErr;
1638 
1639   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
1640 
1641     // Find which track matches this Id
1642     if (trk->id().value() == _ctvmq.list[nt]) {
1643 
1644       // Position of track nt get offset for track nt
1645       toff = _ctvmq.toff[nt];
1646 
1647       // Curvature error
1648       pt      = sqrt(_ctvmq.trkp4[0][nt]*_ctvmq.trkp4[0][nt] +
1649                      _ctvmq.trkp4[1][nt]*_ctvmq.trkp4[1][nt]);
1650       curv    = _ctvmq.pscale/pt;
1651       curvErr = sqrt(_ctvmfr.vmat[toff+0][toff+0]);
1652       ptErr   = _ctvmq.pscale/curv/curv*curvErr;
1653       return ptErr;
1654     }
1655   }
1656 
1657   return 0;
1658 }
1659 
1660 /*=======================================================================
1661                          VertexFit::getPosMomErr
1662   =======================================================================*/
1663 void VertexFit::getPosMomErr(HepMatrix& errors) const
1664 {
1665   // A c++ rewrite of the FORTRAN MASSAG function
1666   // The result of this function is an error matrix in position-momentum basis.
1667   // A 7x7 matrix of errors where the rows/columns are x, y, z, px, py, pz, e.
1668 
1669   double cosph[_maxtrk];
1670   double sinph[_maxtrk];
1671   double cosdph[_maxtrk];
1672   double sindph[_maxtrk];
1673   Hep3Vector mom3[_maxtrk];
1674   HepLorentzVector pmom[_maxtrk];
1675   HepLorentzVector total_mom;
1676 
1677   for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
1678     for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
1679       if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
1680       cosph[ltrk] = cos(_ctvmq.par[ltrk][1]);
1681       sinph[ltrk] = sin(_ctvmq.par[ltrk][1]);
1682       double dphi = 0;
1683       sindph[ltrk] = 2 * _ctvmq.par[ltrk][0] * 
1684        (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
1685         _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk]);
1686       cosdph[ltrk] = sqrt(1.0 - sindph[ltrk] * sindph[ltrk]);
1687       if(fabs(sindph[ltrk]) <= 1.0){
1688         dphi = asin(sindph[ltrk]);
1689       }
1690       double pt = _ctvmq.pscale * fabs(1./_ctvmq.par[ltrk][0]);
1691       mom3[ltrk].setX(pt * cos(_ctvmq.par[ltrk][1] + dphi));
1692       mom3[ltrk].setY(pt * sin(_ctvmq.par[ltrk][1] + dphi));
1693       mom3[ltrk].setZ(pt * _ctvmq.par[ltrk][2]);
1694       double e  = sqrt(_ctvmq.tmass[ltrk] * _ctvmq.tmass[ltrk] 
1695                        + mom3[ltrk].mag2());
1696       pmom[ltrk].setVect(mom3[ltrk]);
1697       pmom[ltrk].setE(e);
1698 
1699       total_mom += pmom[ltrk];
1700     }
1701   }
1702 
1703 // Easy so far, but now it gets ugly.
1704 // Fill dp_dpar with the derivatives of the position and momentum with respect
1705 // to the parameters.
1706 
1707   int ctvmft_dim = 3 * (_ctvmq.nvertx + _ctvmq.ntrack);
1708   HepMatrix deriv(ctvmft_dim, 7, 0);
1709   HepMatrix dp_dpar[_maxvtx] = {deriv, deriv, deriv};
1710 
1711 // Fill the x, y, z rows: 
1712 
1713   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
1714     for(int lcomp = 0; lcomp < 3; lcomp++){
1715       dp_dpar[nvtx][(3 * nvtx) + lcomp][lcomp] = 1.0;
1716     }
1717 
1718 // Fill the px, py, pz, e rows:
1719 
1720     for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
1721       for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
1722         if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
1723 
1724 // Find the derivatives of dphi with respect to x, y, curvature, and phi0:
1725 
1726         double dphi_dx = 2.0 * _ctvmq.par[ltrk][0] * cosph[ltrk]/cosdph[ltrk];
1727         double dphi_dy = 2.0 * _ctvmq.par[ltrk][0] * sinph[ltrk]/cosdph[ltrk];
1728         double dphi_dc = 2.0 * 
1729           (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
1730            _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk])/cosdph[ltrk];
1731         double dphi_dp = 2.0 * _ctvmq.par[ltrk][0] *
1732           (-_ctvmq.xyzvrt[lvtx + 1][0] * sinph[ltrk] + 
1733             _ctvmq.xyzvrt[lvtx + 1][1] * cosph[ltrk])/cosdph[ltrk];
1734 
1735 // Now load the derivative matrix
1736 
1737         int lvele = 3 * lvtx;
1738 // dPx/dx:
1739         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dx;
1740 // dPy/dx:
1741         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dx;
1742 // dPz/dx:
1743         dp_dpar[nvtx][lvele][5]  = 0.; 
1744 // dPy/dx:
1745         dp_dpar[nvtx][lvele][6]  = 0.; 
1746 
1747         lvele++;
1748 // dPx/dy:
1749         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dy;
1750 // dPy/dy:
1751         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dy;
1752 // dPz/dy:
1753         dp_dpar[nvtx][lvele][5]  = 0.; 
1754 // dE/dy:
1755         dp_dpar[nvtx][lvele][6]  = 0.; 
1756 
1757         lvele++;
1758 // dPx/dz:
1759         dp_dpar[nvtx][lvele][3]  = 0.;
1760 // dPy/dz:
1761         dp_dpar[nvtx][lvele][4]  = 0.;
1762 // dPz/dz:
1763         dp_dpar[nvtx][lvele][5]  = 0.; 
1764 // dE/dz:
1765         dp_dpar[nvtx][lvele][6]  = 0.; 
1766 
1767         lvele = 3 * (ltrk + _ctvmq.nvertx);
1768 // dPx/dcurv[ltrk]:
1769         dp_dpar[nvtx][lvele][3]  = -(pmom[ltrk].x()/_ctvmq.par[ltrk][0])
1770                                    - pmom[ltrk].y() * dphi_dc;
1771 // dPy/dcurv[ltrk]:
1772         dp_dpar[nvtx][lvele][4]  = -(pmom[ltrk].y()/_ctvmq.par[ltrk][0])
1773                                    + pmom[ltrk].x() * dphi_dc;
1774 // dPz/dcurv[ltrk]:
1775         dp_dpar[nvtx][lvele][5]  = -(pmom[ltrk].z()/_ctvmq.par[ltrk][0]); 
1776 // dE/dcurv[ltrk]:
1777         dp_dpar[nvtx][lvele][6]  = 
1778           -mom3[ltrk].mag2()/(_ctvmq.par[ltrk][0] * pmom[ltrk].e()); 
1779 
1780         lvele++;
1781 // dPx/dphi[ltrk]:
1782         dp_dpar[nvtx][lvele][3]  = -pmom[ltrk].y() * (1.0 + dphi_dp);
1783 // dPy/dphi[ltrk]:
1784         dp_dpar[nvtx][lvele][4]  =  pmom[ltrk].x() * (1.0 + dphi_dp);
1785 // dPz/dphi[ltrk]:
1786         dp_dpar[nvtx][lvele][5]  = 0.;
1787 // dE/dphi[ltrk]:
1788         dp_dpar[nvtx][lvele][6]  = 0.;
1789 
1790         lvele++;
1791 // dPx/dcot[ltrk]:
1792         dp_dpar[nvtx][lvele][3]  = 0;
1793 // dPy/dcot[ltrk]:
1794         dp_dpar[nvtx][lvele][4]  = 0;
1795 // dPz/dcot[ltrk]:
1796         dp_dpar[nvtx][lvele][5]  = pmom[ltrk].perp();
1797 // dE/dcot[ltrk]:
1798         dp_dpar[nvtx][lvele][6]  = 
1799          pmom[ltrk].perp2() * _ctvmq.par[ltrk][2] / pmom[ltrk].e();
1800       }
1801     }
1802   }
1803 // Now calculate the new error matrix
1804 
1805   // Extract the interesting bits from VMAT.
1806 
1807   HepMatrix vmat(ctvmft_dim,ctvmft_dim,0);
1808   for(int lpar = 0; lpar < ctvmft_dim; lpar++){
1809     int l = lpar%3;
1810     int lvele = 0;
1811     if(lpar < 3  * _ctvmq.nvertx){
1812       int lvtx = lpar/3;
1813       lvele = _ctvmq.voff[lvtx] + l;
1814     }
1815     else{
1816       int ltrk = (lpar - 3 * _ctvmq.nvertx)/3;
1817       lvele = _ctvmq.toff[ltrk] + l;
1818     }
1819     for(int kpar = 0; kpar < ctvmft_dim; kpar++){ 
1820       int k = kpar%3;
1821       int kvele = 0;
1822       if(kpar < 3  * _ctvmq.nvertx){
1823         int kvtx = kpar/3;
1824         kvele = _ctvmq.voff[kvtx] + k;
1825       }
1826       else{
1827         int ktrk = (kpar - 3 * _ctvmq.nvertx)/3;
1828         kvele = _ctvmq.toff[ktrk] + k;
1829       }
1830       vmat[kpar][lpar] = _ctvmfr.vmat[kvele][lvele];
1831     }
1832   }
1833   //  Just about done
1834   HepMatrix ans(7,7,0);
1835   HepMatrix answer[_maxvtx] = {ans, ans, ans};
1836   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
1837     answer[nvtx] = (dp_dpar[nvtx].T() * vmat) * dp_dpar[nvtx];
1838   }
1839   errors = answer[0];
1840 }
1841 
1842 /*=======================================================================
1843                          VertexFit::vOff
1844   =======================================================================*/
1845 int VertexFit::vOff(vertexNumber jv) const
1846 { if(jv < VERTEX_1 || jv > _maxvtx) {return -999;}
1847   else {return _ctvmq.voff[jv-1];}
1848 }
1849 
1850 /*=======================================================================
1851                          VertexFit::tOff
1852   =======================================================================*/
1853 int VertexFit::tOff(const CdfTrack* trk) const
1854 { for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
1855     if(trk->id().value()==_ctvmq.list[kt]) return _ctvmq.toff[kt];
1856   }
1857   return -999;
1858 }
1859 
1860 /*=======================================================================
1861                          VertexFit::pOff
1862   =======================================================================*/
1863 int VertexFit::pOff(int lp) const
1864 { if(lp < 1) {
1865     return -999;
1866   }
1867   else return _ctvmq.poff[lp-1];
1868 }
1869 
1870 /*=======================================================================
1871                          VertexFit::cOff
1872   =======================================================================*/
1873 int VertexFit::cOff(int lc) const
1874 { if(lc < 1) {
1875     return -999;
1876   }
1877   else return _ctvmq.coff[lc-1];
1878 }
1879 
1880 /*=======================================================================
1881                          VertexFit::mOff
1882   =======================================================================*/
1883 int VertexFit::mOff() const
1884 { return _ctvmq.moff; }
1885 
1886 /*=======================================================================
1887                          VertexFit::VMat
1888   =======================================================================*/
1889 double VertexFit::VMat(int i, int j) const
1890 { if(i <0 || j < 0) {return -999.;}
1891   else return _ctvmfr.vmat[i][j];
1892 }
1893 
1894 /*=======================================================================
1895                          VertexFit::allocateVertexNumber
1896   =======================================================================*/
1897 
1898 // Facilitates CandNode recursion.  CandNodes have no way of deciding
1899 // which vertex they are, and these trivial functions help them do that.
1900 VertexFit::vertexNumber VertexFit::allocateVertexNumber()
1901 {
1902   if ((_currentAllocatedVertexNumber < PRIMARY_VERTEX) ||
1903       (_currentAllocatedVertexNumber > _maxvtx)) {
1904     std::cout << "VertexFit::allocateVertexNumber: out of range!" << std::endl;
1905     return PRIMARY_VERTEX;
1906   }
1907   return vertexNumber(++_currentAllocatedVertexNumber);
1908 }
1909 
1910 void VertexFit::resetAllocatedVertexNumber()
1911 {
1912   _currentAllocatedVertexNumber = 0;
1913 }
1914 
1915 /*=======================================================================
1916                          VertexFit::restoreFromCommons()
1917   =======================================================================*/
1918 
1919 void VertexFit::restoreFromCommons() {
1920         stat=0;
1921         _ctvmq_com  = (CTVMQ*) ctvmq_address_();
1922         _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
1923         _fiddle_com = (FIDDLE*) fiddle_address_();
1924         _trkprm_com = (TRKPRM*) trkprm_address_();
1925         _ctvmq  = *_ctvmq_com;
1926         _ctvmfr = *_ctvmfr_com;
1927         _fiddle = *_fiddle_com;
1928         _trkprm = *_trkprm_com;
1929 }
1930 
1931 
1932 #endif  /* VERTEXFIT_CC */
1933
[ source navigation ] [ diff markup ] [ identifier search ] [ freetext search ] [ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Valid HTML 4.01!
Send problems or questions to cdfcode@fnal.gov