---

w_YFMVTR.cpp

Go to the documentation of this file.

#include "w_YFMVTR.h"
#include <iostream.h>
//=======================================================================
// w_YFMVTR methods
// 11/2001 Kay Huettmann
//=======================================================================

//-----------------------------------------------------------------------
// getResult(..)
// YFMVTR wrapper (YTOP vertex fit)
// Parameters are passed via a common block, since cfortan is limited to 27 parameters
//-----------------------------------------------------------------------
bool w_YFMVTR::getResult(const vector<Helix>& vHelix,
                           HepLorentzVector& track4V, HepSymMatrix& covTrack, 
                           Hep3Vector& vertex3V, HepSymMatrix& covVertex, 
                           HepMatrix& corrVertexTrack,
                           float& chi2) const {
 
   // set dummy vertex constraint 
  Vertex theVertexC;
  
  // initialize output parameters
  HepLorentzVector temp4V(0.,0.,0.,0.);
  track4V = temp4V;
 
  HepSymMatrix tempCov4V(4,0);
  covTrack = tempCov4V;

  Hep3Vector temp3V(0.,0.,0.);
  vertex3V = temp3V;

  HepSymMatrix tempCov3V(3,0);
  covVertex = tempCov3V;

  HepMatrix tempCorr(3,3,0);
  corrVertexTrack = tempCorr;

  float _chi2 = -1.;
  chi2 = _chi2;


  if (vHelix.size() > 10) 
    cerr << "__w_YFMVTR::getResult: Too many tracks (" << vHelix.size() 
         << "). Surplus tracks will be ignored in the fit !!" << endl;

  // Count charged and neutral tracks
  int n_charged = 0;
  int n_neutral = 0;
  vector<const Helix>::iterator itH = vHelix.begin();
  for ( ; itH < vHelix.end(); itH++ ) {
    if (itH->charge() != 0) {
      n_charged++;
    } else {
      n_neutral++;
    }
  }

  // fill YMFVTR common block
  FCB.nvx = 0;
  FCB.nhx = 0;
  FCB.neu = 0;
  FCB.lvapc = 1;
  FCB.lpsum = 1;
  FCB.lmcalc = 1;
  FCB.lvfix = 0;
  FCB.vxin[0] = 0.;
  FCB.vvxin[0] = 0.;

  // Check whether there is a vertex constraint
  if ( _pVertexC ) {
    FCB.nvx = 1;
    FCB.lvapc = 0;
    // vertex coordinates ...
    FCB.vxin[0] = _pVertexC->v3V().x();
    FCB.vxin[1] = _pVertexC->v3V().y();
    FCB.vxin[2] = _pVertexC->v3V().z();
    // covariance matrix 
    int j = 0;
    for (int m = 1; m <= 3; m++) {
      for (int n = 1; n <= m; n++) {
        FCB.vvxin[j++] = _pVertexC->cov_matrix(m,n) ;
      }
    } 
  }
  
  int ntrk = 0; // total number of tracks
  int nchg = 0;
  int nneu = 0;
  int indf;
  int indv;
  
  float cmass[10];
  float nmass[10];
  
 
  for (itH = vHelix.begin() ; itH < vHelix.end(); itH++ ) {

    if ((itH - vHelix.begin()) > 10) break;   // too many tracks
   
    ntrk++; 
    indf = 5*(ntrk -1);
    indv = 15*(ntrk -1);

    //
    // ---  charged track  ---  
    //
    if (itH->charge() != 0) {
      nchg++;
      FCB.nhx++;  
      FCB.ixhx[nchg -1] = ntrk;    // ntrk is a FORTRAN array index, so it has to start at 1
      cmass[nchg -1] = itH->mass();
    
      // track helix parameters
      FCB.hxin[indf   ] = itH->IR();
      FCB.hxin[indf+ 1] = itH->TL();
      FCB.hxin[indf+ 2] = itH->P0();
      FCB.hxin[indf+ 3] = itH->D0();
      FCB.hxin[indf+ 4] = itH->Z0();

      // track covariance matrix
      int j = 0;
      for (int m = 1; m <= 5; m++) {
        for (int n = 1; n <= m; n++) { 
          j++;
          FCB.vhxin[indv +j -1] = itH->cov_matrix(m,n);
        }
      } 
     


    } else {
      //
      // ---  neutral track  ---
      //
      nneu++;
      FCB.neu++;
      FCB.ixnu[nneu -1] = ntrk;
      nmass[nneu -1] = itH->mass();
    
      // track parameters
      FCB.hxin[indf   ] = itH->IR();
      FCB.hxin[indf+ 1] = itH->TL();
      FCB.hxin[indf+ 2] = itH->P0();
      FCB.hxin[indf+ 3] = itH->D0();
      FCB.hxin[indf+ 4] = itH->Z0();

      // track covariance matrix
      int j = 0;
      for (int m = 1; m <= 5; m++) {
        for (int n = 1; n <= m; n++) { 
          j++;
          FCB.vtnui[indv +j -1] = itH->cov_matrix(m,n);
        }
      } 



    }
  
  }

  int j = 0;
  int i;
  for (i = 0; i < nchg; i++) {
    FCB.ampc[j++] = cmass[i];
  }
  for (i = 0; i < nneu; i++) {
    FCB.ampc[j++] = nmass[i];
  }
  for (i = 0; i < ntrk*5; i++) {
    FCB.tnui[i] = FCB.hxin[i];
  }

  FCB.npidc = 1;
  

  
  

  // call the Fortran routine
  int ifail = W_YFMVTR();

  if (! ifail) {

    float ptrk[5];
    float vptrk[10];
   
    ptrk[0] = FCB.psum[0];
    ptrk[1] = FCB.psum[1];
    ptrk[2] = FCB.psum[2];
    ptrk[3] = sqrt(FCB.amass*FCB.amass + ptrk[0]*ptrk[0] +
                   ptrk[1]*ptrk[1] + ptrk[2]*ptrk[2]);
    ptrk[4] = FCB.amass;


   // calculate chi2/d.o.f.
   int ndf = ( _pVertexC ) ?  ntrk*2 : ntrk*2 -3;
   _chi2 = FCB.chisq/((float) ndf);

   // protection against crazy values from YTOP fit (taken from QYFIT)
   if (_chi2 < 500.) {
     if (ptrk[4] > 500. || ptrk[0] > 1000. || ptrk[1] > 1000. ||
         ptrk[2] > 1000. ) {
       _chi2 = -1.;
     }
   } else {
     _chi2 = -1.;
   }
 
   // calculate correct 4-momentum covariance matrix
   if (_chi2 > 0.) {
     
     float r1 = ptrk[0]/ptrk[3];
     float r2 = ptrk[1]/ptrk[3];
     float r3 = ptrk[2]/ptrk[3];
     float r4 = ptrk[4]/ptrk[3]; 
     
     for (int i = 0; i <= 6; i++) {
       vptrk[i] = FCB.vpsum[i];
     }

     vptrk[9] = r1*r1*vptrk[0] + r2*r2*vptrk[2] +
       r3*r3*vptrk[5] + r4*r4*FCB.dmass*FCB.dmass +
       2.*r1*(r2*vptrk[1] + r3*vptrk[3]) +
       2.*r2*r3*vptrk[4] +
       2.*r4*(r1*FCB.vmps[0] + r2*FCB.vmps[1] + r3*FCB.vmps[2]);
     
     vptrk[6] = r1*vptrk[0] + r2*vptrk[1] + r3*vptrk[3]
       +  r4*FCB.vmps[0];
   
     vptrk[7] = r1*vptrk[1] + r2*vptrk[2] + r3*vptrk[4]
       + r4*FCB.vmps[1];

     vptrk[8] = r1*vptrk[3] + r2*vptrk[4] + r3*vptrk[5]
       + r4*FCB.vmps[2];

     // prepare return values

     //  mother track
     track4V.setPx(ptrk[0]);
     track4V.setPy(ptrk[1]);
     track4V.setPz(ptrk[2]);
     track4V.setE(ptrk[3]);

     // vertex
      vertex3V.setX(FCB.vxout[0]);
      vertex3V.setY(FCB.vxout[1]);
      vertex3V.setZ(FCB.vxout[2]);

     // 4-momentum covariance matrix
      int j = 0;
      int m;
      int n;
      for (m = 1 ; m <= 4; m++) {
        for (n = 1; n <= m; n++) {
          covTrack(m,n)= vptrk[j++];
        }
      }

     // vertex covariance matrix
      j = 0;
      for (m = 1; m <= 3; m++) { 
        for (n = 1; n <= m; n++) {
          covVertex(m,n)= FCB.vvxou[j++];
        }
      }
 
      // correlation between momentum annd vertex
      j = 0;
      for (m = 1; m <= 3; m++) {
        for (n = 1; n <= 3; n++) {
          // NOTE: You must take the transposed of the Fortran array.
          corrVertexTrack(n,m)= FCB.vpsvx[j++];
        }
      }

     // chi2
     chi2 = _chi2; 

   }
 }



 return (ifail != 0);
}

---