Version: SMASH-2.0
rootoutput.cc
Go to the documentation of this file.
1 /*
2  *
3  * Copyright (c) 2014-2020
4  * SMASH Team
5  *
6  * GNU General Public License (GPLv3 or later)
7  *
8  */
9 
10 #include "smash/rootoutput.h"
11 #include "TFile.h"
12 #include "TTree.h"
13 #include "smash/action.h"
14 #include "smash/clock.h"
15 #include "smash/forwarddeclarations.h"
16 #include "smash/particles.h"
17 
18 namespace smash {
19 static constexpr int LHyperSurfaceCrossing = LogArea::HyperSurfaceCrossing::id;
20 static constexpr int LOutput = LogArea::Output::id;
21 
22 const int RootOutput::max_buffer_size_ = 500000;
23 
225 RootOutput::RootOutput(const bf::path &path, const std::string &name,
226  const OutputParameters &out_par)
227  : OutputInterface(name),
228  filename_(path / (name + ".root")),
229  write_collisions_(name == "Collisions" || name == "Dileptons" ||
230  name == "Photons"),
231  write_particles_(name == "Particles"),
232  write_initial_conditions_(name == "SMASH_IC"),
233  particles_only_final_(out_par.part_only_final),
234  autosave_frequency_(1000),
235  part_extended_(out_par.part_extended),
236  coll_extended_(out_par.coll_extended),
237  ic_extended_(out_par.ic_extended) {
238  filename_unfinished_ = filename_;
239  filename_unfinished_ += ".unfinished";
240  root_out_file_ =
241  make_unique<TFile>(filename_unfinished_.native().c_str(), "NEW");
242  init_trees();
243 }
244 
245 void RootOutput::init_trees() {
246  if (write_particles_ || write_initial_conditions_) {
247  particles_tree_ = new TTree("particles", "particles");
248 
249  particles_tree_->Branch("ev", &ev_, "ev/I");
250  particles_tree_->Branch("tcounter", &tcounter_, "tcounter/I");
251  particles_tree_->Branch("npart", &npart_, "npart/I");
252  particles_tree_->Branch("test_p", &test_p_, "test_p/I");
253  particles_tree_->Branch("modus_l", &modus_l_, "modus_l/D");
254  particles_tree_->Branch("current_t", &current_t_, "current_t/D");
255  particles_tree_->Branch("impact_b", &impact_b_, "impact_b/D");
256  particles_tree_->Branch("empty_event", &empty_event_, "empty_event/O");
257 
258  particles_tree_->Branch("pdgcode", &pdgcode_[0], "pdgcode[npart]/I");
259  particles_tree_->Branch("charge", &charge_[0], "charge[npart]/I");
260 
261  particles_tree_->Branch("p0", &p0_[0], "p0[npart]/D");
262  particles_tree_->Branch("px", &px_[0], "px[npart]/D");
263  particles_tree_->Branch("py", &py_[0], "py[npart]/D");
264  particles_tree_->Branch("pz", &pz_[0], "pz[npart]/D");
265 
266  particles_tree_->Branch("t", &t_[0], "t[npart]/D");
267  particles_tree_->Branch("x", &x_[0], "x[npart]/D");
268  particles_tree_->Branch("y", &y_[0], "y[npart]/D");
269  particles_tree_->Branch("z", &z_[0], "z[npart]/D");
270 
271  particles_tree_->Branch("E_kinetic_tot", &E_kinetic_tot_,
272  "E_kinetic_tot/D");
273  particles_tree_->Branch("E_fields_tot", &E_fields_tot_, "E_fields_tot/D");
274  particles_tree_->Branch("E_tot", &E_tot_, "E_tot/D");
275 
276  if (part_extended_ || ic_extended_) {
277  particles_tree_->Branch("ncoll", &coll_per_part_[0], "ncoll[npart]/I");
278  particles_tree_->Branch("form_time", &formation_time_[0],
279  "form_time[npart]/D");
280  particles_tree_->Branch("xsecfac", &xsec_factor_[0], "xsecfac[npart]/D");
281  particles_tree_->Branch("proc_id_origin", &proc_id_origin_[0],
282  "proc_id_origin[npart]/I");
283  particles_tree_->Branch("proc_type_origin", &proc_type_origin_[0],
284  "proc_type_origin[npart]/I");
285  particles_tree_->Branch("time_last_coll", &time_last_coll_[0],
286  "time_last_coll[npart]/D");
287  particles_tree_->Branch("pdg_mother1", &pdg_mother1_[0],
288  "pdg_mother1[npart]/I");
289  particles_tree_->Branch("pdg_mother2", &pdg_mother2_[0],
290  "pdg_mother2[npart]/I");
291  }
292  }
293 
294  if (write_collisions_) {
295  collisions_tree_ = new TTree("collisions", "collisions");
296 
297  collisions_tree_->Branch("nin", &nin_, "nin/I");
298  collisions_tree_->Branch("nout", &nout_, "nout/I");
299  collisions_tree_->Branch("npart", &npart_, "npart/I");
300  collisions_tree_->Branch("ev", &ev_, "ev/I");
301  collisions_tree_->Branch("weight", &wgt_, "weight/D");
302  collisions_tree_->Branch("partial_weight", &par_wgt_, "partial_weight/D");
303 
304  collisions_tree_->Branch("pdgcode", &pdgcode_[0], "pdgcode[npart]/I");
305  collisions_tree_->Branch("charge", &charge_[0], "charge[npart]/I");
306 
307  collisions_tree_->Branch("p0", &p0_[0], "p0[npart]/D");
308  collisions_tree_->Branch("px", &px_[0], "px[npart]/D");
309  collisions_tree_->Branch("py", &py_[0], "py[npart]/D");
310  collisions_tree_->Branch("pz", &pz_[0], "pz[npart]/D");
311 
312  collisions_tree_->Branch("t", &t_[0], "t[npart]/D");
313  collisions_tree_->Branch("x", &x_[0], "x[npart]/D");
314  collisions_tree_->Branch("y", &y_[0], "y[npart]/D");
315  collisions_tree_->Branch("z", &z_[0], "z[npart]/D");
316 
317  if (coll_extended_) {
318  collisions_tree_->Branch("ncoll", &coll_per_part_[0], "ncoll[npart]/I");
319  collisions_tree_->Branch("form_time", &formation_time_[0],
320  "form_time[npart]/D");
321  collisions_tree_->Branch("xsecfac", &xsec_factor_[0], "xsecfac[npart]/D");
322  collisions_tree_->Branch("proc_id_origin", &proc_id_origin_[0],
323  "proc_id_origin[npart]/I");
324  collisions_tree_->Branch("proc_type_origin", &proc_type_origin_[0],
325  "proc_type_origin[npart]/I");
326  collisions_tree_->Branch("time_last_coll", &time_last_coll_[0],
327  "time_last_coll[npart]/D");
328  collisions_tree_->Branch("pdg_mother1", &pdg_mother1_[0],
329  "pdg_mother1[npart]/I");
330  collisions_tree_->Branch("pdg_mother2", &pdg_mother2_[0],
331  "pdg_mother2[npart]/I");
332  }
333  }
334 }
335 
340 RootOutput::~RootOutput() {
341  // kOverwrite option prevents from writing extra TKey objects into root file
342  root_out_file_->Write("", TObject::kOverwrite);
343  root_out_file_->Close();
344  bf::rename(filename_unfinished_, filename_);
345 }
346 
347 void RootOutput::at_eventstart(const Particles &particles,
348  const int event_number, const EventInfo &event) {
349  // save event number
350  current_event_ = event_number;
351 
352  modus_l_ = event.modus_length;
353  test_p_ = event.test_particles;
354  current_t_ = event.current_time;
355  E_kinetic_tot_ = event.total_kinetic_energy;
356  E_fields_tot_ = event.total_mean_field_energy;
357  E_tot_ = event.total_energy;
358 
359  if (write_particles_ && particles_only_final_ == OutputOnlyFinal::No) {
360  output_counter_ = 0;
361  // This is to have only one output of positive impact parameter per event
362  impact_b_ = -1.0;
363  empty_event_ = false;
364  particles_to_tree(particles);
365  output_counter_++;
366  }
367 }
368 
369 void RootOutput::at_intermediate_time(const Particles &particles,
370  const std::unique_ptr<Clock> &,
371  const DensityParameters &,
372  const EventInfo &event) {
373  modus_l_ = event.modus_length;
374  test_p_ = event.test_particles;
375  current_t_ = event.current_time;
376  E_kinetic_tot_ = event.total_kinetic_energy;
377  E_fields_tot_ = event.total_mean_field_energy;
378  E_tot_ = event.total_energy;
379 
380  if (write_particles_ && particles_only_final_ == OutputOnlyFinal::No) {
381  particles_to_tree(particles);
382  output_counter_++;
383  }
384 }
385 
386 void RootOutput::at_eventend(const Particles &particles,
387  const int /*event_number*/,
388  const EventInfo &event) {
389  modus_l_ = event.modus_length;
390  test_p_ = event.test_particles;
391  current_t_ = event.current_time;
392  E_kinetic_tot_ = event.total_kinetic_energy;
393  E_fields_tot_ = event.total_mean_field_energy;
394  E_tot_ = event.total_energy;
395 
396  impact_b_ = event.impact_parameter;
397  empty_event_ = event.empty_event;
398  if (write_particles_ &&
399  !(event.empty_event &&
400  particles_only_final_ == OutputOnlyFinal::IfNotEmpty)) {
401  particles_to_tree(particles);
402  }
403  /* Forced regular dump from operational memory to disk. Very demanding!
404  * If program crashes written data will NOT be lost. */
405  if (current_event_ > 0 && current_event_ % autosave_frequency_ == 0) {
406  if (write_particles_ || write_initial_conditions_) {
407  particles_tree_->AutoSave("SaveSelf");
408  }
409  if (write_collisions_) {
410  collisions_tree_->AutoSave("SaveSelf");
411  }
412  }
413 
414  if (write_initial_conditions_) {
415  // If the runtime is too short some particles might not yet have
416  // reached the hypersurface. Warning is printed.
417  if (particles.size() != 0) {
418  logg[LHyperSurfaceCrossing].warn(
419  "End time might be too small for initial conditions output. "
420  "Hypersurface has not yet been crossed by ",
421  particles.size(), " particle(s).");
422  }
423  }
424 }
425 
426 void RootOutput::at_interaction(const Action &action,
427  const double /*density*/) {
428  if (write_collisions_) {
429  collisions_to_tree(action.incoming_particles(), action.outgoing_particles(),
430  action.get_total_weight(), action.get_partial_weight());
431  }
432 
433  if (write_initial_conditions_ &&
434  action.get_type() == ProcessType::HyperSurfaceCrossing) {
435  particles_to_tree(action.incoming_particles());
436  }
437 }
438 
439 template <typename T>
440 void RootOutput::particles_to_tree(T &particles) {
441  int i = 0;
442 
443  ev_ = current_event_;
444  tcounter_ = output_counter_;
445  bool exceeded_buffer_message = true;
446 
447  for (const auto &p : particles) {
448  // Buffer full - flush to tree, else fill with particles
449  if (i >= max_buffer_size_) {
450  if (exceeded_buffer_message) {
451  logg[LOutput].warn()
452  << "\nThe number of particles N = " << particles.size()
453  << " exceeds the maximum buffer size B = " << max_buffer_size_
454  << ".\nceil(N/B) = "
455  << std::ceil(particles.size() /
456  static_cast<double>(max_buffer_size_))
457  << " separate ROOT Tree entries will be created at this output."
458  << "\nMaximum buffer size (max_buffer_size_) can be changed in "
459  << "rootoutput.h\n\n";
460  exceeded_buffer_message = false;
461  }
462  npart_ = max_buffer_size_;
463  i = 0;
464  particles_tree_->Fill();
465  } else {
466  pdgcode_[i] = p.pdgcode().get_decimal();
467  charge_[i] = p.type().charge();
468 
469  p0_[i] = p.momentum().x0();
470  px_[i] = p.momentum().x1();
471  py_[i] = p.momentum().x2();
472  pz_[i] = p.momentum().x3();
473 
474  t_[i] = p.position().x0();
475  x_[i] = p.position().x1();
476  y_[i] = p.position().x2();
477  z_[i] = p.position().x3();
478 
479  if (part_extended_ || ic_extended_) {
480  const auto h = p.get_history();
481  formation_time_[i] = p.formation_time();
482  xsec_factor_[i] = p.xsec_scaling_factor();
483  time_last_coll_[i] = h.time_last_collision;
484  coll_per_part_[i] = h.collisions_per_particle;
485  proc_id_origin_[i] = h.id_process;
486  proc_type_origin_[i] = static_cast<int>(h.process_type);
487  pdg_mother1_[i] = h.p1.get_decimal();
488  pdg_mother2_[i] = h.p2.get_decimal();
489  }
490 
491  i++;
492  }
493  }
494  // Flush rest to tree
495  if (i > 0) {
496  npart_ = i;
497  particles_tree_->Fill();
498  }
499 }
500 
501 void RootOutput::collisions_to_tree(const ParticleList &incoming,
502  const ParticleList &outgoing,
503  const double weight,
504  const double partial_weight) {
505  ev_ = current_event_;
506  nin_ = incoming.size();
507  nout_ = outgoing.size();
508  npart_ = nin_ + nout_;
509  wgt_ = weight;
510  par_wgt_ = partial_weight;
511 
512  int i = 0;
513 
514  /* It is assumed that nin + nout < max_buffer_size_
515  * This is true for any possible reaction for current buffer size: 10000
516  * But if one wants initial/final particles written to collisions
517  * then implementation should be updated. */
518 
519  for (const ParticleList &plist : {incoming, outgoing}) {
520  for (const auto &p : plist) {
521  pdgcode_[i] = p.pdgcode().get_decimal();
522  charge_[i] = p.type().charge();
523 
524  p0_[i] = p.momentum().x0();
525  px_[i] = p.momentum().x1();
526  py_[i] = p.momentum().x2();
527  pz_[i] = p.momentum().x3();
528 
529  t_[i] = p.position().x0();
530  x_[i] = p.position().x1();
531  y_[i] = p.position().x2();
532  z_[i] = p.position().x3();
533 
534  if (coll_extended_) {
535  const auto h = p.get_history();
536  formation_time_[i] = p.formation_time();
537  xsec_factor_[i] = p.xsec_scaling_factor();
538  time_last_coll_[i] = h.time_last_collision;
539  coll_per_part_[i] = h.collisions_per_particle;
540  proc_id_origin_[i] = h.id_process;
541  proc_type_origin_[i] = static_cast<int>(h.process_type);
542  pdg_mother1_[i] = h.p1.get_decimal();
543  pdg_mother2_[i] = h.p2.get_decimal();
544  }
545 
546  i++;
547  }
548  }
549 
550  collisions_tree_->Fill();
551 }
552 } // namespace smash
smash::RootOutput::~RootOutput
~RootOutput()
Destructor.
Definition: rootoutput.cc:340
smash::Action::get_type
virtual ProcessType get_type() const
Get the process type.
Definition: action.h:131
smash
Definition: action.h:24
smash::RootOutput::wgt_
double wgt_
Property that is written to ROOT output.
Definition: rootoutput.h:207
smash::RootOutput::filename_unfinished_
bf::path filename_unfinished_
Filename of output as long as simulation is still running.
Definition: rootoutput.h:135
smash::RootOutput::charge_
std::vector< int > charge_
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::E_tot_
double E_tot_
Property that is written to ROOT output.
Definition: rootoutput.h:208
smash::RootOutput::pdg_mother1_
std::vector< int > pdg_mother1_
Property that is written to ROOT output.
Definition: rootoutput.h:204
smash::RootOutput::coll_extended_
const bool coll_extended_
Whether extended collisions output is on.
Definition: rootoutput.h:238
smash::LOutput
static constexpr int LOutput
Definition: oscaroutput.cc:23
smash::DensityParameters
A class to pre-calculate and store parameters relevant for density calculation.
Definition: density.h:108
smash::RootOutput::y_
std::vector< double > y_
Property that is written to ROOT output.
Definition: rootoutput.h:192
smash::RootOutput::z_
std::vector< double > z_
Property that is written to ROOT output.
Definition: rootoutput.h:193
smash::Particles::size
size_t size() const
Definition: particles.h:87
smash::RootOutput::proc_id_origin_
std::vector< int > proc_id_origin_
Property that is written to ROOT output.
Definition: rootoutput.h:202
smash::RootOutput::pdg_mother2_
std::vector< int > pdg_mother2_
Property that is written to ROOT output.
Definition: rootoutput.h:205
smash::RootOutput::xsec_factor_
std::vector< double > xsec_factor_
Property that is written to ROOT output.
Definition: rootoutput.h:196
smash::RootOutput::p0_
std::vector< double > p0_
Property that is written to ROOT output.
Definition: rootoutput.h:186
smash::RootOutput::at_eventend
void at_eventend(const Particles &particles, const int event_number, const EventInfo &event) override
update event number and impact parameter, and writes intermediate particles to a tree.
Definition: rootoutput.cc:386
smash::RootOutput::ev_
int ev_
Property that is written to ROOT output.
Definition: rootoutput.h:206
smash::Action::get_partial_weight
virtual double get_partial_weight() const =0
Return the specific weight for the chosen outgoing channel, which is mainly used for the partial weig...
smash::RootOutput::py_
std::vector< double > py_
Property that is written to ROOT output.
Definition: rootoutput.h:188
smash::RootOutput::RootOutput
RootOutput(const bf::path &path, const std::string &name, const OutputParameters &out_par)
Construct ROOT output.
Definition: rootoutput.cc:225
smash::RootOutput::filename_
const bf::path filename_
Filename of output.
Definition: rootoutput.h:133
smash::RootOutput::E_fields_tot_
double E_fields_tot_
Property that is written to ROOT output.
Definition: rootoutput.h:208
smash::RootOutput::write_initial_conditions_
bool write_initial_conditions_
Option to write particles tree for initial conditions.
Definition: rootoutput.h:219
smash::RootOutput::pz_
std::vector< double > pz_
Property that is written to ROOT output.
Definition: rootoutput.h:189
smash::RootOutput::formation_time_
std::vector< double > formation_time_
Property that is written to ROOT output.
Definition: rootoutput.h:194
action.h
smash::RootOutput::test_p_
int test_p_
Property that is written to ROOT output.
Definition: rootoutput.h:206
smash::EventInfo
Structure to contain custom data for output.
Definition: outputinterface.h:35
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
clock.h
smash::Action::outgoing_particles
const ParticleList & outgoing_particles() const
Get the list of particles that resulted from the action.
Definition: action.h:245
smash::RootOutput::init_trees
void init_trees()
Basic initialization routine, creating the TTree objects for particles and collisions.
Definition: rootoutput.cc:245
smash::RootOutput::current_t_
double current_t_
Property that is written to ROOT output.
Definition: rootoutput.h:207
smash::RootOutput::at_interaction
void at_interaction(const Action &action, const double density) override
Writes collisions to a tree defined by treename.
Definition: rootoutput.cc:426
forwarddeclarations.h
smash::RootOutput::output_counter_
int output_counter_
Number of output in a given event.
Definition: rootoutput.h:169
smash::RootOutput::autosave_frequency_
int autosave_frequency_
Root file cannot be read if it was not properly closed and finalized.
Definition: rootoutput.h:233
smash::RootOutput::empty_event_
bool empty_event_
Property that is written to ROOT output.
Definition: rootoutput.h:209
OutputOnlyFinal::IfNotEmpty
Print only final-state particles, and those only if the event is not empty.
smash::OutputParameters
Helper structure for Experiment to hold output options and parameters.
Definition: outputparameters.h:25
smash::RootOutput::nout_
int nout_
Property that is written to ROOT output.
Definition: rootoutput.h:206
rootoutput.h
smash::RootOutput::write_particles_
bool write_particles_
Option to write particles tree.
Definition: rootoutput.h:216
smash::RootOutput::coll_per_part_
std::vector< int > coll_per_part_
Property that is written to ROOT output.
Definition: rootoutput.h:201
smash::OutputInterface
Abstraction of generic output.
Definition: outputinterface.h:65
smash::RootOutput::at_eventstart
void at_eventstart(const Particles &particles, const int event_number, const EventInfo &event) override
update event number and writes intermediate particles to a tree.
Definition: rootoutput.cc:347
smash::RootOutput::pdgcode_
std::vector< int > pdgcode_
Property that is written to ROOT output.
Definition: rootoutput.h:199
smash::Action::get_total_weight
virtual double get_total_weight() const =0
Return the total weight value, which is mainly used for the weight output entry.
smash::RootOutput::npart_
int npart_
Property that is written to ROOT output.
Definition: rootoutput.h:206
smash::LHyperSurfaceCrossing
static constexpr int LHyperSurfaceCrossing
Definition: hypersurfacecrossingaction.cc:16
smash::RootOutput::modus_l_
double modus_l_
Property that is written to ROOT output.
Definition: rootoutput.h:207
smash::RootOutput::collisions_tree_
TTree * collisions_tree_
TTree for collision output.
Definition: rootoutput.h:151
smash::RootOutput::impact_b_
double impact_b_
Property that is written to ROOT output.
Definition: rootoutput.h:207
smash::RootOutput::write_collisions_
bool write_collisions_
Option to write collisions tree.
Definition: rootoutput.h:213
smash::RootOutput::current_event_
int current_event_
Number of current event.
Definition: rootoutput.h:171
smash::RootOutput::px_
std::vector< double > px_
Property that is written to ROOT output.
Definition: rootoutput.h:187
particles.h
smash::RootOutput::t_
std::vector< double > t_
Property that is written to ROOT output.
Definition: rootoutput.h:190
smash::RootOutput::time_last_coll_
std::vector< double > time_last_coll_
Property that is written to ROOT output.
Definition: rootoutput.h:197
smash::RootOutput::root_out_file_
std::unique_ptr< TFile > root_out_file_
Pointer to root output file.
Definition: rootoutput.h:137
smash::Particles
Definition: particles.h:33
smash::EventInfo::empty_event
bool empty_event
True if no collisions happened.
Definition: outputinterface.h:51
smash::RootOutput::particles_tree_
TTree * particles_tree_
TTree for particles output.
Definition: rootoutput.h:144
smash::RootOutput::collisions_to_tree
void collisions_to_tree(const ParticleList &incoming, const ParticleList &outgoing, const double weight, const double partial_weight)
Writes collisions to a tree defined by treename.
Definition: rootoutput.cc:501
smash::RootOutput::part_extended_
const bool part_extended_
Whether extended particle output is on.
Definition: rootoutput.h:236
smash::Action
Definition: action.h:35
smash::RootOutput::nin_
int nin_
Property that is written to ROOT output.
Definition: rootoutput.h:206
smash::RootOutput::at_intermediate_time
void at_intermediate_time(const Particles &particles, const std::unique_ptr< Clock > &clock, const DensityParameters &dens_param, const EventInfo &event) override
Writes intermediate particles to a tree defined by treename, if it is allowed (i.e....
Definition: rootoutput.cc:369
smash::RootOutput::par_wgt_
double par_wgt_
Property that is written to ROOT output.
Definition: rootoutput.h:207
smash::RootOutput::max_buffer_size_
static const int max_buffer_size_
Maximal buffer size.
Definition: rootoutput.h:179
smash::ProcessType::HyperSurfaceCrossing
Hypersurface crossing Particles are removed from the evolution and printed to a separate output to se...
smash::RootOutput::particles_to_tree
void particles_to_tree(T &particles)
Writes particles to a tree defined by treename.
Definition: rootoutput.cc:440
smash::pdg::p
constexpr int p
Proton.
Definition: pdgcode_constants.h:28
smash::RootOutput::particles_only_final_
OutputOnlyFinal particles_only_final_
Print only final particles in the event, no intermediate output.
Definition: rootoutput.h:222
smash::Action::incoming_particles
const ParticleList & incoming_particles() const
Get the list of particles that go into the action.
Definition: action.cc:57
smash::NeedsToWrap::No
smash::RootOutput::E_kinetic_tot_
double E_kinetic_tot_
Property that is written to ROOT output.
Definition: rootoutput.h:208
smash::RootOutput::tcounter_
int tcounter_
Property that is written to ROOT output.
Definition: rootoutput.h:206
smash::RootOutput::ic_extended_
const bool ic_extended_
Whether extended ic output is on.
Definition: rootoutput.h:240
smash::RootOutput::x_
std::vector< double > x_
Property that is written to ROOT output.
Definition: rootoutput.h:191
smash::RootOutput::proc_type_origin_
std::vector< int > proc_type_origin_
Property that is written to ROOT output.
Definition: rootoutput.h:203