OPALX (Object Oriented Parallel Accelerator Library for Exascal) MINIorX
OPALX
StatWriter.cpp
Go to the documentation of this file.
1//
2// Class StatWriter
3// This class writes bunch statistics (*.stat).
4//
5// Copyright (c) 2019, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
6// Christof Metzger-Kraus, Open Sourcerer
7// All rights reserved
8//
9// This file is part of OPAL.
10//
11// OPAL is free software: you can redistribute it and/or modify
12// it under the terms of the GNU General Public License as published by
13// the Free Software Foundation, either version 3 of the License, or
14// (at your option) any later version.
15//
16// You should have received a copy of the GNU General Public License
17// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
18//
19#include "StatWriter.h"
20
22#include "PartBunch/PartBunch.h"
23#include "Physics/Units.h"
24#include "Utilities/Timer.h"
25
26#include <sstream>
27
28StatWriter::StatWriter(const std::string& fname, bool restart) : StatBaseWriter(fname, restart) {
29}
30
31void StatWriter::fillHeader(const losses_t& losses) {
32 if (this->hasColumns()) {
33 return;
34 }
35
36 columns_m.addColumn("t", "double", "ns", "Time");
37 columns_m.addColumn("s", "double", "m", "Path length");
38 columns_m.addColumn("numParticles", "long", "1", "Number of Macro Particles");
39 columns_m.addColumn("charge", "double", "1", "Bunch Charge");
40 columns_m.addColumn("energy", "double", "MeV", "Mean Bunch Energy");
41
42 columns_m.addColumn("rms_x", "double", "m", "RMS Beamsize in x");
43 columns_m.addColumn("rms_y", "double", "m", "RMS Beamsize in y");
44 columns_m.addColumn("rms_s", "double", "m", "RMS Beamsize in s");
45
46 columns_m.addColumn("rms_px", "double", "1", "RMS Normalized Momenta in x");
47 columns_m.addColumn("rms_py", "double", "1", "RMS Normalized Momenta in y");
48 columns_m.addColumn("rms_ps", "double", "1", "RMS Normalized Momenta in s");
49
50 columns_m.addColumn("emit_x", "double", "m", "Normalized Emittance x");
51 columns_m.addColumn("emit_y", "double", "m", "Normalized Emittance y");
52 columns_m.addColumn("emit_s", "double", "m", "Normalized Emittance s");
53
54 columns_m.addColumn("mean_x", "double", "m", "Mean Beam Position in x");
55 columns_m.addColumn("mean_y", "double", "m", "Mean Beam Position in y");
56 columns_m.addColumn("mean_s", "double", "m", "Mean Beam Position in s");
57
58 columns_m.addColumn("ref_x", "double", "m", "x coordinate of reference particle in lab cs");
59 columns_m.addColumn("ref_y", "double", "m", "y coordinate of reference particle in lab cs");
60 columns_m.addColumn("ref_z", "double", "m", "z coordinate of reference particle in lab cs");
61
62 columns_m.addColumn("ref_px", "double", "1", "x momentum of reference particle in lab cs");
63 columns_m.addColumn("ref_py", "double", "1", "y momentum of reference particle in lab cs");
64 columns_m.addColumn("ref_pz", "double", "1", "z momentum of reference particle in lab cs");
65
66 columns_m.addColumn("max_x", "double", "m", "Max Beamsize in x");
67 columns_m.addColumn("max_y", "double", "m", "Max Beamsize in y");
68 columns_m.addColumn("max_s", "double", "m", "Max Beamsize in s");
69
70 columns_m.addColumn("xpx", "double", "1", "Correlation xpx");
71 columns_m.addColumn("ypy", "double", "1", "Correlation ypy");
72 columns_m.addColumn("zpz", "double", "1", "Correlation zpz");
73
74 columns_m.addColumn("Dx", "double", "m", "Dispersion in x");
75 columns_m.addColumn("DDx", "double", "1", "Derivative of dispersion in x");
76 columns_m.addColumn("Dy", "double", "m", "Dispersion in y");
77 columns_m.addColumn("DDy", "double", "1", "Derivative of dispersion in y");
78
79 columns_m.addColumn("Bx_ref", "double", "T", "Bx-Field component of ref particle");
80 columns_m.addColumn("By_ref", "double", "T", "By-Field component of ref particle");
81 columns_m.addColumn("Bz_ref", "double", "T", "Bz-Field component of ref particle");
82
83 columns_m.addColumn("Ex_ref", "double", "MV/m", "Ex-Field component of ref particle");
84 columns_m.addColumn("Ey_ref", "double", "MV/m", "Ey-Field component of ref particle");
85 columns_m.addColumn("Ez_ref", "double", "MV/m", "Ez-Field component of ref particle");
86
87 columns_m.addColumn("dE", "double", "MeV", "energy spread of the beam");
88 columns_m.addColumn("dt", "double", "ns", "time step size");
89 columns_m.addColumn("partsOutside", "double", "1", "outside n*sigma of the beam");
90
91 columns_m.addColumn("DebyeLength", "double", "m", "Debye length in the boosted frame");
92 columns_m.addColumn("plasmaParameter", "double", "1", "Plasma parameter that gives no. of particles in a Debye sphere");
93 columns_m.addColumn("temperature", "double", "K", "Temperature of the beam");
94 columns_m.addColumn("rmsDensity", "double", "1", "RMS number density of the beam");
95
96
98 /*
99 if (Options::computePercentiles) {
100 columns_m.addColumn("68_Percentile_x", "double", "m",
101 "68.27 percentile (1 sigma of normal distribution) of x-component of
102 position"); columns_m.addColumn("68_Percentile_y", "double", "m", "68.27 percentile (1 sigma of
103 normal distribution) of y-component of position"); columns_m.addColumn("68_Percentile_z",
104 "double", "m", "68.27 percentile (1 sigma of normal distribution) of z-component of position");
105
106 columns_m.addColumn("95_Percentile_x", "double", "m",
107 "95.45 percentile (2 sigma of normal distribution) of x-component of
108 position"); columns_m.addColumn("95_Percentile_y", "double", "m", "95.45 percentile (2 sigma of
109 normal distribution) of y-component of position"); columns_m.addColumn("95_Percentile_z",
110 "double", "m", "95.45 percentile (2 sigma of normal distribution) of z-component of position");
111
112 columns_m.addColumn("99_Percentile_x", "double", "m",
113 "99.73 percentile (3 sigma of normal distribution) of x-component of
114 position"); columns_m.addColumn("99_Percentile_y", "double", "m", "99.73 percentile (3 sigma of
115 normal distribution) of y-component of position"); columns_m.addColumn("99_Percentile_z",
116 "double", "m", "99.73 percentile (3 sigma of normal distribution) of z-component of position");
117
118 columns_m.addColumn("99_99_Percentile_x", "double", "m",
119 "99.994 percentile (4 sigma of normal distribution) of x-component of
120 position"); columns_m.addColumn("99_99_Percentile_y", "double", "m", "99.994 percentile (4 sigma
121 of normal distribution) of y-component of position"); columns_m.addColumn("99_99_Percentile_z",
122 "double", "m", "99.994 percentile (4 sigma of normal distribution) of z-component of position");
123
124 columns_m.addColumn("normalizedEps68Percentile_x", "double", "m",
125 "x-component of normalized emittance at 68 percentile (1 sigma of normal
126 distribution)"); columns_m.addColumn("normalizedEps68Percentile_y", "double", "m", "y-component
127 of normalized emittance at 68 percentile (1 sigma of normal distribution)");
128 columns_m.addColumn("normalizedEps68Percentile_z", "double", "m",
129 "z-component of normalized emittance at 68 percentile (1 sigma of normal
130 distribution)");
131
132 columns_m.addColumn("normalizedEps95Percentile_x", "double", "m",
133 "x-component of normalized emittance at 95 percentile (2 sigma of normal
134 distribution)"); columns_m.addColumn("normalizedEps95Percentile_y", "double", "m", "y-component
135 of normalized emittance at 95 percentile (2 sigma of normal distribution)");
136 columns_m.addColumn("normalizedEps95Percentile_z", "double", "m",
137 "z-component of normalized emittance at 95 percentile (2 sigma of normal
138 distribution)");
139
140 columns_m.addColumn("normalizedEps99Percentile_x", "double", "m",
141 "x-component of normalized emittance at 99 percentile (3 sigma of normal
142 distribution)"); columns_m.addColumn("normalizedEps99Percentile_y", "double", "m", "y-component
143 of normalized emittance at 99 percentile (3 sigma of normal distribution)");
144 columns_m.addColumn("normalizedEps99Percentile_z", "double", "m",
145 "z-component of normalized emittance at 99 percentile (3 sigma of normal
146 distribution)");
147
148 columns_m.addColumn("normalizedEps99_99Percentile_x", "double", "m",
149 "x-component of normalized emittance at 99.99 percentile (4 sigma of
150 normal distribution)"); columns_m.addColumn("normalizedEps99_99Percentile_y", "double", "m",
151 "y-component of normalized emittance at 99.99 percentile (4 sigma of
152 normal distribution)"); columns_m.addColumn("normalizedEps99_99Percentile_z", "double", "m",
153 "z-component of normalized emittance at 99.99 percentile (4 sigma of
154 normal distribution)");
155 }
156 */
157 if (OpalData::getInstance()->isInOPALCyclMode() && ippl::Comm->size() == 1) {
158 columns_m.addColumn("R0_x", "double", "m", "R0 Particle position in x");
159 columns_m.addColumn("R0_y", "double", "m", "R0 Particle position in y");
160 columns_m.addColumn("R0_s", "double", "m", "R0 Particle position in z");
161
162 columns_m.addColumn("P0_x", "double", "1", "R0 Particle momentum in x");
163 columns_m.addColumn("P0_y", "double", "1", "R0 Particle momentum in y");
164 columns_m.addColumn("P0_s", "double", "1", "R0 Particle momentum in z");
165 }
166
167 if (OpalData::getInstance()->isInOPALCyclMode()) {
168 columns_m.addColumn("halo_x", "double", "1", "Halo in x");
169 columns_m.addColumn("halo_y", "double", "1", "Halo in y");
170 columns_m.addColumn("halo_z", "double", "1", "Halo in z");
171
172 columns_m.addColumn("azimuth", "double", "deg", "Azimuth in global coordinates");
173 }
174
175 for (size_t i = 0; i < losses.size(); ++i) {
176 columns_m.addColumn(losses[i].first, "long", "1", "Number of lost particles in element");
177 }
178
179 if (mode_m == std::ios::app)
180 return;
181
182 OPALTimer::Timer simtimer;
183 std::string dateStr(simtimer.date());
184 std::string timeStr(simtimer.time());
185
186 std::stringstream ss;
187 ss << "Statistics data '" << OpalData::getInstance()->getInputFn() << "' " << dateStr << " "
188 << timeStr;
189
190 this->addDescription(ss.str(), "stat parameters");
191
192 this->addDefaultParameters();
193
194 this->addInfo("ascii", 1);
195}
196
198 PartBunch_t* beam, Vector_t<double, 3> FDext[], const losses_t& losses, const double& azimuth,
199 const size_t npOutside) {
201 std::shared_ptr<ParticleContainer_t> pc = beam->getParticleContainer();
202 double Ekin = pc->getMeanKineticEnergy();
203
204 double pathLength = beam->get_sPos();
205
208
209 double Q = beam->getCharge();
210
211 if (ippl::Comm->rank() != 0) {
212 return;
213 }
214
215 fillHeader(losses);
216
217 this->open();
218
219 this->writeHeader();
220
221 columns_m.addColumnValue("t", beam->getT() * Units::s2ns); // 1
222 columns_m.addColumnValue("s", pathLength); // 2
223 columns_m.addColumnValue("numParticles", beam->getTotalNum()); // 3
224 columns_m.addColumnValue("charge", Q); // 4
225 columns_m.addColumnValue("energy", Ekin); // 5
226
227 columns_m.addColumnValue("rms_x", pc->getRmsR()(0)); // 6
228 columns_m.addColumnValue("rms_y", pc->getRmsR()(1)); // 7
229 columns_m.addColumnValue("rms_s", pc->getRmsR()(2)); // 8
230
231 columns_m.addColumnValue("rms_px", pc->getRmsP()(0)); // 9
232 columns_m.addColumnValue("rms_py", pc->getRmsP()(1)); // 10
233 columns_m.addColumnValue("rms_ps", pc->getRmsP()(2)); // 11
234
235 columns_m.addColumnValue("emit_x", beam->get_norm_emit()(0)); // 12
236 columns_m.addColumnValue("emit_y", beam->get_norm_emit()(1)); // 13
237 columns_m.addColumnValue("emit_s", beam->get_norm_emit()(2)); // 14
238
239 columns_m.addColumnValue("mean_x", pc->getMeanR()(0)); // 15
240 columns_m.addColumnValue("mean_y", pc->getMeanR()(1)); // 16
241 columns_m.addColumnValue("mean_s", pc->getMeanR()(2)); // 17
242
243 columns_m.addColumnValue("ref_x", beam->RefPartR_m(0)); // 18
244 columns_m.addColumnValue("ref_y", beam->RefPartR_m(1)); // 19
245 columns_m.addColumnValue("ref_z", beam->RefPartR_m(2)); // 20
246
247 columns_m.addColumnValue("ref_px", beam->RefPartP_m(0)); // 21
248 columns_m.addColumnValue("ref_py", beam->RefPartP_m(1)); // 22
249 columns_m.addColumnValue("ref_pz", beam->RefPartP_m(2)); // 23
250
251 columns_m.addColumnValue("max_x", pc->getMaxR()(0)); // 24
252 columns_m.addColumnValue("max_y", pc->getMaxR()(1)); // 25
253 columns_m.addColumnValue("max_s", pc->getMaxR()(2)); // 26
254
255 // Write out Courant Snyder parameters.
256 columns_m.addColumnValue("xpx", beam->get_rprms()(0)); // 27
257 columns_m.addColumnValue("ypy", beam->get_rprms()(1)); // 28
258 columns_m.addColumnValue("zpz", beam->get_rprms()(2)); // 29
259
260 // Write out dispersion.
261 columns_m.addColumnValue("Dx", beam->get_Dx()); // 30
262 columns_m.addColumnValue("DDx", beam->get_DDx()); // 31
263 columns_m.addColumnValue("Dy", beam->get_Dy()); // 32
264 columns_m.addColumnValue("DDy", beam->get_DDy()); // 33
265
266 // Write head/reference particle/tail field information.
267 columns_m.addColumnValue("Bx_ref", FDext[0](0)); // 34 B-ref x
268 columns_m.addColumnValue("By_ref", FDext[0](1)); // 35 B-ref y
269 columns_m.addColumnValue("Bz_ref", FDext[0](2)); // 36 B-ref z
270
271 columns_m.addColumnValue("Ex_ref", FDext[1](0)); // 37 E-ref x
272 columns_m.addColumnValue("Ey_ref", FDext[1](1)); // 38 E-ref y
273 columns_m.addColumnValue("Ez_ref", FDext[1](2)); // 39 E-ref z
274
275 columns_m.addColumnValue("dE", beam->getdE()); // 40 dE energy spread
276 columns_m.addColumnValue("dt", beam->getdT() * Units::s2ns); // 41 dt time step size
277 columns_m.addColumnValue("partsOutside", npOutside); // 42 number of particles outside n*sigma
278
279 columns_m.addColumnValue("DebyeLength", beam->get_debyeLength()); // 43 Debye length in the boosted frame
280 columns_m.addColumnValue("plasmaParameter", beam->get_plasmaParameter()); // 43 plasma parameter
281 columns_m.addColumnValue("temperature", beam->get_temperature()); // 44 Temperature
282 columns_m.addColumnValue("rmsDensity", beam->get_rmsDensity()); // 45 RMS number density
283
284 /*
285 if (Options::computePercentiles) {
286 columns_m.addColumnValue("68_Percentile_x", beam->get_68Percentile()[0]);
287 columns_m.addColumnValue("68_Percentile_y", beam->get_68Percentile()[1]);
288 columns_m.addColumnValue("68_Percentile_z", beam->get_68Percentile()[2]);
289
290 columns_m.addColumnValue("95_Percentile_x", beam->get_95Percentile()[0]);
291 columns_m.addColumnValue("95_Percentile_y", beam->get_95Percentile()[1]);
292 columns_m.addColumnValue("95_Percentile_z", beam->get_95Percentile()[2]);
293
294 columns_m.addColumnValue("99_Percentile_x", beam->get_99Percentile()[0]);
295 columns_m.addColumnValue("99_Percentile_y", beam->get_99Percentile()[1]);
296 columns_m.addColumnValue("99_Percentile_z", beam->get_99Percentile()[2]);
297
298 columns_m.addColumnValue("99_99_Percentile_x", beam->get_99_99Percentile()[0]);
299 columns_m.addColumnValue("99_99_Percentile_y", beam->get_99_99Percentile()[1]);
300 columns_m.addColumnValue("99_99_Percentile_z", beam->get_99_99Percentile()[2]);
301
302 columns_m.addColumnValue(
303 "normalizedEps68Percentile_x", beam->get_normalizedEps_68Percentile()[0]);
304 columns_m.addColumnValue(
305 "normalizedEps68Percentile_y", beam->get_normalizedEps_68Percentile()[1]);
306 columns_m.addColumnValue(
307 "normalizedEps68Percentile_z", beam->get_normalizedEps_68Percentile()[2]);
308
309 columns_m.addColumnValue(
310 "normalizedEps95Percentile_x", beam->get_normalizedEps_95Percentile()[0]);
311 columns_m.addColumnValue(
312 "normalizedEps95Percentile_y", beam->get_normalizedEps_95Percentile()[1]);
313 columns_m.addColumnValue(
314 "normalizedEps95Percentile_z", beam->get_normalizedEps_95Percentile()[2]);
315
316 columns_m.addColumnValue(
317 "normalizedEps99Percentile_x", beam->get_normalizedEps_99Percentile()[0]);
318 columns_m.addColumnValue(
319 "normalizedEps99Percentile_y", beam->get_normalizedEps_99Percentile()[1]);
320 columns_m.addColumnValue(
321 "normalizedEps99Percentile_z", beam->get_normalizedEps_99Percentile()[2]);
322
323 columns_m.addColumnValue(
324 "normalizedEps99_99Percentile_x", beam->get_normalizedEps_99_99Percentile()[0]);
325 columns_m.addColumnValue(
326 "normalizedEps99_99Percentile_y", beam->get_normalizedEps_99_99Percentile()[1]);
327 columns_m.addColumnValue(
328 "normalizedEps99_99Percentile_z", beam->get_normalizedEps_99_99Percentile()[2]);
329 }
330 */
331 if (OpalData::getInstance()->isInOPALCyclMode()) {
332 if (ippl::Comm->size() == 1) {
333 if (beam->getLocalNum() > 0) {
334 columns_m.addColumnValue("R0_x", beam->R(0)[0]);
335 columns_m.addColumnValue("R0_y", beam->R(0)[1]);
336 columns_m.addColumnValue("R0_s", beam->R(0)[2]);
337 columns_m.addColumnValue("P0_x", beam->P(0)[0]);
338 columns_m.addColumnValue("P0_y", beam->P(0)[1]);
339 columns_m.addColumnValue("P0_s", beam->P(0)[2]);
340 } else {
341 columns_m.addColumnValue("R0_x", 0.0);
342 columns_m.addColumnValue("R0_y", 0.0);
343 columns_m.addColumnValue("R0_s", 0.0);
344 columns_m.addColumnValue("P0_x", 0.0);
345 columns_m.addColumnValue("P0_y", 0.0);
346 columns_m.addColumnValue("P0_s", 0.0);
347 }
348 }
349 Vector_t<double, 3> halo = beam->get_halo();
350 columns_m.addColumnValue("halo_x", halo(0));
351 columns_m.addColumnValue("halo_y", halo(1));
352 columns_m.addColumnValue("halo_z", halo(2));
353
354 columns_m.addColumnValue("azimuth", azimuth);
355 }
356
357 for (size_t i = 0; i < losses.size(); ++i) {
358 long unsigned int loss = losses[i].second;
359 columns_m.addColumnValue(losses[i].first, loss);
360 }
361
362 this->writeRow();
363
364 this->close();
365}
ParticleContainer< double, 3 > ParticleContainer_t
Definition Component.h:30
PartBunch< PLayout_t< double, 3 >, double, 3 > PartBunch_t
ippl::Vector< T, Dim > Vector_t
constexpr double s2ns
Definition Units.h:44
std::string getInputFn()
get opals input filename
Definition OpalData.cpp:681
static OpalData * getInstance()
Definition OpalData.cpp:195
double get_Dy() const
double get_sPos() const
double getdT() const
double get_plasmaParameter() const
Definition PartBunch.h:726
double get_Dx() const
double get_rmsDensity() const
Definition PartBunch.h:730
double getdE() const
double getT() const
double getCharge() const
get the total charge per simulation particle
Vector_t< T, Dim > get_halo() const
double get_debyeLength() const
Definition PartBunch.h:722
Vector_t< T, Dim > RefPartR_m
std::shared_ptr< ParticleContainer_t > getParticleContainer()
Definition PartBunch.h:248
size_t getLocalNum() const
double get_DDy() const
Vector_t< double, Dim > R(size_t i)
Definition PartBunch.h:276
ParticleAttrib< Vector_t< T, Dim > > P
double get_DDx() const
double get_temperature() const
Definition PartBunch.h:714
size_t getTotalNum() const
Vector_t< T, Dim > RefPartP_m
Vector_t< T, Dim > get_rprms() const
Vector_t< T, Dim > get_norm_emit() const
SDDSColumnSet columns_m
Definition SDDSWriter.h:110
bool hasColumns() const
Definition SDDSWriter.h:171
void addDefaultParameters()
void addDescription(const std::string &text, const std::string &content)
Definition SDDSWriter.h:142
void writeHeader()
Write SDDS header.
std::ios_base::openmode mode_m
First write to the statistics output file.
Definition SDDSWriter.h:108
void writeRow()
Definition SDDSWriter.h:159
void addInfo(const std::string &mode, const size_t &no_row_counts)
Definition SDDSWriter.h:155
StatBaseWriter(const std::string &fname, bool restart)
std::vector< std::pair< std::string, unsigned int > > losses_t
Definition StatWriter.h:26
void write(PartBunch_t *beam, Vector_t< double, 3 > FDext[], const losses_t &losses=losses_t(), const double &azimuth=-1, const size_t npOutside=0)
Write statistical data.
void fillHeader(const losses_t &losses=losses_t())
StatWriter(const std::string &fname, bool restart)
std::string time() const
Return time.
Definition Timer.cpp:42
std::string date() const
Return date.
Definition Timer.cpp:35