IAEAphsp/src/G4IAEAphspWriter.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Author: M.A. Cortes-Giraldo, Universidad de Sevilla
0027 //
0028 // History changelog prior creation of this example:
0029 // - 17/10/2009: inheritance removed (not needed)
0030 // - 17/10/2009: version 1.0
0031 // - 07/10/2024: version 2.0 for Geant4 example (MT compliant)
0032 // - 18/10/2025: version 3.0
0033 //   - Creation of IAEASourceIdRegistry for thread-safe source_id assignation
0034 //
0035
0036
0037 #include "G4IAEAphspWriter.hh"
0038
0039 #include "globals.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4Run.hh"
0042 #include "G4Step.hh"
0043 #include "G4Track.hh"
0044
0045 #include "iaea_phsp.h"
0046 #include "IAEASourceIdRegistry.hh"
0047 #include "G4IAEAphspWriterStack.hh"
0048
0049 #include <map>
0050 #include <sstream>
0051 #include <vector>
0052
0053 //==============================================================================
0054
0055 G4IAEAphspWriter::G4IAEAphspWriter(const G4String filename)
0056 {
0057   fFileName = filename;
0058   fZphspVec = new std::vector<G4double>;
0059   fConstVariables = new std::map<G4int, G4double>;
0060   fOrigHistories = new std::vector<G4int>;
0061   fIAEASourcesOpen = 0;
0062   G4cout << "G4IAEAphspWriter object constructed." << G4endl;
0063 }
0064
0065
0066 G4IAEAphspWriter::~G4IAEAphspWriter()
0067 {
0068   if (fZphspVec) delete fZphspVec;
0069   if (fConstVariables) delete fConstVariables;
0070   if (fOrigHistories) delete fOrigHistories;
0071 }
0072
0073
0074 //==============================================================================
0075
0076 void G4IAEAphspWriter::AddZphsp(const G4double zphsp)
0077 {
0078   fZphspVec->push_back(zphsp);
0079   G4cout << "G4IAEAphspWriter: Registered phase-space plane at z = "
0080      << zphsp/cm << " cm." << G4endl;
0081
0082   // fOrigHistories must incorporate a new element (equal size as fZphspVec)
0083   fOrigHistories->push_back(0);
0084 }
0085
0086
0087 //==============================================================================
0088
0089 void G4IAEAphspWriter::SetDataFromIAEAStack(const G4IAEAphspWriterStack* stack)
0090 {
0091   if (!stack) {
0092     G4ExceptionDescription msg;
0093     msg <<  "No G4IAEAphspWriterStack has been constructed!" << G4endl;
0094     G4Exception("G4IAEAphspWriter::SetDataFromIAEAStack()",
0095         "IAEAphspWriter001", FatalException, msg );
0096     return;
0097   }
0098   else {
0099     fFileName = stack->GetFileName();
0100
0101     if (stack->GetZphspVec()->size() > 0) {
0102       (*fZphspVec) = *(stack->GetZphspVec()); // copy objects, not pointers
0103
0104       G4cout << "G4IAEAphspWriter::fFileName = " << fFileName << G4endl;
0105       G4cout << "G4IAEAphspWriter::fZphspVec->size() = "
0106          << fZphspVec->size() << G4endl;
0107
0108       // fOrigHistories must have as many element as fZphspVec
0109       fOrigHistories->assign(fZphspVec->size(), 0);
0110     }
0111     else {
0112       G4ExceptionDescription msg;
0113       msg << "No phsp plane z-coordinate has been defined!" << G4endl;
0114       G4Exception("G4IAEAphspWriter::SetDataFromIAEAStack()",
0115           "IAEAphspWriter002", FatalErrorInArgument, msg );
0116       return;
0117     }
0118   }
0119 }
0120
0121
0122 //==============================================================================
0123
0124 void G4IAEAphspWriter::SetConstVariable(const G4int idx, const G4double val)
0125 {
0126   if (idx < 0 || idx > 6) {
0127     G4cout << "No constant variable applies for index " << idx
0128        << ". Doing nothing." << G4endl;
0129     return;
0130   }
0131
0132   fConstVariables->insert( std::pair<G4int, G4double>(idx, val) );
0133   switch (idx) {
0134   case 0:
0135     G4cout << "Variable 'x' set to constant value " << val << " cm" << G4endl;
0136     break;
0137   case 1:
0138     G4cout << "Variable 'y' set to constant value " << val << " cm" << G4endl;
0139     break;
0140   case 2:
0141     G4cout << "Variable 'z' set to constant value " << val << " cm" << G4endl;
0142     break;
0143   case 3:
0144     G4cout << "Variable 'u' set to constant value " << val << G4endl;
0145     break;
0146   case 4:
0147     G4cout << "Variable 'v' set to constant value " << val << G4endl;
0148     break;
0149   case 5:
0150     G4cout << "Variable 'w' set to constant value " << val << G4endl;
0151     break;
0152   case 6:
0153     G4cout << "Variable 'wt' set to constant value " << val << G4endl;
0154   }
0155 }
0156
0157
0158 //==============================================================================
0159
0160 void G4IAEAphspWriter::WriteIAEAParticle(const size_t idx, const G4int incHist,
0161                      const G4int pdg, const G4double kinE,
0162                      const G4double wt,
0163                      const G4ThreeVector pos,
0164                      const G4ThreeVector momDir)
0165 {
0166   IAEA_I32 partType;
0167   switch(pdg) {
0168   case 22:
0169     partType = 1;  // gamma
0170     break;
0171   case 11:
0172     partType = 2;  // electron
0173     break;
0174   case -11:
0175     partType = 3;  // positron
0176     break;
0177   case 2112:
0178     partType = 4;  // neutron
0179     break;
0180   case 2212:
0181     partType = 5;  // proton
0182     break;
0183   default:
0184     G4ExceptionDescription msg;
0185     msg << "PDG " << pdg
0186     << " is not supported by IAEAphsp format and will not be recorded."
0187     << G4endl;
0188     G4Exception("G4IAEAphspWriter::WriteIAEAParticle()",
0189         "IAEAphspWriter003", JustWarning, msg);
0190     return;
0191   }
0192
0193   IAEA_I32 sourceID = static_cast<IAEA_I32>(idx+fIAEASourcesOpen);
0194
0195   IAEA_I32 nStat = static_cast<IAEA_I32>(incHist);
0196   IAEA_Float energy = static_cast<IAEA_Float>(kinE/MeV);
0197   IAEA_Float weight = static_cast<IAEA_Float>(wt);
0198   IAEA_Float x = static_cast<IAEA_Float>( pos.x()/cm );
0199   IAEA_Float y = static_cast<IAEA_Float>( pos.y()/cm );
0200   IAEA_Float z = static_cast<IAEA_Float>( pos.z()/cm );
0201   IAEA_Float u = static_cast<IAEA_Float>( momDir.x() );
0202   IAEA_Float v = static_cast<IAEA_Float>( momDir.y() );
0203   IAEA_Float w = static_cast<IAEA_Float>( momDir.z() );
0204
0205   // Extra variables
0206   IAEA_Float extraFloat = -1; // no extra floats stored
0207   IAEA_I32 extraInt = nStat;
0208
0209   // And finally store the particle following the IAEA routines
0210   iaea_write_particle(&sourceID, &nStat, &partType,
0211               &energy, &weight,
0212               &x, &y, &z, &u, &v, &w, &extraFloat, &extraInt);
0213 }
0214
0215
0216
0217 //==============================================================================
0218
0219 void G4IAEAphspWriter::WriteIAEAParticle(const G4Step* aStep,
0220                      const G4int zStopIdx)
0221 {
0222   IAEA_I32 sourceID =
0223     static_cast<IAEA_I32>(zStopIdx+fIAEASourcesOpen); // beware
0224   G4double zStop = (*fZphspVec)[zStopIdx];
0225
0226   // The particle type and kinetic energy
0227   // -------------------------------------
0228   const G4Track* aTrack = aStep->GetTrack();
0229   G4int PDGCode = aTrack->GetDefinition()->GetPDGEncoding();
0230   IAEA_I32 partType;
0231   G4double postE = aStep->GetPostStepPoint()->GetKineticEnergy();
0232   G4double preE = aStep->GetPreStepPoint()->GetKineticEnergy();
0233   IAEA_Float kinEnergyMeV;
0234   G4ThreeVector postR = aStep->GetPostStepPoint()->GetPosition();
0235   G4ThreeVector preR = aStep->GetPreStepPoint()->GetPosition();
0236   G4double postZ = postR.z();
0237   G4double preZ = preR.z();
0238
0239   switch(PDGCode) {
0240   case 22:
0241     partType = 1;  // gamma
0242     kinEnergyMeV = static_cast<IAEA_Float>(preE/MeV);
0243     break;
0244   case 11:
0245     partType = 2;  // electron
0246     kinEnergyMeV =
0247       static_cast<IAEA_Float>( (preE+
0248                 (postE-preE)*(zStop-preZ)/(postZ-preZ))/MeV );
0249     break;
0250   case -11:
0251     partType = 3;  // positron
0252     kinEnergyMeV =
0253       static_cast<IAEA_Float>( (preE+
0254                 (postE-preE)*(zStop-preZ)/(postZ-preZ))/MeV );
0255     break;
0256   case 2112:
0257     partType = 4;  // neutron
0258     kinEnergyMeV = static_cast<IAEA_Float>(preE/MeV);
0259     break;
0260   case 2212:
0261     partType = 5;  // proton
0262     kinEnergyMeV =
0263       static_cast<IAEA_Float>( (preE+
0264                 (postE-preE)*(zStop-preZ)/(postZ-preZ))/MeV );
0265     break;
0266   default:
0267     G4String pname = aTrack->GetDefinition()->GetParticleName();
0268     G4String errmsg = "'" + pname + "' is not supported by IAEAphsp format"
0269       + " and will not be recorded.";
0270     G4Exception("G4IAEAphspWriter::WriteIAEAParticle()",
0271         "IAEAphspWriter004", JustWarning, errmsg.c_str() );
0272     return;
0273   }
0274
0275   // Track weight
0276   IAEA_Float wt = static_cast<IAEA_Float>(aTrack->GetWeight());
0277
0278   // Position
0279   G4double postX = postR.x();
0280   G4double preX = preR.x();
0281   G4double postY = postR.y();
0282   G4double preY = preR.y();
0283   IAEA_Float x =
0284     static_cast<IAEA_Float>( (preX+
0285                   (postX-preX)*(zStop-preZ)/(postZ-preZ))/cm );
0286   IAEA_Float y =
0287     static_cast<IAEA_Float>( (preY+
0288                   (postY-preY)*(zStop-preZ)/(postZ-preZ))/cm );
0289   IAEA_Float z =
0290     static_cast<IAEA_Float>( zStop/cm );
0291
0292   // Momentum direction
0293   G4ThreeVector momDir = aStep->GetPreStepPoint()->GetMomentumDirection();
0294   IAEA_Float u = static_cast<IAEA_Float>(momDir.x());
0295   IAEA_Float v = static_cast<IAEA_Float>(momDir.y());
0296   IAEA_Float w = static_cast<IAEA_Float>(momDir.z());
0297
0298   // Extra variables
0299   IAEA_Float extraFloat = -1; // no extra floats stored
0300   // IAEA_I32 extraInt = static_cast<IAEA_I32>((*fIncrNumberVector)[zStopIdx]);
0301   // IAEA_I32 nStat = extraInt;
0302   IAEA_I32 extraInt = -1;
0303   IAEA_I32 nStat = 0;  //MACG Admitted values, this MUST change
0304
0305   // And finally store the particle following the IAEA routines
0306   iaea_write_particle(&sourceID, &nStat, &partType,
0307               &kinEnergyMeV, &wt,
0308               &x, &y, &z, &u, &v, &w, &extraFloat, &extraInt);
0309 }
0310
0311
0312
0313 //==============================================================================
0314
0315 void G4IAEAphspWriter::OpenIAEAphspOutFiles(const G4Run* aRun)
0316 {
0317   // Open all the files intended to store
0318   // the phase spaces following the IAEA format.
0319
0320   const IAEA_I32 accessWrite = 2;  // 2 = Writing mode in IAEA routines
0321
0322   size_t nZphsps = fZphspVec->size();
0323   for (size_t ii = 0; ii < nZphsps; ii++) {
0324     // Set the source ID and file name in a unique way
0325     std::stringstream sstr;
0326     sstr << ((*fZphspVec)[ii]/cm);
0327     G4String zphsp(sstr.str());
0328     G4String fullName = fFileName + "_" + zphsp + "cm";
0329
0330     // This part is only added when running several runs
0331     // during the simulation.
0332     G4int runID = aRun->GetRunID();
0333     if (runID > 0) {
0334       std::stringstream sstr2;
0335       sstr2 << runID;
0336       G4String runIDStr(sstr2.str());
0337       fullName += "_";
0338       fullName += runIDStr;
0339     }
0340
0341
0342     // Create the file to store the IAEA phase space
0343
0344     // Reserve a global ID and request it explicitly
0345     G4int reserved = IAEASourceIdRegistry::Instance().ReserveNextLowest();
0346     if (reserved < 0) {
0347       G4ExceptionDescription ed;
0348       ed << "No free IAEA source IDs available for writer" << G4endl;
0349       G4Exception("G4IAEAphspWriter::OpenIAEAphspOutFiles",
0350           "IAEAphspWriter005",FatalException, ed);
0351     }
0352     IAEA_I32 sourceWrite = static_cast<IAEA_I32>(reserved);
0353     char* filename = const_cast<char*>(fullName.data());
0354     IAEA_I32 result = 0;
0355
0356     iaea_new_source( &sourceWrite, filename, &accessWrite,
0357              &result, fullName.size()+1 );
0358
0359     if (result < 0 || sourceWrite < 0) {
0360       IAEASourceIdRegistry::Instance().Release(reserved);
0361       G4ExceptionDescription ed;
0362       ed << "IAEAphsp output file opening operation failed!" << G4endl;
0363       G4Exception("G4IAEAphspWriter::OpenIAEAphspOutFiles()",
0364           "IAEAphspWriter006", FatalException, ed);
0365     }
0366
0367     // This difference tells about the number of IAEA files already open.
0368     fIAEASourcesOpen = sourceWrite - ii;
0369     G4cout << "G4IAEAphspWriter::OpenOutputIAEAphspFiles() ==> "
0370        << "\"" << fullName << "\"   IAEAphsp id = " << sourceWrite << "."
0371        << G4endl;
0372
0373
0374     // Set the global information and options.
0375
0376     // Set constant variables
0377     std::map<G4int, G4double>::iterator itmap;
0378     for (itmap = fConstVariables->begin();
0379      itmap != fConstVariables->end(); itmap++) {
0380       IAEA_I32 varIdx = static_cast<IAEA_I32>( (*itmap).first );
0381       IAEA_Float varValue = static_cast<IAEA_I32>( (*itmap).second );
0382       iaea_set_constant_variable(&sourceWrite, &varIdx, &varValue);
0383     }
0384     //MACG
0385     // Set constant Z
0386     // IAEA_I32 varIdx = 2;  // 0=x, 1=y, 2=z, 3=u, 4=v, 5=w, 6=wt
0387     // IAEA_Float varValue = static_cast<IAEA_Float>((*fZphspVec)[ii]/cm);
0388     // iaea_set_constant_variable(&sourceWrite, &varIdx, &varValue);
0389
0390     // Extra variables
0391     IAEA_I32 extraFloats = 0;
0392     IAEA_I32 extraInts = 1;
0393     iaea_set_extra_numbers(&sourceWrite, &extraFloats, &extraInts);
0394
0395     // Extra variables types
0396     IAEA_I32 longIdx = 0;
0397     IAEA_I32 longType = 1; // incremental history number
0398     iaea_set_type_extralong_variable(&sourceWrite, &longIdx, &longType);
0399   }
0400 }
0401
0402
0403 //==============================================================================
0404
0405 void G4IAEAphspWriter::CloseIAEAphspOutFiles()
0406 {
0407   // Close the IAEA files
0408   G4int nZphsps = fZphspVec->size();
0409   for (G4int ii = 0; ii < nZphsps; ii++) {
0410     const IAEA_I32 sourceID = static_cast<IAEA_I32>(ii+fIAEASourcesOpen);
0411     IAEA_I64 nEvts = static_cast<IAEA_I64>( fOrigHistories->at(ii) );
0412     iaea_set_total_original_particles(&sourceID, &nEvts);
0413
0414     IAEA_I32 result = 0;
0415     iaea_print_header(&sourceID, &result);
0416     if (result < 0) {
0417       G4Exception("G4IAEAphspWriter::EndOfRunAction()",
0418           "IAEAphspWriter007", JustWarning,
0419           "IAEA phsp source not found");
0420     }
0421
0422     iaea_destroy_source(&sourceID, &result);
0423     if (result > 0) {
0424       IAEASourceIdRegistry::Instance().Release(static_cast<G4int>(sourceID));
0425       G4cout << "Phase-space file at z_phsp = " << (*fZphspVec)[ii]/cm
0426          << " cm (IAEA source id #" << sourceID << ") closed successfully!"
0427          << G4endl << G4endl;
0428     }
0429     else {
0430       G4Exception("G4IAEAphspWriter::EndOfRunAction()",
0431           "IAEAphspWriter008", JustWarning,
0432           "IAEA file not closed properly");
0433     }
0434   }
0435 }