Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFParticle_nTuple.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFParticle_nTuple.cc
1 #include "KFParticle_nTuple.h"
2 
3 #include "KFParticle_Tools.h"
4 
5 #include <ffaobjects/EventHeaderv1.h>
6 
7 #include <phool/PHNodeIterator.h> // for PHNodeIterator
8 #include <phool/getClass.h>
9 
10 #include <KFParticle.h>
11 
12 #include <Rtypes.h>
13 #include <TString.h> // for TString, operator+
14 #include <TTree.h>
15 
16 #include <algorithm> // for max
17 #include <cmath>
18 #include <cstdlib> // for abs, size_t
19 #include <map> // for map, _Rb_tree_iterator, map<>:...
20 #include <utility> // for pair
21 
22 class PHCompositeNode;
23 class PHNode;
24 
26 KFParticle_Tools kfpTupleTools;
27 
28 KFParticle_nTuple::KFParticle_nTuple()
29  : m_has_intermediates_nTuple(false)
30  , m_extrapolateTracksToSV_nTuple(true)
31  , m_constrain_to_vertex_nTuple(false)
32  , m_get_all_PVs(false)
33  , m_truth_matching(false)
34  , m_detector_info(false)
35  , m_calo_info(true)
36  , m_use_intermediate_name(true)
37  , m_get_charge_conjugate_nTuple(true)
38  , m_tree(nullptr)
39 {
40 } // Constructor
41 
42 void KFParticle_nTuple::initializeVariables()
43 {
44  // m_calculated_mother_cov = -99;
45 }
46 
47 void KFParticle_nTuple::initializeBranches()
48 {
49  delete m_tree;
50  m_tree = new TTree("DecayTree", "DecayTree");
51  m_tree->OptimizeBaskets();
52  m_tree->SetAutoSave(-5e6); // Save the output file every 5MB
53 
54  std::string mother_name;
55  if (m_mother_name.empty())
56  {
57  mother_name = "mother";
58  }
59  else
60  {
61  mother_name = m_mother_name;
62  }
63 
64  size_t pos;
65  std::string undrscr = "_";
66  std::string nothing = "";
67  std::map<std::string, std::string> forbiddenStrings;
68  forbiddenStrings["/"] = undrscr;
69  forbiddenStrings["("] = undrscr;
70  forbiddenStrings[")"] = nothing;
71  forbiddenStrings["+"] = "plus";
72  forbiddenStrings["-"] = "minus";
73  forbiddenStrings["*"] = "star";
74  for (auto const& [badString, goodString] : forbiddenStrings)
75  {
76  while ((pos = mother_name.find(badString)) != std::string::npos)
77  {
78  mother_name.replace(pos, 1, goodString);
79  }
80  }
81 
82  m_tree->Branch(TString(mother_name) + "_mass", &m_calculated_mother_mass, TString(mother_name) + "_mass/F");
83  m_tree->Branch(TString(mother_name) + "_massErr", &m_calculated_mother_mass_err, TString(mother_name) + "_massErr/F");
84  if (m_constrain_to_vertex_nTuple)
85  {
86  m_tree->Branch(TString(mother_name) + "_decayTime", &m_calculated_mother_decaytime, TString(mother_name) + "_decayTime/F");
87  m_tree->Branch(TString(mother_name) + "_decayTimeErr", &m_calculated_mother_decaytime_err, TString(mother_name) + "_decayTimeErr/F");
88  m_tree->Branch(TString(mother_name) + "_decayLength", &m_calculated_mother_decaylength, TString(mother_name) + "_decayLength/F");
89  m_tree->Branch(TString(mother_name) + "_decayLengthErr", &m_calculated_mother_decaylength_err, TString(mother_name) + "_decayLengthErr/F");
90  m_tree->Branch(TString(mother_name) + "_DIRA", &m_calculated_mother_dira, TString(mother_name) + "_DIRA/F");
91  m_tree->Branch(TString(mother_name) + "_FDchi2", &m_calculated_mother_fdchi2, TString(mother_name) + "_FDchi2/F");
92  m_tree->Branch(TString(mother_name) + "_IP", &m_calculated_mother_ip, TString(mother_name) + "_IP/F");
93  m_tree->Branch(TString(mother_name) + "_IPchi2", &m_calculated_mother_ipchi2, TString(mother_name) + "_IPchi2/F");
94  m_tree->Branch(TString(mother_name) + "_IPErr", &m_calculated_mother_ip_err, TString(mother_name) + "_IPErr/F");
95  m_tree->Branch(TString(mother_name) + "_IP_xy", &m_calculated_mother_ip_xy, TString(mother_name) + "_IP_xy/F");
96  }
97  if (m_get_all_PVs)
98  {
99  m_tree->Branch(TString(mother_name) + "_IP_allPV", &allPV_mother_IP);
100  m_tree->Branch(TString(mother_name) + "_IPchi2_allPV", &allPV_mother_IPchi2);
101  }
102  m_tree->Branch(TString(mother_name) + "_x", &m_calculated_mother_x, TString(mother_name) + "_x/F");
103  m_tree->Branch(TString(mother_name) + "_y", &m_calculated_mother_y, TString(mother_name) + "_y/F");
104  m_tree->Branch(TString(mother_name) + "_z", &m_calculated_mother_z, TString(mother_name) + "_z/F");
105  m_tree->Branch(TString(mother_name) + "_px", &m_calculated_mother_px, TString(mother_name) + "_px/F");
106  m_tree->Branch(TString(mother_name) + "_py", &m_calculated_mother_py, TString(mother_name) + "_py/F");
107  m_tree->Branch(TString(mother_name) + "_pz", &m_calculated_mother_pz, TString(mother_name) + "_pz/F");
108  m_tree->Branch(TString(mother_name) + "_pE", &m_calculated_mother_pe, TString(mother_name) + "_pE/F");
109  m_tree->Branch(TString(mother_name) + "_p", &m_calculated_mother_p, TString(mother_name) + "_p/F");
110  m_tree->Branch(TString(mother_name) + "_pErr", &m_calculated_mother_p_err, TString(mother_name) + "_pErr/F");
111  m_tree->Branch(TString(mother_name) + "_pT", &m_calculated_mother_pt, TString(mother_name) + "_pT/F");
112  m_tree->Branch(TString(mother_name) + "_pTErr", &m_calculated_mother_pt_err, TString(mother_name) + "_pTErr/F");
113  m_tree->Branch(TString(mother_name) + "_charge", &m_calculated_mother_q, TString(mother_name) + "_charge/B");
114  m_tree->Branch(TString(mother_name) + "_pseudorapidity", &m_calculated_mother_eta, TString(mother_name) + "_pseudorapidity/F");
115  m_tree->Branch(TString(mother_name) + "_rapidity", &m_calculated_mother_rapidity, TString(mother_name) + "_rapidity/F");
116  m_tree->Branch(TString(mother_name) + "_theta", &m_calculated_mother_theta, TString(mother_name) + "_theta/F");
117  m_tree->Branch(TString(mother_name) + "_phi", &m_calculated_mother_phi, TString(mother_name) + "_phi/F");
118  m_tree->Branch(TString(mother_name) + "_vertex_volume", &m_calculated_mother_v, TString(mother_name) + "_vertex_volume/F");
119  m_tree->Branch(TString(mother_name) + "_chi2", &m_calculated_mother_chi2, TString(mother_name) + "_chi2/F");
120  m_tree->Branch(TString(mother_name) + "_nDoF", &m_calculated_mother_ndof, TString(mother_name) + "_nDoF/I");
121  m_tree->Branch(TString(mother_name) + "_PDG_ID", &m_calculated_mother_pdgID, TString(mother_name) + "_PDG_ID/I");
122  m_tree->Branch(TString(mother_name) + "_Covariance", &m_calculated_mother_cov, TString(mother_name) + "_Covariance[21]/F", 21);
123 
124  std::vector<std::string> intermediateNameMapping; // What intermediate is associate to what track
125  if (m_has_intermediates_nTuple)
126  {
127  for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
128  {
129  std::string intermediate_name;
130  if (!m_use_intermediate_name)
131  {
132  intermediate_name = "intermediate_" + std::to_string(i + 1);
133  }
134  else
135  {
136  intermediate_name = m_intermediate_name_ntuple[i];
137  }
138 
139  // 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
140  for (auto const& [badString, goodString] : forbiddenStrings)
141  {
142  while ((pos = intermediate_name.find(badString)) != std::string::npos)
143  {
144  intermediate_name.replace(pos, 1, goodString);
145  }
146  }
147 
148  m_tree->Branch(TString(intermediate_name) + "_mass", &m_calculated_intermediate_mass[i], TString(intermediate_name) + "_mass/F");
149  m_tree->Branch(TString(intermediate_name) + "_massErr", &m_calculated_intermediate_mass_err[i], TString(intermediate_name) + "_massErr/F");
150  m_tree->Branch(TString(intermediate_name) + "_decayTime", &m_calculated_intermediate_decaytime[i], TString(intermediate_name) + "_decayTime/F");
151  m_tree->Branch(TString(intermediate_name) + "_decayTimeErr", &m_calculated_intermediate_decaytime_err[i], TString(intermediate_name) + "_decayTimeErr/F");
152  m_tree->Branch(TString(intermediate_name) + "_decayLength", &m_calculated_intermediate_decaylength[i], TString(intermediate_name) + "_decayLength/F");
153  m_tree->Branch(TString(intermediate_name) + "_decayLengthErr", &m_calculated_intermediate_decaylength_err[i], TString(intermediate_name) + "_decayLengthErr/F");
154  m_tree->Branch(TString(intermediate_name) + "_DIRA", &m_calculated_intermediate_dira[i], TString(intermediate_name) + "_DIRA/F");
155  m_tree->Branch(TString(intermediate_name) + "_FDchi2", &m_calculated_intermediate_fdchi2[i], TString(intermediate_name) + "_FDchi2/F");
156  if (m_constrain_to_vertex_nTuple)
157  {
158  m_tree->Branch(TString(intermediate_name) + "_IP", &m_calculated_intermediate_ip[i], TString(intermediate_name) + "_IP/F");
159  m_tree->Branch(TString(intermediate_name) + "_IPchi2", &m_calculated_intermediate_ipchi2[i], TString(intermediate_name) + "_IPchi2/F");
160  m_tree->Branch(TString(intermediate_name) + "_IPErr", &m_calculated_intermediate_ip_err[i], TString(intermediate_name) + "_IPErr/F");
161  m_tree->Branch(TString(intermediate_name) + "_IP_xy", &m_calculated_intermediate_ip_xy[i], TString(intermediate_name) + "_IP_xy/F");
162  }
163  if (m_get_all_PVs)
164  {
165  m_tree->Branch(TString(intermediate_name) + "_IP_allPV", &allPV_intermediates_IP[i]);
166  m_tree->Branch(TString(intermediate_name) + "_IPchi2_allPV", &allPV_intermediates_IPchi2[i]);
167  }
168  m_tree->Branch(TString(intermediate_name) + "_x", &m_calculated_intermediate_x[i], TString(intermediate_name) + "_x/F");
169  m_tree->Branch(TString(intermediate_name) + "_y", &m_calculated_intermediate_y[i], TString(intermediate_name) + "_y/F");
170  m_tree->Branch(TString(intermediate_name) + "_z", &m_calculated_intermediate_z[i], TString(intermediate_name) + "_z/F");
171  m_tree->Branch(TString(intermediate_name) + "_px", &m_calculated_intermediate_px[i], TString(intermediate_name) + "_px/F");
172  m_tree->Branch(TString(intermediate_name) + "_py", &m_calculated_intermediate_py[i], TString(intermediate_name) + "_py/F");
173  m_tree->Branch(TString(intermediate_name) + "_pz", &m_calculated_intermediate_pz[i], TString(intermediate_name) + "_pz/F");
174  m_tree->Branch(TString(intermediate_name) + "_pE", &m_calculated_intermediate_pe[i], TString(intermediate_name) + "_pE/F");
175  m_tree->Branch(TString(intermediate_name) + "_p", &m_calculated_intermediate_p[i], TString(intermediate_name) + "_p/F");
176  m_tree->Branch(TString(intermediate_name) + "_pErr", &m_calculated_intermediate_p_err[i], TString(intermediate_name) + "_pErr/F");
177  m_tree->Branch(TString(intermediate_name) + "_pT", &m_calculated_intermediate_pt[i], TString(intermediate_name) + "_pT/F");
178  m_tree->Branch(TString(intermediate_name) + "_pTErr", &m_calculated_intermediate_pt_err[i], TString(intermediate_name) + "_pTErr/F");
179  m_tree->Branch(TString(intermediate_name) + "_charge", &m_calculated_intermediate_q[i], TString(intermediate_name) + "_charge/B");
180  m_tree->Branch(TString(intermediate_name) + "_pseudorapidity", &m_calculated_intermediate_eta[i], TString(intermediate_name) + "_pseudorapidity/F");
181  m_tree->Branch(TString(intermediate_name) + "_rapidity", &m_calculated_intermediate_rapidity[i], TString(intermediate_name) + "_rapidity/F");
182  m_tree->Branch(TString(intermediate_name) + "_theta", &m_calculated_intermediate_theta[i], TString(intermediate_name) + "_theta/F");
183  m_tree->Branch(TString(intermediate_name) + "_phi", &m_calculated_intermediate_phi[i], TString(intermediate_name) + "_phi/F");
184  m_tree->Branch(TString(intermediate_name) + "_vertex_volume", &m_calculated_intermediate_v[i], TString(intermediate_name) + "_vertex_volume/F");
185  m_tree->Branch(TString(intermediate_name) + "_chi2", &m_calculated_intermediate_chi2[i], TString(intermediate_name) + "_chi2/F");
186  m_tree->Branch(TString(intermediate_name) + "_nDoF", &m_calculated_intermediate_ndof[i], TString(intermediate_name) + "_nDoF/I");
187  m_tree->Branch(TString(intermediate_name) + "_PDG_ID", &m_calculated_intermediate_pdgID[i], TString(intermediate_name) + "_PDG_ID/I");
188  m_tree->Branch(TString(intermediate_name) + "_Covariance", &m_calculated_intermediate_cov[i], TString(intermediate_name) + "_Covariance[21]/F", 21);
189 
190  for (int j = 0; j < m_num_tracks_from_intermediate_nTuple[i]; ++j)
191  {
192  intermediateNameMapping.push_back(intermediate_name + "_");
193  }
194  }
195  }
196 
197  int num_intermediate_tracks = 0;
198  for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
199  {
200  num_intermediate_tracks += m_num_tracks_from_intermediate_nTuple[i];
201  }
202 
203  for (int i = 0; i < m_num_tracks_nTuple; ++i)
204  {
205  std::string daughter_number = "track_" + std::to_string(i + 1);
206 
207  if (m_has_intermediates_nTuple && i < num_intermediate_tracks)
208  {
209  daughter_number.insert(0, intermediateNameMapping[i]);
210  }
211 
212  m_tree->Branch(TString(daughter_number) + "_mass", &m_calculated_daughter_mass[i], TString(daughter_number) + "_mass/F");
213  if (m_constrain_to_vertex_nTuple)
214  {
215  m_tree->Branch(TString(daughter_number) + "_IP", &m_calculated_daughter_ip[i], TString(daughter_number) + "_IP/F");
216  m_tree->Branch(TString(daughter_number) + "_IPchi2", &m_calculated_daughter_ipchi2[i], TString(daughter_number) + "_IPchi2/F");
217  m_tree->Branch(TString(daughter_number) + "_IPErr", &m_calculated_daughter_ip_err[i], TString(daughter_number) + "_IPErr/F");
218  m_tree->Branch(TString(daughter_number) + "_IP_xy", &m_calculated_daughter_ip_xy[i], TString(daughter_number) + "_IP_xy/F");
219  }
220  if (m_get_all_PVs)
221  {
222  m_tree->Branch(TString(daughter_number) + "_IP_allPV", &allPV_daughter_IP[i]);
223  m_tree->Branch(TString(daughter_number) + "_IPchi2_allPV", &allPV_daughter_IPchi2[i]);
224  }
225  m_tree->Branch(TString(daughter_number) + "_x", &m_calculated_daughter_x[i], TString(daughter_number) + "_x/F");
226  m_tree->Branch(TString(daughter_number) + "_y", &m_calculated_daughter_y[i], TString(daughter_number) + "_y/F");
227  m_tree->Branch(TString(daughter_number) + "_z", &m_calculated_daughter_z[i], TString(daughter_number) + "_z/F");
228  m_tree->Branch(TString(daughter_number) + "_px", &m_calculated_daughter_px[i], TString(daughter_number) + "_px/F");
229  m_tree->Branch(TString(daughter_number) + "_py", &m_calculated_daughter_py[i], TString(daughter_number) + "_py/F");
230  m_tree->Branch(TString(daughter_number) + "_pz", &m_calculated_daughter_pz[i], TString(daughter_number) + "_pz/F");
231  m_tree->Branch(TString(daughter_number) + "_pE", &m_calculated_daughter_pe[i], TString(daughter_number) + "_pE/F");
232  m_tree->Branch(TString(daughter_number) + "_p", &m_calculated_daughter_p[i], TString(daughter_number) + "_p/F");
233  m_tree->Branch(TString(daughter_number) + "_pErr", &m_calculated_daughter_p_err[i], TString(daughter_number) + "_pErr/F");
234  m_tree->Branch(TString(daughter_number) + "_pT", &m_calculated_daughter_pt[i], TString(daughter_number) + "_pT/F");
235  m_tree->Branch(TString(daughter_number) + "_pTErr", &m_calculated_daughter_pt_err[i], TString(daughter_number) + "_pTErr/F");
236  m_tree->Branch(TString(daughter_number) + "_jT", &m_calculated_daughter_jt[i], TString(daughter_number) + "_jT/F");
237  m_tree->Branch(TString(daughter_number) + "_charge", &m_calculated_daughter_q[i], TString(daughter_number) + "_charge/B");
238  m_tree->Branch(TString(daughter_number) + "_pseudorapidity", &m_calculated_daughter_eta[i], TString(daughter_number) + "_pseudorapidity/F");
239  m_tree->Branch(TString(daughter_number) + "_rapidity", &m_calculated_daughter_rapidity[i], TString(daughter_number) + "_rapidity/F");
240  m_tree->Branch(TString(daughter_number) + "_theta", &m_calculated_daughter_theta[i], TString(daughter_number) + "_theta/F");
241  m_tree->Branch(TString(daughter_number) + "_phi", &m_calculated_daughter_phi[i], TString(daughter_number) + "_phi/F");
242  m_tree->Branch(TString(daughter_number) + "_chi2", &m_calculated_daughter_chi2[i], TString(daughter_number) + "_chi2/F");
243  m_tree->Branch(TString(daughter_number) + "_nDoF", &m_calculated_daughter_ndof[i], TString(daughter_number) + "_nDoF/I");
244  m_tree->Branch(TString(daughter_number) + "_track_ID", &m_calculated_daughter_trid[i], TString(daughter_number) + "_track_ID/I");
245  m_tree->Branch(TString(daughter_number) + "_PDG_ID", &m_calculated_daughter_pdgID[i], TString(daughter_number) + "_PDG_ID/I");
246  m_tree->Branch(TString(daughter_number) + "_Covariance", &m_calculated_daughter_cov[i], TString(daughter_number) + "_Covariance[21]/F", 21);
247 
248  if (m_calo_info)
249  {
250  initializeCaloBranches(m_tree, i, daughter_number);
251  }
252  if (m_truth_matching)
253  {
254  initializeTruthBranches(m_tree, i, daughter_number, m_constrain_to_vertex_nTuple);
255  }
256  if (m_detector_info)
257  {
258  initializeDetectorBranches(m_tree, i, daughter_number);
259  }
260  }
261 
262  int iter = 0;
263  for (int i = 0; i < m_num_tracks_nTuple; ++i)
264  {
265  for (int j = 0; j < m_num_tracks_nTuple; ++j)
266  {
267  if (i < j)
268  {
269  std::string dca_branch_name = "track_" + std::to_string(i + 1) + "_track_" + std::to_string(j + 1) + "_DCA";
270  std::string dca_leaf_name = dca_branch_name + "/F";
271  m_tree->Branch(dca_branch_name.c_str(), &m_daughter_dca[iter], dca_leaf_name.c_str());
272  ++iter;
273  }
274  }
275  }
276 
277  if (m_constrain_to_vertex_nTuple)
278  {
279  m_tree->Branch("primary_vertex_x", &m_calculated_vertex_x, "primary_vertex_x/F");
280  m_tree->Branch("primary_vertex_y", &m_calculated_vertex_y, "primary_vertex_y/F");
281  m_tree->Branch("primary_vertex_z", &m_calculated_vertex_z, "primary_vertex_z/F");
282  m_tree->Branch("primary_vertex_nTracks", &m_calculated_vertex_nTracks, "primary_vertex_nTracks/I");
283  m_tree->Branch("primary_vertex_volume", &m_calculated_vertex_v, "primary_vertex_volume/F");
284  m_tree->Branch("primary_vertex_chi2", &m_calculated_vertex_chi2, "primary_vertex_chi2/F");
285  m_tree->Branch("primary_vertex_nDoF", &m_calculated_vertex_ndof, "primary_vertex_nDoF/I");
286  // m_tree->Branch( "primary_vertex_Covariance", m_calculated_vertex_cov, "primary_vertex_Covariance[6]/F", 6 );
287  m_tree->Branch("primary_vertex_Covariance", &m_calculated_vertex_cov, "primary_vertex_Covariance[6]/F", 6);
288  }
289  if (m_get_all_PVs)
290  {
291  m_tree->Branch("all_primary_vertex_x", &allPV_x);
292  m_tree->Branch("all_primary_vertex_y", &allPV_y);
293  m_tree->Branch("all_primary_vertex_z", &allPV_z);
294  }
295 
296  m_tree->Branch("secondary_vertex_mass_pionPID", &m_sv_mass, "secondary_vertex_mass_pionPID/F");
297 
298  m_tree->Branch("nPrimaryVertices", &m_nPVs, "nPrimaryVertices/I");
299  m_tree->Branch("nEventTracks", &m_multiplicity, "nEventTracks/I");
300 
301  m_tree->Branch("runNumber", &m_runNumber, "runNumber/I");
302  m_tree->Branch("eventNumber", &m_evtNumber, "eventNumber/I");
303 }
304 
305 void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode,
306  KFParticle motherParticle,
307  const KFParticle& vertex,
308  std::vector<KFParticle> daughters,
309  std::vector<KFParticle> intermediates,
310  int nPVs, int multiplicity)
311 {
312  const float speedOfLight = 2.99792458e-2;
313 
314  KFParticle temp;
315  KFParticle* daughterArray = &daughters[0];
316 
317  bool switchTrackPosition;
318 
319  int num_tracks_used_by_intermediates = 0;
320  for (int k = 0; k < m_num_intermediate_states_nTuple; ++k) // Rearrange tracks from intermediate states
321  {
322  for (int i = 0; i < m_num_tracks_from_intermediate_nTuple[k]; ++i)
323  {
324  for (int j = i + 1; j < m_num_tracks_from_intermediate_nTuple[k]; ++j)
325  {
326  int particleAElement = i + num_tracks_used_by_intermediates;
327  int particleBElement = j + num_tracks_used_by_intermediates;
328  int particleAPID = daughterArray[particleAElement].GetPDG();
329  int particleBPID = daughterArray[particleBElement].GetPDG();
330 
331  if (m_get_charge_conjugate_nTuple)
332  {
333  float daughterA_mass = kfpTupleTools.getParticleMass(particleAPID);
334  float daughterB_mass = kfpTupleTools.getParticleMass(particleBPID);
335  switchTrackPosition = daughterA_mass > daughterB_mass;
336  }
337  else
338  {
339  switchTrackPosition = particleAPID > particleBPID;
340  }
341  if (switchTrackPosition)
342  {
343  temp = daughterArray[particleAElement];
344  daughterArray[particleAElement] = daughterArray[particleBElement];
345  daughterArray[particleBElement] = temp;
346  }
347  }
348  }
349  num_tracks_used_by_intermediates += m_num_tracks_from_intermediate_nTuple[k];
350  }
351 
352  int num_remaining_tracks = m_num_tracks_nTuple - num_tracks_used_by_intermediates;
353 
354  for (int i = 0; i < num_remaining_tracks; i++)
355  {
356  for (int j = i + 1; j < num_remaining_tracks; j++)
357  {
358  int particleAElement = i + num_tracks_used_by_intermediates;
359  int particleBElement = j + num_tracks_used_by_intermediates;
360  int particleAPID = daughterArray[particleAElement].GetPDG();
361  int particleBPID = daughterArray[particleBElement].GetPDG();
362 
363  if (m_get_charge_conjugate_nTuple)
364  {
365  float daughterA_mass = kfpTupleTools.getParticleMass(particleAPID);
366  float daughterB_mass = kfpTupleTools.getParticleMass(particleBPID);
367  switchTrackPosition = daughterA_mass > daughterB_mass;
368  }
369  else
370  {
371  switchTrackPosition = particleAPID > particleBPID;
372  }
373  if (switchTrackPosition)
374  {
375  temp = daughterArray[particleAElement];
376  daughterArray[particleAElement] = daughterArray[particleBElement];
377  daughterArray[particleBElement] = temp;
378  }
379  }
380  }
381 
382  if (m_extrapolateTracksToSV_nTuple)
383  {
384  for (unsigned int i = 0; i < daughters.size(); ++i)
385  {
386  daughterArray[i].SetProductionVertex(motherParticle);
387  }
388  }
389 
390  if (m_constrain_to_vertex_nTuple)
391  {
392  m_calculated_mother_dira = kfpTupleTools.eventDIRA(motherParticle, vertex);
393  m_calculated_mother_fdchi2 = kfpTupleTools.flightDistanceChi2(motherParticle, vertex);
394  m_calculated_mother_ip = motherParticle.GetDistanceFromVertex(vertex);
395  m_calculated_mother_ipchi2 = motherParticle.GetDeviationFromVertex(vertex);
396  m_calculated_mother_ip_err = m_calculated_mother_ip / std::sqrt(m_calculated_mother_ipchi2);
397  m_calculated_mother_ip_xy = motherParticle.GetDistanceFromVertexXY(vertex);
398  }
399  m_calculated_mother_x = motherParticle.GetX();
400  m_calculated_mother_y = motherParticle.GetY();
401  m_calculated_mother_z = motherParticle.GetZ();
402  m_calculated_mother_px = motherParticle.GetPx();
403  m_calculated_mother_py = motherParticle.GetPy();
404  m_calculated_mother_pz = motherParticle.GetPz();
405  m_calculated_mother_pe = motherParticle.GetE();
406  m_calculated_mother_p = motherParticle.GetP();
407  m_calculated_mother_p_err = motherParticle.GetErrP();
408  m_calculated_mother_pt = motherParticle.GetPt();
409  m_calculated_mother_pt_err = motherParticle.GetErrPt();
410  m_calculated_mother_q = motherParticle.Q();
411  m_calculated_mother_eta = motherParticle.GetEta();
412  m_calculated_mother_rapidity = motherParticle.GetRapidity();
413  m_calculated_mother_theta = motherParticle.GetTheta();
414  m_calculated_mother_phi = motherParticle.GetPhi();
415  m_calculated_mother_v = kfpTupleTools.calculateEllipsoidVolume(motherParticle);
416  m_calculated_mother_chi2 = motherParticle.GetChi2();
417  m_calculated_mother_ndof = motherParticle.GetNDF();
418  m_calculated_mother_pdgID = motherParticle.GetPDG();
419  // m_calculated_mother_cov = &motherParticle.CovarianceMatrix()[0];
420  for (int j = 0; j < 21; ++j)
421  {
422  m_calculated_mother_cov[j] = motherParticle.GetCovariance(j);
423  }
424 
425  motherParticle.GetMass(m_calculated_mother_mass, m_calculated_mother_mass_err);
426 
427  KFParticle* intermediateArray = &intermediates[0];
428  if (m_has_intermediates_nTuple)
429  {
430  for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
431  {
432  m_calculated_intermediate_dira[i] = kfpTupleTools.eventDIRA(intermediateArray[i], motherParticle);
433  m_calculated_intermediate_fdchi2[i] = kfpTupleTools.flightDistanceChi2(intermediateArray[i], motherParticle);
434  if (m_constrain_to_vertex_nTuple)
435  {
436  m_calculated_intermediate_ip[i] = intermediateArray[i].GetDistanceFromVertex(vertex);
437  m_calculated_intermediate_ipchi2[i] = intermediateArray[i].GetDeviationFromVertex(vertex);
438  m_calculated_intermediate_ip_err[i] = m_calculated_intermediate_ip[i] / std::sqrt(m_calculated_intermediate_ipchi2[i]);
439  m_calculated_intermediate_ip_xy[i] = intermediateArray[i].GetDistanceFromVertexXY(vertex);
440  }
441  m_calculated_intermediate_x[i] = intermediateArray[i].GetX();
442  m_calculated_intermediate_y[i] = intermediateArray[i].GetY();
443  m_calculated_intermediate_z[i] = intermediateArray[i].GetZ();
444  m_calculated_intermediate_px[i] = intermediateArray[i].GetPx();
445  m_calculated_intermediate_py[i] = intermediateArray[i].GetPy();
446  m_calculated_intermediate_pz[i] = intermediateArray[i].GetPz();
447  m_calculated_intermediate_pe[i] = intermediateArray[i].GetE();
448  m_calculated_intermediate_p[i] = intermediateArray[i].GetP();
449  m_calculated_intermediate_p_err[i] = intermediateArray[i].GetErrP();
450  m_calculated_intermediate_pt[i] = intermediateArray[i].GetPt();
451  m_calculated_intermediate_pt_err[i] = intermediateArray[i].GetErrPt();
452  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
453  m_calculated_intermediate_eta[i] = intermediateArray[i].GetEta();
454  m_calculated_intermediate_rapidity[i] = intermediateArray[i].GetRapidity();
455  m_calculated_intermediate_theta[i] = intermediateArray[i].GetTheta();
456  m_calculated_intermediate_phi[i] = intermediateArray[i].GetPhi();
457  m_calculated_intermediate_v[i] = kfpTupleTools.calculateEllipsoidVolume(intermediateArray[i]);
458  m_calculated_intermediate_chi2[i] = intermediateArray[i].GetChi2();
459  m_calculated_intermediate_ndof[i] = intermediateArray[i].GetNDF();
460  m_calculated_intermediate_pdgID[i] = intermediateArray[i].GetPDG();
461  // m_calculated_intermediate_cov[i] = &intermediateArray[i].CovarianceMatrix()[0];
462  for (int j = 0; j < 21; ++j)
463  {
464  m_calculated_intermediate_cov[i][j] = intermediateArray[i].GetCovariance(j);
465  }
466 
467  intermediateArray[i].GetMass(m_calculated_intermediate_mass[i], m_calculated_intermediate_mass_err[i]);
468  intermediateArray[i].SetProductionVertex(motherParticle);
469  intermediateArray[i].GetLifeTime(m_calculated_intermediate_decaytime[i], m_calculated_intermediate_decaytime_err[i]);
470  intermediateArray[i].GetDecayLength(m_calculated_intermediate_decaylength[i], m_calculated_intermediate_decaylength_err[i]);
471 
472  m_calculated_intermediate_decaytime[i] /= speedOfLight;
473  m_calculated_intermediate_decaytime_err[i] /= speedOfLight;
474  }
475  }
476 
477  for (int i = 0; i < m_num_tracks_nTuple; ++i)
478  {
479  m_calculated_daughter_mass[i] = daughterArray[i].GetMass();
480  if (m_constrain_to_vertex_nTuple)
481  {
482  m_calculated_daughter_ip[i] = daughterArray[i].GetDistanceFromVertex(vertex);
483  m_calculated_daughter_ipchi2[i] = daughterArray[i].GetDeviationFromVertex(vertex);
484  m_calculated_daughter_ip_err[i] = m_calculated_daughter_ip[i] / std::sqrt(m_calculated_daughter_ipchi2[i]);
485  m_calculated_daughter_ip_xy[i] = daughterArray[i].GetDistanceFromVertexXY(vertex);
486  }
487  m_calculated_daughter_x[i] = daughterArray[i].GetX();
488  m_calculated_daughter_y[i] = daughterArray[i].GetY();
489  m_calculated_daughter_z[i] = daughterArray[i].GetZ();
490  m_calculated_daughter_px[i] = daughterArray[i].GetPx();
491  m_calculated_daughter_py[i] = daughterArray[i].GetPy();
492  m_calculated_daughter_pz[i] = daughterArray[i].GetPz();
493  m_calculated_daughter_pe[i] = daughterArray[i].GetE();
494  m_calculated_daughter_p[i] = daughterArray[i].GetP();
495  m_calculated_daughter_p_err[i] = daughterArray[i].GetErrP();
496  m_calculated_daughter_pt[i] = daughterArray[i].GetPt();
497  m_calculated_daughter_pt_err[i] = daughterArray[i].GetErrPt();
498  m_calculated_daughter_q[i] = daughterArray[i].Q();
499  m_calculated_daughter_eta[i] = daughterArray[i].GetEta();
500  m_calculated_daughter_rapidity[i] = daughterArray[i].GetRapidity();
501  m_calculated_daughter_theta[i] = daughterArray[i].GetTheta();
502  m_calculated_daughter_phi[i] = daughterArray[i].GetPhi();
503  m_calculated_daughter_chi2[i] = daughterArray[i].GetChi2();
504  m_calculated_daughter_ndof[i] = daughterArray[i].GetNDF();
505  m_calculated_daughter_trid[i] = daughterArray[i].Id();
506  m_calculated_daughter_pdgID[i] = daughterArray[i].GetPDG();
507  // m_calculated_daughter_cov[i] = &daughterArray[i].CovarianceMatrix()[0];
508  for (int j = 0; j < 21; ++j)
509  {
510  m_calculated_daughter_cov[i][j] = daughterArray[i].GetCovariance(j);
511  }
512 
513  if (m_calo_info)
514  {
515  fillCaloBranch(topNode, m_tree, daughterArray[i], i);
516  }
517  if (m_truth_matching)
518  {
519  fillTruthBranch(topNode, m_tree, daughterArray[i], i, vertex, m_constrain_to_vertex_nTuple);
520  }
521  if (m_truth_matching)
522  {
523  getHepMCInfo(topNode, m_tree, daughterArray[i], i);
524  }
525  if (m_detector_info)
526  {
527  fillDetectorBranch(topNode, m_tree, daughterArray[i], i);
528  }
529  }
530 
531  int iter = 0;
532  // Calcualte jT wrt their own mother, not grandmother
533  for (int k = 0; k < m_num_intermediate_states_nTuple; ++k)
534  {
535  for (int j = 0; j < m_num_tracks_from_intermediate_nTuple[k]; ++j)
536  {
537  m_calculated_daughter_jt[iter] = kfpTupleTools.calculateJT(intermediateArray[k], daughterArray[iter]);
538  ++iter;
539  }
540  }
541  for (int k = 0; k < num_remaining_tracks; k++)
542  {
543  m_calculated_daughter_jt[iter] = kfpTupleTools.calculateJT(motherParticle, daughterArray[iter]);
544  ++iter;
545  }
546 
547  iter = 0;
548  for (int i = 0; i < m_num_tracks_nTuple; ++i)
549  {
550  for (int j = 0; j < m_num_tracks_nTuple; ++j)
551  {
552  if (i < j)
553  {
554  m_daughter_dca[iter] = daughterArray[i].GetDistanceFromParticle(daughterArray[j]);
555  ++iter;
556  }
557  }
558  }
559 
560  if (m_get_all_PVs)
561  {
562  std::vector<KFParticle> sortedDaughterVector;
563  sortedDaughterVector.reserve(m_num_tracks_nTuple);
564  for (int i = 0; i < m_num_tracks_nTuple; ++i)
565  {
566  sortedDaughterVector.push_back(daughterArray[i]);
567  }
568  allPVInfo(topNode, m_tree, motherParticle, sortedDaughterVector, intermediates);
569  }
570 
571  if (m_constrain_to_vertex_nTuple)
572  {
573  motherParticle.SetProductionVertex(vertex);
574  motherParticle.GetLifeTime(m_calculated_mother_decaytime, m_calculated_mother_decaytime_err);
575  motherParticle.GetDecayLength(m_calculated_mother_decaylength, m_calculated_mother_decaylength_err);
576 
577  m_calculated_mother_decaytime /= speedOfLight;
578  m_calculated_mother_decaytime_err /= speedOfLight;
579 
580  m_calculated_vertex_x = vertex.GetX();
581  m_calculated_vertex_y = vertex.GetY();
582  m_calculated_vertex_z = vertex.GetZ();
583  m_calculated_vertex_v = kfpTupleTools.calculateEllipsoidVolume(vertex);
584  m_calculated_vertex_chi2 = vertex.GetChi2();
585  m_calculated_vertex_ndof = vertex.GetNDF();
586  // m_calculated_vertex_cov = &vertex.CovarianceMatrix()[0];
587  for (int j = 0; j < 6; ++j)
588  {
589  m_calculated_vertex_cov[j] = vertex.GetCovariance(j);
590  }
591  m_calculated_vertex_nTracks = kfpTupleTools.getTracksFromVertex(topNode, vertex, m_vtx_map_node_name_nTuple);
592  }
593 
594  m_sv_mass = calc_secondary_vertex_mass_noPID(daughters);
595 
596  m_nPVs = nPVs;
597  m_multiplicity = multiplicity;
598 
599  PHNodeIterator nodeIter(topNode);
600 
601  PHNode* evtNode = dynamic_cast<PHNode*>(nodeIter.findFirst("EventHeader"));
602 
603  if (evtNode)
604  {
605  EventHeaderv1* evtHeader = findNode::getClass<EventHeaderv1>(topNode, "EventHeader");
606  m_runNumber = evtHeader->get_RunNumber();
607  m_evtNumber = evtHeader->get_EvtSequence();
608  }
609  else
610  {
611  m_runNumber = m_evtNumber = -1;
612  }
613 
614  m_tree->Fill();
615 
616  if (m_truth_matching || m_detector_info)
617  {
618  clearVectors();
619  }
620 }
621 
622 float KFParticle_nTuple::calc_secondary_vertex_mass_noPID(std::vector<KFParticle> kfp_daughters)
623 {
624  KFParticle mother_noPID;
625  KFParticle* daughterArray = &kfp_daughters[0];
626 
627  for (int i = 0; i < m_num_tracks_nTuple; ++i)
628  {
629  KFParticle daughter_noPID;
630  float f_trackParameters[6], f_trackCovariance[21];
631  for (int j = 0; j < 6; ++j)
632  {
633  f_trackParameters[j] = daughterArray[i].GetParameter(j);
634  }
635  for (int j = 0; j < 21; ++j)
636  {
637  f_trackCovariance[j] = daughterArray[i].GetCovariance(j);
638  }
639  daughter_noPID.Create(f_trackParameters, f_trackCovariance, daughterArray[i].Q(), -1);
640  mother_noPID.AddDaughter(daughter_noPID);
641  }
642 
643  return mother_noPID.GetMass();
644 }