utility.cxx

Go to the documentation of this file.

#include "EvtRecEvent/EvtRecVeeVertex.h"
#include "DTagAlg/LocalLambdaSelector.h"
#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);
  }
 
  vector<int> trkIdxCol;
  trkIdxCol.clear();
  for (int i = 0; i < nvalid; i++)
    trkIdxCol.push_back(i);
  vtxfit->AddVertex(0, vxpar, trkIdxCol);
  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);
 
}
 
HepLorentzVector utility::vfit(string channel, vector<int> kaonid, vector<int> pionid, vector<int> protonid,  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()+protonid.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);
    }
 
    for(int i=0; i<protonid.size();++i){
            EvtRecTrackIterator itTrk=charged_begin + protonid[i];
            RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
            pold+=getp4(mdcKalTrk, 4);
            WTrackParameter wtrk(xmass[4],mdcKalTrk->getZHelixP(),mdcKalTrk->getZErrorP() );
            vtxfit->AddTrack(kaonid.size()+pionid.size()+i, wtrk);
    }
 
    vector<int> trkIdxCol;
    trkIdxCol.clear();
    for (int i = 0; i < nvalid; i++)
            trkIdxCol.push_back(i);
    vtxfit->AddVertex(0, vxpar, trkIdxCol);
    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);
 
}
 
 
 
vector<double> utility::SecondaryVFit(EvtRecVeeVertex* ks, IVertexDbSvc* vtxsvc){
 
  // by default: chi2 = -999, length = -999, error = 999
  double vfitchi2 = -999;
  double vfitlength = -999;
  double vfiterror = 999;
 
  vector<double> results;
  results.push_back(vfitchi2);
  results.push_back(vfitlength);
  results.push_back(vfiterror);
 
  // --------------------------------------------------
  // Do a secondary vertex fit and check the flight significance
  // --------------------------------------------------
 
  EvtRecTrack* veeTrack1 = ks->pairDaughters().first;
  RecMdcKalTrack* veeKalTrack1 = veeTrack1->mdcKalTrack();
  WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
 
  EvtRecTrack* veeTrack2 = ks->pairDaughters().second;
  RecMdcKalTrack* veeKalTrack2 = veeTrack2->mdcKalTrack();
  WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
 
  VertexParameter wideVertex;
  HepPoint3D vWideVertex(0., 0., 0.);
  HepSymMatrix evWideVertex(3, 0);
 
  evWideVertex[0][0] = 1.0e12;
  evWideVertex[1][1] = 1.0e12;
  evWideVertex[2][2] = 1.0e12;
 
  wideVertex.setVx(vWideVertex);
  wideVertex.setEvx(evWideVertex);
 
  // First, perform a vertex fit
  VertexFit* vtxfit = VertexFit::instance();
  vtxfit->init();
 
  // add the daughters
  vtxfit->AddTrack(0,veeInitialWTrack1);
  vtxfit->AddTrack(1,veeInitialWTrack2);
  vtxfit->AddVertex(0,wideVertex,0,1);
 
  // do the fit
  vtxfit->Fit(0);
  vtxfit->Swim(0);
  vtxfit->BuildVirtualParticle(0);
 
  // Now perform the secondary vertex fit
  SecondVertexFit* svtxfit = SecondVertexFit::instance();
  svtxfit->init();
 
  // add the primary vertex
  VertexParameter beamSpot;
  HepPoint3D vBeamSpot(0., 0., 0.);
  HepSymMatrix evBeamSpot(3, 0);
 
  if(vtxsvc->isVertexValid()){
    double* dbv = vtxsvc->PrimaryVertex();
    double*  vv = vtxsvc->SigmaPrimaryVertex();
    for (unsigned int ivx = 0; ivx < 3; ivx++){
      vBeamSpot[ivx] = dbv[ivx];
      evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
    }
  }
  else{
    cout << "KSSELECTOR ERROR:  Could not find a vertex from VertexDbSvc" << endl;
    return results;
  }
 
  beamSpot.setVx(vBeamSpot);
  beamSpot.setEvx(evBeamSpot);
 
  VertexParameter primaryVertex(beamSpot);
  svtxfit->setPrimaryVertex(primaryVertex);
 
  // add the secondary vertex
  svtxfit->setVpar(vtxfit->vpar(0));
 
  // add the Ks or lambda
  svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
 
  // do the second vertex fit
  // if the fit fails, the default values will not be changed
  if( !svtxfit->Fit() ) return results;
 
  // save the new ks parameters
  vfitlength = svtxfit->decayLength();
  vfiterror = svtxfit->decayLengthError();
  vfitchi2 = svtxfit->chisq();
 
  results.clear();
  results.push_back(vfitchi2);
  results.push_back(vfitlength);
  results.push_back(vfiterror);
 
  return results;
}
 
vector<double> utility::SecondaryVFit_Lambda(EvtRecVeeVertex* lambda, IVertexDbSvc* vtxsvc){
 
    // by default: chi2 = -999, length = -999, error = 999
    double vfitchi2 = -999;
    double vfitlength = -999;
    double vfiterror = 999;
 
    vector<double> results;
    results.push_back(vfitchi2);
    results.push_back(vfitlength);
    results.push_back(vfiterror);
 
    int index[2]  = {0, 1};
    if (lambda->vertexId() < 0){
            index[0] = 1;  
            index[1] = 0;
    } 
 
    EvtRecTrack* recTrk_proton = lambda->daughter(index[0]);
    EvtRecTrack* recTrk_pion = lambda->daughter(index[1]);
 
    // --------------------------------------------------
    // Do a secondary vertex fit and check the flight significance
    // --------------------------------------------------
 
    RecMdcKalTrack* mdcKalTrk_proton = recTrk_proton->mdcKalTrack();
    mdcKalTrk_proton->setPidType(RecMdcKalTrack::proton);
    WTrackParameter WTrk_proton;
    if(lambdaSelector.IfProtonPID()) WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelixP(), mdcKalTrk_proton->getZErrorP());
    else  WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelix() , mdcKalTrk_proton->getZError());
 
    RecMdcKalTrack* mdcKalTrk_pion = recTrk_pion->mdcKalTrack();
    mdcKalTrk_pion->setPidType(RecMdcKalTrack::pion);
    WTrackParameter WTrk_pion = WTrackParameter(0.13957018, mdcKalTrk_pion->getZHelix(), mdcKalTrk_pion->getZError());
 
    VertexParameter wideVertex;
    HepPoint3D vWideVertex(0., 0., 0.);
    HepSymMatrix evWideVertex(3, 0);
 
    evWideVertex[0][0] = 1.0e12;
    evWideVertex[1][1] = 1.0e12;
    evWideVertex[2][2] = 1.0e12;
 
    wideVertex.setVx(vWideVertex);
    wideVertex.setEvx(evWideVertex);
 
    // First, perform a vertex fit
    VertexFit* vtxfit = VertexFit::instance();
    vtxfit->init();
 
    // add the daughters
    vtxfit->AddTrack(0,WTrk_proton);
    vtxfit->AddTrack(1,WTrk_pion);
    vtxfit->AddVertex(0,wideVertex,0,1);
 
    // do the fit
    vtxfit->Fit(0);
    vtxfit->Swim(0);
    vtxfit->BuildVirtualParticle(0);
 
    // Now perform the secondary vertex fit
    SecondVertexFit* svtxfit = SecondVertexFit::instance();
    svtxfit->init();
 
    // add the primary vertex
    VertexParameter beamSpot;
    HepPoint3D vBeamSpot(0., 0., 0.);
    HepSymMatrix evBeamSpot(3, 0);
 
    if(vtxsvc->isVertexValid()){
            double* dbv = vtxsvc->PrimaryVertex();
            double*  vv = vtxsvc->SigmaPrimaryVertex();
            for (unsigned int ivx = 0; ivx < 3; ivx++){
                    vBeamSpot[ivx] = dbv[ivx];
                    evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
            }
    }
    else{
            cout << "LambdaSELECTOR ERROR:  Could not find a vertex from VertexDbSvc" << endl;
            return results;
    }
 
    beamSpot.setVx(vBeamSpot);
    beamSpot.setEvx(evBeamSpot);
 
    VertexParameter primaryVertex(beamSpot);
    svtxfit->setPrimaryVertex(primaryVertex);
 
    // add the secondary vertex
    svtxfit->setVpar(vtxfit->vpar(0));
 
    // add the Ks or lambda
    svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
 
    // do the second vertex fit
    // if the fit fails, the default values will not be changed
    if( !svtxfit->Fit() ) return results;
 
    // save the new ks parameters
    vfitlength = svtxfit->decayLength();
    vfiterror = svtxfit->decayLengthError();
    vfitchi2 = svtxfit->chisq();
 
    results.clear();
    results.push_back(vfitchi2);
    results.push_back(vfitlength);
    results.push_back(vfiterror);
 
    return results;
}