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