IAEAphsp/src/G4IAEAphspWriterStack.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
0027 //
0028 // 2025-08-27: Its objects work in local runs, their mission is to store the
0029 //   info of the particles to be written in the IAEAphsp output files at
0030 //   the end of the run. In other words, they constitute a stack for
0031 //   IAEAphsp particles until the run finishes, when the IAEAphsp file is
0032 //   actually written.
0033 //
0034
0035
0036 #include "G4IAEAphspWriterStack.hh"
0037
0038 #include "globals.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4ThreeVector.hh"
0041 #include "G4Step.hh"
0042 #include "G4IAEAphspWriter.hh"
0043
0044 #include <set>
0045 #include <vector>
0046
0047
0048 //==============================================================================
0049
0050 G4IAEAphspWriterStack::G4IAEAphspWriterStack(const G4String filename)
0051 {
0052   fFileName = filename;
0053   fZphspVec = new std::vector<G4double>;
0054   fIncrNumberVec = new std::vector<G4int>;
0055   fPassingTracksVec = new std::vector< std::set<G4int>* >;
0056   fPDGMtrx = new std::vector< std::vector<G4int>* >;
0057   fPosMtrx = new std::vector< std::vector<G4ThreeVector>* >;
0058   fMomMtrx = new std::vector< std::vector<G4ThreeVector>* >;
0059   fEneMtrx = new std::vector< std::vector<G4double>* >;
0060   fWtMtrx = new std::vector< std::vector<G4double>* >;
0061   fNstatMtrx = new std::vector< std::vector<G4int>* >;
0062
0063   G4cout << "G4IAEAphspWriterStack object constructed for files \""
0064      << fFileName << "_<zphsp>cm.IAEA*\"" << G4endl;
0065 }
0066
0067
0068 //==============================================================================
0069
0070 G4IAEAphspWriterStack::~G4IAEAphspWriterStack()
0071 {
0072   if (fZphspVec) delete fZphspVec;
0073   if (fIncrNumberVec)    delete fIncrNumberVec;
0074   if (fPassingTracksVec) delete fPassingTracksVec;
0075   if (fPDGMtrx)    delete fPDGMtrx;
0076   if (fPosMtrx)    delete fPosMtrx;
0077   if (fMomMtrx)    delete fMomMtrx;
0078   if (fEneMtrx)    delete fEneMtrx;
0079   if (fWtMtrx)     delete fWtMtrx;
0080   if (fNstatMtrx)  delete fNstatMtrx;
0081 }
0082
0083
0084 //==============================================================================
0085
0086 void G4IAEAphspWriterStack::AddZphsp(const G4double zphsp)
0087 {
0088   fZphspVec->push_back(zphsp);
0089   G4cout << "G4IAEAphspWriterStack: Registered phase-space plane at z = "
0090      << zphsp/cm << " cm." << G4endl;
0091 }
0092
0093
0094 //==============================================================================
0095
0096 void G4IAEAphspWriterStack::ClearZphspVec()
0097 {
0098   G4cout << "G4IAEAphspWriterStack: Removing all registered phase-space planes!"
0099      << G4endl;
0100   fZphspVec->clear();
0101 }
0102
0103
0104 //==============================================================================
0105
0106 void G4IAEAphspWriterStack::SetDataFromWriter(const G4IAEAphspWriter* writer)
0107 {
0108   if (!writer) {
0109     G4ExceptionDescription msg;
0110     msg <<  "No G4IAEAphspWriter has been constructed!" << G4endl;
0111     G4Exception("G4IAEAphspWriterStack::SetDataFromWriter()",
0112         "IAEAphspWriterStack001", FatalException, msg );
0113     return;
0114   }
0115   else {
0116     fFileName = writer->GetFileName();
0117
0118     if (writer->GetZphspVec()->size() > 0) {
0119       (*fZphspVec) = *(writer->GetZphspVec()); // copy objects, not pointers
0120
0121       G4cout << "G4IAEAphspWriterStack::fFileName = " << fFileName << G4endl;
0122       G4cout << "G4IAEAphspWriterStack::fZphspVec->size() = "
0123          << fZphspVec->size() << G4endl;
0124     }
0125     else {
0126       G4ExceptionDescription msg;
0127       msg << "No phsp plane z-coordinate has been defined!" << G4endl;
0128       G4Exception("G4IAEAphspWriterStack::SetDataFromWriter()",
0129           "IAEAphspWriterStack002", FatalErrorInArgument, msg );
0130       return;
0131     }
0132   }
0133 }
0134
0135
0136 //==============================================================================
0137
0138 void G4IAEAphspWriterStack::PrepareRun()
0139 {
0140   size_t nZphsps = fZphspVec->size();
0141
0142   fIncrNumberVec->reserve(nZphsps);
0143   fPassingTracksVec->reserve(nZphsps);
0144   fPDGMtrx->reserve(nZphsps);
0145   fPosMtrx->reserve(nZphsps);
0146   fMomMtrx->reserve(nZphsps);
0147   fEneMtrx->reserve(nZphsps);
0148   fWtMtrx->reserve(nZphsps);
0149   fNstatMtrx->reserve(nZphsps);
0150
0151   for (size_t ii = 0; ii < nZphsps; ii++) {
0152     fIncrNumberVec->push_back(0);
0153
0154     auto aSet = new std::set<G4int>;
0155     fPassingTracksVec->push_back(aSet);
0156     auto pdgVec = new std::vector<G4int>;
0157     fPDGMtrx->push_back(pdgVec);
0158     auto posVec = new std::vector<G4ThreeVector>;
0159     fPosMtrx->push_back(posVec);
0160     auto momVec = new std::vector<G4ThreeVector>;
0161     fMomMtrx->push_back(momVec);
0162     auto eneVec = new std::vector<G4double>;
0163     fEneMtrx->push_back(eneVec);
0164     auto wtVec = new std::vector<G4double>;
0165     fWtMtrx->push_back(wtVec);
0166     auto nstatVec = new std::vector<G4int>;
0167     fNstatMtrx->push_back(nstatVec);
0168   }
0169   G4cout << "G4IAEAphspWriterStack::PrepareRun() done!" << G4endl;
0170 }
0171
0172
0173 //==============================================================================
0174
0175 void G4IAEAphspWriterStack::PrepareNextEvent()
0176 {
0177   // Update all the incremental history numbers.
0178   for ( auto& ii : (*fIncrNumberVec) )
0179     ii++;
0180
0181   // Remove the track ID's stored during this event.
0182   for ( auto& trackIDs : (*fPassingTracksVec) )
0183     trackIDs->clear();
0184
0185   // -- DEBUG!!
0186   // G4cout << "G4IAEAphspWriterStack ready for the next event!" << G4endl;
0187 }
0188
0189
0190
0191 //==============================================================================
0192
0193 void G4IAEAphspWriterStack::StoreParticleIfEligible(const G4Step* aStep)
0194 {
0195   const G4ThreeVector postR = aStep->GetPostStepPoint()->GetPosition();
0196   const G4ThreeVector preR = aStep->GetPreStepPoint()->GetPosition();
0197   const G4double postZ = postR.z();
0198   const G4double preZ = preR.z();
0199
0200   // Check what phsp planes are being crossed
0201   size_t phspIdx = 0;
0202   for (const auto& phspZ : (*fZphspVec) ) {
0203     if ( (postZ-phspZ)*(preZ-phspZ) < 0 ) {
0204       // Get trackID and check if it is already in fPassingTracksVec[ii]
0205       const G4int trackID = aStep->GetTrack()->GetTrackID();
0206       std::set<G4int>::iterator is;
0207       is = (*fPassingTracksVec)[phspIdx]->find(trackID);
0208
0209       if ( is == (*fPassingTracksVec)[phspIdx]->end() ) {
0210     // This particle has not crossed this phsp plane before.
0211     // Then, put it on the stack if it is of a type foreseen
0212     // by the IAEAphsp format
0213     const G4int pdgCode =
0214       aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0215     if (pdgCode == 22 || pdgCode == 11 || pdgCode == -11 ||
0216         pdgCode == 2112 || pdgCode == 2212)
0217       StoreIAEAParticle(aStep, phspIdx, pdgCode);
0218       }
0219     }
0220     phspIdx++;
0221   }
0222 }
0223
0224
0225 //==============================================================================
0226
0227 void G4IAEAphspWriterStack::StoreIAEAParticle(const G4Step* aStep,
0228                           const G4int phspIndex,
0229                           const G4int pdgCode)
0230 {
0231   const G4double zStop = (*fZphspVec)[phspIndex];
0232   const G4Track* aTrack = aStep->GetTrack();
0233
0234   // Get step info
0235   // --------------------------------
0236   const G4ThreeVector postR = aStep->GetPostStepPoint()->GetPosition();
0237   const G4ThreeVector preR = aStep->GetPreStepPoint()->GetPosition();
0238   const G4double postZ = postR.z();
0239   const G4double preZ = preR.z();
0240
0241   // Set kinetic energy
0242   G4double kinEnergy;
0243   if (pdgCode == 22 || pdgCode == 2112) {  // gamma or neutron
0244     kinEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
0245   }
0246   else if (pdgCode == 11 || pdgCode == -11 ||
0247        pdgCode == 2212 ) {  // electron, positron or proton
0248     const G4double postE = aStep->GetPostStepPoint()->GetKineticEnergy();
0249     const G4double preE = aStep->GetPreStepPoint()->GetKineticEnergy();
0250     kinEnergy = preE + (postE-preE)*(zStop-preZ)/(postZ-preZ);
0251   }
0252   else { // not a particle for the IAEA format
0253     G4ExceptionDescription ED;
0254     ED << "\"" << aTrack->GetDefinition()->GetParticleName()
0255        << "\" is not supported by the IAEAphsp format; not recorded.";
0256     G4Exception("G4IAEAphspWriterStack::StoreIAEAParticle()",
0257         "IAEAphspWriterStack003", JustWarning, ED );
0258     return;
0259   }
0260
0261   // Position
0262   const G4ThreeVector phspPos = preR + (postR-preR)*(zStop-preZ)/(postZ-preZ);
0263
0264   // Momentum direction
0265   const G4ThreeVector phspMomDir =
0266     aStep->GetPreStepPoint()->GetMomentumDirection();
0267
0268   // Track weight
0269   const G4double wt = aTrack->GetWeight();
0270
0271   // n_stat value
0272   const G4int nStat = (*fIncrNumberVec)[phspIndex];
0273
0274   // Store info in stacking vectors
0275   // ------------------------------
0276   ((*fPDGMtrx)[phspIndex])->push_back(pdgCode);
0277   ((*fNstatMtrx)[phspIndex])->push_back(nStat);
0278   ((*fPosMtrx)[phspIndex])->push_back(phspPos);
0279   ((*fMomMtrx)[phspIndex])->push_back(phspMomDir);
0280   ((*fEneMtrx)[phspIndex])->push_back(kinEnergy);
0281   ((*fWtMtrx)[phspIndex])->push_back(wt);
0282
0283   // Once stored, reset the incremental history number (n_stat = 0)
0284   (*fIncrNumberVec)[phspIndex] = 0;
0285
0286   // And now register this trackID to protect against multiple crossers
0287   (*fPassingTracksVec)[phspIndex]->insert( aTrack->GetTrackID() );
0288
0289   // -- DEBUG!!
0290   // G4cout << "G4IAEAphspWriterStack: Particle stored in phsp plane ["
0291   //     << phspIndex << "] at place #" << ((*fPDGMtrx)[phspIndex])->size()
0292   //     << " with the following values:" << G4endl;
0293   // G4cout << "\tPDG = " << pdgCode << "   nStat = " << nStat
0294   //     << "   phspPos = " << phspPos << "   phspMomDir = " << phspMomDir
0295   //     << "   kinEnergy = " << kinEnergy << "   weight = " << wt
0296   //     << G4endl;
0297
0298 }
0299
0300
0301 //==============================================================================
0302
0303 void G4IAEAphspWriterStack::ClearRunVectors()
0304 {
0305   // Clear the run vectors at run termination
0306
0307   fIncrNumberVec->clear();
0308
0309   for (auto& trackIDs : (*fPassingTracksVec) ) delete trackIDs;
0310   fPassingTracksVec->clear();
0311
0312   for (auto& pdgVec : (*fPDGMtrx) ) delete pdgVec;
0313   fPDGMtrx->clear();
0314
0315   for (auto& posVec : (*fPosMtrx) ) delete posVec;
0316   fPosMtrx->clear();
0317
0318   for (auto& momVec : (*fMomMtrx) ) delete momVec;
0319   fMomMtrx->clear();
0320
0321   for (auto& eneVec : (*fEneMtrx) ) delete eneVec;
0322   fEneMtrx->clear();
0323
0324   for (auto& wtVec : (*fWtMtrx) ) delete wtVec;
0325   fWtMtrx->clear();
0326
0327   for (auto& nStatVec : (*fNstatMtrx) ) delete nStatVec;
0328   fNstatMtrx->clear();
0329
0330   G4cout << "G4IAEAphspWriterStack run vectors cleaned!" << G4endl;
0331 }