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