utility.cxx

Go to the documentation of this file.

#include "DTagAlg/utility.h"
 
HepLorentzVector utility::getp4(RecMdcKalTrack* mdcKalTrack, int pid){
  
  HepVector zhelix;
  double mass=0;
 
  if(pid==0){
    zhelix=mdcKalTrack->getZHelixE();
    mass=0.000511;
  }
  else if(pid==1){
    zhelix=mdcKalTrack->getZHelixMu();
    mass= 0.105658;
  }
  else if(pid==2){
    zhelix=mdcKalTrack->getZHelix();
    mass=0.139570;
  }
  else if(pid==3){
    zhelix=mdcKalTrack->getZHelixK();
    mass= 0.493677;
  }
  else{
    zhelix=mdcKalTrack->getZHelixP();
    mass= 0.938272;
  }
 
  double dr(0),phi0(0),kappa(0),dz(0),tanl(0);
  dr=zhelix[0];
  phi0=zhelix[1];
  kappa=zhelix[2];
  dz=zhelix[3];
  tanl=zhelix[4];
 
  int charge=0;
  if (kappa > 0.0000000001)
    charge = 1;
  else if (kappa < -0.0000000001)
    charge = -1;
 
  double pxy=0;
  if(kappa!=0) pxy = 1.0/fabs(kappa);
 
  double px = pxy * (-sin(phi0));
  double py = pxy * cos(phi0);
  double pz = pxy * tanl;
 
  double e  = sqrt( pxy*pxy + pz*pz + mass*mass );
 
  return HepLorentzVector(px, py, pz, e);
 
 
}
 
 
HepLorentzVector utility::vfit(string channel, vector<int> kaonid, vector<int> pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
  //use charged tracks only, except pions from Ks
 
  HepLorentzVector pchange(0,0,0,0);
 
  int nvalid= kaonid.size()+pionid.size();
  if(nvalid<=1)
    return pchange;
 
  HepSymMatrix Evx(3, 0);
  double bx = 1E+6; Evx[0][0] = bx*bx;
  double by = 1E+6; Evx[1][1] = by*by;
  double bz = 1E+6; Evx[2][2] = bz*bz;
  VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
  double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
 
 
  VertexFit* vtxfit = VertexFit::instance();
  vtxfit->init();
 
 
  HepLorentzVector pold(0,0,0,0);
 
  for(int i=0; i<kaonid.size();++i){
    EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
    RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
    pold+=getp4(mdcKalTrk, 3);
    WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
    vtxfit->AddTrack(i, wtrk);
  }
 
  for(int i=0; i<pionid.size();++i){
    EvtRecTrackIterator itTrk=charged_begin + pionid[i];
    RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
    pold+=getp4(mdcKalTrk, 2);
    WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
    vtxfit->AddTrack(kaonid.size()+i, wtrk);
  }
 
 
  vtxfit->AddVertex(0, vxpar, 0, 1);
  if(!vtxfit->Fit(0))
    return pchange;
 
  vtxfit->Swim(0);
 
  HepLorentzVector pnew(0,0,0,0);
 
  for(int i=0; i<nvalid;++i){
    WTrackParameter wtrk = vtxfit->wtrk(i);
    HepVector trk_val = HepVector(7,0);
    trk_val = wtrk.w();
    HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
    pnew+=P_trk;
  }
 
  return (pnew-pold);
 
}