diff --git a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc index 4797d957ec..449f16d732 100644 --- a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc +++ b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc @@ -215,6 +215,14 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) m_true_mother_p = std::sqrt(std::pow(mother->momentum().px(), 2) + std::pow(mother->momentum().py(), 2) + std::pow(mother->momentum().pz(), 2)); // Must have an old HepMC build, no mag function m_true_mother_eta = mother->momentum().eta(); m_true_mother_phi = mother->momentum().phi(); + if (mother->momentum().e() > abs(mother->momentum().pz())) + { + m_true_mother_rapidity = 0.5 * log((mother->momentum().e() + mother->momentum().pz())/(mother->momentum().e() - mother->momentum().pz())); + } + else + { + m_true_mother_rapidity = -999.; + } HepMC::GenVertex *thisVtx = mother->production_vertex(); m_primary_vtx_x = thisVtx->point3d().x(); @@ -241,6 +249,7 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) std::abs(decay[i].second)) != std::end(trackableParticles)) { ++index; + if (theEvent && decay[i].first.second > -1) { HepMC::GenParticle *daughterHepMC = theEvent->barcode_to_particle(decay[i].first.second); @@ -297,8 +306,75 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) continue; } - if (motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second) + if (motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second && m_nDaughters == 2) { + if (Verbosity() >= VERBOSITY_MORE || true) // fix later + { + daughterG4->identify(); + } + + m_is_primary = m_truthInfo->is_sPHENIX_primary(motherG4); + + CLHEP::Hep3Vector *mother3Vector = new CLHEP::Hep3Vector(motherG4->get_px(), motherG4->get_py(), motherG4->get_pz()); + motherTrueLV->setVectM((*mother3Vector), getParticleMass(decay[0].second)); + m_true_mother_pT = motherTrueLV->perp(); + m_true_mother_p = mother3Vector->mag(); + m_true_mother_eta = motherTrueLV->pseudoRapidity(); + m_true_mother_phi = motherTrueLV->phi(); + m_true_mother_rapidity = motherTrueLV->rapidity(); + + PHG4VtxPoint *thisVtx = m_truthInfo->GetVtx(motherG4->get_vtx_id()); + m_primary_vtx_x = thisVtx->get_x(); + m_primary_vtx_y = thisVtx->get_y(); + m_primary_vtx_z = thisVtx->get_z(); + + daughterTrueLV->setVectM(CLHEP::Hep3Vector(daughterG4->get_px(), daughterG4->get_py(), daughterG4->get_pz()), getParticleMass(decay[i].second)); + daughterSumTrueLV += *daughterTrueLV; + + // Now get the decay vertex position + thisVtx = m_truthInfo->GetVtx(daughterG4->get_vtx_id()); + m_secondary_vtx_x = thisVtx->get_x(); + m_secondary_vtx_y = thisVtx->get_y(); + m_secondary_vtx_z = thisVtx->get_z(); + + m_true_track_PID[index] = daughterG4->get_pid(); + truth_ID = daughterG4->get_track_id(); + + delete mother3Vector; + } + else if (m_nDaughters == 3) + { + if (i != 4 && motherG4->get_pid() == decay[3].second && motherG4->get_barcode() == decay[3].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second) + { + PHG4Particle *motherG4_temp = nullptr; + if (motherG4->get_parent_id() != 0) + { + motherG4_temp = m_truthInfo->GetParticle(motherG4->get_parent_id()); + } + else + { + continue; + } + + if (motherG4_temp->get_pid() == decay[0].second && motherG4_temp->get_barcode() == decay[0].first.second) + { + motherG4 = motherG4_temp; + } + else + { + continue; + } + } + else if (i != 4) + { + continue; + } + + if (i==4 && !(motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second)) + { + continue; + } + if (Verbosity() >= VERBOSITY_MORE) { daughterG4->identify(); @@ -312,6 +388,7 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) m_true_mother_p = mother3Vector->mag(); m_true_mother_eta = motherTrueLV->pseudoRapidity(); m_true_mother_phi = motherTrueLV->phi(); + m_true_mother_rapidity = motherTrueLV->rapidity(); PHG4VtxPoint *thisVtx = m_truthInfo->GetVtx(motherG4->get_vtx_id()); m_primary_vtx_x = thisVtx->get_x(); @@ -327,16 +404,18 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) m_secondary_vtx_y = thisVtx->get_y(); m_secondary_vtx_z = thisVtx->get_z(); - m_true_track_PID[index] = daughterG4->get_pid(); + m_true_track_PID[index] = daughterG4->get_pid(); truth_ID = daughterG4->get_track_id(); delete mother3Vector; + break; } } } m_true_track_pT[index] = (float) daughterTrueLV->perp(); m_true_track_eta[index] = (float) daughterTrueLV->pseudoRapidity(); + m_true_track_rapidity[index] = (float) daughterTrueLV->rapidity(); m_true_track_phi[index] = (float) daughterTrueLV->phi(); m_min_true_track_pT = std::min(m_true_track_pT[index], m_min_true_track_pT); m_max_true_track_pT = std::max(m_true_track_pT[index], m_max_true_track_pT); @@ -480,6 +559,7 @@ void HFTrackEfficiency::initializeBranches() m_tree->Branch("reco_mother_pT", &m_reco_mother_pT, "reco_mother_pT/F"); m_tree->Branch("true_mother_p", &m_true_mother_p, "true_mother_p/F"); m_tree->Branch("true_mother_eta", &m_true_mother_eta, "true_mother_eta/F"); + m_tree->Branch("true_mother_rapidity", &m_true_mother_rapidity, "true_mother_rapidity/F"); m_tree->Branch("true_mother_phi", &m_true_mother_phi, "true_mother_phi/F"); m_tree->Branch("min_true_track_pT", &m_min_true_track_pT, "min_true_track_pT/F"); m_tree->Branch("min_reco_track_pT", &m_min_reco_track_pT, "min_reco_track_pT/F"); @@ -495,6 +575,7 @@ void HFTrackEfficiency::initializeBranches() m_tree->Branch("reco_" + TString(daughter_number) + "_pT", &m_reco_track_pT[iTrack], "reco_" + TString(daughter_number) + "_pT/F"); m_tree->Branch("true_" + TString(daughter_number) + "_eta", &m_true_track_eta[iTrack], "true_" + TString(daughter_number) + "_eta/F"); m_tree->Branch("reco_" + TString(daughter_number) + "_eta", &m_reco_track_eta[iTrack], "reco_" + TString(daughter_number) + "_eta/F"); + m_tree->Branch("true_" + TString(daughter_number) + "_rapidity", &m_true_track_rapidity[iTrack], "true_" + TString(daughter_number) + "_rapidity/F"); m_tree->Branch("true_" + TString(daughter_number) + "_phi", &m_true_track_phi[iTrack], "true_" + TString(daughter_number) + "_phi/F"); m_tree->Branch("reco_" + TString(daughter_number) + "_phi", &m_reco_track_phi[iTrack], "reco_" + TString(daughter_number) + "_phi/F"); m_tree->Branch("true_" + TString(daughter_number) + "_PID", &m_true_track_PID[iTrack], "true_" + TString(daughter_number) + "_PID/F"); @@ -521,6 +602,7 @@ void HFTrackEfficiency::resetBranches() m_reco_mother_pT = std::numeric_limits::quiet_NaN(); m_true_mother_p = std::numeric_limits::quiet_NaN(); m_true_mother_eta = std::numeric_limits::quiet_NaN(); + m_true_mother_rapidity = std::numeric_limits::quiet_NaN(); m_min_true_track_pT = std::numeric_limits::max(); m_min_reco_track_pT = std::numeric_limits::max(); m_max_true_track_pT = -1 * std::numeric_limits::max(); @@ -533,6 +615,7 @@ void HFTrackEfficiency::resetBranches() m_reco_track_pT[iTrack] = std::numeric_limits::quiet_NaN(); m_true_track_eta[iTrack] = std::numeric_limits::quiet_NaN(); m_reco_track_eta[iTrack] = std::numeric_limits::quiet_NaN(); + m_true_track_rapidity[iTrack] = std::numeric_limits::quiet_NaN(); m_true_track_phi[iTrack] = std::numeric_limits::quiet_NaN(); m_reco_track_phi[iTrack] = std::numeric_limits::quiet_NaN(); m_true_track_PID[iTrack] = std::numeric_limits::quiet_NaN(); diff --git a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.h b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.h index f667c3bc27..c5645240b0 100644 --- a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.h +++ b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.h @@ -96,6 +96,7 @@ class HFTrackEfficiency : public SubsysReco float m_reco_mother_pT{std::numeric_limits::quiet_NaN()}; float m_true_mother_p{std::numeric_limits::quiet_NaN()}; float m_true_mother_eta{std::numeric_limits::quiet_NaN()}; + float m_true_mother_rapidity{std::numeric_limits::quiet_NaN()}; float m_true_mother_phi{std::numeric_limits::quiet_NaN()}; float m_min_true_track_pT{std::numeric_limits::max()}; float m_min_reco_track_pT{std::numeric_limits::max()}; @@ -107,6 +108,7 @@ class HFTrackEfficiency : public SubsysReco float m_reco_track_pT[m_maxTracks]{std::numeric_limits::quiet_NaN()}; float m_true_track_eta[m_maxTracks]{std::numeric_limits::quiet_NaN()}; float m_reco_track_eta[m_maxTracks]{std::numeric_limits::quiet_NaN()}; + float m_true_track_rapidity[m_maxTracks]{std::numeric_limits::quiet_NaN()}; float m_true_track_phi[m_maxTracks]{std::numeric_limits::quiet_NaN()}; float m_reco_track_phi[m_maxTracks]{std::numeric_limits::quiet_NaN()}; float m_true_track_PID[m_maxTracks]{std::numeric_limits::quiet_NaN()}; diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc index bc9f0b4c3a..392e0c4267 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc @@ -205,6 +205,8 @@ void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode) m_tree->Branch(TString(daughter_number) + "_IPchi2", &m_calculated_daughter_ipchi2[i], TString(daughter_number) + "_IPchi2/F"); m_tree->Branch(TString(daughter_number) + "_IPErr", &m_calculated_daughter_ip_err[i], TString(daughter_number) + "_IPErr/F"); m_tree->Branch(TString(daughter_number) + "_IP_xy", &m_calculated_daughter_ip_xy[i], TString(daughter_number) + "_IP_xy/F"); + m_tree->Branch(TString(daughter_number) + "_DCA_sig", &m_calculated_daughter_PV_dca_sig[i], TString(daughter_number) + "_DCA_sig/F"); + m_tree->Branch(TString(daughter_number) + "_DCA_sig_xy", &m_calculated_daughter_PV_dca_xy_sig[i], TString(daughter_number) + "_DCA_sig_xy/F"); } if (m_get_all_PVs) { @@ -270,6 +272,14 @@ void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode) std::string dca_branch_name_xy = dca_branch_name + "_xy"; std::string dca_leaf_name_xy = dca_branch_name_xy + "/F"; m_tree->Branch(dca_branch_name_xy.c_str(), &m_daughter_dca_xy[iter], dca_leaf_name_xy.c_str()); + + std::string dca_sig_branch_name = "track_" + std::to_string(i + 1) + "_track_" + std::to_string(j + 1) + "_DCA_sig"; + std::string dca_sig_leaf_name = dca_sig_branch_name + "/F"; + m_tree->Branch(dca_sig_branch_name.c_str(), &m_daughter_dca_sig[iter], dca_sig_leaf_name.c_str()); + + std::string dca_sig_branch_name_xy = dca_sig_branch_name + "_xy"; + std::string dca_sig_leaf_name_xy = dca_sig_branch_name_xy + "/F"; + m_tree->Branch(dca_sig_branch_name_xy.c_str(), &m_daughter_dca_sig_xy[iter], dca_sig_leaf_name_xy.c_str()); ++iter; } @@ -496,6 +506,8 @@ void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode, m_calculated_daughter_ipchi2[i] = daughterArray[i].GetDeviationFromVertex(vertex_fillbranch); m_calculated_daughter_ip_err[i] = m_calculated_daughter_ip[i] / std::sqrt(m_calculated_daughter_ipchi2[i]); m_calculated_daughter_ip_xy[i] = daughterArray[i].GetDistanceFromVertexXY(vertex_fillbranch); + m_calculated_daughter_PV_dca_sig[i] = daughterArray[i].GetDeviationFromVertex(vertex_fillbranch); + m_calculated_daughter_PV_dca_xy_sig[i] = daughterArray[i].GetDeviationFromVertexXY(vertex_fillbranch); } m_calculated_daughter_x[i] = daughterArray[i].GetX(); m_calculated_daughter_y[i] = daughterArray[i].GetY(); @@ -596,6 +608,8 @@ void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode, { m_daughter_dca[iter] = daughterArray[i].GetDistanceFromParticle(daughterArray[j]); m_daughter_dca_xy[iter] = daughterArray[i].GetDistanceFromParticleXY(daughterArray[j]); + m_daughter_dca_sig[iter] = daughterArray[i].GetDeviationFromParticle(daughterArray[j]); + m_daughter_dca_sig_xy[iter] = daughterArray[i].GetDeviationFromParticleXY(daughterArray[j]); ++iter; } } diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h index 45256fc1cc..fd01a27a2f 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h @@ -174,6 +174,8 @@ class KFParticle_nTuple : public KFParticle_truthAndDetTools, public KFParticle_ float m_calculated_daughter_mass[max_tracks]{0}; float m_calculated_daughter_ip[max_tracks]{0}; float m_calculated_daughter_ip_xy[max_tracks]{0}; + float m_calculated_daughter_PV_dca_sig[max_tracks]{0}; + float m_calculated_daughter_PV_dca_xy_sig[max_tracks]{0}; float m_calculated_daughter_ipchi2[max_tracks]{0}; float m_calculated_daughter_ip_err[max_tracks]{0}; float m_calculated_daughter_x[max_tracks]{0}; @@ -207,6 +209,8 @@ class KFParticle_nTuple : public KFParticle_truthAndDetTools, public KFParticle_ float m_daughter_dca[99]{0}; float m_daughter_dca_xy[99]{0}; + float m_daughter_dca_sig[99]{0}; + float m_daughter_dca_sig_xy[99]{0}; float m_calculated_vertex_x{-1}; float m_calculated_vertex_y{-1};