Line data Source code
1 : #include "crpropa/module/TextOutput.h"
2 : #include "crpropa/module/ParticleCollector.h"
3 : #include "crpropa/Units.h"
4 : #include "crpropa/Version.h"
5 : #include "crpropa/Random.h"
6 : #include "crpropa/base64.h"
7 :
8 : #include "kiss/string.h"
9 :
10 : #ifdef CRPROPA_HAVE_ZLIB
11 : #include "zlib/izstream.hpp"
12 : #include "zlib/ozstream.hpp"
13 : #endif
14 :
15 : #include <sstream>
16 : #include <cstdio>
17 : #include <cinttypes>
18 : #include <stdexcept>
19 : #include <iostream>
20 :
21 :
22 :
23 : namespace crpropa {
24 :
25 0 : TextOutput::TextOutput() : Output(), out(&std::cout), storeRandomSeeds(false) {
26 0 : }
27 :
28 7 : TextOutput::TextOutput(OutputType outputtype) : Output(outputtype), out(&std::cout), storeRandomSeeds(false) {
29 7 : }
30 :
31 0 : TextOutput::TextOutput(std::ostream &out) : Output(), out(&out), storeRandomSeeds(false) {
32 0 : }
33 :
34 0 : TextOutput::TextOutput(std::ostream &out,
35 0 : OutputType outputtype) : Output(outputtype), out(&out), storeRandomSeeds(false) {
36 0 : }
37 :
38 4 : TextOutput::TextOutput(const std::string &filename) : Output(), outfile(filename.c_str(),
39 2 : std::ios::binary), out(&outfile), filename(
40 4 : filename), storeRandomSeeds(false) {
41 2 : if (!outfile.is_open())
42 2 : throw std::runtime_error(std::string("Cannot create file: ") + filename);
43 3 : if (kiss::ends_with(filename, ".gz"))
44 0 : gzip();
45 3 : }
46 :
47 1 : TextOutput::TextOutput(const std::string &filename,
48 2 : OutputType outputtype) : Output(outputtype), outfile(filename.c_str(),
49 1 : std::ios::binary), out(&outfile), filename(
50 2 : filename), storeRandomSeeds(false) {
51 1 : if (!outfile.is_open())
52 0 : throw std::runtime_error(std::string("Cannot create file: ") + filename);
53 2 : if (kiss::ends_with(filename, ".gz"))
54 0 : gzip();
55 1 : }
56 :
57 9 : void TextOutput::printHeader() const {
58 9 : *out << "#";
59 9 : if (fields.test(TrajectoryLengthColumn))
60 6 : *out << "\tD";
61 9 : if (fields.test(TimeColumn))
62 3 : *out << "\ttime";
63 9 : if (fields.test(RedshiftColumn))
64 2 : *out << "\tz";
65 9 : if (fields.test(SerialNumberColumn))
66 3 : *out << "\tSN";
67 9 : if (fields.test(CurrentIdColumn))
68 8 : *out << "\tID";
69 9 : if (fields.test(CurrentEnergyColumn))
70 8 : *out << "\tE";
71 9 : if (fields.test(CurrentPositionColumn) && oneDimensional)
72 2 : *out << "\tX";
73 9 : if (fields.test(CurrentPositionColumn) && not oneDimensional)
74 3 : *out << "\tX\tY\tZ";
75 9 : if (fields.test(CurrentDirectionColumn) && not oneDimensional)
76 3 : *out << "\tPx\tPy\tPz";
77 9 : if (fields.test(SerialNumberColumn))
78 3 : *out << "\tSN0";
79 9 : if (fields.test(SourceIdColumn))
80 6 : *out << "\tID0";
81 9 : if (fields.test(SourceEnergyColumn))
82 6 : *out << "\tE0";
83 9 : if (fields.test(SourcePositionColumn) && oneDimensional)
84 1 : *out << "\tX0";
85 9 : if (fields.test(SourcePositionColumn) && not oneDimensional)
86 2 : *out << "\tX0\tY0\tZ0";
87 9 : if (fields.test(SourceDirectionColumn) && not oneDimensional)
88 2 : *out << "\tP0x\tP0y\tP0z";
89 9 : if (fields.test(SerialNumberColumn))
90 3 : *out << "\tSN1";
91 9 : if (fields.test(CreatedIdColumn))
92 2 : *out << "\tID1";
93 9 : if (fields.test(CreatedEnergyColumn))
94 2 : *out << "\tE1";
95 9 : if (fields.test(CreatedPositionColumn) && oneDimensional)
96 1 : *out << "\tX1";
97 9 : if (fields.test(CreatedPositionColumn) && not oneDimensional)
98 1 : *out << "\tX1\tY1\tZ1";
99 9 : if (fields.test(CreatedDirectionColumn) && not oneDimensional)
100 1 : *out << "\tP1x\tP1y\tP1z";
101 9 : if (fields.test(WeightColumn))
102 2 : *out << "\tW";
103 9 : if (fields.test(CandidateTagColumn))
104 3 : *out << "\ttag";
105 9 : for(std::vector<Property>::const_iterator iter = properties.begin();
106 10 : iter != properties.end(); ++iter)
107 : {
108 1 : *out << "\t" << (*iter).name;
109 : }
110 :
111 9 : *out << "\n#\n";
112 9 : if (fields.test(TrajectoryLengthColumn))
113 6 : *out << "# D Trajectory length [" << lengthScale / Mpc
114 6 : << " Mpc]\n";
115 9 : if (fields.test(TimeColumn))
116 3 : *out << "# time Time [" << timeScale / Myr
117 3 : << " Myr]\n";
118 9 : if (fields.test(RedshiftColumn))
119 2 : *out << "# z Redshift\n";
120 9 : if (fields.test(SerialNumberColumn))
121 3 : *out << "# SN/SN0/SN1 Serial number. Unique (within this run) id of the particle.\n";
122 1 : if (fields.test(CurrentIdColumn) || fields.test(CreatedIdColumn)
123 10 : || fields.test(SourceIdColumn))
124 8 : *out << "# ID/ID0/ID1 Particle type (PDG MC numbering scheme)\n";
125 1 : if (fields.test(CurrentEnergyColumn) || fields.test(CreatedEnergyColumn)
126 10 : || fields.test(SourceEnergyColumn))
127 8 : *out << "# E/E0/E1 Energy [" << energyScale / EeV << " EeV]\n";
128 4 : if (fields.test(CurrentPositionColumn) || fields.test(CreatedPositionColumn)
129 13 : || fields.test(SourcePositionColumn))
130 5 : *out << "# X/X0/X1... Position [" << lengthScale / Mpc << " Mpc]\n";
131 : if (fields.test(CurrentDirectionColumn)
132 5 : || fields.test(CreatedDirectionColumn)
133 14 : || fields.test(SourceDirectionColumn))
134 4 : *out << "# Px/P0x/P1x... Heading (unit vector of momentum)\n";
135 9 : if (fields.test(WeightColumn))
136 2 : *out << "# W Weights" << " \n";
137 9 : if (fields.test(CandidateTagColumn)) {
138 3 : *out << "# tag Candidate tag can be given by the source feature (user defined tag) or by the following interaction process \n";
139 3 : *out << "#\tES \tElasticScattering \n" << "#\tEPP \tElectronPairProduction \n" << "#\tEMPP\tEMPairProduction\n"
140 : << "#\tEMDP\tEMDoublePairProduction\n" << "#\tEMTP\tEMTripletPairProduction \n" << "#\tEMIC\tEMInverseComptonScattering\n"
141 : << "#\tND \tNuclearDecay\n" << "#\tPD \tPhotoDisintegration\n" << "#\tPPP \tPhotoPionProduction\n" << "#\tSYN \tSynchrotronRadiation\n"
142 3 : << "#\tPRIM/SEC\t primary / secondary particle\n";
143 : }
144 9 : for(std::vector<Property>::const_iterator iter = properties.begin();
145 10 : iter != properties.end(); ++iter)
146 : {
147 2 : *out << "# " << (*iter).name << " " << (*iter).comment << "\n";
148 : }
149 :
150 9 : *out << "# no index = current, 0 = at source, 1 = at point of creation\n#\n";
151 18 : *out << "# CRPropa version: " << g_GIT_DESC << "\n#\n";
152 :
153 9 : if (storeRandomSeeds)
154 : {
155 0 : *out << "# Random seeds:\n";
156 0 : std::vector< std::vector<uint32_t> > seeds = Random::getSeedThreads();
157 :
158 0 : for (size_t i =0; i < seeds.size(); i++)
159 : {
160 0 : std::string encoded_data = Base64::encode((unsigned char*) &seeds[i][0], sizeof(seeds[i][0]) * seeds[i].size() / sizeof(unsigned char));
161 0 : *out << "# Thread " << i << ": ";
162 0 : *out << encoded_data;
163 0 : *out << "\n";
164 : }
165 0 : }
166 9 : }
167 :
168 19 : void TextOutput::process(Candidate *c) const {
169 19 : if (fields.none() && properties.empty())
170 0 : return;
171 :
172 : size_t buffersize = 2048;
173 19 : char buffer[buffersize];
174 : size_t p = 0;
175 :
176 19 : std::locale old_locale = std::locale::global(std::locale::classic());
177 :
178 19 : if (fields.test(TrajectoryLengthColumn))
179 16 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
180 16 : c->getTrajectoryLength() / lengthScale);
181 19 : if (fields.test(TimeColumn))
182 13 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
183 13 : c->getTime() / timeScale);
184 :
185 19 : if (fields.test(RedshiftColumn))
186 12 : p += std::snprintf(buffer + p, buffersize - p, "%1.5E\t", c->getRedshift());
187 :
188 19 : if (fields.test(SerialNumberColumn))
189 13 : p += std::snprintf(buffer + p, buffersize - p, "%10" PRIu64 "\t",
190 : c->getSerialNumber());
191 19 : if (fields.test(CurrentIdColumn))
192 18 : p += std::snprintf(buffer + p, buffersize - p, "%10i\t", c->current.getId());
193 19 : if (fields.test(CurrentEnergyColumn))
194 18 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
195 18 : c->current.getEnergy() / energyScale);
196 19 : if (fields.test(CurrentPositionColumn)) {
197 15 : if (oneDimensional) {
198 2 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
199 2 : c->current.getPosition().x / lengthScale);
200 : } else {
201 13 : const Vector3d pos = c->current.getPosition() / lengthScale;
202 13 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
203 : pos.z);
204 : }
205 : }
206 19 : if (fields.test(CurrentDirectionColumn)) {
207 14 : if (not oneDimensional) {
208 13 : const Vector3d pos = c->current.getDirection();
209 13 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
210 : pos.z);
211 : }
212 : }
213 :
214 19 : if (fields.test(SerialNumberColumn))
215 13 : p += std::snprintf(buffer + p, buffersize - p, "%10" PRIu64 "\t", c->getSourceSerialNumber());
216 19 : if (fields.test(SourceIdColumn))
217 16 : p += std::snprintf(buffer + p, buffersize - p, "%10i\t", c->source.getId());
218 19 : if (fields.test(SourceEnergyColumn))
219 16 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
220 16 : c->source.getEnergy() / energyScale);
221 19 : if (fields.test(SourcePositionColumn)) {
222 13 : if (oneDimensional) {
223 1 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
224 1 : c->source.getPosition().x / lengthScale);
225 : } else {
226 12 : const Vector3d pos = c->source.getPosition() / lengthScale;
227 12 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
228 : pos.z);
229 : }
230 : }
231 19 : if (fields.test(SourceDirectionColumn)) {
232 13 : if (not oneDimensional) {
233 12 : const Vector3d pos = c->source.getDirection();
234 12 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
235 : pos.z);
236 : }
237 :
238 : }
239 :
240 19 : if (fields.test(SerialNumberColumn))
241 13 : p += std::snprintf(buffer + p, buffersize - p, "%10" PRIu64 "\t",
242 : c->getCreatedSerialNumber());
243 19 : if (fields.test(CreatedIdColumn))
244 12 : p += std::snprintf(buffer + p, buffersize - p, "%10i\t", c->created.getId());
245 19 : if (fields.test(CreatedEnergyColumn))
246 12 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
247 12 : c->created.getEnergy() / energyScale);
248 19 : if (fields.test(CreatedPositionColumn)) {
249 12 : if (oneDimensional) {
250 1 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t",
251 1 : c->created.getPosition().x / lengthScale);
252 : } else {
253 11 : const Vector3d pos = c->created.getPosition() / lengthScale;
254 11 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
255 : pos.z);
256 : }
257 : }
258 19 : if (fields.test(CreatedDirectionColumn)) {
259 12 : if (not oneDimensional) {
260 11 : const Vector3d pos = c->created.getDirection();
261 11 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y,
262 : pos.z);
263 : }
264 : }
265 19 : if (fields.test(WeightColumn)) {
266 12 : p += std::snprintf(buffer + p, buffersize - p, "%8.5E\t", c->getWeight());
267 : }
268 19 : if (fields.test(CandidateTagColumn)) {
269 13 : p += std::snprintf(buffer + p, buffersize - p, "%s\t", c->getTagOrigin().c_str());
270 : }
271 :
272 19 : for(std::vector<Output::Property>::const_iterator iter = properties.begin();
273 20 : iter != properties.end(); ++iter) {
274 1 : Variant v;
275 1 : if (c->hasProperty((*iter).name)) {
276 0 : v = c->getProperty((*iter).name);
277 : } else {
278 1 : v = (*iter).defaultValue;
279 : }
280 1 : p += std::snprintf(buffer + p, buffersize - p, "%s", v.toString("\t").c_str());
281 1 : p += std::snprintf(buffer + p, buffersize - p, "\t");
282 1 : }
283 19 : buffer[p - 1] = '\n';
284 :
285 19 : std::locale::global(old_locale);
286 :
287 38 : #pragma omp critical(FileOutput)
288 : {
289 19 : if (count == 0)
290 9 : printHeader();
291 19 : Output::process(c);
292 19 : out->write(buffer, p);
293 : }
294 :
295 38 : }
296 :
297 1 : void TextOutput::load(const std::string &filename, ParticleCollector *collector){
298 :
299 : std::string line;
300 : std::istream *in;
301 1 : std::ifstream infile(filename.c_str());
302 :
303 1 : Output output;
304 1 : double lengthScale = output.getLengthScale();
305 1 : double timeScale = output.getTimeScale();
306 1 : double energyScale = output.getEnergyScale();
307 :
308 1 : if (!infile.good())
309 0 : throw std::runtime_error("crpropa::TextOutput: could not open file " + filename);
310 : in = &infile;
311 :
312 2 : if (kiss::ends_with(filename, ".gz")){
313 : #ifdef CRPROPA_HAVE_ZLIB
314 0 : in = new zstream::igzstream(*in);
315 : #else
316 : throw std::runtime_error("CRPropa was built without Zlib compression!");
317 : #endif
318 : }
319 :
320 39 : while (std::getline(*in, line)) {
321 38 : std::stringstream stream(line);
322 38 : if (stream.peek() == '#')
323 : continue;
324 :
325 11 : ref_ptr<Candidate> c = new Candidate();
326 : double val_d; int val_i;
327 : double x, y, z;
328 : stream >> val_d;
329 11 : c->setTrajectoryLength(val_d * lengthScale); // D
330 : stream >> val_d;
331 11 : c->setTime(val_d * timeScale); // time
332 : stream >> val_d;
333 11 : c->setRedshift(val_d); // z
334 11 : stream >> val_i;
335 11 : c->setSerialNumber(val_i); // SN
336 11 : stream >> val_i;
337 11 : c->current.setId(val_i); // ID
338 : stream >> val_d;
339 11 : c->current.setEnergy(val_d * energyScale); // E
340 : stream >> x >> y >> z;
341 11 : c->current.setPosition(Vector3d(x, y, z) * lengthScale); // X, Y, Z
342 : stream >> x >> y >> z;
343 11 : c->current.setDirection(Vector3d(x, y, z) * lengthScale); // Px, Py, Pz
344 11 : stream >> val_i; // SN0 (TODO: Reconstruct the parent-child relationship)
345 11 : stream >> val_i;
346 11 : c->source.setId(val_i); // ID0
347 : stream >> val_d;
348 11 : c->source.setEnergy(val_d * energyScale); // E0
349 : stream >> x >> y >> z;
350 11 : c->source.setPosition(Vector3d(x, y, z) * lengthScale); // X0, Y0, Z0
351 : stream >> x >> y >> z;
352 11 : c->source.setDirection(Vector3d(x, y, z) * lengthScale); // P0x, P0y, P0z
353 11 : stream >> val_i; // SN1
354 11 : stream >> val_i;
355 11 : c->created.setId(val_i); // ID1
356 : stream >> val_d;
357 11 : c->created.setEnergy(val_d * energyScale); // E1
358 : stream >> x >> y >> z;
359 11 : c->created.setPosition(Vector3d(x, y, z) * lengthScale); // X1, Y1, Z1
360 : stream >> x >> y >> z;
361 11 : c->created.setDirection(Vector3d(x, y, z) * lengthScale); // P1x, P1y, P1z
362 : stream >> val_d;
363 11 : c->setWeight(val_d); // W
364 :
365 22 : collector->process(c);
366 38 : }
367 1 : infile.close();
368 2 : }
369 :
370 0 : std::string TextOutput::getDescription() const {
371 0 : return "TextOutput";
372 : }
373 :
374 10 : void TextOutput::close() {
375 : #ifdef CRPROPA_HAVE_ZLIB
376 10 : zstream::ogzstream *zs = dynamic_cast<zstream::ogzstream *>(out);
377 10 : if (zs) {
378 0 : zs->close();
379 0 : delete out;
380 0 : out = 0;
381 : }
382 : #endif
383 10 : outfile.flush();
384 10 : }
385 :
386 10 : TextOutput::~TextOutput() {
387 9 : close();
388 10 : }
389 :
390 0 : void TextOutput::gzip() {
391 : #ifdef CRPROPA_HAVE_ZLIB
392 0 : out = new zstream::ogzstream(*out);
393 : #else
394 : throw std::runtime_error("CRPropa was built without Zlib compression!");
395 : #endif
396 0 : }
397 :
398 0 : void TextOutput::dumpIndexList(std::vector<int> indices) {
399 0 : #pragma omp critical(FileOutput)
400 : {
401 0 : std::stringstream ss;
402 0 : ss << "#" << "\t";
403 0 : for (int i = 0; i < indices.size(); i++)
404 0 : ss << indices[i] << "\t";
405 :
406 : const std::string cstr = ss.str();
407 0 : out-> write(cstr.c_str(), cstr.length());
408 0 : }
409 0 : }
410 :
411 : } // namespace crpropa
|
