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 "TrackingUserHL/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     }
660   }
661 // If there was a problem with any track, fail the fit.

662   if(!mstat) {
663     stat = 1;
664     return false;
665   }
666 
667 // First copy information into CTVMFT common blocks

668    *_ctvmq_com=_ctvmq;
669    *_ctvmfr_com=_ctvmfr;
670    *_fiddle_com=_fiddle;
671    *_trkprm_com=_trkprm;
672 // Do the vertex fit.

673    int print=0;
674    int level=0;
675    ctvmft_(print,level,stat);
676 // Now copy information from CTVMFT common blocks to local storage

677    _ctvmq=*_ctvmq_com;
678    _ctvmfr=*_ctvmfr_com;
679    _fiddle=*_fiddle_com;
680    _trkprm=*_trkprm_com;
681    
682    return stat==0;
683 }
684 
685 /*=======================================================================

686                          VertexFit::print

687   =======================================================================*/
688 
689 void VertexFit::print() const
690 {
691    print(std::cout);
692 }
693 
694 void VertexFit::print(std::ostream& os) const
695 {
696    os << "****************************** VertexFit "
697       << "******************************" << std::endl;
698    os << "Number of tracks: " << _ctvmq.ntrack << std::endl;
699    os << "   Tracks: ";
700    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
701       if(jt != 0) os << ",  ";
702       CdfTrack::Id trkId(_ctvmq.list[jt]);
703       os << trkId;
704    }
705    os << std::endl;
706    os << "Number of vertices: " << _ctvmq.nvertx << std::endl;
707    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
708       os << "   Vertex " << jv+1 << " tracks: ";
709       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
710          if(_ctvmq.trkvtx[jv][jt]) {
711             os << " " << _ctvmq.list[jt];
712          }
713       }
714       os << std::endl;
715    }
716    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
717       if(_ctvmq.vtxpnt[0][jv]==0) {
718          os << "   Vertex " << jv+1 << " points to the primary vertex ";
719       }
720       else if(_ctvmq.vtxpnt[0][jv]>0) {
721          os << "   Vertex " << jv+1 << " points to vertex "
722               << _ctvmq.vtxpnt[0][jv];
723       }
724       if(_ctvmq.vtxpnt[1][jv]==1) {
725          os << " in 2 dimensions" << std::endl;
726       }
727       else if(_ctvmq.vtxpnt[1][jv]==2) {
728          os << " in 3 dimensions" << std::endl;
729       }
730       else if(_ctvmq.vtxpnt[1][jv]==3) {
731          os << ", a single track vertex" << std::endl;
732       }
733       if(_ctvmq.cvtx[jv]>0) {
734          os << "   Vertex " << jv+1 << " is a conversion" << std::endl;
735       }
736    }
737    os << "Number of mass constraints: " << _ctvmq.nmassc << std::endl;
738    for(int jmc=0; jmc<_ctvmq.nmassc; ++jmc) {
739       os << "   Tracks ";
740       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
741          if(_ctvmq.trkmcn[jmc][jt]) {
742             os << " " << _ctvmq.list[jt];
743          }
744       }
745       os << " constrained to mass " << _ctvmq.cmass[jmc]
746          << " Gev/c^2" << std::endl;
747    }
748    if(stat==-999) {
749       os << "No fit has been done." << std::endl;
750    }
751    else {
752       os << "***** Results of Fit *****" << std::endl;
753       printErr(os);
754       os << "   Status = " << stat << std::endl;
755       os << "   Chi-square = " << _ctvmq.chisqr[0]
756          << " for " << _ctvmq.ndof << " degrees of freedom." << std::endl;
757       os << "   => probability = " << prob() << std::endl;
758       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
759          os << "Vertex " << jv+1 << " position: "
760             << _ctvmq.xyzvrt[jv+1][0] << " "
761             << _ctvmq.xyzvrt[jv+1][1] << " "
762             << _ctvmq.xyzvrt[jv+1][2] << std::endl;
763       }
764       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
765          os << "Track " << _ctvmq.list[jt] << " P4: " 
766             << _ctvmq.trkp4[0][jt] << " "
767             << _ctvmq.trkp4[1][jt] << " "
768             << _ctvmq.trkp4[2][jt] << " "
769             << _ctvmq.trkp4[3][jt] << " " << std::endl;
770       }
771    }
772    os << "****************************************"
773       << "**************************" << std::endl;
774 
775    return;
776 }
777 
778 /*=======================================================================

779                          VertexFit::printErr

780   =======================================================================*/
781 
782 void VertexFit::printErr() const
783 {
784    printErr(std::cout);
785 }
786 
787 void VertexFit::printErr(std::ostream& os) const
788 {
789    os << "VertexFit: IJKERR = " << _ctvmq.ijkerr[0] << ", "
790                         << _ctvmq.ijkerr[1] << ", "
791                         << _ctvmq.ijkerr[2] << std::endl;
792    if(status()==0 && _ctvmq.ijkerr[0]==0) return;
793    if(_ctvmq.ijkerr[0] == -1) {
794       os << "   Problem with GETTRK:  track requested is not in list."
795          << std::endl
796          << "   This should not happen - Contact VertexFit expert "
797          << _expert << "." <<std::endl;
798    }
799    else if(_ctvmq.ijkerr[0]==1) {
800       os << "   Problem in CTVM00:" << std::endl;
801       if(_ctvmq.ijkerr[1]==1) {
802          os << "      Number of tracks is " << _ctvmq.ntrack
803             << "." << std::endl;
804          if(_ctvmq.ntrack < 2) {
805             os << ", which is too few (must be at least 2)." << std::endl;
806          }
807          else if(_ctvmq.ntrack > _maxtrk) {
808             os << ", which is too many (maximum is " << _maxtrk
809                << ")." << std::endl;
810          }
811          else {
812             os << "      Problem with number of tracks"
813                << " for unknown reasons." << std::endl;
814          }
815       }
816       else if(_ctvmq.ijkerr[1]==2) {
817          os << "      Number of vertices is " << _ctvmq.nvertx
818             << "." << std::endl;
819          if(_ctvmq.nvertx < 1) {
820             os << ", which is too few (must be at least 1)." << std::endl;
821          }
822          else if(_ctvmq.nvertx > _maxvtx) {
823             os << ", which is too many (maximum is " << _maxvtx
824                << ")." << std::endl;
825          }
826          else {
827             os << std::endl << "      Problem with number of vertices"
828                << " for unknown reasons." << std::endl;
829          }
830       }
831       else if(_ctvmq.ijkerr[1]==3) {
832          os << "      Number of mass constraints is " << _ctvmq.nmassc
833             << "." << std::endl;
834          if(_ctvmq.nmassc < 0) {
835             os << ", which is negative." << std::endl;
836          }
837          else if(_ctvmq.nmassc > _maxmcn) {
838             os << ", which is too many (maximum is " << _maxmcn
839                << ")." << std::endl;
840          }
841          else {
842             os << std::endl << "      Problem with number of mass"
843                << " constraints for unknown reasons." << std::endl;
844          }
845       }
846       else if(_ctvmq.ijkerr[1]==11) {
847          os << "      Vertex " << _ctvmq.ijkerr[2]
848             << " has less than one track." << std::endl;
849       }
850       else if(_ctvmq.ijkerr[1]==12) {
851          os << "      Vertex " << _ctvmq.ijkerr[2]
852             << " is a conversion vertex with a number of tracks"
853             << " different than two." << std::endl;
854       }
855       else if(_ctvmq.ijkerr[1]==13) {
856          os << "      Vertex " << _ctvmq.ijkerr[2]
857             << " is a one track vertex that has no multi-track"
858             << " descendents." << std::endl;
859       }
860       else if(_ctvmq.ijkerr[1]==14) {
861          os << "      Vertex " << _ctvmq.ijkerr[2]
862             << " does not point at a vertex with a lower number."
863             << std::endl;
864       }
865       else if(_ctvmq.ijkerr[1]==15) {
866          os << "      Vertex " << _ctvmq.ijkerr[2]
867             << " has a parent vertex that is a conversion." << std::endl;
868       }
869       else if(_ctvmq.ijkerr[1]==16) {
870          os << "      Vertex " << _ctvmq.ijkerr[2]
871             << " does 1 track pointing to a vertex with"
872             << " more than 1 track." << std::endl;
873       }
874       else if(_ctvmq.ijkerr[1]==17) {
875          os << "      Vertex " << _ctvmq.ijkerr[2]
876             << " does 0 track pointing to a vertex with"
877             << " more than 0 track (?)." << std::endl;
878       }
879       else if(_ctvmq.ijkerr[1]==19) {
880          os << "      Primary vertex error matrix is singular." << std::endl;
881       }
882       else if(_ctvmq.ijkerr[1]==21) {
883          os << "      Track with Id " << _ctvmq.ijkerr[2]
884             << "is not in any vertex." << std::endl;
885       }
886       else if(_ctvmq.ijkerr[1]==22) {
887          os << "      Track with Id " << _ctvmq.ijkerr[2]
888             << "is in multiple vertices." << std::endl;
889       }
890       else if(_ctvmq.ijkerr[1]==23) {
891          os << "      Track with Id " << _ctvmq.ijkerr[2]
892             << "occurs more than once." << std::endl;
893       }
894       else if(_ctvmq.ijkerr[1]==31) {
895          os << "      A mass constraint has less than 2 tracks." << std::endl;
896       }
897       else if(_ctvmq.ijkerr[1]==32) {
898          os << "      The sum masses of the tracks in a mass constraint"
899             << " exceeds the constraint mass." << std::endl;
900       }
901       else if(_ctvmq.ijkerr[1]==33) {
902          os << "      Beamline constraint. Beam covariance not set properly."
903             << " Negative diagonal elements." << std::endl;
904       }
905       else if(_ctvmq.ijkerr[1]==34) {
906          os << "      Beamline constraint. Beam covariance not set properly."
907             << " Off-diagonal elements not zero." << std::endl;
908       }
909       else if(_ctvmq.ijkerr[1]==36) {
910          os << "      Beamline constraint. Number of vertices = " 
911             << _ctvmq.nvertx << " Should be 1." << std::endl;
912       }
913    }
914    else if(_ctvmq.ijkerr[0] == 2) {
915       if(_ctvmq.ijkerr[1] == 20) {
916          os << "   Problem in CTVM00: " << std::endl;
917          os << "      Track has negative Id = "
918             << _ctvmq.list[_ctvmq.ijkerr[2]-1] << "." << std::endl;
919       }
920       else {
921          os << "   Problem in CTVMFA with vertex "
922             << _ctvmq.ijkerr[2] << ": " << std::endl;
923          os << "      Failure in vertex first approximation." << std::endl;
924          if(_ctvmq.ijkerr[1] == 1) {
925             os << "      Tracks are concentric circles." << std::endl;
926          }
927          if(_ctvmq.ijkerr[1] == 2) {
928             os << "      Conversion vertex has widely separated"
929                << " exterior circles at midpoint." << std::endl;
930          }
931          if(_ctvmq.ijkerr[1] == 3) {
932             os << "      Conversion vertex has widely separated"
933                << " interior circles at midpoint." << std::endl;
934          }
935          if(_ctvmq.ijkerr[1] == 4) {
936             os << "      Vertex has widely separated"
937                << " exterior circles at approximate vertex." << std::endl;
938          }
939          if(_ctvmq.ijkerr[1] == 5) {
940             os << "      Vertex has widely separated"
941                << " interior circles at approximate vertex." << std::endl;
942          }
943          if(_ctvmq.ijkerr[1] == 6) {
944             os << "      Rv is too large at the chosen"
945                << " intersection point." << std::endl;
946          }
947          if(_ctvmq.ijkerr[1] == 7) {
948             os << "      Delta z is too large at the chosen"
949                << " intersection point." << std::endl;
950          }
951          if(_ctvmq.ijkerr[1] == 8) {
952             os << "      A track's turning to the chosen vertex"
953                << " is too large." << std::endl;
954          }
955          if(_ctvmq.ijkerr[1] == 9) {
956             os << "      There is no solution with an adequately"
957                << " positive arc length." << std::endl;
958          }
959          if(_ctvmq.ijkerr[1] == 21) {
960             os << "      zero-track vertexing: either/both vertex "
961                << " momenta are too small (<0.01 MeV)." << std::endl;
962          }
963          if(_ctvmq.ijkerr[1] == 22) {
964             os << "      zero-track vertexing: Two lines (tracks) are "
965                << " parallel/antiparallel." << std::endl;
966          }
967 
968       }
969    }
970    else if(_ctvmq.ijkerr[0] == 3) {
971       os << "   Problem in CTVM01 with track with Id = "
972          << _ctvmq.list[_ctvmq.ijkerr[2]-1] << ": " << std::endl;
973       if(_ctvmq.ijkerr[1] == 1) {
974          os << "      GETTRK cannot find Id in list." << std::endl;
975       }
976       if(_ctvmq.ijkerr[1] == 2) {
977          os << "      Covariance matrix could not be inverted."  << endl 
978             << "      Offending track number (in order addded) is "
979             << _ctvmq.ijkerr[2] << "." << std::endl;
980       }
981       if(_ctvmq.ijkerr[1] == 3) {
982          os << "      Track turns through too large an angle"
983             << " to the vertex." << std::endl;
984       }
985       if(_ctvmq.ijkerr[1] == 4) {
986          os << "      Track moves too far backward to vertex." << std::endl;
987       }
988       if(_ctvmq.ijkerr[1] == 5) {
989          os << "      Track with curvature > 0.01." << std::endl
990             << "      Offending track number is "
991             << _ctvmq.ijkerr[2] << "." << std::endl;
992       }
993    }
994    else if(status() == 9) {
995       os << "   General fit problem: " << std::endl;
996       if(_ctvmq.ijkerr[1] == 1) {
997          os << "      Singular solution matrix." << std::endl;
998       }
999       if(_ctvmq.ijkerr[1] == 2 || _ctvmq.ijkerr[1] == 3) {
1000          os << "      Too many iterations ( "
1001             << _ctvmq.ijkerr[2] << "(." << std::endl;
1002       }
1003       if(_ctvmq.ijkerr[1] == 4) {
1004          os << "      Convergence failure." << std::endl;
1005       }
1006       if(_ctvmq.ijkerr[1] == 5) {
1007          os << "      Bad convergence." << std::endl;
1008       }
1009       if(_ctvmq.ijkerr[1] == 9) {
1010          os << "      Ill-formed  covariance matrix." << std::endl;
1011       }
1012    }
1013    else {
1014       os << "   The error codes above are not recognized." << std::endl
1015          << "   Contact VertexFit expert " << _expert << "." << std::endl;
1016    }
1017    return;
1018 }
1019 
1020 
1021 /*=======================================================================

1022                          VertexFit::getIJKErr

1023   =======================================================================*/
1024 
1025 void VertexFit::getIJKErr(int& err0, int& err1, int& err2) const
1026 {
1027   err0 = _ctvmq.ijkerr[0];
1028   err1 = _ctvmq.ijkerr[1];
1029   err2 = _ctvmq.ijkerr[2];
1030   return;
1031 }
1032 
1033 int VertexFit::getIJKErr0() const
1034 {
1035   return _ctvmq.ijkerr[0];
1036 }
1037 int VertexFit::getIJKErr1() const
1038 {
1039   return _ctvmq.ijkerr[1];
1040 }
1041 int VertexFit::getIJKErr2() const
1042 {
1043   return _ctvmq.ijkerr[2];
1044 }
1045 
1046 /*=======================================================================

1047                        VertexFit::getErrTrackId()

1048   =======================================================================*/
1049 
1050 int VertexFit::getErrTrackId() const
1051  {
1052    if (status() == 0) return 0;
1053    int trkId = 0;
1054    // Problems with track in CTVM01 or track has negative id in CTVM00

1055    // See PrintErr() for a more detailed list of error codes.

1056    if ((_ctvmq.ijkerr[0] == 2 && _ctvmq.ijkerr[1] == 20) ||
1057        _ctvmq.ijkerr[0] == 3) {
1058      trkId = _ctvmq.list[_ctvmq.ijkerr[2]-1];
1059    }
1060    
1061    return trkId;
1062  }
1063 
1064 /*=======================================================================

1065                          VertexFit::expert

1066   =======================================================================*/
1067 
1068 string VertexFit::expert() const
1069 {
1070    return _expert;
1071 }
1072 
1073 /*=======================================================================

1074                          VertexFit::status

1075   =======================================================================*/
1076 
1077 int VertexFit::status() const
1078 {
1079 
1080    return stat;
1081 }
1082 
1083 /*=======================================================================

1084                          VertexFit::chisq

1085   =======================================================================*/
1086 
1087 float VertexFit::chisq() const
1088 {
1089 // Chi-square of fit

1090       return _ctvmq.chisqr[0];
1091 }
1092 
1093 /*=======================================================================

1094                          VertexFit::ndof

1095   =======================================================================*/
1096 
1097 int VertexFit::ndof() const
1098 {
1099 // Number of degrees of freedom of fit.

1100    if(_ctvmq.chisqr[0] >= 0) {
1101       return _ctvmq.ndof; }
1102    else {
1103       return 0;}
1104 }
1105 
1106 /*=======================================================================

1107                          VertexFit::prob

1108   =======================================================================*/
1109 
1110 float VertexFit::prob() const
1111 {
1112 // Probability of chi-square of fit

1113    if(_ctvmq.chisqr[0]>=0.) {
1114       float chisq=_ctvmq.chisqr[0];
1115       int nd=_ctvmq.ndof;
1116       return prob_(chisq,nd);
1117    }
1118    else {
1119       return -999.;
1120    }
1121 }
1122 
1123 /*=======================================================================

1124                          VertexFit::chisq(CdfTrack*)

1125   =======================================================================*/
1126 
1127 float VertexFit::chisq(const CdfTrack* trk) const
1128 {
1129 // This method returns the chisquare contribution for one track 

1130 // If fit not successful or not done yet, return -1.

1131    if(_ctvmq.chisqr[0] < 0) return -1.;
1132 // Look for this track in the list of tracks.

1133    for(int jt = 0; jt < _ctvmq.ntrack; ++jt) {
1134      if(trk->id().value() == _ctvmq.list[jt]) {
1135 // Found the track, so return its chisquare contribution.

1136        return _ctvmq.chit[jt];
1137      }
1138    }
1139 // If track is not in list, return -1.

1140    return -1.;
1141 }
1142 
1143 /*=======================================================================

1144                          VertexFit::chisq_rphi()

1145   =======================================================================*/
1146 
1147 float VertexFit::chisq_rphi() const
1148 {
1149 // This method returns the chisquare contribution in the r-phi plane.

1150    int index[3] = {0,1,3};
1151 // If fit not successful or not done yet, return -1.

1152    if(_ctvmq.chisqr[0] < 0) return -1.;
1153 // Loop over the tracks in the event.

1154    float chisq = 0.;
1155    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1156 // Double loop over the relevant parameter indices.

1157      for(int k1=0; k1<3; ++k1) {
1158        for(int k2=0; k2<3; ++k2) {
1159 // Add contribution to chisquare.

1160          chisq += _ctvmq.pardif[jt][index[k1]] *
1161                   _ctvmq.g[jt][index[k1]][index[k2]] *
1162                   _ctvmq.pardif[jt][index[k2]];
1163        }
1164      }
1165    }
1166 // Return the chisquare.

1167    return chisq;
1168 }
1169 
1170 
1171 /*=======================================================================

1172                          VertexFit::chisq_z()

1173   =======================================================================*/
1174 
1175 float VertexFit::chisq_z() const
1176 {
1177 // This method returns the chisquare contribution in the z direction.

1178    int index[2] = {2,4};
1179 // If fit not successful or not done yet, return -1.

1180    if(_ctvmq.chisqr[0] < 0) return -1.;
1181 // Loop over the tracks in the event.

1182    float chisq = 0.;
1183    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1184 // Double loop over the relevant parameter indices.

1185      for(int k1=0; k1<2; ++k1) {
1186        for(int k2=0; k2<2; ++k2) {
1187 // Add contribution to chisquare.

1188          chisq += _ctvmq.pardif[jt][index[k1]] *
1189                   _ctvmq.g[jt][index[k1]][index[k2]] *
1190                   _ctvmq.pardif[jt][index[k2]];
1191        }
1192      }
1193    }
1194 // Return the chisquare.

1195    return chisq;
1196 }
1197 
1198 
1199 /*=======================================================================

1200                          VertexFit::chisq_rphiz()

1201   =======================================================================*/
1202 
1203 float VertexFit::chisq_rphiz() const
1204 {
1205 // This method returns the chisquare contribution of the cross

1206 // terms in the r-phi and z directions.

1207    int index1[2] = {2,4};
1208    int index2[3] = {0,1,3};
1209 // If fit not successful or not done yet, return -1.

1210    if(_ctvmq.chisqr[0] < 0) return -1.;
1211 // Loop over the tracks in the event.

1212    float chisq = 0.;
1213    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1214 // Double loop over the relevant parameter indices.

1215      for(int k1=0; k1<2; ++k1) {
1216        for(int k2=0; k2<3; ++k2) {
1217 // Add contribution to chisquare.

1218          chisq += _ctvmq.pardif[jt][index1[k1]] *
1219                   _ctvmq.g[jt][index1[k1]][index2[k2]] *
1220                   _ctvmq.pardif[jt][index2[k2]];
1221        }
1222      }
1223    }
1224 // Return the chisquare.

1225    return 2.*chisq;
1226 }
1227 
1228 /*=======================================================================

1229                          VertexFit::chisq_rphi(CdfTrack*)

1230   =======================================================================*/
1231 
1232 float VertexFit::chisq_rphi(const CdfTrack* trk) const
1233 {
1234 // This method returns the chisquare contribution in the r-phi plane.

1235    int index[3] = {0,1,3};
1236 // If fit not successful or not done yet, return -1.

1237    if(_ctvmq.chisqr[0] < 0) return -1.;
1238 // Loop over the tracks in the event, looking for the one we want

1239    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1240      if(trk->id().value() == _ctvmq.list[jt]) {
1241 // Found the track, so calculate its chisquare contribution.

1242        float chisq = 0.;
1243 // Double loop over the relevant parameter indices.

1244        for(int k1=0; k1<3; ++k1) {
1245          for(int k2=0; k2<3; ++k2) {
1246 // Add contribution to chisquare.

1247            chisq += _ctvmq.pardif[jt][index[k1]] *
1248                     _ctvmq.g[jt][index[k1]][index[k2]] *
1249                     _ctvmq.pardif[jt][index[k2]];
1250          }
1251        }
1252        return chisq;
1253      }
1254    }
1255 // Track not found, return -1.

1256    return -1.;
1257 }
1258 
1259 
1260 /*=======================================================================

1261                          VertexFit::chisq_z(CdfTrack*)

1262   =======================================================================*/
1263 
1264 float VertexFit::chisq_z(const CdfTrack* trk) const
1265 {
1266 // This method returns the chisquare contribution in the z direction.

1267    int index[2] = {2,4};
1268 // If fit not successful or not done yet, return -1.

1269    if(_ctvmq.chisqr[0] < 0) return -1.;
1270 // Loop over the tracks in the event, looking for the one we want.

1271    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1272      if(trk->id().value() == _ctvmq.list[jt]) {
1273 // Found the track, so calculate its chisquare contribution.

1274        float chisq = 0.;
1275 // Double loop over the relevant parameter indices.

1276        for(int k1=0; k1<2; ++k1) {
1277          for(int k2=0; k2<2; ++k2) {
1278 // Add contribution to chisquare.

1279            chisq += _ctvmq.pardif[jt][index[k1]] *
1280                     _ctvmq.g[jt][index[k1]][index[k2]] *
1281                     _ctvmq.pardif[jt][index[k2]];
1282          }
1283        }
1284        return chisq;
1285      }
1286    }
1287 // Track not found - return -1.

1288    return -1.;
1289 }
1290 
1291 
1292 /*=======================================================================

1293                          VertexFit::chisq_rphiz(CdfTrack*)

1294   =======================================================================*/
1295 
1296 float VertexFit::chisq_rphiz(const CdfTrack* trk) const
1297 {
1298 // This method returns the chisquare contribution of the cross

1299 // terms in the r-phi and z directions.

1300    int index1[2] = {2,4};
1301    int index2[3] = {0,1,3};
1302 // If fit not successful or not done yet, return -1.

1303    if(_ctvmq.chisqr[0] < 0) return -1.;
1304 // Loop over the tracks in the event.

1305    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1306      if(trk->id().value() == _ctvmq.list[jt]) {
1307 // Found the track, so calculate its chisquare contribution.

1308        float chisq = 0.;
1309 // Double loop over the relevant parameter indices.

1310        for(int k1=0; k1<2; ++k1) {
1311          for(int k2=0; k2<3; ++k2) {
1312 // Add contribution to chisquare.

1313            chisq += _ctvmq.pardif[jt][index1[k1]] *
1314                     _ctvmq.g[jt][index1[k1]][index2[k2]] *
1315                     _ctvmq.pardif[jt][index2[k2]];
1316          }
1317        }
1318        return 2.*chisq;
1319      }
1320    }
1321 // Track not found so return -1.

1322    return -1.;
1323 }
1324 
1325 /*=======================================================================

1326                          VertexFit::getTrackP4

1327   =======================================================================*/
1328 
1329 HepLorentzVector VertexFit::getTrackP4(const CdfTrack* trk) const
1330 {
1331    if(stat != 0) return HepLorentzVector(0,0,0,0);
1332 // return four momentum of fit track 

1333    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1334 //   Find which track matches this Id.

1335       if(trk->id().value()==_ctvmq.list[jt]) {
1336          HepLorentzVector p4((double)_ctvmq.trkp4[0][jt],
1337                              (double)_ctvmq.trkp4[1][jt],
1338                              (double)_ctvmq.trkp4[2][jt],
1339                              (double)_ctvmq.trkp4[3][jt]);
1340          return p4;
1341       }
1342    }
1343    return HepLorentzVector(0,0,0,0);
1344 }
1345 
1346 /*=======================================================================

1347                          VertexFit::getHelix

1348   =======================================================================*/
1349 
1350 Helix VertexFit::getHelix(const CdfTrack* trk) const
1351 {
1352    if(stat != 0) return Helix(0,0,0,0,0);
1353 // return helix parameters of fit track 

1354    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
1355 //   Find which track matches this Id.

1356       if(trk->id().value()==_ctvmq.list[jt]) {
1357          Helix trhelix((double)_ctvmq.par[jt][2],
1358                        (double)_ctvmq.par[jt][0],
1359                        (double)_ctvmq.par[jt][4],
1360                        (double)_ctvmq.par[jt][3],
1361                        (double)_ctvmq.par[jt][1]);
1362          return trhelix;
1363       }
1364    }
1365    return Helix(0,0,0,0,0);
1366 }
1367 
1368 /*=======================================================================

1369                          VertexFit::getMass

1370   =======================================================================*/
1371 
1372 float VertexFit::getMass(const CdfTrack* trk1,
1373                          const CdfTrack* trk2, float& dmass) const
1374 {
1375 
1376 //

1377 // Vertex Fitting is prone to floating point errors.  Here, we

1378 // turn those off during the duration of this routine: Joe and Farrukh

1379 //

1380 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1381   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1382   sigaction(SIGFPE, &myaction, &oldaction);
1383   if (setjmp(env)!=0) {
1384     sigaction(SIGFPE, &oldaction,0);
1385     return -999;
1386   }
1387 #endif
1388 
1389    dmass=-999.;
1390    if(stat!=0) return -999.;
1391 // Get the fit invariant mass of tracks trk1 and trk2.

1392 // dmass is the error on the mass.

1393    int ntrk=2;
1394    const CdfTrack* trks[2];
1395    trks[0]=trk1;
1396    trks[1]=trk2;
1397 
1398   double result = getMass(ntrk,trks,dmass);
1399    //return something for FPE Protection-> Joe and Farrukh

1400 
1401 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1402   sigaction(SIGFPE, &oldaction,0);
1403 #endif
1404 
1405   return result;
1406 }
1407 
1408 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
1409                                 const CdfTrack* trk3,float& dmass) const
1410 {
1411 
1412 //

1413 // Vertex Fitting is prone to floating point errors.  Here, we

1414 // turn those off during the duration of this routine: Joe and Farrukh

1415 //

1416 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1417   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1418   sigaction(SIGFPE, &myaction, &oldaction);
1419   if (setjmp(env)!=0) {
1420     sigaction(SIGFPE, &oldaction,0);
1421     return -999;
1422   }
1423 #endif
1424 
1425 
1426 
1427    dmass=-999.;
1428    if(stat!=0) return -999.;
1429 // Get the fit invariant mass of tracks trk1, trk2, and trk3.

1430 // dmass is the error on the mass.

1431    int ntrk=3;
1432    const CdfTrack* trks[3];
1433    trks[0]=trk1;
1434    trks[1]=trk2;
1435    trks[2]=trk3;
1436 
1437 
1438    double result = getMass(ntrk,trks,dmass);
1439 
1440 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1441   sigaction(SIGFPE, &oldaction,0);
1442 #endif
1443 
1444   return result;
1445 }
1446 
1447 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
1448      const CdfTrack* trk3, const CdfTrack* trk4, float& dmass) const
1449 {
1450    dmass=-999.;
1451    if(stat!=0) return -999.;
1452 // Get the fit invariant mass of tracks trk1, trk2, trk3, and trk4.

1453 // dmass is the error on the mass.

1454    int ntrk=4;
1455    const CdfTrack* trks[4];
1456    trks[0]=trk1;
1457    trks[1]=trk2;
1458    trks[2]=trk3;
1459    trks[3]=trk4;
1460 //

1461 // Vertex Fitting is prone to floating point errors.  Here, we

1462 // turn those off during the duration of this routine: Joe and Farrukh

1463 //

1464 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1465   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1466   sigaction(SIGFPE, &myaction, &oldaction);
1467   if (setjmp(env)!=0) {
1468     sigaction(SIGFPE, &oldaction,0);
1469     return -999;
1470   }
1471 #endif
1472 
1473   double result = getMass(ntrk,trks,dmass);
1474 
1475    //return something for FPE Protection-> Joe and Farrukh

1476 
1477 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1478   sigaction(SIGFPE, &oldaction,0);
1479 #endif
1480 
1481   return result;
1482 }
1483 
1484 float VertexFit::getMass(int ntrk, const CdfTrack* tracks[],
1485                                                 float& dmass) const
1486 {
1487 
1488 // Vertex Fitting is prone to floating point errors.  Here, we

1489 // turn those off during the duration of this routine: Joe and Farrukh

1490 //

1491 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1492   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
1493   sigaction(SIGFPE, &myaction, &oldaction);
1494   if (setjmp(env)!=0) {
1495     sigaction(SIGFPE, &oldaction,0);
1496     return -999;
1497   }
1498 #endif
1499 
1500    dmass=-999.;
1501    if(stat!=0) return -999.;
1502 // Get fit invariant mass of ntrk tracks listed in array tracks.

1503 // dmass is the error on the mass.

1504    dmass=-999.;
1505    if(ntrk <= 0) return 0;
1506    int jtrks[_maxtrk];
1507    for(int jt=0; jt<ntrk; ++jt) {
1508       bool found=false;
1509       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
1510          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
1511             found=true;
1512             jtrks[jt]=kt+1;
1513          }
1514       }
1515       if(!found) return 0;
1516    }
1517 // Copy information into CTVMFT common blocks

1518    *_ctvmq_com=_ctvmq;
1519    *_ctvmfr_com=_ctvmfr;
1520    int ntr=ntrk;
1521    float mass;
1522    double p4[4];
1523    mcalc_(ntr, jtrks, mass, dmass, p4);
1524 
1525    //return something for FPE-> Joe and Farrukh

1526 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
1527   sigaction(SIGFPE, &oldaction,0);
1528 #endif
1529 
1530    return mass;
1531 }
1532 
1533 /*=======================================================================

1534                          VertexFit::getDecayLength

1535   =======================================================================*/
1536 
1537 float VertexFit::getDecayLength(vertexNumber nv, vertexNumber mv, 
1538                                 const Hep3Vector& dir, float& dlerr) const
1539 {
1540    dlerr=-999.;
1541    if(stat!=0) return -999.;
1542 // Get the signed decay length from vertex nv to vertex mv

1543 // along the x-y direction defined by the 3 vector dir.

1544 // dlerr is the error on the decay length including the

1545 // full error matrix.

1546 // Check that vertices are within range.

1547    if(nv<0 || nv>=_ctvmq.nvertx) return -999.;
1548    if(mv<1 || mv>_ctvmq.nvertx) return -999.;
1549    if(nv>=mv) return -999.;
1550    float dir_t=sqrt(dir.x()*dir.x()+dir.y()*dir.y());
1551    if(dir_t <= 0.) return -999.;
1552    Hep3Vector dv=getVertex(mv)-getVertex(nv);
1553    float dl=(dv.x()*dir.x()+dv.y()*dir.y())/dir_t;
1554 // Set up the column matrix of derivatives

1555    HepMatrix A(4,1);
1556    A(1,1)=dir.x()/dir_t;
1557    A(2,1)=dir.y()/dir_t;
1558    A(3,1)=-dir.x()/dir_t;
1559    A(4,1)=-dir.y()/dir_t;
1560 // Check if first vertex (nv) is primary vertex.  If it is,

1561 // check if it was used in the primary vertex.  If not,

1562 // all of the corresponding error matrix elements are those

1563 // supplied for the primary vertex.

1564    int nvf=0;
1565    if(nv==0) {
1566       nvf=-1;
1567       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
1568          if(_ctvmq.vtxpnt[0][jv]==0) nvf=0;
1569       }
1570    }
1571 // Get the relevant submatrix of the full error matrix.

1572    HepMatrix V(4,4,0);
1573    if(nvf < 0) {
1574       V(1,1)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+1);
1575       V(1,2)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2);
1576       V(2,1)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+1);
1577       V(2,2)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+2);
1578       V(3,3)=_ctvmq.exyzpv[0][0];
1579       V(3,4)=_ctvmq.exyzpv[0][1];
1580       V(4,3)=_ctvmq.exyzpv[1][0];
1581       V(4,4)=_ctvmq.exyzpv[1][1];
1582    }
1583    else {
1584 //    Get the indices into the error matrix vmat

1585       int index[4]={_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2,0,0};
1586       if(nv==0) {
1587          index[2]=1;
1588          index[3]=2;
1589       }
1590       else {
1591         index[2] = _ctvmq.voff[nv-1]+1;
1592         index[3] = _ctvmq.voff[nv-1]+2;
1593       }
1594       for(int j=0; j<4; ++j) {
1595          for(int k=0; k<4; ++k) {
1596             V[j][k]=getErrorMatrix(index[j],index[k]);
1597          }
1598       }
1599    }
1600 // Calculate square of dlerr

1601    dlerr=(A.T()*V*A)(1,1);
1602    if(dlerr >= 0.) {
1603       dlerr=sqrt(dlerr);
1604    }
1605    else {
1606       dlerr=-sqrt(-dlerr);
1607    }
1608 
1609    return dl;
1610 }
1611 
1612 /*=======================================================================

1613                          VertexFit::get_dr

1614   =======================================================================*/
1615 
1616 float VertexFit::get_dr(vertexNumber nv, vertexNumber mv,
1617                                            float& drerr) const
1618 {
1619    drerr=-999.;
1620    if(stat!=0) return -999.;
1621 // Get the transvese distance between vertices nv and mv

1622 // and return it as the function value.  drerr is the

1623 // uncertainty on the transverse distance, calculated from the

1624 // full error matrix including correlations.

1625    float dxyz[3];
1626    float dr;
1627    float dz;
1628    float dl[3];
1629 
1630    int mvert=mv;
1631    int nvert=nv;
1632 // Copy information into CTVMFT common blocks

1633    *_ctvmq_com=_ctvmq;
1634    *_ctvmfr_com=_ctvmfr;
1635 // Do calculation

1636    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
1637    drerr=dl[0];
1638    return dr;
1639 }
1640 
1641 /*=======================================================================

1642                          VertexFit::get_dz

1643   =======================================================================*/
1644 
1645 float VertexFit::get_dz(vertexNumber nv, vertexNumber mv,
1646                                              float& dzerr) const
1647 {
1648    dzerr=-999.;
1649    if(stat!=0) return -999.;
1650 // Get the longitudinal distance between vertices nv and mv

1651 // and return it as the function value.  dzerr is the

1652 // uncertainty on the longitudinal distance, calculated from the

1653 // full error matrix including correlations.

1654    float dxyz[3];
1655    float dr;
1656    float dz;
1657    float dl[3];
1658 
1659    int mvert=mv;
1660    int nvert=nv;
1661 // Copy information into CTVMFT common blocks

1662    *_ctvmq_com=_ctvmq;
1663    *_ctvmfr_com=_ctvmfr;
1664 // Do calculation

1665    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
1666    dzerr=dl[1];
1667    return dz;
1668 }
1669 
1670 /*=======================================================================

1671                          VertexFit::getVertex

1672   =======================================================================*/
1673 
1674 Hep3Vector VertexFit::getVertex(vertexNumber nv) const
1675 {
1676    if(stat!=0) return -999.;
1677 // Return x,y,z position of vertex nv.

1678    if(nv<0 || nv>_ctvmq.nvertx) return Hep3Vector(0,0,0);
1679 // Check if first vertex (nv) is primary vertex.  If it is,

1680 // check if it was used in the primary vertex.  If not,

1681 // all of the corresponding error matrix elements are those

1682 // supplied for the primary vertex.

1683    int nvf=0;
1684    if(nv==0) {
1685       nvf=-1;
1686       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
1687          if(_ctvmq.vtxpnt[0][jv]==0) nvf=0;
1688       }
1689    }
1690    Hep3Vector vertex;
1691 // If primary vertex was not part of fit, take vertex

1692 // location as supplied.

1693    if(nvf < 0) {
1694       vertex.setX(_ctvmq.xyzpv0[0]);
1695       vertex.setY(_ctvmq.xyzpv0[1]);
1696       vertex.setZ(_ctvmq.xyzpv0[2]);
1697    }
1698    else{
1699       vertex.setX(_ctvmq.xyzvrt[nv][0]);
1700       vertex.setY(_ctvmq.xyzvrt[nv][1]);
1701       vertex.setZ(_ctvmq.xyzvrt[nv][2]);
1702    }
1703 // If we have a different reference point, need to add it back in

1704    vertex += _referencePoint;
1705    return vertex;
1706 }
1707 
1708 /*=======================================================================

1709                          VertexFit::getErrorMatrix

1710   =======================================================================*/
1711 
1712 double VertexFit::getErrorMatrix(int j, int k) const
1713 {
1714    if(stat!=0) return -999.;
1715 // Return the j,k element of the full error matrix VMAT.

1716 // Since the CTVMFT documentation assumes the indices

1717 // start from 1 (ala Fortran), we will also assume this

1718 // and convert the C++ indices.  Note also the that order of

1719 // Fortran and C++ indices is different.  We assume that

1720 // j and k are in the Fortran order.

1721    if(j<1 || k<1 || j>_maxdim+1 || k>_maxdim) {
1722       return -999.;
1723    }
1724    return _ctvmfr.vmat[k-1][j-1];
1725 }
1726 
1727 HepSymMatrix VertexFit::getErrorMatrix(VertexFit::vertexNumber nv) const 
1728 {  
1729 // Return errors for vertex nv.

1730 
1731   HepSymMatrix cov(3,0);
1732   
1733 // If this is the primary vertex, return the error matrix the

1734 // user supplied.

1735   if(nv==PRIMARY_VERTEX) {
1736     for(int j=0; j<3; j++)
1737       for(int k=0; k<3; k++) 
1738         cov[j][k] = _ctvmq.exyzpv[j][k];
1739     return cov;
1740   }
1741   
1742   if(stat!=0) return cov;
1743 // Return x,y,z position of vertex nv.

1744   if(nv<VERTEX_1 || nv>_ctvmq.nvertx) return cov;
1745 // get offset for vertex nv

1746   int voff = _ctvmq.voff[nv-1];
1747 // fill matrix

1748   for(int i = 0 ; i < 3 ; ++i)
1749     for(int j = i ; j < 3 ; ++j)
1750       cov[i][j] = _ctvmfr.vmat[voff+i][voff+j];
1751   return cov;
1752 }
1753 
1754 HepSymMatrix VertexFit::getErrorMatrix(const CdfTrack* trk) const 
1755 {
1756   HepSymMatrix cov(3,0);
1757   if (stat != 0) return cov;
1758 
1759   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
1760 
1761     // Find which track matches this Id

1762     if (trk->id().value() == _ctvmq.list[nt]) {
1763 
1764       // Position of track nt get offset for track nt

1765       int toff = _ctvmq.toff[nt];
1766 
1767       // Fill matrix -- Crv,Phi,Ctg

1768       for(int i = 0 ; i < 3 ; ++i)
1769         for(int j = i ; j < 3 ; ++j)
1770           cov[i][j] = _ctvmfr.vmat[toff+i][toff+j];
1771     }
1772   }
1773   return cov;
1774 }
1775 
1776 float VertexFit::getPtError(const CdfTrack* trk) const
1777 {  
1778   if (stat != 0) return 0;
1779   
1780   int   toff;
1781   float pt,curv,curvErr,ptErr;
1782 
1783   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
1784 
1785     // Find which track matches this Id

1786     if (trk->id().value() == _ctvmq.list[nt]) {
1787 
1788       // Position of track nt get offset for track nt

1789       toff = _ctvmq.toff[nt];
1790 
1791       // Curvature error

1792       pt      = sqrt(_ctvmq.trkp4[0][nt]*_ctvmq.trkp4[0][nt] +
1793                      _ctvmq.trkp4[1][nt]*_ctvmq.trkp4[1][nt]);
1794       curv    = _ctvmq.pscale/pt;
1795       curvErr = sqrt(_ctvmfr.vmat[toff+0][toff+0]);
1796       ptErr   = _ctvmq.pscale/curv/curv*curvErr;
1797       return ptErr;
1798     }
1799   }
1800 
1801   return 0;
1802 }
1803 
1804 /*=======================================================================

1805                          VertexFit::getPosMomErr

1806   =======================================================================*/
1807 void VertexFit::getPosMomErr(HepMatrix& errors) const
1808 {
1809   // A c++ rewrite of the FORTRAN MASSAG function

1810   // The result of this function is an error matrix in position-momentum basis.

1811   // A 7x7 matrix of errors where the rows/columns are x, y, z, px, py, pz, e.

1812 
1813   double cosph[_maxtrk];
1814   double sinph[_maxtrk];
1815   double cosdph[_maxtrk];
1816   double sindph[_maxtrk];
1817   Hep3Vector mom3[_maxtrk];
1818   HepLorentzVector pmom[_maxtrk];
1819   HepLorentzVector total_mom;
1820 
1821   for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
1822     for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
1823       if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
1824       cosph[ltrk] = cos(_ctvmq.par[ltrk][1]);
1825       sinph[ltrk] = sin(_ctvmq.par[ltrk][1]);
1826       double dphi = 0;
1827       sindph[ltrk] = 2 * _ctvmq.par[ltrk][0] * 
1828        (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
1829         _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk]);
1830       cosdph[ltrk] = sqrt(1.0 - sindph[ltrk] * sindph[ltrk]);
1831       if(fabs(sindph[ltrk]) <= 1.0){
1832         dphi = asin(sindph[ltrk]);
1833       }
1834       double pt = _ctvmq.pscale * fabs(1./_ctvmq.par[ltrk][0]);
1835       mom3[ltrk].setX(pt * cos(_ctvmq.par[ltrk][1] + dphi));
1836       mom3[ltrk].setY(pt * sin(_ctvmq.par[ltrk][1] + dphi));
1837       mom3[ltrk].setZ(pt * _ctvmq.par[ltrk][2]);
1838       double e  = sqrt(_ctvmq.tmass[ltrk] * _ctvmq.tmass[ltrk] 
1839                        + mom3[ltrk].mag2());
1840       pmom[ltrk].setVect(mom3[ltrk]);
1841       pmom[ltrk].setE(e);
1842 
1843       total_mom += pmom[ltrk];
1844     }
1845   }
1846 
1847 // Easy so far, but now it gets ugly.

1848 // Fill dp_dpar with the derivatives of the position and momentum with respect

1849 // to the parameters.

1850 
1851   int ctvmft_dim = 3 * (_ctvmq.nvertx + _ctvmq.ntrack);
1852   HepMatrix deriv(ctvmft_dim, 7, 0);
1853   HepMatrix dp_dpar[_maxvtx] = {deriv, deriv, deriv};
1854 
1855 // Fill the x, y, z rows: 

1856 
1857   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
1858     for(int lcomp = 0; lcomp < 3; lcomp++){
1859       dp_dpar[nvtx][(3 * nvtx) + lcomp][lcomp] = 1.0;
1860     }
1861 
1862 // Fill the px, py, pz, e rows:

1863 
1864     for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
1865       for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
1866         if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
1867 
1868 // Find the derivatives of dphi with respect to x, y, curvature, and phi0:

1869 
1870         double dphi_dx = 2.0 * _ctvmq.par[ltrk][0] * cosph[ltrk]/cosdph[ltrk];
1871         double dphi_dy = 2.0 * _ctvmq.par[ltrk][0] * sinph[ltrk]/cosdph[ltrk];
1872         double dphi_dc = 2.0 * 
1873           (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
1874            _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk])/cosdph[ltrk];
1875         double dphi_dp = 2.0 * _ctvmq.par[ltrk][0] *
1876           (-_ctvmq.xyzvrt[lvtx + 1][0] * sinph[ltrk] + 
1877             _ctvmq.xyzvrt[lvtx + 1][1] * cosph[ltrk])/cosdph[ltrk];
1878 
1879 // Now load the derivative matrix

1880 
1881         int lvele = 3 * lvtx;
1882 // dPx/dx:

1883         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dx;
1884 // dPy/dx:

1885         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dx;
1886 // dPz/dx:

1887         dp_dpar[nvtx][lvele][5]  = 0.; 
1888 // dPy/dx:

1889         dp_dpar[nvtx][lvele][6]  = 0.; 
1890 
1891         lvele++;
1892 // dPx/dy:

1893         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dy;
1894 // dPy/dy:

1895         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dy;
1896 // dPz/dy:

1897         dp_dpar[nvtx][lvele][5]  = 0.; 
1898 // dE/dy:

1899         dp_dpar[nvtx][lvele][6]  = 0.; 
1900 
1901         lvele++;
1902 // dPx/dz:

1903         dp_dpar[nvtx][lvele][3]  = 0.;
1904 // dPy/dz:

1905         dp_dpar[nvtx][lvele][4]  = 0.;
1906 // dPz/dz:

1907         dp_dpar[nvtx][lvele][5]  = 0.; 
1908 // dE/dz:

1909         dp_dpar[nvtx][lvele][6]  = 0.; 
1910 
1911         lvele = 3 * (ltrk + _ctvmq.nvertx);
1912 // dPx/dcurv[ltrk]:

1913         dp_dpar[nvtx][lvele][3]  = -(pmom[ltrk].x()/_ctvmq.par[ltrk][0])
1914                                    - pmom[ltrk].y() * dphi_dc;
1915 // dPy/dcurv[ltrk]:

1916         dp_dpar[nvtx][lvele][4]  = -(pmom[ltrk].y()/_ctvmq.par[ltrk][0])
1917                                    + pmom[ltrk].x() * dphi_dc;
1918 // dPz/dcurv[ltrk]:

1919         dp_dpar[nvtx][lvele][5]  = -(pmom[ltrk].z()/_ctvmq.par[ltrk][0]); 
1920 // dE/dcurv[ltrk]:

1921         dp_dpar[nvtx][lvele][6]  = 
1922           -mom3[ltrk].mag2()/(_ctvmq.par[ltrk][0] * pmom[ltrk].e()); 
1923 
1924         lvele++;
1925 // dPx/dphi[ltrk]:

1926         dp_dpar[nvtx][lvele][3]  = -pmom[ltrk].y() * (1.0 + dphi_dp);
1927 // dPy/dphi[ltrk]:

1928         dp_dpar[nvtx][lvele][4]  =  pmom[ltrk].x() * (1.0 + dphi_dp);
1929 // dPz/dphi[ltrk]:

1930         dp_dpar[nvtx][lvele][5]  = 0.;
1931 // dE/dphi[ltrk]:

1932         dp_dpar[nvtx][lvele][6]  = 0.;
1933 
1934         lvele++;
1935 // dPx/dcot[ltrk]:

1936         dp_dpar[nvtx][lvele][3]  = 0;
1937 // dPy/dcot[ltrk]:

1938         dp_dpar[nvtx][lvele][4]  = 0;
1939 // dPz/dcot[ltrk]:

1940         dp_dpar[nvtx][lvele][5]  = pmom[ltrk].perp();
1941 // dE/dcot[ltrk]:

1942         dp_dpar[nvtx][lvele][6]  = 
1943          pmom[ltrk].perp2() * _ctvmq.par[ltrk][2] / pmom[ltrk].e();
1944       }
1945     }
1946   }
1947 // Now calculate the new error matrix

1948 
1949   // Extract the interesting bits from VMAT.

1950 
1951   HepMatrix vmat(ctvmft_dim,ctvmft_dim,0);
1952   for(int lpar = 0; lpar < ctvmft_dim; lpar++){
1953     int l = lpar%3;
1954     int lvele = 0;
1955     if(lpar < 3  * _ctvmq.nvertx){
1956       int lvtx = lpar/3;
1957       lvele = _ctvmq.voff[lvtx] + l;
1958     }
1959     else{
1960       int ltrk = (lpar - 3 * _ctvmq.nvertx)/3;
1961       lvele = _ctvmq.toff[ltrk] + l;
1962     }
1963     for(int kpar = 0; kpar < ctvmft_dim; kpar++){ 
1964       int k = kpar%3;
1965       int kvele = 0;
1966       if(kpar < 3  * _ctvmq.nvertx){
1967         int kvtx = kpar/3;
1968         kvele = _ctvmq.voff[kvtx] + k;
1969       }
1970       else{
1971         int ktrk = (kpar - 3 * _ctvmq.nvertx)/3;
1972         kvele = _ctvmq.toff[ktrk] + k;
1973       }
1974       vmat[kpar][lpar] = _ctvmfr.vmat[kvele][lvele];
1975     }
1976   }
1977   //  Just about done

1978   HepMatrix ans(7,7,0);
1979   HepMatrix answer[_maxvtx] = {ans, ans, ans};
1980   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
1981     answer[nvtx] = (dp_dpar[nvtx].T() * vmat) * dp_dpar[nvtx];
1982   }
1983   errors = answer[0];
1984 }
1985 
1986 /*=======================================================================

1987                          VertexFit::vOff

1988   =======================================================================*/
1989 int VertexFit::vOff(vertexNumber jv) const
1990 { if(jv < VERTEX_1 || jv > _maxvtx) {return -999;}
1991   else {return _ctvmq.voff[jv-1];}
1992 }
1993 
1994 /*=======================================================================

1995                          VertexFit::tOff

1996   =======================================================================*/
1997 int VertexFit::tOff(const CdfTrack* trk) const
1998 { for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
1999     if(trk->id().value()==_ctvmq.list[kt]) return _ctvmq.toff[kt];
2000   }
2001   return -999;
2002 }
2003 
2004 /*=======================================================================

2005                          VertexFit::pOff

2006   =======================================================================*/
2007 int VertexFit::pOff(int lp) const
2008 { if(lp < 1) {
2009     return -999;
2010   }
2011   else return _ctvmq.poff[lp-1];
2012 }
2013 
2014 /*=======================================================================

2015                          VertexFit::cOff

2016   =======================================================================*/
2017 int VertexFit::cOff(int lc) const
2018 { if(lc < 1) {
2019     return -999;
2020   }
2021   else return _ctvmq.coff[lc-1];
2022 }
2023 
2024 /*=======================================================================

2025                          VertexFit::mOff

2026   =======================================================================*/
2027 int VertexFit::mOff() const
2028 { return _ctvmq.moff; }
2029 
2030 /*=======================================================================

2031                          VertexFit::VMat

2032   =======================================================================*/
2033 double VertexFit::VMat(int i, int j) const
2034 { if(i <0 || j < 0) {return -999.;}
2035   else return _ctvmfr.vmat[i][j];
2036 }
2037 
2038 /*=======================================================================

2039                          VertexFit::allocateVertexNumber

2040   =======================================================================*/
2041 
2042 // Facilitates CandNode recursion.  CandNodes have no way of deciding

2043 // which vertex they are, and these trivial functions help them do that.

2044 VertexFit::vertexNumber VertexFit::allocateVertexNumber()
2045 {
2046   if ((_currentAllocatedVertexNumber < PRIMARY_VERTEX) ||
2047       (_currentAllocatedVertexNumber > _maxvtx)) {
2048     std::cout << "VertexFit::allocateVertexNumber: out of range!" << std::endl;
2049     return PRIMARY_VERTEX;
2050   }
2051   return vertexNumber(++_currentAllocatedVertexNumber);
2052 }
2053 
2054 void VertexFit::resetAllocatedVertexNumber()
2055 {
2056   _currentAllocatedVertexNumber = 0;
2057 }
2058 
2059 /*=======================================================================

2060                          VertexFit::restoreFromCommons()

2061   =======================================================================*/
2062 
2063 void VertexFit::restoreFromCommons() {
2064         stat=0;
2065         _ctvmq_com  = (CTVMQ*) ctvmq_address_();
2066         _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
2067         _fiddle_com = (FIDDLE*) fiddle_address_();
2068         _trkprm_com = (TRKPRM*) trkprm_address_();
2069         _ctvmq  = *_ctvmq_com;
2070         _ctvmfr = *_ctvmfr_com;
2071         _fiddle = *_fiddle_com;
2072         _trkprm = *_trkprm_com;
2073 }
2074 
2075 
2076 /*=======================================================================

2077                          VertexFit::setTrackReferencePoint

2078   =======================================================================*/
2079 void VertexFit::setTrackReferencePoint( const Hep3Vector & ref )
2080 {
2081   // Set extrapolation flag

2082   _extrapolateTrackErrors = true; 
2083   // Keep new reference point

2084   _referencePoint = ref;
2085   // Keep track of old primary vertex in CDF coordinates

2086   _cdfPrimaryVertex = _primaryVertex;
2087   // Set new primary vertex relative to new reference point

2088   setPrimaryVertex( _cdfPrimaryVertex - _referencePoint);
2089 }
2090 
2091 /*=======================================================================

2092                          VertexFit::moveReferencePoint()

2093   =======================================================================*/
2094 
2095 void VertexFit::moveReferencePoint( HepVector & v,
2096                                     HepSymMatrix & m ) {
2097 
2098   // Move the reference point of the track to _referencePoint

2099   if ( !_extrapolateTrackErrors ) 
2100     return;
2101   
2102   TrackExtrapolator e;
2103 
2104   // Modifies v and m

2105   e.moveReferencePoint( v, m, _referencePoint );
2106 }
2107 
2108 
2109 #endif  /* VERTEXFIT_CC */
2110 
2111
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Send problems or questions to cdfcode@fnal.gov