FullJetFinder.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FullJetFinder.cc

#include "FullJetFinder.h"

#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/PHTFileServer.h>

#include <phool/PHCompositeNode.h>
#include <phool/getClass.h>

#include <jetbase/JetContainer.h>
#include <jetbase/JetMap.h>
#include <jetbase/Jet.h>

#include <centrality/CentralityInfo.h>

#include <globalvertex/GlobalVertex.h>
#include <globalvertex/GlobalVertexMap.h>

#include <particleflowreco/ParticleFlowElementv1.h>
#include <particleflowreco/ParticleFlowElementContainer.h>

#define HomogeneousField
#include <KFParticle.h>

#include <g4main/PHG4TruthInfoContainer.h>
#include <g4main/PHG4Particle.h>

#include <trackbase_historic/SvtxTrackMap.h>
#include <trackbase_historic/TrackAnalysisUtils.h>

#include <trackbase/MvtxDefs.h>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#include <HepMC/GenEvent.h>
#include <HepMC/GenVertex.h>
#pragma GCC diagnostic pop

#include <phhepmc/PHHepMCGenEvent.h>
#include <phhepmc/PHHepMCGenEventMap.h>

#include <TDatabasePDG.h>

#include <TTree.h>
#include <TH1.h>

#include <cmath>



//____________________________________________________________________________..
FullJetFinder::FullJetFinder(const std::string& outputfilename):
 SubsysReco("FullJetFinder")
 , m_recoJetName()
 , m_truthJetName()
 , m_outputFileName(outputfilename)
 , m_ptRangeReco(0, 100)
 , m_ptRangeTruth(0, 100)
 , m_doTruthJets(0)
 , m_T()
{
  std::cout << "FullJetFinder::FullJetFinder(const std::string &name) Calling ctor" << std::endl;
}

//____________________________________________________________________________..
FullJetFinder::~FullJetFinder()
{
  std::cout << "FullJetFinder::~FullJetFinder() Calling dtor" << std::endl;
}

//_________________________________________________________
void FullJetFinder::Container::Reset()
{
  reco_jet_n = -1;
  truth_jet_n = -1;
  centrality = -1;
  impactparam = -1;
  recojets.clear();
  truthjets.clear();
  primaryVertex.clear();
}

//____________________________________________________________________________..
int FullJetFinder::Init(PHCompositeNode *topNode)
{
  std::cout << "FullJetFinder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
  PHTFileServer::get().open(m_outputFileName, "RECREATE");
  std::cout << "FullJetFinder::Init - Output to " << m_outputFileName << std::endl;

  if (m_inputs == 0 || m_inputs > 10){
    std::cout << "MyJetAnalysis::Init - Error number of inputs outside (0,10> " << std::endl;
    exit(-1);
  }

  // configure Tree
  for(int i = 0 ; i < m_inputs; i++){
    //m_T[i] = new TTree(m_TreeNameCollection.at(i).data(), m_TreeNameCollection.at(i).data());
    m_T[i] = new TTree("Data", "Data");
    m_container[i] = new Container;
    m_T[i]->Branch( "JetContainer", &m_container[i] );
  
    TString tmp = "";
    m_chi2ndf[i] = new TH1D(tmp +"chi2ndf_" + m_TreeNameCollection.at(i).data(),tmp +"chi2ndf_" + m_TreeNameCollection.at(i).data(),1000,0,100);
    m_mvtxcl[i] = new TH1I(tmp +"mvtxcl_" + m_TreeNameCollection.at(i).data(),tmp +"mvtxcl_" + m_TreeNameCollection.at(i).data(),6,0,6);
    m_inttcl[i] = new TH1I(tmp +"inttcl_" + m_TreeNameCollection.at(i).data(),tmp +"inttcl_" + m_TreeNameCollection.at(i).data(),6,0,6);
    m_mtpccl[i]= new TH1I(tmp +"tpccl_" + m_TreeNameCollection.at(i).data(),tmp +"tpccl_" + m_TreeNameCollection.at(i).data(),100,0,100);
  }

  m_stat = new TH1I("event_stat","event_stat",10,-0.5,9.5);
  m_stat->GetXaxis()->SetBinLabel(1,"n_events");
  m_stat->GetXaxis()->SetBinLabel(2,"GV_exists");
  m_stat->GetXaxis()->SetBinLabel(3,"GV_notempty");



  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::InitRun(PHCompositeNode *topNode)
{
  std::cout << "FullJetFinder::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::process_event(PHCompositeNode *topNode)
{  
  m_stat->Fill(0);
  //centrality
  CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
  if (!cent_node){
    std::cout << "MyJetAnalysis::process_event - Error can not find centrality node " << std::endl;
    exit(-1);
  }

  // interface to reco jets
  if ((m_recoJetName.size() != m_truthJetName.size()) || m_recoJetName.empty() || m_truthJetName.empty()){
    std::cout << "MyJetAnalysis::process_event - Error in m_recoJetName size !=  m_truthJetName.size() or empty" << std::endl;
    exit(-1);
  }

  if (m_recoJetName.size() != static_cast<long unsigned int>(m_inputs)){
    std::cout << "MyJetAnalysis::process_event - Error number of AddInput() calls does not match with number of inputs specified in constructor" << std::endl;
    exit(-1);
  }



  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
  if (!vertexmap){
    std::cout << "MyJetAnalysis::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
    assert(vertexmap);  // force quit
    exit(-1);
  }

  m_stat->Fill(1);

  if (vertexmap->empty()){
    std::cout << "MyJetAnalysis::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
    return Fun4AllReturnCodes::ABORTEVENT;
  }

  m_stat->Fill(2);


//std::cout<<"processing event "<<std::endl;
/*
new event
A -0.000921115 0.00567926 2.26759 0 0
source 400
vtx -0.000921115 0.00567926 2.26759

A 0 0 -11.2419 1 1
source 300
vtx 0 0 -11.2419

A -0.000921115 0.00567926 2.26759 2 2
source 400
vtx -0.000921115 0.00567926 2.26759*/



  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
  if (!truthinfo)
  {
    std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
     exit(-1);
  }

  //loop over multiple JetReco->add_algo()
  for(long unsigned int input = 0; input < m_recoJetName.size(); input++){
    //reset per event container
    m_container[input]->Reset();
    
    //get the event centrality/impact parameter from HIJING
    m_container[input]->centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);
    m_container[input]->impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);

    

    // interface to reco jets
    JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName.at(input));
    if (!jets){
      std::cout << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node " << m_recoJetName.at(input) << std::endl;
      exit(-1);
    }

    //interface to truth jets
    //JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName.at(input));
    JetContainer* jetsMC = findNode::getClass<JetContainer>(topNode, m_truthJetName.at(input));
    if (!jetsMC && m_doTruthJets){
      std::cout << "MyJetAnalysis::process_event - Error can not find DST Truth JetContainer node "     << m_truthJetName.at(input) << std::endl;
      exit(-1);
    }

    //get associated hepMCGenEvent
    HepMC::GenEvent *hepMCGenEvent = nullptr;
    if(m_doTruthJets){
      PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");

      if (!hepmceventmap){
        std::cout << PHWHERE
            << "HEPMC event map node is missing, can't collected HEPMC truth particles"
            << std::endl;
        return 0;
      }

      PHHepMCGenEvent *hepmcevent = hepmceventmap->get(1);

      if (!hepmcevent)
      {
        std::cout << PHWHERE
            << "PHHepMCGenEvent node is missing, can't collected HEPMC truth particles"
            << std::endl;
        return 0;
      }

      hepMCGenEvent = hepmcevent->getEvent();

      if(!hepMCGenEvent) return Fun4AllReturnCodes::ABORTEVENT;
    }

    //get Particle Flow container
    ParticleFlowElementContainer *pflowContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");

    if(!pflowContainer){
      std::cout << PHWHERE
          << "ParticleFlowElements node is missing, can't collect particles"
          << std::endl;
      return -1;
    }

    for(auto vertex : *vertexmap){
    PrimaryVertex primary;
    std::vector<GlobalVertex::VTXTYPE> source;
    source.clear();
    primary.id = vertex.second->get_id();
    primary.x = vertex.second->get_x();
    primary.y = vertex.second->get_y();
    primary.z = vertex.second->get_z();
    primary.x_unc = sqrt(vertex.second->get_error(0,0));
    primary.y_unc = sqrt(vertex.second->get_error(1,1));
    primary.z_unc = sqrt(vertex.second->get_error(2,2));
    primary.chisq = vertex.second->get_chisq();
    primary.ndf = vertex.second->get_ndof();
    
     //std::cout<<std::endl<<"A "<<vertex.second->get_x()<<" "<<vertex.second->get_y()<<" "<<vertex.second->get_z()<<" "<<vertex.second->get_id()<<" "<<vertex.first<<std::endl;
    for (auto iter = vertex.second->begin_vertexes();   iter != vertex.second->end_vertexes();   ++iter){
      source.push_back(iter->first);
      GlobalVertex::VertexVector vtx = iter->second;
      //std::cout<<"source "<<iter->first<<std::endl;
     // for(auto vx :vtx){
        //std::cout<<"vtx "<<vx->get_x()<<" "<<vx->get_y()<<" "<<vx->get_z()<<" "<<vx->get_chisq()<<std::endl;
      //}
    }
    if((std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::SVTX) != source.end()) && (std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::MBD) != source.end())) primary.vtxtype = GlobalVertex::VTXTYPE::SVTX_MBD;
    else if(std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::SVTX) != source.end()) primary.vtxtype = GlobalVertex::VTXTYPE::SVTX;
    else if(std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::MBD) != source.end()) primary.vtxtype = GlobalVertex::VTXTYPE::MBD;
    m_container[input]->primaryVertex.push_back(primary);
  }

    int nrecojet = 0;

    Jet::PROPERTY recojet_area_index = jets->property_index(Jet::PROPERTY::prop_area);

    for (Jet* jet : *jets){
      if(jet->get_pt() < m_ptRangeReco.first || jet->get_pt() > m_ptRangeReco.second) continue;
      if (not (std::abs(jet->get_eta()) <= 1.1 - (doFiducial?jetR.at(input):0))) continue;

      //std::cout<<jet->get_pt()<<std::endl;
      nrecojet++;

      int nChtracks = 0;

      RecoJets recojet;

      //loop over jet constituents
      for (const auto& comp : jet->get_comp_vec()){
        ParticleFlowElement *pflow = pflowContainer->getParticleFlowElement(comp.second);
        //if charged track
        if(pflow->get_track()){
          chConstituent track_properties;

          //get vertex associated to tarck
          int id = pflow->get_track()->get_vertex_id();
          GlobalVertex *vtx = vertexmap->get(id);
          //skip if track without vertex
          if (vtx == nullptr){
            std::cout << "MyJetAnalysis::process_event - Fatal Error - Gvtx null." << std::endl;
            continue;
          }
          nChtracks++;

          float DCA_xy, DCA_xy_unc;
          GetDistanceFromVertexXY(pflow->get_track(),vtx,DCA_xy,DCA_xy_unc); //rerutn DCA_XY vector, val ||DCA|| + unc DCA
          double dot = (pflow->get_track()->get_x() - vtx->get_x()) * jet->get_px() + (pflow->get_track()->get_y() - vtx->get_y()) * jet->get_py();
          double sign = int(dot/std::abs(dot));

          //tarcking team way - no uncertainty of vertex in DCA calculation
          /*Acts::Vector3 vtxActs(vtx->get_x(), vtx->get_y(), vtx->get_z());
          std::pair<std::pair<float, float>, std::pair<float, float>> DCApair;
          DCApair = TrackAnalysisUtils::get_dca(pflow->get_track(),vtxActs);
          double dot2 = (pflow->get_track()->get_x() - vtx->get_x()) * jet->get_px() + (pflow->get_track()->get_y() - vtx->get_y()) * jet->get_py();
          double sign2 = int(dot2/std::abs(dot2));*/
      
          //get some track quality values
          int n_mvtx_hits = 0;
          int n_intt_hits = 0;
          int n_tpc_hits = 0;

          for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(pflow->get_track())){
            switch (TrkrDefs::getTrkrId(ckey)){
              case TrkrDefs::mvtxId:
                n_mvtx_hits++;
                break;
              case TrkrDefs::inttId:
                n_intt_hits++;
                break;
              case TrkrDefs::tpcId:
                n_tpc_hits++;
                break;
            }
          }


          /*m_chi2ndf[input]->Fill(pflow->get_track()->get_chisq()/pflow->get_track()->get_ndf());
          m_mvtxcl[input]->Fill(m_nmaps);
          m_inttcl[input]->Fill(m_nintt);
          m_mtpccl[input]->Fill(m_ntpc);*/
          track_properties.pflowtype = pflow->get_type();
          track_properties.vtx_id = id;
          track_properties.e = pflow->get_e();
          track_properties.eta = pflow->get_eta();
          track_properties.phi = pflow->get_phi();
          track_properties.pt = pflow->get_pt();
          track_properties.DCA_xy = DCA_xy;
          track_properties.DCA_xy_unc = DCA_xy_unc;
          track_properties.sDCA_xy = sign*std::abs(DCA_xy/DCA_xy_unc);
          track_properties.n_mvtx = n_mvtx_hits;
          track_properties.n_intt = n_intt_hits;
          track_properties.n_tpc = n_tpc_hits;
          track_properties.chisq = pflow->get_track()->get_chisq();
          track_properties.ndf = pflow->get_track()->get_ndf();

          recojet.chConstituents.push_back(track_properties);
        } // end if(pflow->get_track())
        else{ //neutral tracl
          neConstituent neutral_properties;
          neutral_properties.pflowtype = pflow->get_type();
          neutral_properties.e = pflow->get_e();
          neutral_properties.eta = pflow->get_eta();
          neutral_properties.phi = pflow->get_phi();
          recojet.neConstituents.push_back(neutral_properties);
        }
      } // end for (const auto& comp : jet->get_comp_vec())

      recojet.id = jet->get_id();
      recojet.area = jet->get_property(recojet_area_index);
      recojet.num_Constituents = static_cast<int>(jet->get_comp_vec().size());
      recojet.num_ChConstituents = nChtracks;
      recojet.px = jet->get_px();
      recojet.py = jet->get_py();
      recojet.pz = jet->get_pz();
      recojet.pt = jet->get_pt();
      recojet.eta = jet->get_eta();
      recojet.phi = jet->get_phi();
      recojet.m = jet->get_mass();
      recojet.e = jet->get_e();

      m_container[input]->recojets.push_back(recojet);
    } // end for (Jet* jet : *jets)
    m_container[input]->reco_jet_n = static_cast<int>(nrecojet);


  
    //get truth jets
    if(m_doTruthJets){
      int ntruthjet = 0;
      TruthJets mtruthjet;
      Jet::PROPERTY truthjet_area_index = jetsMC->property_index(Jet::PROPERTY::prop_area);

      //Jet::PROPERTY recojet_area_index = jetsMC->property_index(Jet::PROPERTY::prop_area);
      
      //truth jet loop
      //for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter){   
      for (Jet* truthjet : *jetsMC){
        //Jet* truthjet = iter->second;
        if(truthjet->get_pt() < m_ptRangeTruth.first || truthjet->get_pt() > m_ptRangeTruth.second) continue;
        if (not (std::abs(truthjet->get_eta()) <= 1.1 - (doFiducial?jetR.at(input):0))) continue;
        if(truthjet->get_pt() < 10.0) continue;
        ntruthjet++;


        mtruthjet.id = truthjet->get_id();
        mtruthjet.area = truthjet->get_property(truthjet_area_index);
        mtruthjet.num_Constituents = truthjet->size_comp();
        mtruthjet.px = truthjet->get_px();
        mtruthjet.py = truthjet->get_py();
        mtruthjet.pz = truthjet->get_pz();
        mtruthjet.pt = truthjet->get_pt();
        mtruthjet.eta = truthjet->get_eta();
        mtruthjet.phi = truthjet->get_phi();
        mtruthjet.m = truthjet->get_mass();
        mtruthjet.e = truthjet->get_e();


        //jet-tarck PDG identification
        TString taggedPDGIDs[] = {"B-Meson","CharmedMeson","CharmedBaryon", "B-Baryon"};
        int forbiddenPDGIDs[] = {1,2,3,4,5,6,7,8,21,22}; 
        int nChTrackstruth = 0;

        //loop over jet constituents
        //for(auto citer = truthjet->begin_comp(); citer != truthjet->end_comp(); ++citer){
        for (const auto& comp : truthjet->get_comp_vec()){
          //PHG4Particle *g4part = truthinfo->GetPrimaryParticle(citer->second);
          PHG4Particle *g4part = truthinfo->GetPrimaryParticle(comp.second);       
                HepMC::GenParticle* constituent = hepMCGenEvent->barcode_to_particle(g4part->get_barcode());
          if(std::abs((TDatabasePDG::Instance()->GetParticle((constituent)->pdg_id()))->Charge()) > 10e-4) nChTrackstruth++;

          //mother track tagged, stop search
          if (std::find(std::begin(taggedPDGIDs), std::end(taggedPDGIDs), TString((TDatabasePDG::Instance()->GetParticle((constituent)->pdg_id()))->ParticleClass())) != std::end(taggedPDGIDs)){
              mtruthjet.constituents_PDG_ID.push_back((constituent)->pdg_id());
              continue;
          }

          //not yet tagged, begin history search
            while (!constituent->is_beam()){
            bool breakOut = false;
            bool taggedHF = false;
            for (HepMC::GenVertex::particle_iterator mother = constituent->production_vertex()->particles_begin(HepMC::parents); mother != constituent->production_vertex()->particles_end(HepMC::parents); ++mother){
              //found HF in parent tree
              if (std::find(std::begin(taggedPDGIDs), std::end(taggedPDGIDs), TString((TDatabasePDG::Instance()->GetParticle((*mother)->pdg_id()))->ParticleClass())) != std::end(taggedPDGIDs)){
                taggedHF = true;
                mtruthjet.constituents_PDG_ID.push_back((*mother)->pdg_id());
                break;
              }
              //reached partonic level, break search
              if (std::find(std::begin(forbiddenPDGIDs), std::end(forbiddenPDGIDs), abs((*mother)->pdg_id())) != std::end(forbiddenPDGIDs)){
                breakOut = true;
                break;
              }
              constituent = *mother;
            }
            if (breakOut || taggedHF) break;
          }
        }// end loop over constituents

        mtruthjet.num_ChConstituents = nChTrackstruth ;
        m_container[input]->truthjets.push_back(mtruthjet);
      }// end jet loop
      m_container[input]->truth_jet_n = static_cast<int>(ntruthjet);
    }//end if do truth
    //fill tree
    m_T[input]->Fill();
  }//end input loop

  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::ResetEvent(PHCompositeNode *topNode)
{
  //std::cout << "FullJetFinder::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::EndRun(const int runnumber)
{
  std::cout << "FullJetFinder::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::End(PHCompositeNode *topNode)
{
  std::cout << "FullJetFinder::End - Output to " << m_outputFileName << std::endl;
  
   /*CutSelection cutselection;
    cutselection.jetptmin = 5;
    cutselection.jetptmin = 30;*/

  if(PHTFileServer::get().cd(m_outputFileName)){
    m_stat->Write();
  }

  for(int i = 0 ; i < m_inputs; i++){
    if(m_T[i] ){
      if(PHTFileServer::get().cd(m_outputFileName)){
        TDirectory *cdtof = gDirectory->mkdir(m_TreeNameCollection.at(i).data());
        cdtof->cd();
        //cdtof-><CutSelection>WriteObject(cutselection, "Selection");
        //cutselection.Write("Selection");
        m_T[i]->Write();
      }
    }
  }
  std::cout << "FullJetFinder::End(PHCompositeNode *topNode) This is the End..." << std::endl;
  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
int FullJetFinder::Reset(PHCompositeNode *topNode)
{
 std::cout << "FullJetFinder::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________..
void FullJetFinder::Print(const std::string &what) const
{
  std::cout << "FullJetFinder::Print(const std::string &what) const Printing info for " << what << std::endl;
}

void FullJetFinder::GetDistanceFromVertexXY(SvtxTrack *m_dst_track, GlobalVertex *m_dst_vertex,float &val, float &err){
 
  KFParticle::SetField(1.4e1);

  KFParticle kfp_particle;

  float f_trackParameters[6] = {m_dst_track->get_x(),
                                m_dst_track->get_y(),
                                m_dst_track->get_z(),
                                m_dst_track->get_px(),
                                m_dst_track->get_py(),
                                m_dst_track->get_pz()};

  float f_trackCovariance[21];
  unsigned int iterate = 0;
  for (unsigned int i = 0; i < 6; ++i)
    for (unsigned int j = 0; j <= i; ++j)
    {
      f_trackCovariance[iterate] = m_dst_track->get_error(i, j);
      ++iterate;
    }

  kfp_particle.Create(f_trackParameters, f_trackCovariance, (Int_t) m_dst_track->get_charge(), -1);
  kfp_particle.NDF() = m_dst_track->get_ndf();
  kfp_particle.Chi2() = m_dst_track->get_chisq();
  kfp_particle.SetId(m_dst_track->get_id());

  float f_vertexParameters[6] = {m_dst_vertex->get_x(),
                                 m_dst_vertex->get_y(),
                                 m_dst_vertex->get_z(), 0, 0, 0};

  float f_vertexCovariance[21];
  unsigned int iterate2 = 0;
  for (unsigned int i = 0; i < 3; ++i)
    for (unsigned int j = 0; j <= i; ++j)
    {
      f_vertexCovariance[iterate2] = m_dst_vertex->get_error(i, j);
      ++iterate2;
    }

  KFParticle kfp_vertex;
  kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
  kfp_vertex.NDF() = m_dst_vertex->get_ndof();
  kfp_vertex.Chi2() = m_dst_vertex->get_chisq();
  kfp_particle.GetDistanceFromVertexXY(kfp_vertex,val,err);
}