---

w_YFMVTC.cpp

Go to the documentation of this file.

#include "w_YFMVTC.h"
#include <iostream.h>

//-----------------------------------------------------------------------
// getResult(..)
// YFMVTC wrapper (YTOP vertex fit with mass constraint)
// Parameters are passed via a common block (cfortan is limited to 27 parameters)
//-----------------------------------------------------------------------
bool w_YFMVTC::getResult(const vector<Helix>& vHelix,
                         HepLorentzVector& track4V, HepSymMatrix& covTrack, 
                         Hep3Vector& vertex3V, HepSymMatrix& covVertex, 
                         HepMatrix& corrVertexTrack,
                         float& chi2) const {
 
  // 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
  VFMCB.nvx = 0;
  VFMCB.nhx = 0;
  VFMCB.neu = 0;
  VFMCB.lvapc = 1;
  VFMCB.lpsum = 1;
  VFMCB.lmcalc = 1;
  VFMCB.lvfix = 0;
  VFMCB.vxin[0] = 0.;
  VFMCB.vvxin[0] = 0.;
  VFMCB.lvapc = 0;

  
  // Check whether there is a mass constraint
  if ( _pMassC ) {
    cout << "w_YFMVTC::getResult() mass constraint = " << *_pMassC << endl;
    VFMCB.lvapc = 1;
    VFMCB.conmas = *_pMassC;
  }

  VFMCB.nvx = 0;

  
  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++;
      VFMCB.nhx++;  
      VFMCB.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
      VFMCB.hxin[indf   ] = itH->IR();
      VFMCB.hxin[indf+ 1] = itH->TL();
      VFMCB.hxin[indf+ 2] = itH->P0();
      VFMCB.hxin[indf+ 3] = itH->D0();
      VFMCB.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++;
          VFMCB.vhxin[indv +j -1] = itH->cov_matrix(m,n);
        }
      } 
     


    } else {
      //
      // ---  neutral track  ---
      //
      nneu++;
      VFMCB.neu++;
      VFMCB.ixnu[nneu -1] = ntrk;
      nmass[nneu -1] = itH->mass();
    
      // track parameters
      VFMCB.hxin[indf   ] = itH->IR();
      VFMCB.hxin[indf+ 1] = itH->TL();
      VFMCB.hxin[indf+ 2] = itH->P0();
      VFMCB.hxin[indf+ 3] = itH->D0();
      VFMCB.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++;
          VFMCB.vtnui[indv +j -1] = itH->cov_matrix(m,n);
        }
      } 



    }
  
  }

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

  VFMCB.npidc = 1;

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

  if (! ifail) {

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


   // calculate chi2/d.o.f.
    int ndf = ( _pMassC ) ?  ntrk*2 -3 +1 : ntrk*2 -3 ; 
   _chi2 = VFMCB.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] = VFMCB.vpsum[i];
     }

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

     vptrk[8] = r1*vptrk[3] + r2*vptrk[4] + r3*vptrk[5]
       + r4*VFMCB.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(VFMCB.vxout[0]);
      vertex3V.setY(VFMCB.vxout[1]);
      vertex3V.setZ(VFMCB.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)= VFMCB.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)= VFMCB.vpsvx[j++];
        }
      }

     // chi2
     chi2 = _chi2; 

   }
 }



 return (ifail != 0);

}

---