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