Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:48 

0001 #include "KFParticle_nTuple.h"
0002 
0003 #include "KFParticle_Tools.h"
0004 
0005 #include <ffaobjects/EventHeader.h>
0006 #include <ffarawobjects/Gl1Packet.h>
0007 
0008 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0009 #include <phool/getClass.h>
0010 #include <trackbase_historic/SvtxTrack.h>   
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 
0013 #include <KFParticle.h>
0014 #include <KFVertex.h>
0015 
0016 #include <Rtypes.h>
0017 #include <TString.h>  // for TString, operator+
0018 #include <TTree.h>
0019 
0020 #include <algorithm>  // for max
0021 #include <cmath>
0022 #include <cstdlib>  // for abs, size_t
0023 #include <map>      // for map, _Rb_tree_iterator, map<>:...
0024 
0025 class PHCompositeNode;
0026 class PHNode;
0027 
0028 /// Create necessary objects
0029 KFParticle_Tools kfpTupleTools;
0030 float TempError;
0031 
0032 void KFParticle_nTuple::initializeVariables()
0033 {
0034   // m_calculated_mother_cov = -99;
0035 }
0036 
0037 void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode)
0038 {
0039   //kfpTupleTools.init_dEdx_fits(); //Cant do this! Two trees open at once!
0040 
0041   delete m_tree;
0042   m_tree = new TTree("DecayTree", "DecayTree");
0043   m_tree->OptimizeBaskets();
0044   m_tree->SetAutoSave(-5e6);  // Save the output file every 5MB
0045 
0046   std::string mother_name;
0047   if (m_mother_name.empty())
0048   {
0049     mother_name = "mother";
0050   }
0051   else
0052   {
0053     mother_name = m_mother_name;
0054   }
0055 
0056   size_t pos;
0057   std::string undrscr = "_";
0058   std::string nothing = "";
0059   std::map<std::string, std::string> forbiddenStrings;
0060   forbiddenStrings["/"] = undrscr;
0061   forbiddenStrings["("] = undrscr;
0062   forbiddenStrings[")"] = nothing;
0063   forbiddenStrings["+"] = "plus";
0064   forbiddenStrings["-"] = "minus";
0065   forbiddenStrings["*"] = "star";
0066   for (auto const& [badString, goodString] : forbiddenStrings)
0067   {
0068     while ((pos = mother_name.find(badString)) != std::string::npos)
0069     {
0070       mother_name.replace(pos, 1, goodString);
0071     }
0072   }
0073 
0074   m_tree->Branch(TString(mother_name) + "_mass", &m_calculated_mother_mass, TString(mother_name) + "_mass/F");
0075   m_tree->Branch(TString(mother_name) + "_massErr", &m_calculated_mother_mass_err, TString(mother_name) + "_massErr/F");
0076   if (m_constrain_to_vertex_nTuple)
0077   {
0078     m_tree->Branch(TString(mother_name) + "_decayTime", &m_calculated_mother_decaytime, TString(mother_name) + "_decayTime/F");
0079     m_tree->Branch(TString(mother_name) + "_decayTimeErr", &m_calculated_mother_decaytime_err, TString(mother_name) + "_decayTimeErr/F");
0080     m_tree->Branch(TString(mother_name) + "_decayLength", &m_calculated_mother_decaylength, TString(mother_name) + "_decayLength/F");
0081     m_tree->Branch(TString(mother_name) + "_decayLengthErr", &m_calculated_mother_decaylength_err, TString(mother_name) + "_decayLengthErr/F");
0082     m_tree->Branch(TString(mother_name) + "_decayLength_xy", &m_calculated_mother_decaylength_xy, TString(mother_name) + "_decayLength_xy/F");
0083     m_tree->Branch(TString(mother_name) + "_decayLengthErr_xy", &m_calculated_mother_decaylength_xy_err, TString(mother_name) + "_decayLengthErr_xy/F");
0084     m_tree->Branch(TString(mother_name) + "_DIRA", &m_calculated_mother_dira, TString(mother_name) + "_DIRA/F");
0085     m_tree->Branch(TString(mother_name) + "_DIRA_xy", &m_calculated_mother_dira_xy, TString(mother_name) + "_DIRA_xy/F");
0086     m_tree->Branch(TString(mother_name) + "_FDchi2", &m_calculated_mother_fdchi2, TString(mother_name) + "_FDchi2/F");
0087     m_tree->Branch(TString(mother_name) + "_IP", &m_calculated_mother_ip, TString(mother_name) + "_IP/F");
0088     m_tree->Branch(TString(mother_name) + "_IPchi2", &m_calculated_mother_ipchi2, TString(mother_name) + "_IPchi2/F");
0089     m_tree->Branch(TString(mother_name) + "_IPErr", &m_calculated_mother_ip_err, TString(mother_name) + "_IPErr/F");
0090     m_tree->Branch(TString(mother_name) + "_IP_xy", &m_calculated_mother_ip_xy, TString(mother_name) + "_IP_xy/F");
0091   }
0092   if (m_get_all_PVs)
0093   {
0094     m_tree->Branch(TString(mother_name) + "_IP_allPV", &allPV_mother_IP);
0095     m_tree->Branch(TString(mother_name) + "_IPchi2_allPV", &allPV_mother_IPchi2);
0096   }
0097   m_tree->Branch(TString(mother_name) + "_x", &m_calculated_mother_x, TString(mother_name) + "_x/F");
0098   m_tree->Branch(TString(mother_name) + "_y", &m_calculated_mother_y, TString(mother_name) + "_y/F");
0099   m_tree->Branch(TString(mother_name) + "_z", &m_calculated_mother_z, TString(mother_name) + "_z/F");
0100   m_tree->Branch(TString(mother_name) + "_px", &m_calculated_mother_px, TString(mother_name) + "_px/F");
0101   m_tree->Branch(TString(mother_name) + "_py", &m_calculated_mother_py, TString(mother_name) + "_py/F");
0102   m_tree->Branch(TString(mother_name) + "_pz", &m_calculated_mother_pz, TString(mother_name) + "_pz/F");
0103   m_tree->Branch(TString(mother_name) + "_pE", &m_calculated_mother_pe, TString(mother_name) + "_pE/F");
0104   m_tree->Branch(TString(mother_name) + "_p", &m_calculated_mother_p, TString(mother_name) + "_p/F");
0105   m_tree->Branch(TString(mother_name) + "_pErr", &m_calculated_mother_p_err, TString(mother_name) + "_pErr/F");
0106   m_tree->Branch(TString(mother_name) + "_pT", &m_calculated_mother_pt, TString(mother_name) + "_pT/F");
0107   m_tree->Branch(TString(mother_name) + "_pTErr", &m_calculated_mother_pt_err, TString(mother_name) + "_pTErr/F");
0108   m_tree->Branch(TString(mother_name) + "_charge", &m_calculated_mother_q, TString(mother_name) + "_charge/B");
0109   m_tree->Branch(TString(mother_name) + "_pseudorapidity", &m_calculated_mother_eta, TString(mother_name) + "_pseudorapidity/F");
0110   m_tree->Branch(TString(mother_name) + "_rapidity", &m_calculated_mother_rapidity, TString(mother_name) + "_rapidity/F");
0111   m_tree->Branch(TString(mother_name) + "_theta", &m_calculated_mother_theta, TString(mother_name) + "_theta/F");
0112   m_tree->Branch(TString(mother_name) + "_phi", &m_calculated_mother_phi, TString(mother_name) + "_phi/F");
0113   m_tree->Branch(TString(mother_name) + "_vertex_volume", &m_calculated_mother_v, TString(mother_name) + "_vertex_volume/F");
0114   m_tree->Branch(TString(mother_name) + "_chi2", &m_calculated_mother_chi2, TString(mother_name) + "_chi2/F");
0115   m_tree->Branch(TString(mother_name) + "_nDoF", &m_calculated_mother_ndof, TString(mother_name) + "_nDoF/i");
0116   m_tree->Branch(TString(mother_name) + "_PDG_ID", &m_calculated_mother_pdgID, TString(mother_name) + "_PDG_ID/I");
0117   m_tree->Branch(TString(mother_name) + "_Covariance", &m_calculated_mother_cov, TString(mother_name) + "_Covariance[21]/F", 21);
0118 
0119   std::vector<std::string> intermediateNameMapping;  // What intermediate is associate to what track
0120   if (m_has_intermediates_nTuple)
0121   {
0122     for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
0123     {
0124       std::string intermediate_name = m_intermediate_name_ntuple[i];
0125 
0126       // Note, TBranch will not allow the leaf to contain a forward slash as it is used to define the branch type. Causes problems with J/psi
0127       for (auto const& [badString, goodString] : forbiddenStrings)
0128       {
0129         while ((pos = intermediate_name.find(badString)) != std::string::npos)
0130         {
0131           intermediate_name.replace(pos, 1, goodString);
0132         }
0133       }
0134 
0135       m_tree->Branch(TString(intermediate_name) + "_mass", &m_calculated_intermediate_mass[i], TString(intermediate_name) + "_mass/F");
0136       m_tree->Branch(TString(intermediate_name) + "_massErr", &m_calculated_intermediate_mass_err[i], TString(intermediate_name) + "_massErr/F");
0137       m_tree->Branch(TString(intermediate_name) + "_decayTime", &m_calculated_intermediate_decaytime[i], TString(intermediate_name) + "_decayTime/F");
0138       m_tree->Branch(TString(intermediate_name) + "_decayTimeErr", &m_calculated_intermediate_decaytime_err[i], TString(intermediate_name) + "_decayTimeErr/F");
0139       m_tree->Branch(TString(intermediate_name) + "_decayLength", &m_calculated_intermediate_decaylength[i], TString(intermediate_name) + "_decayLength/F");
0140       m_tree->Branch(TString(intermediate_name) + "_decayLengthErr", &m_calculated_intermediate_decaylength_err[i], TString(intermediate_name) + "_decayLengthErr/F");
0141       m_tree->Branch(TString(intermediate_name) + "_decayLength_xy", &m_calculated_intermediate_decaylength_xy[i], TString(intermediate_name) + "_decayLength_xy/F");
0142       m_tree->Branch(TString(intermediate_name) + "_decayLengthErr_xy", &m_calculated_intermediate_decaylength_xy_err[i], TString(intermediate_name) + "_decayLengthErr_xy/F");
0143       m_tree->Branch(TString(intermediate_name) + "_DIRA", &m_calculated_intermediate_dira[i], TString(intermediate_name) + "_DIRA/F");
0144       m_tree->Branch(TString(intermediate_name) + "_FDchi2", &m_calculated_intermediate_fdchi2[i], TString(intermediate_name) + "_FDchi2/F");
0145       if (m_constrain_to_vertex_nTuple)
0146       {
0147         m_tree->Branch(TString(intermediate_name) + "_IP", &m_calculated_intermediate_ip[i], TString(intermediate_name) + "_IP/F");
0148         m_tree->Branch(TString(intermediate_name) + "_IPchi2", &m_calculated_intermediate_ipchi2[i], TString(intermediate_name) + "_IPchi2/F");
0149         m_tree->Branch(TString(intermediate_name) + "_IPErr", &m_calculated_intermediate_ip_err[i], TString(intermediate_name) + "_IPErr/F");
0150         m_tree->Branch(TString(intermediate_name) + "_IP_xy", &m_calculated_intermediate_ip_xy[i], TString(intermediate_name) + "_IP_xy/F");
0151       }
0152       if (m_get_all_PVs)
0153       {
0154         m_tree->Branch(TString(intermediate_name) + "_IP_allPV", &allPV_intermediates_IP[i]);
0155         m_tree->Branch(TString(intermediate_name) + "_IPchi2_allPV", &allPV_intermediates_IPchi2[i]);
0156       }
0157       m_tree->Branch(TString(intermediate_name) + "_x", &m_calculated_intermediate_x[i], TString(intermediate_name) + "_x/F");
0158       m_tree->Branch(TString(intermediate_name) + "_y", &m_calculated_intermediate_y[i], TString(intermediate_name) + "_y/F");
0159       m_tree->Branch(TString(intermediate_name) + "_z", &m_calculated_intermediate_z[i], TString(intermediate_name) + "_z/F");
0160       m_tree->Branch(TString(intermediate_name) + "_px", &m_calculated_intermediate_px[i], TString(intermediate_name) + "_px/F");
0161       m_tree->Branch(TString(intermediate_name) + "_py", &m_calculated_intermediate_py[i], TString(intermediate_name) + "_py/F");
0162       m_tree->Branch(TString(intermediate_name) + "_pz", &m_calculated_intermediate_pz[i], TString(intermediate_name) + "_pz/F");
0163       m_tree->Branch(TString(intermediate_name) + "_pE", &m_calculated_intermediate_pe[i], TString(intermediate_name) + "_pE/F");
0164       m_tree->Branch(TString(intermediate_name) + "_p", &m_calculated_intermediate_p[i], TString(intermediate_name) + "_p/F");
0165       m_tree->Branch(TString(intermediate_name) + "_pErr", &m_calculated_intermediate_p_err[i], TString(intermediate_name) + "_pErr/F");
0166       m_tree->Branch(TString(intermediate_name) + "_pT", &m_calculated_intermediate_pt[i], TString(intermediate_name) + "_pT/F");
0167       m_tree->Branch(TString(intermediate_name) + "_pTErr", &m_calculated_intermediate_pt_err[i], TString(intermediate_name) + "_pTErr/F");
0168       m_tree->Branch(TString(intermediate_name) + "_charge", &m_calculated_intermediate_q[i], TString(intermediate_name) + "_charge/B");
0169       m_tree->Branch(TString(intermediate_name) + "_pseudorapidity", &m_calculated_intermediate_eta[i], TString(intermediate_name) + "_pseudorapidity/F");
0170       m_tree->Branch(TString(intermediate_name) + "_rapidity", &m_calculated_intermediate_rapidity[i], TString(intermediate_name) + "_rapidity/F");
0171       m_tree->Branch(TString(intermediate_name) + "_theta", &m_calculated_intermediate_theta[i], TString(intermediate_name) + "_theta/F");
0172       m_tree->Branch(TString(intermediate_name) + "_phi", &m_calculated_intermediate_phi[i], TString(intermediate_name) + "_phi/F");
0173       m_tree->Branch(TString(intermediate_name) + "_vertex_volume", &m_calculated_intermediate_v[i], TString(intermediate_name) + "_vertex_volume/F");
0174       m_tree->Branch(TString(intermediate_name) + "_chi2", &m_calculated_intermediate_chi2[i], TString(intermediate_name) + "_chi2/F");
0175       m_tree->Branch(TString(intermediate_name) + "_nDoF", &m_calculated_intermediate_ndof[i], TString(intermediate_name) + "_nDoF/i");
0176       m_tree->Branch(TString(intermediate_name) + "_PDG_ID", &m_calculated_intermediate_pdgID[i], TString(intermediate_name) + "_PDG_ID/I");
0177       m_tree->Branch(TString(intermediate_name) + "_Covariance", &m_calculated_intermediate_cov[i], TString(intermediate_name) + "_Covariance[21]/F", 21);
0178 
0179       for (int j = 0; j < m_num_tracks_from_intermediate_nTuple[i]; ++j)
0180       {
0181         intermediateNameMapping.push_back(intermediate_name + "_");
0182       }
0183     }
0184   }
0185 
0186   int num_intermediate_tracks = 0;
0187   for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
0188   {
0189     num_intermediate_tracks += m_num_tracks_from_intermediate_nTuple[i];
0190   }
0191 
0192   for (int i = 0; i < m_num_tracks_nTuple; ++i)
0193   {
0194     std::string daughter_number = "track_" + std::to_string(i + 1);
0195 
0196     if (m_has_intermediates_nTuple && i < num_intermediate_tracks)
0197     {
0198       daughter_number.insert(0, intermediateNameMapping[i]);
0199     }
0200 
0201     m_tree->Branch(TString(daughter_number) + "_mass", &m_calculated_daughter_mass[i], TString(daughter_number) + "_mass/F");
0202     if (m_constrain_to_vertex_nTuple)
0203     {
0204       m_tree->Branch(TString(daughter_number) + "_IP", &m_calculated_daughter_ip[i], TString(daughter_number) + "_IP/F");
0205       m_tree->Branch(TString(daughter_number) + "_IPchi2", &m_calculated_daughter_ipchi2[i], TString(daughter_number) + "_IPchi2/F");
0206       m_tree->Branch(TString(daughter_number) + "_IPErr", &m_calculated_daughter_ip_err[i], TString(daughter_number) + "_IPErr/F");
0207       m_tree->Branch(TString(daughter_number) + "_IP_xy", &m_calculated_daughter_ip_xy[i], TString(daughter_number) + "_IP_xy/F");
0208     }
0209     if (m_get_all_PVs)
0210     {
0211       m_tree->Branch(TString(daughter_number) + "_IP_allPV", &allPV_daughter_IP[i]);
0212       m_tree->Branch(TString(daughter_number) + "_IPchi2_allPV", &allPV_daughter_IPchi2[i]);
0213     }
0214     m_tree->Branch(TString(daughter_number) + "_x", &m_calculated_daughter_x[i], TString(daughter_number) + "_x/F");
0215     m_tree->Branch(TString(daughter_number) + "_y", &m_calculated_daughter_y[i], TString(daughter_number) + "_y/F");
0216     m_tree->Branch(TString(daughter_number) + "_z", &m_calculated_daughter_z[i], TString(daughter_number) + "_z/F");
0217     m_tree->Branch(TString(daughter_number) + "_px", &m_calculated_daughter_px[i], TString(daughter_number) + "_px/F");
0218     m_tree->Branch(TString(daughter_number) + "_py", &m_calculated_daughter_py[i], TString(daughter_number) + "_py/F");
0219     m_tree->Branch(TString(daughter_number) + "_pz", &m_calculated_daughter_pz[i], TString(daughter_number) + "_pz/F");
0220     m_tree->Branch(TString(daughter_number) + "_pE", &m_calculated_daughter_pe[i], TString(daughter_number) + "_pE/F");
0221     m_tree->Branch(TString(daughter_number) + "_p", &m_calculated_daughter_p[i], TString(daughter_number) + "_p/F");
0222     m_tree->Branch(TString(daughter_number) + "_pErr", &m_calculated_daughter_p_err[i], TString(daughter_number) + "_pErr/F");
0223     m_tree->Branch(TString(daughter_number) + "_pT", &m_calculated_daughter_pt[i], TString(daughter_number) + "_pT/F");
0224     m_tree->Branch(TString(daughter_number) + "_pTErr", &m_calculated_daughter_pt_err[i], TString(daughter_number) + "_pTErr/F");
0225     m_tree->Branch(TString(daughter_number) + "_jT", &m_calculated_daughter_jt[i], TString(daughter_number) + "_jT/F");
0226     m_tree->Branch(TString(daughter_number) + "_charge", &m_calculated_daughter_q[i], TString(daughter_number) + "_charge/B");
0227     m_tree->Branch(TString(daughter_number) + "_bunch_crossing", &m_calculated_daughter_bunch_crossing[i], TString(daughter_number) + "_bunch_crossing/I");
0228     m_tree->Branch(TString(daughter_number) + "_pseudorapidity", &m_calculated_daughter_eta[i], TString(daughter_number) + "_pseudorapidity/F");
0229     m_tree->Branch(TString(daughter_number) + "_rapidity", &m_calculated_daughter_rapidity[i], TString(daughter_number) + "_rapidity/F");
0230     m_tree->Branch(TString(daughter_number) + "_theta", &m_calculated_daughter_theta[i], TString(daughter_number) + "_theta/F");
0231     m_tree->Branch(TString(daughter_number) + "_phi", &m_calculated_daughter_phi[i], TString(daughter_number) + "_phi/F");
0232     m_tree->Branch(TString(daughter_number) + "_chi2", &m_calculated_daughter_chi2[i], TString(daughter_number) + "_chi2/F");
0233     m_tree->Branch(TString(daughter_number) + "_nDoF", &m_calculated_daughter_ndof[i], TString(daughter_number) + "_nDoF/i");
0234     m_tree->Branch(TString(daughter_number) + "_track_ID", &m_calculated_daughter_trid[i], TString(daughter_number) + "_track_ID/I");
0235     m_tree->Branch(TString(daughter_number) + "_PDG_ID", &m_calculated_daughter_pdgID[i], TString(daughter_number) + "_PDG_ID/I");
0236     m_tree->Branch(TString(daughter_number) + "_Covariance", &m_calculated_daughter_cov[i], TString(daughter_number) + "_Covariance[21]/F", 21);
0237     m_tree->Branch(TString(daughter_number) + "_calculated_dEdx", &m_calculated_daughter_dedx[i], TString(daughter_number) + "_calculated_dEdx/F");
0238     //m_tree->Branch(TString(daughter_number) + "_expected_pion_dEdx", &m_calculated_daughter_expected_dedx_pion[i], TString(daughter_number) + "_expected_pion_dEdx/F");
0239     //m_tree->Branch(TString(daughter_number) + "_expected_kaon_dEdx", &m_calculated_daughter_expected_dedx_kaon[i], TString(daughter_number) + "_expected_kaon_dEdx/F");
0240     //m_tree->Branch(TString(daughter_number) + "_expected_proton_dEdx", &m_calculated_daughter_expected_dedx_proton[i], TString(daughter_number) + "_expected_proton_dEdx/F");
0241 
0242     if (m_calo_info)
0243     {
0244       initializeCaloBranches(m_tree, i, daughter_number);
0245     }
0246     if (m_truth_matching)
0247     {
0248       initializeTruthBranches(m_tree, i, daughter_number, m_constrain_to_vertex_nTuple);
0249     }
0250     if (m_detector_info)
0251     {
0252       initializeDetectorBranches(m_tree, i, daughter_number);
0253     }
0254   }
0255 
0256   int iter = 0;
0257   for (int i = 0; i < m_num_tracks_nTuple; ++i)
0258   {
0259     for (int j = 0; j < m_num_tracks_nTuple; ++j)
0260     {
0261       if (i < j)
0262       {
0263         std::string dca_branch_name = "track_" + std::to_string(i + 1) + "_track_" + std::to_string(j + 1) + "_DCA";
0264         std::string dca_leaf_name = dca_branch_name + "/F";
0265         m_tree->Branch(dca_branch_name.c_str(), &m_daughter_dca[iter], dca_leaf_name.c_str());
0266 
0267         std::string dca_branch_name_xy = dca_branch_name + "_xy";
0268         std::string dca_leaf_name_xy = dca_branch_name_xy + "/F";
0269         m_tree->Branch(dca_branch_name_xy.c_str(), &m_daughter_dca_xy[iter], dca_leaf_name_xy.c_str());
0270 
0271         ++iter;
0272       }
0273     }
0274   }
0275 
0276   if (m_constrain_to_vertex_nTuple)
0277   {
0278     m_tree->Branch("primary_vertex_x", &m_calculated_vertex_x, "primary_vertex_x/F");
0279     m_tree->Branch("primary_vertex_y", &m_calculated_vertex_y, "primary_vertex_y/F");
0280     m_tree->Branch("primary_vertex_z", &m_calculated_vertex_z, "primary_vertex_z/F");
0281     m_tree->Branch("primary_vertex_volume", &m_calculated_vertex_v, "primary_vertex_volume/F");
0282     m_tree->Branch("primary_vertex_chi2", &m_calculated_vertex_chi2, "primary_vertex_chi2/F");
0283     m_tree->Branch("primary_vertex_nDoF", &m_calculated_vertex_ndof, "primary_vertex_nDoF/i");
0284     m_tree->Branch("primary_vertex_ID", &m_calculated_vertex_ID, "primary_vertex_ID/I");
0285     // m_tree->Branch( "primary_vertex_Covariance",   m_calculated_vertex_cov, "primary_vertex_Covariance[6]/F", 6 );
0286     m_tree->Branch("primary_vertex_Covariance", &m_calculated_vertex_cov, "primary_vertex_Covariance[6]/F", 6);
0287   }
0288   if (m_get_all_PVs)
0289   {
0290     m_tree->Branch("all_primary_vertex_x", &allPV_x);
0291     m_tree->Branch("all_primary_vertex_y", &allPV_y);
0292     m_tree->Branch("all_primary_vertex_z", &allPV_z);
0293   }
0294 
0295   m_tree->Branch("secondary_vertex_mass_pionPID", &m_sv_mass, "secondary_vertex_mass_pionPID/F");
0296 
0297   m_tree->Branch("nPrimaryVerticesOfBC", &m_nPVs, "nPrimaryVerticesOfBC/I");
0298   m_tree->Branch("nTracksOfBC", &m_multiplicity, "nTracksOfBC/I");
0299   m_tree->Branch("nTracksOfVertex", &m_nTracksOfVertex, "nTracksOfVertex/I");
0300 
0301   m_tree->Branch("runNumber", &m_runNumber, "runNumber/I");
0302   m_tree->Branch("eventNumber", &m_evtNumber, "eventNumber/I");
0303   m_tree->Branch("BCO", &m_bco, "BCO/L");
0304 
0305   if (m_get_trigger_info)
0306   {
0307     m_trigger_info_available = buildTriggerBranches(topNode, m_tree);
0308   }
0309 }
0310 
0311 void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode,
0312                                    KFParticle motherParticle,
0313                                    const KFParticle& vertex_fillbranch,
0314                                    std::vector<KFParticle> daughters,
0315                                    std::vector<KFParticle> intermediates)
0316 {
0317   const float speedOfLight = 2.99792458e-2;
0318 
0319   KFParticle temp;
0320   KFParticle* daughterArray = &daughters[0];
0321 
0322   bool switchTrackPosition;
0323 
0324   int num_tracks_used_by_intermediates = 0;
0325   for (int k = 0; k < m_num_intermediate_states_nTuple; ++k)  // Rearrange tracks from intermediate states
0326   {
0327     for (int i = 0; i < m_num_tracks_from_intermediate_nTuple[k]; ++i)
0328     {
0329       for (int j = i + 1; j < m_num_tracks_from_intermediate_nTuple[k]; ++j)
0330       {
0331         int particleAElement = i + num_tracks_used_by_intermediates;
0332         int particleBElement = j + num_tracks_used_by_intermediates;
0333         int particleAPID = daughterArray[particleAElement].GetPDG();
0334         int particleBPID = daughterArray[particleBElement].GetPDG();
0335 
0336         if (m_get_charge_conjugate_nTuple)
0337         {
0338           float daughterA_mass = kfpTupleTools.getParticleMass(particleAPID);
0339           float daughterB_mass = kfpTupleTools.getParticleMass(particleBPID);
0340           switchTrackPosition = daughterA_mass > daughterB_mass;
0341         }
0342         else
0343         {
0344           switchTrackPosition = particleAPID > particleBPID;
0345         }
0346         if (switchTrackPosition)
0347         {
0348           temp = daughterArray[particleAElement];
0349           daughterArray[particleAElement] = daughterArray[particleBElement];
0350           daughterArray[particleBElement] = temp;
0351         }
0352       }
0353     }
0354     num_tracks_used_by_intermediates += m_num_tracks_from_intermediate_nTuple[k];
0355   }
0356 
0357   int num_remaining_tracks = m_num_tracks_nTuple - num_tracks_used_by_intermediates;
0358 
0359   for (int i = 0; i < num_remaining_tracks; i++)
0360   {
0361     for (int j = i + 1; j < num_remaining_tracks; j++)
0362     {
0363       int particleAElement = i + num_tracks_used_by_intermediates;
0364       int particleBElement = j + num_tracks_used_by_intermediates;
0365       int particleAPID = daughterArray[particleAElement].GetPDG();
0366       int particleBPID = daughterArray[particleBElement].GetPDG();
0367 
0368       if (m_get_charge_conjugate_nTuple)
0369       {
0370         float daughterA_mass = kfpTupleTools.getParticleMass(particleAPID);
0371         float daughterB_mass = kfpTupleTools.getParticleMass(particleBPID);
0372         switchTrackPosition = daughterA_mass > daughterB_mass;
0373       }
0374       else
0375       {
0376         switchTrackPosition = particleAPID > particleBPID;
0377       }
0378       if (switchTrackPosition)
0379       {
0380         temp = daughterArray[particleAElement];
0381         daughterArray[particleAElement] = daughterArray[particleBElement];
0382         daughterArray[particleBElement] = temp;
0383       }
0384     }
0385   }
0386 
0387   if (m_extrapolateTracksToSV_nTuple)
0388   {
0389     for (unsigned int i = 0; i < daughters.size(); ++i)
0390     {
0391       daughterArray[i].SetProductionVertex(motherParticle);
0392     }
0393   }
0394 
0395   if (m_constrain_to_vertex_nTuple)
0396   {
0397     m_calculated_mother_dira = kfpTupleTools.eventDIRA(motherParticle, vertex_fillbranch);
0398     m_calculated_mother_dira_xy = kfpTupleTools.eventDIRA(motherParticle, vertex_fillbranch, false);
0399     m_calculated_mother_fdchi2 = kfpTupleTools.flightDistanceChi2(motherParticle, vertex_fillbranch);
0400     m_calculated_mother_ip = motherParticle.GetDistanceFromVertex(vertex_fillbranch);
0401     m_calculated_mother_ipchi2 = motherParticle.GetDeviationFromVertex(vertex_fillbranch);
0402     m_calculated_mother_ip_err = m_calculated_mother_ip / std::sqrt(m_calculated_mother_ipchi2);
0403     m_calculated_mother_ip_xy = motherParticle.GetDistanceFromVertexXY(vertex_fillbranch);
0404   }
0405   m_calculated_mother_x = motherParticle.GetX();
0406   m_calculated_mother_y = motherParticle.GetY();
0407   m_calculated_mother_z = motherParticle.GetZ();
0408   m_calculated_mother_px = motherParticle.GetPx();
0409   m_calculated_mother_py = motherParticle.GetPy();
0410   m_calculated_mother_pz = motherParticle.GetPz();
0411   m_calculated_mother_pe = motherParticle.GetE();
0412   m_calculated_mother_p = motherParticle.GetP();
0413   m_calculated_mother_p_err = motherParticle.GetErrP();
0414   m_calculated_mother_pt = motherParticle.GetPt();
0415   m_calculated_mother_pt_err = motherParticle.GetErrPt();
0416   m_calculated_mother_q = motherParticle.Q();
0417   //m_calculated_mother_eta = motherParticle.GetEta();
0418   motherParticle.GetEta(m_calculated_mother_eta,TempError);
0419   m_calculated_mother_rapidity = motherParticle.GetRapidity();
0420   m_calculated_mother_theta = motherParticle.GetTheta();
0421   //m_calculated_mother_phi = motherParticle.GetPhi();
0422   motherParticle.GetPhi(m_calculated_mother_phi, TempError);
0423   m_calculated_mother_v = kfpTupleTools.calculateEllipsoidVolume(motherParticle);
0424   m_calculated_mother_pdgID = motherParticle.GetPDG();
0425   // m_calculated_mother_cov          = &motherParticle.CovarianceMatrix()[0];
0426   for (int j = 0; j < 21; ++j)
0427   {
0428     m_calculated_mother_cov[j] = motherParticle.GetCovariance(j);
0429   }
0430 
0431   motherParticle.GetMass(m_calculated_mother_mass, m_calculated_mother_mass_err);
0432 
0433   KFParticle* intermediateArray = &intermediates[0];
0434   if (m_has_intermediates_nTuple)
0435   {
0436     for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
0437     {
0438       m_calculated_intermediate_dira[i] = kfpTupleTools.eventDIRA(intermediateArray[i], motherParticle);
0439       m_calculated_intermediate_fdchi2[i] = kfpTupleTools.flightDistanceChi2(intermediateArray[i], motherParticle);
0440       if (m_constrain_to_vertex_nTuple)
0441       {
0442         m_calculated_intermediate_ip[i] = intermediateArray[i].GetDistanceFromVertex(vertex_fillbranch);
0443         m_calculated_intermediate_ipchi2[i] = intermediateArray[i].GetDeviationFromVertex(vertex_fillbranch);
0444         m_calculated_intermediate_ip_err[i] = m_calculated_intermediate_ip[i] / std::sqrt(m_calculated_intermediate_ipchi2[i]);
0445         m_calculated_intermediate_ip_xy[i] = intermediateArray[i].GetDistanceFromVertexXY(vertex_fillbranch);
0446       }
0447       m_calculated_intermediate_x[i] = intermediateArray[i].GetX();
0448       m_calculated_intermediate_y[i] = intermediateArray[i].GetY();
0449       m_calculated_intermediate_z[i] = intermediateArray[i].GetZ();
0450       m_calculated_intermediate_px[i] = intermediateArray[i].GetPx();
0451       m_calculated_intermediate_py[i] = intermediateArray[i].GetPy();
0452       m_calculated_intermediate_pz[i] = intermediateArray[i].GetPz();
0453       m_calculated_intermediate_pe[i] = intermediateArray[i].GetE();
0454       m_calculated_intermediate_p[i] = intermediateArray[i].GetP();
0455       m_calculated_intermediate_p_err[i] = intermediateArray[i].GetErrP();
0456       m_calculated_intermediate_pt[i] = intermediateArray[i].GetPt();
0457       m_calculated_intermediate_pt_err[i] = intermediateArray[i].GetErrPt();
0458       m_calculated_intermediate_q[i] = intermediateArray[i].Q();  // I used to cast this as an int. clang-tidy want Uchar_t to Char_t to int
0459       //m_calculated_intermediate_eta[i] = intermediateArray[i].GetEta();
0460       intermediateArray[i].GetEta(m_calculated_intermediate_eta[i], TempError);
0461       m_calculated_intermediate_rapidity[i] = intermediateArray[i].GetRapidity();
0462       m_calculated_intermediate_theta[i] = intermediateArray[i].GetTheta();
0463       //m_calculated_intermediate_phi[i] = intermediateArray[i].GetPhi();
0464       intermediateArray[i].GetPhi(m_calculated_intermediate_phi[i], TempError);
0465       m_calculated_intermediate_v[i] = kfpTupleTools.calculateEllipsoidVolume(intermediateArray[i]);
0466       m_calculated_intermediate_pdgID[i] = intermediateArray[i].GetPDG();
0467       // m_calculated_intermediate_cov[i]          = &intermediateArray[i].CovarianceMatrix()[0];
0468       for (int j = 0; j < 21; ++j)
0469       {
0470         m_calculated_intermediate_cov[i][j] = intermediateArray[i].GetCovariance(j);
0471       }
0472 
0473       intermediateArray[i].GetMass(m_calculated_intermediate_mass[i], m_calculated_intermediate_mass_err[i]);
0474       intermediateArray[i].SetProductionVertex(motherParticle);
0475       intermediateArray[i].GetLifeTime(m_calculated_intermediate_decaytime[i], m_calculated_intermediate_decaytime_err[i]);
0476       intermediateArray[i].GetDecayLength(m_calculated_intermediate_decaylength[i], m_calculated_intermediate_decaylength_err[i]);
0477       intermediateArray[i].GetDecayLengthXY(m_calculated_intermediate_decaylength_xy[i], m_calculated_intermediate_decaylength_xy_err[i]);
0478 
0479       m_calculated_intermediate_decaytime[i] /= speedOfLight;
0480       m_calculated_intermediate_decaytime_err[i] /= speedOfLight;
0481     }
0482   }
0483 
0484   isTrackEMCalmatch = true;
0485   for (int i = 0; i < m_num_tracks_nTuple; ++i)
0486   {
0487     m_calculated_daughter_mass[i] = daughterArray[i].GetMass();
0488     if (m_constrain_to_vertex_nTuple)
0489     {
0490       m_calculated_daughter_ip[i] = daughterArray[i].GetDistanceFromVertex(vertex_fillbranch);
0491       m_calculated_daughter_ipchi2[i] = daughterArray[i].GetDeviationFromVertex(vertex_fillbranch);
0492       m_calculated_daughter_ip_err[i] = m_calculated_daughter_ip[i] / std::sqrt(m_calculated_daughter_ipchi2[i]);
0493       m_calculated_daughter_ip_xy[i] = daughterArray[i].GetDistanceFromVertexXY(vertex_fillbranch);
0494     }
0495     m_calculated_daughter_x[i] = daughterArray[i].GetX();
0496     m_calculated_daughter_y[i] = daughterArray[i].GetY();
0497     m_calculated_daughter_z[i] = daughterArray[i].GetZ();
0498     m_calculated_daughter_px[i] = daughterArray[i].GetPx();
0499     m_calculated_daughter_py[i] = daughterArray[i].GetPy();
0500     m_calculated_daughter_pz[i] = daughterArray[i].GetPz();
0501     m_calculated_daughter_pe[i] = daughterArray[i].GetE();
0502     m_calculated_daughter_p[i] = daughterArray[i].GetP();
0503     m_calculated_daughter_p_err[i] = daughterArray[i].GetErrP();
0504     m_calculated_daughter_pt[i] = daughterArray[i].GetPt();
0505     m_calculated_daughter_pt_err[i] = daughterArray[i].GetErrPt();
0506     m_calculated_daughter_q[i] = daughterArray[i].Q();
0507     //m_calculated_daughter_eta[i] = daughterArray[i].GetEta();
0508     daughterArray[i].GetEta(m_calculated_daughter_eta[i], TempError);
0509     m_calculated_daughter_rapidity[i] = daughterArray[i].GetRapidity();
0510     m_calculated_daughter_theta[i] = daughterArray[i].GetTheta();
0511     //m_calculated_daughter_phi[i] = daughterArray[i].GetPhi();
0512     daughterArray[i].GetPhi(m_calculated_daughter_phi[i], TempError);
0513     m_calculated_daughter_chi2[i] = daughterArray[i].GetChi2();
0514     m_calculated_daughter_ndof[i] = daughterArray[i].GetNDF();
0515     m_calculated_daughter_trid[i] = daughterArray[i].Id();
0516     m_calculated_daughter_pdgID[i] = daughterArray[i].GetPDG();
0517     // m_calculated_daughter_cov[i]          = &daughterArray[i].CovarianceMatrix()[0];
0518     for (int j = 0; j < 21; ++j)
0519     {
0520       m_calculated_daughter_cov[i][j] = daughterArray[i].GetCovariance(j);
0521     }
0522 
0523     //Now get bunch crossing number for the daughter particle
0524     SvtxTrackMap *thisTrackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple.c_str());
0525     SvtxTrack *thisTrack = getTrack(daughterArray[i].Id(), thisTrackMap);
0526     m_calculated_daughter_bunch_crossing[i] = thisTrack->get_crossing(); 
0527     m_calculated_daughter_dedx[i] = kfpTupleTools.get_dEdx(topNode, daughterArray[i]);
0528     //m_calculated_daughter_expected_dedx_pion[i] = kfpTupleTools.get_dEdx_fitValue((Int_t) daughterArray[i].GetQ() * daughterArray[i].GetP(), 211);
0529     //m_calculated_daughter_expected_dedx_kaon[i] = kfpTupleTools.get_dEdx_fitValue((Int_t) daughterArray[i].GetQ() * daughterArray[i].GetP(), 321);
0530     //m_calculated_daughter_expected_dedx_proton[i] = kfpTupleTools.get_dEdx_fitValue((Int_t) daughterArray[i].GetQ() * daughterArray[i].GetP(), 2212);
0531 
0532     bool tempEMCalmatch = false;
0533     if (m_calo_info)
0534     {
0535       fillCaloBranch(topNode, m_tree, daughterArray[i], i, tempEMCalmatch);
0536       if (m_require_track_emcal_match && !tempEMCalmatch)
0537       {
0538         isTrackEMCalmatch = false;
0539       }
0540     }
0541     if (m_truth_matching)
0542     {
0543       fillTruthBranch(topNode, m_tree, daughterArray[i], i, vertex_fillbranch, m_constrain_to_vertex_nTuple);
0544     }
0545     if (m_truth_matching)
0546     {
0547       getHepMCInfo(topNode, m_tree, daughterArray[i], i);
0548     }
0549     if (m_detector_info)
0550     {
0551       fillDetectorBranch(topNode, m_tree, daughterArray[i], i);
0552     }
0553   }
0554 
0555   KFVertex motherDecayVertex;
0556 
0557   int iter = 0;
0558   // Calcualte jT wrt their own mother, not grandmother
0559   for (int k = 0; k < m_num_intermediate_states_nTuple; ++k)
0560   {
0561     KFVertex intermediateDecayVertex;
0562     for (int j = 0; j < m_num_tracks_from_intermediate_nTuple[k]; ++j)
0563     {
0564       m_calculated_daughter_jt[iter] = kfpTupleTools.calculateJT(intermediateArray[k], daughterArray[iter]);
0565       intermediateDecayVertex += daughterArray[iter];
0566       ++iter;
0567     }
0568     m_calculated_intermediate_chi2[k] = intermediateDecayVertex.GetChi2();
0569     m_calculated_intermediate_ndof[k] = intermediateDecayVertex.GetNDF();
0570     motherDecayVertex += intermediateArray[k];
0571   }
0572   for (int k = 0; k < num_remaining_tracks; k++)
0573   {
0574     m_calculated_daughter_jt[iter] = kfpTupleTools.calculateJT(motherParticle, daughterArray[iter]);
0575     motherDecayVertex += daughterArray[iter];
0576     ++iter;
0577   }
0578 
0579   m_calculated_mother_chi2 = motherDecayVertex.GetChi2();
0580   m_calculated_mother_ndof = motherDecayVertex.GetNDF();
0581 
0582   iter = 0;
0583   for (int i = 0; i < m_num_tracks_nTuple; ++i)
0584   {
0585     for (int j = 0; j < m_num_tracks_nTuple; ++j)
0586     {
0587       if (i < j)
0588       {
0589         m_daughter_dca[iter] = daughterArray[i].GetDistanceFromParticle(daughterArray[j]);
0590         m_daughter_dca_xy[iter] = daughterArray[i].GetDistanceFromParticleXY(daughterArray[j]);
0591         ++iter;
0592       }
0593     }
0594   }
0595 
0596   if (m_get_all_PVs)
0597   {
0598     std::vector<KFParticle> sortedDaughterVector;
0599     sortedDaughterVector.reserve(m_num_tracks_nTuple);
0600     for (int i = 0; i < m_num_tracks_nTuple; ++i)
0601     {
0602       sortedDaughterVector.push_back(daughterArray[i]);
0603     }
0604     allPVInfo(topNode, m_tree, motherParticle, sortedDaughterVector, intermediates);
0605   }
0606 
0607   if (m_constrain_to_vertex_nTuple)
0608   {
0609     motherParticle.SetProductionVertex(vertex_fillbranch);
0610     motherParticle.GetLifeTime(m_calculated_mother_decaytime, m_calculated_mother_decaytime_err);
0611     motherParticle.GetDecayLength(m_calculated_mother_decaylength, m_calculated_mother_decaylength_err);
0612     motherParticle.GetDecayLengthXY(m_calculated_mother_decaylength_xy, m_calculated_mother_decaylength_xy_err);
0613 
0614     m_calculated_mother_decaytime /= speedOfLight;
0615     m_calculated_mother_decaytime_err /= speedOfLight;
0616 
0617     m_calculated_vertex_x = vertex_fillbranch.GetX();
0618     m_calculated_vertex_y = vertex_fillbranch.GetY();
0619     m_calculated_vertex_z = vertex_fillbranch.GetZ();
0620     m_calculated_vertex_v = kfpTupleTools.calculateEllipsoidVolume(vertex_fillbranch);
0621     m_calculated_vertex_chi2 = vertex_fillbranch.GetChi2();
0622     m_calculated_vertex_ndof = vertex_fillbranch.GetNDF();
0623 
0624     // it only makes sense to calculate PVID for non-fake vertex
0625     // (this otherwise crashes if m_use_fake_pv_nTuple is true in an event with no real vertices)
0626     if(m_use_fake_pv_nTuple) m_calculated_vertex_ID = -100; // error value returned by getPVID
0627     else m_calculated_vertex_ID = getPVID(topNode, vertex_fillbranch);
0628     // m_calculated_vertex_cov          = &vertex_fillbranch.CovarianceMatrix()[0];
0629     for (int j = 0; j < 6; ++j)
0630     {
0631       m_calculated_vertex_cov[j] = vertex_fillbranch.GetCovariance(j);
0632     }
0633   }
0634 
0635   m_sv_mass = calc_secondary_vertex_mass_noPID(daughters);
0636 
0637   kfpTupleTools.getTracksFromBC(topNode, m_calculated_daughter_bunch_crossing[0], m_vtx_map_node_name_nTuple, m_multiplicity, m_nPVs);  
0638   // cannot retrieve vertex map info from fake PV, hence the second condition
0639   if (m_constrain_to_vertex_nTuple && !m_use_fake_pv_nTuple)
0640   {
0641     m_nTracksOfVertex = kfpTupleTools.getTracksFromVertex(topNode, vertex_fillbranch, m_vtx_map_node_name_nTuple);
0642   }
0643   else
0644   {
0645     m_nTracksOfVertex = 0;
0646   }
0647 
0648   PHNodeIterator nodeIter(topNode);
0649 
0650   PHNode* evtNode = dynamic_cast<PHNode*>(nodeIter.findFirst("EventHeader"));
0651 
0652   if (evtNode)
0653   {
0654     EventHeader* evtHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0655     m_runNumber = evtHeader->get_RunNumber();
0656     m_evtNumber = evtHeader->get_EvtSequence();
0657 
0658     auto gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0659     if (!gl1packet)
0660     {
0661       gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0662     }
0663     m_bco = m_trigger_info_available ? gl1packet->lValue(0, "BCO") + m_calculated_daughter_bunch_crossing[0] : 0;
0664   }
0665   else
0666   {
0667     m_runNumber = m_evtNumber = m_bco = -1;
0668   }
0669 
0670   if (m_trigger_info_available)
0671   {
0672     fillTriggerBranches(topNode);
0673   }
0674 
0675   if (fillConditionMet())
0676   {
0677     m_tree->Fill();
0678   }
0679 
0680   if (m_truth_matching || m_detector_info)
0681   {
0682     clearVectors();
0683   }
0684 
0685   if (m_trigger_info_available)
0686   {
0687     resetTriggerBranches();
0688   }
0689 }
0690 
0691 float KFParticle_nTuple::calc_secondary_vertex_mass_noPID(std::vector<KFParticle> kfp_daughters)
0692 {
0693   KFParticle mother_noPID;
0694   KFParticle* daughterArray = &kfp_daughters[0];
0695 
0696   for (int i = 0; i < m_num_tracks_nTuple; ++i)
0697   {
0698     KFParticle daughter_noPID;
0699     float f_trackParameters[6], f_trackCovariance[21];
0700     for (int j = 0; j < 6; ++j)
0701     {
0702       f_trackParameters[j] = daughterArray[i].GetParameter(j);
0703     }
0704     for (int j = 0; j < 21; ++j)
0705     {
0706       f_trackCovariance[j] = daughterArray[i].GetCovariance(j);
0707     }
0708     daughter_noPID.Create(f_trackParameters, f_trackCovariance, daughterArray[i].Q(), -1);
0709     mother_noPID.AddDaughter(daughter_noPID);
0710   }
0711 
0712   return mother_noPID.GetMass();
0713 }
0714 
0715 bool KFParticle_nTuple::fillConditionMet()
0716 {
0717   // return true if this is a track-only analysis
0718   if (!m_calo_info)
0719   {
0720     return true;
0721   }
0722 
0723   // return true if do not require track-calo matching
0724   if (!m_require_track_emcal_match)
0725   {
0726     return true;
0727   }
0728 
0729   // if requiring track-calo matching, the match result is returned
0730   return isTrackEMCalmatch;
0731 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!