Version: SMASH-1.8
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 
87 RootOutput::RootOutput(const bf::path &path, const std::string &name,
88  const OutputParameters &out_par)
89  : OutputInterface(name),
90  filename_(path / (name + ".root")),
91  write_collisions_(name == "Collisions" || name == "Dileptons" ||
92  name == "Photons"),
93  write_particles_(name == "Particles"),
94  write_initial_conditions_(name == "SMASH_IC"),
95  particles_only_final_(out_par.part_only_final),
96  autosave_frequency_(1000),
97  part_extended_(out_par.part_extended),
98  coll_extended_(out_par.coll_extended),
99  ic_extended_(out_par.ic_extended) {
100  filename_unfinished_ = filename_;
101  filename_unfinished_ += ".unfinished";
102  root_out_file_ =
103  make_unique<TFile>(filename_unfinished_.native().c_str(), "NEW");
104  init_trees();
105 }
106 
107 void RootOutput::init_trees() {
108  if (write_particles_ || write_initial_conditions_) {
109  particles_tree_ = new TTree("particles", "particles");
110 
111  particles_tree_->Branch("npart", &npart, "npart/I");
112  particles_tree_->Branch("impact_b", &impact_b, "impact_b/D");
113  particles_tree_->Branch("empty_event", &empty_event_, "empty_event/O");
114  particles_tree_->Branch("ev", &ev, "ev/I");
115  particles_tree_->Branch("tcounter", &tcounter, "tcounter/I");
116 
117  particles_tree_->Branch("pdgcode", &pdgcode[0], "pdgcode[npart]/I");
118  particles_tree_->Branch("charge", &charge[0], "charge[npart]/I");
119 
120  particles_tree_->Branch("p0", &p0[0], "p0[npart]/D");
121  particles_tree_->Branch("px", &px[0], "px[npart]/D");
122  particles_tree_->Branch("py", &py[0], "py[npart]/D");
123  particles_tree_->Branch("pz", &pz[0], "pz[npart]/D");
124 
125  particles_tree_->Branch("t", &t[0], "t[npart]/D");
126  particles_tree_->Branch("x", &x[0], "x[npart]/D");
127  particles_tree_->Branch("y", &y[0], "y[npart]/D");
128  particles_tree_->Branch("z", &z[0], "z[npart]/D");
129 
130  if (part_extended_ || ic_extended_) {
131  particles_tree_->Branch("ncoll", &coll_per_part_[0], "ncoll[npart]/I");
132  particles_tree_->Branch("form_time", &formation_time_[0],
133  "form_time[npart]/D");
134  particles_tree_->Branch("xsecfac", &xsec_factor_[0], "xsecfac[npart]/D");
135  particles_tree_->Branch("proc_id_origin", &proc_id_origin_[0],
136  "proc_id_origin[npart]/I");
137  particles_tree_->Branch("proc_type_origin", &proc_type_origin_[0],
138  "proc_type_origin[npart]/I");
139  particles_tree_->Branch("time_last_coll", &time_last_coll_[0],
140  "time_last_coll[npart]/D");
141  particles_tree_->Branch("pdg_mother1", &pdg_mother1_[0],
142  "pdg_mother1[npart]/I");
143  particles_tree_->Branch("pdg_mother2", &pdg_mother2_[0],
144  "pdg_mother2[npart]/I");
145  }
146  }
147 
148  if (write_collisions_) {
149  collisions_tree_ = new TTree("collisions", "collisions");
150 
151  collisions_tree_->Branch("nin", &nin, "nin/I");
152  collisions_tree_->Branch("nout", &nout, "nout/I");
153  collisions_tree_->Branch("npart", &npart, "npart/I");
154  collisions_tree_->Branch("ev", &ev, "ev/I");
155  collisions_tree_->Branch("weight", &wgt, "weight/D");
156  collisions_tree_->Branch("partial_weight", &par_wgt, "partial_weight/D");
157 
158  collisions_tree_->Branch("pdgcode", &pdgcode[0], "pdgcode[npart]/I");
159  collisions_tree_->Branch("charge", &charge[0], "charge[npart]/I");
160 
161  collisions_tree_->Branch("p0", &p0[0], "p0[npart]/D");
162  collisions_tree_->Branch("px", &px[0], "px[npart]/D");
163  collisions_tree_->Branch("py", &py[0], "py[npart]/D");
164  collisions_tree_->Branch("pz", &pz[0], "pz[npart]/D");
165 
166  collisions_tree_->Branch("t", &t[0], "t[npart]/D");
167  collisions_tree_->Branch("x", &x[0], "x[npart]/D");
168  collisions_tree_->Branch("y", &y[0], "y[npart]/D");
169  collisions_tree_->Branch("z", &z[0], "z[npart]/D");
170 
171  if (coll_extended_) {
172  collisions_tree_->Branch("ncoll", &coll_per_part_[0], "ncoll[npart]/I");
173  collisions_tree_->Branch("form_time", &formation_time_[0],
174  "form_time[npart]/D");
175  collisions_tree_->Branch("xsecfac", &xsec_factor_[0], "xsecfac[npart]/D");
176  collisions_tree_->Branch("proc_id_origin", &proc_id_origin_[0],
177  "proc_id_origin[npart]/I");
178  collisions_tree_->Branch("proc_type_origin", &proc_type_origin_[0],
179  "proc_type_origin[npart]/I");
180  collisions_tree_->Branch("time_last_coll", &time_last_coll_[0],
181  "time_last_coll[npart]/D");
182  collisions_tree_->Branch("pdg_mother1", &pdg_mother1_[0],
183  "pdg_mother1[npart]/I");
184  collisions_tree_->Branch("pdg_mother2", &pdg_mother2_[0],
185  "pdg_mother2[npart]/I");
186  }
187  }
188 }
189 
194 RootOutput::~RootOutput() {
195  // kOverwrite option prevents from writing extra TKey objects into root file
196  root_out_file_->Write("", TObject::kOverwrite);
197  root_out_file_->Close();
198  bf::rename(filename_unfinished_, filename_);
199 }
200 
201 void RootOutput::at_eventstart(const Particles &particles,
202  const int event_number) {
203  // save event number
204  current_event_ = event_number;
205 
206  if (write_particles_ && particles_only_final_ == OutputOnlyFinal::No) {
207  output_counter_ = 0;
208  // This is to have only one output of positive impact parameter per event
209  impact_b = -1.0;
210  empty_event_ = false;
211  particles_to_tree(particles);
212  output_counter_++;
213  }
214 }
215 
216 void RootOutput::at_intermediate_time(const Particles &particles,
217  const std::unique_ptr<Clock> &,
218  const DensityParameters &) {
219  if (write_particles_ && particles_only_final_ == OutputOnlyFinal::No) {
220  particles_to_tree(particles);
221  output_counter_++;
222  }
223 }
224 
225 void RootOutput::at_eventend(const Particles &particles,
226  const int /*event_number*/,
227  double impact_parameter, bool empty_event) {
228  impact_b = impact_parameter;
229  empty_event_ = empty_event;
230  if (write_particles_ &&
231  !(empty_event && particles_only_final_ == OutputOnlyFinal::IfNotEmpty)) {
232  particles_to_tree(particles);
233  }
234  /* Forced regular dump from operational memory to disk. Very demanding!
235  * If program crashes written data will NOT be lost. */
236  if (current_event_ > 0 && current_event_ % autosave_frequency_ == 0) {
237  if (write_particles_ || write_initial_conditions_) {
238  particles_tree_->AutoSave("SaveSelf");
239  }
240  if (write_collisions_) {
241  collisions_tree_->AutoSave("SaveSelf");
242  }
243  }
244 
245  if (write_initial_conditions_) {
246  // If the runtime is too short some particles might not yet have
247  // reached the hypersurface. Warning is printed.
248  if (particles.size() != 0) {
249  logg[LHyperSurfaceCrossing].warn(
250  "End time might be too small for initial conditions output. "
251  "Hypersurface has not yet been crossed by ",
252  particles.size(), " particle(s).");
253  }
254  }
255 }
256 
257 void RootOutput::at_interaction(const Action &action,
258  const double /*density*/) {
259  if (write_collisions_) {
260  collisions_to_tree(action.incoming_particles(), action.outgoing_particles(),
261  action.get_total_weight(), action.get_partial_weight());
262  }
263 
264  if (write_initial_conditions_ &&
265  action.get_type() == ProcessType::HyperSurfaceCrossing) {
266  particles_to_tree(action.incoming_particles());
267  }
268 }
269 
270 template <typename T>
271 void RootOutput::particles_to_tree(T &particles) {
272  int i = 0;
273 
274  tcounter = output_counter_;
275  ev = current_event_;
276 
277  for (const auto &p : particles) {
278  // Buffer full - flush to tree, else fill with particles
279  if (i >= max_buffer_size_) {
280  npart = max_buffer_size_;
281  i = 0;
282  particles_tree_->Fill();
283  } else {
284  t[i] = p.position().x0();
285  x[i] = p.position().x1();
286  y[i] = p.position().x2();
287  z[i] = p.position().x3();
288 
289  p0[i] = p.momentum().x0();
290  px[i] = p.momentum().x1();
291  py[i] = p.momentum().x2();
292  pz[i] = p.momentum().x3();
293 
294  pdgcode[i] = p.pdgcode().get_decimal();
295  charge[i] = p.type().charge();
296 
297  if (part_extended_ || ic_extended_) {
298  const auto h = p.get_history();
299  formation_time_[i] = p.formation_time();
300  xsec_factor_[i] = p.xsec_scaling_factor();
301  time_last_coll_[i] = h.time_last_collision;
302  coll_per_part_[i] = h.collisions_per_particle;
303  proc_id_origin_[i] = h.id_process;
304  proc_type_origin_[i] = static_cast<int>(h.process_type);
305  pdg_mother1_[i] = h.p1.get_decimal();
306  pdg_mother2_[i] = h.p2.get_decimal();
307  }
308 
309  i++;
310  }
311  }
312  // Flush rest to tree
313  if (i > 0) {
314  npart = i;
315  particles_tree_->Fill();
316  }
317 }
318 
319 void RootOutput::collisions_to_tree(const ParticleList &incoming,
320  const ParticleList &outgoing,
321  const double weight,
322  const double partial_weight) {
323  ev = current_event_;
324  nin = incoming.size();
325  nout = outgoing.size();
326  npart = nin + nout;
327  wgt = weight;
328  par_wgt = partial_weight;
329 
330  int i = 0;
331 
332  /* It is assumed that nin + nout < max_buffer_size_
333  * This is true for any possible reaction for current buffer size: 10000
334  * But if one wants initial/final particles written to collisions
335  * then implementation should be updated. */
336 
337  for (const ParticleList &plist : {incoming, outgoing}) {
338  for (const auto &p : plist) {
339  t[i] = p.position().x0();
340  x[i] = p.position().x1();
341  y[i] = p.position().x2();
342  z[i] = p.position().x3();
343 
344  p0[i] = p.momentum().x0();
345  px[i] = p.momentum().x1();
346  py[i] = p.momentum().x2();
347  pz[i] = p.momentum().x3();
348 
349  pdgcode[i] = p.pdgcode().get_decimal();
350  charge[i] = p.type().charge();
351 
352  if (coll_extended_) {
353  const auto h = p.get_history();
354  formation_time_[i] = p.formation_time();
355  xsec_factor_[i] = p.xsec_scaling_factor();
356  time_last_coll_[i] = h.time_last_collision;
357  coll_per_part_[i] = h.collisions_per_particle;
358  proc_id_origin_[i] = h.id_process;
359  proc_type_origin_[i] = static_cast<int>(h.process_type);
360  pdg_mother1_[i] = h.p1.get_decimal();
361  pdg_mother2_[i] = h.p2.get_decimal();
362  }
363 
364  i++;
365  }
366  }
367 
368  collisions_tree_->Fill();
369 }
370 } // namespace smash
smash::RootOutput::~RootOutput
~RootOutput()
Destructor.
Definition: rootoutput.cc:194
smash::Action::get_type
virtual ProcessType get_type() const
Get the process type.
Definition: action.h:131
smash
Definition: action.h:24
smash::RootOutput::npart
int npart
Property that is written to ROOT output.
Definition: rootoutput.h:202
smash::RootOutput::filename_unfinished_
bf::path filename_unfinished_
Filename of output as long as simulation is still running.
Definition: rootoutput.h:152
smash::RootOutput::pdgcode
std::array< int, max_buffer_size_ > pdgcode
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::coll_extended_
const bool coll_extended_
Whether extended collisions output is on.
Definition: rootoutput.h:233
smash::DensityParameters
A class to pre-calculate and store parameters relevant for density calculation.
Definition: density.h:107
smash::RootOutput::impact_b
double impact_b
Property that is written to ROOT output.
Definition: rootoutput.h:203
smash::RootOutput::ev
int ev
Property that is written to ROOT output.
Definition: rootoutput.h:202
smash::Particles::size
size_t size() const
Definition: particles.h:87
smash::RootOutput::nout
int nout
Property that is written to ROOT output.
Definition: rootoutput.h:202
smash::RootOutput::z
std::array< double, max_buffer_size_ > z
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::y
std::array< double, max_buffer_size_ > y
Property that is written to ROOT output.
Definition: rootoutput.h:198
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::formation_time_
std::array< double, max_buffer_size_ > formation_time_
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::RootOutput
RootOutput(const bf::path &path, const std::string &name, const OutputParameters &out_par)
Construct ROOT output.
Definition: rootoutput.cc:87
smash::RootOutput::filename_
const bf::path filename_
Filename of output.
Definition: rootoutput.h:150
smash::RootOutput::write_initial_conditions_
bool write_initial_conditions_
Option to write particles tree for initial conditions.
Definition: rootoutput.h:214
action.h
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:107
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:257
forwarddeclarations.h
smash::RootOutput::output_counter_
int output_counter_
Number of output in a given event.
Definition: rootoutput.h:186
smash::RootOutput::wgt
double wgt
Property that is written to ROOT output.
Definition: rootoutput.h:203
smash::RootOutput::autosave_frequency_
int autosave_frequency_
Root file cannot be read if it was not properly closed and finalized.
Definition: rootoutput.h:228
smash::RootOutput::empty_event_
bool empty_event_
Property that is written to ROOT output.
Definition: rootoutput.h:204
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::coll_per_part_
std::array< int, max_buffer_size_ > coll_per_part_
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::at_eventend
void at_eventend(const Particles &particles, const int event_number, double impact_parameter, bool empty_event) override
update event number and impact parameter, and writes intermediate particles to a tree.
Definition: rootoutput.cc:225
rootoutput.h
smash::RootOutput::write_particles_
bool write_particles_
Option to write particles tree.
Definition: rootoutput.h:211
smash::OutputInterface
Abstraction of generic output.
Definition: outputinterface.h:35
smash::RootOutput::charge
std::array< int, max_buffer_size_ > charge
Property that is written to ROOT output.
Definition: rootoutput.h:200
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::LHyperSurfaceCrossing
static constexpr int LHyperSurfaceCrossing
Definition: hypersurfacecrossingaction.cc:18
smash::RootOutput::x
std::array< double, max_buffer_size_ > x
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::collisions_tree_
TTree * collisions_tree_
TTree for collision output.
Definition: rootoutput.h:168
smash::RootOutput::px
std::array< double, max_buffer_size_ > px
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::nin
int nin
Property that is written to ROOT output.
Definition: rootoutput.h:202
smash::RootOutput::py
std::array< double, max_buffer_size_ > py
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::write_collisions_
bool write_collisions_
Option to write collisions tree.
Definition: rootoutput.h:208
smash::RootOutput::current_event_
int current_event_
Number of current event.
Definition: rootoutput.h:188
particles.h
smash::RootOutput::xsec_factor_
std::array< double, max_buffer_size_ > xsec_factor_
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::par_wgt
double par_wgt
Property that is written to ROOT output.
Definition: rootoutput.h:203
smash::RootOutput::pdg_mother1_
std::array< int, max_buffer_size_ > pdg_mother1_
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::root_out_file_
std::unique_ptr< TFile > root_out_file_
Pointer to root output file.
Definition: rootoutput.h:154
smash::Particles
Definition: particles.h:33
smash::RootOutput::at_intermediate_time
void at_intermediate_time(const Particles &particles, const std::unique_ptr< Clock > &clock, const DensityParameters &dens_param) override
Writes intermediate particles to a tree defined by treename, if it is allowed (i.e....
Definition: rootoutput.cc:216
smash::RootOutput::particles_tree_
TTree * particles_tree_
TTree for particles output.
Definition: rootoutput.h:161
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:319
smash::RootOutput::p0
std::array< double, max_buffer_size_ > p0
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::part_extended_
const bool part_extended_
Whether extended particle output is on.
Definition: rootoutput.h:231
smash::Action
Definition: action.h:35
smash::RootOutput::max_buffer_size_
static const int max_buffer_size_
Maximal buffer size.
Definition: rootoutput.h:191
smash::ProcessType::HyperSurfaceCrossing
Hypersurface crossing Particles are removed from the evolution and printed to a separate output to se...
smash::RootOutput::time_last_coll_
std::array< double, max_buffer_size_ > time_last_coll_
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::particles_to_tree
void particles_to_tree(T &particles)
Writes particles to a tree defined by treename.
Definition: rootoutput.cc:271
smash::RootOutput::tcounter
int tcounter
Property that is written to ROOT output.
Definition: rootoutput.h:202
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:217
smash::RootOutput::at_eventstart
void at_eventstart(const Particles &particles, const int event_number) override
update event number and writes intermediate particles to a tree.
Definition: rootoutput.cc:201
smash::Action::incoming_particles
const ParticleList & incoming_particles() const
Get the list of particles that go into the action.
Definition: action.cc:59
smash::RootOutput::proc_id_origin_
std::array< int, max_buffer_size_ > proc_id_origin_
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::pz
std::array< double, max_buffer_size_ > pz
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::NeedsToWrap::No
smash::RootOutput::ic_extended_
const bool ic_extended_
Whether extended ic output is on.
Definition: rootoutput.h:235
smash::RootOutput::t
std::array< double, max_buffer_size_ > t
Property that is written to ROOT output.
Definition: rootoutput.h:198
smash::RootOutput::proc_type_origin_
std::array< int, max_buffer_size_ > proc_type_origin_
Property that is written to ROOT output.
Definition: rootoutput.h:200
smash::RootOutput::pdg_mother2_
std::array< int, max_buffer_size_ > pdg_mother2_
Property that is written to ROOT output.
Definition: rootoutput.h:200