IAEAphsp/src/IAEAphspRun.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
0027
0028 #include "IAEAphspRun.hh"
0029
0030 #include "globals.hh"
0031 #include "G4Event.hh"
0032 #include "G4Run.hh"
0033
0034 #include <vector>
0035
0036 #include "G4IAEAphspWriter.hh"
0037 #include "G4IAEAphspWriterStack.hh"
0038
0039
0040 //==============================================================================
0041
0042 IAEAphspRun::IAEAphspRun()
0043 :G4Run()
0044 {
0045   G4cout << "Creating default IAEAphspRun object" << G4endl;
0046 }
0047
0048
0049
0050 //==============================================================================
0051
0052 IAEAphspRun::IAEAphspRun(G4IAEAphspWriterStack* iaeaStack)
0053 :G4Run()
0054 {
0055   G4cout << "Creating IAEAphspRun object with IAEAphspWriterStack" << G4endl;
0056   fIAEAphspWriterStack = iaeaStack;
0057   fIAEAphspWriterStack->PrepareRun();
0058 }
0059
0060
0061
0062 //==============================================================================
0063
0064 IAEAphspRun::~IAEAphspRun()
0065 {
0066   G4cout << "Destroying IAEAphspRun object" << G4endl;
0067
0068   if (fIAEAphspWriterStack)
0069     fIAEAphspWriterStack->ClearRunVectors();  // deletion only in RunAction!
0070
0071   if (fIAEAphspWriter) delete fIAEAphspWriter;
0072 }
0073
0074
0075
0076 //==============================================================================
0077 //  RecordEvent() is a method called at end of event after EndOfEventAction().
0078
0079 void IAEAphspRun::RecordEvent(const G4Event* aEvent)
0080 {
0081   G4Run::RecordEvent(aEvent);   // Mandatory to increment 'numberOfEvent'
0082   // G4cout << "IAEAphspRun: numberOfEvent = " << numberOfEvent << G4endl;
0083   // G4cout << "Event ID = " << aEvent->GetEventID() << G4endl;
0084
0085   if (fIAEAphspWriterStack)
0086     fIAEAphspWriterStack->PrepareNextEvent();
0087 }
0088
0089
0090
0091 //==============================================================================
0092 // Merge info from local IAEAphspRun object to the global IAEAphspRun object
0093
0094 void IAEAphspRun::Merge(const G4Run* aRun)
0095 {
0096   G4cout << "IAEAphspRun::Merge() started" << G4endl;
0097
0098   const IAEAphspRun* localRun = static_cast<const IAEAphspRun*>(aRun);
0099   auto localPhspStack = localRun->GetIAEAphspWriterStack();
0100
0101   if (localPhspStack) {     // only if we have IAEAphsp files
0102     DumpToIAEAphspFiles(localPhspStack);
0103
0104     // Update the number of original histories to all files
0105     const G4int histories = localRun->GetNumberOfEvent();
0106     const size_t nPhsp = localPhspStack->GetZphspVec()->size();
0107     for (size_t jj = 0; jj < nPhsp; jj++)
0108       fIAEAphspWriter->SumOrigHistories(jj, histories);
0109   }
0110
0111   G4Run::Merge(aRun);
0112 }
0113
0114
0115 //==============================================================================
0116 // Dump info contained in local IAEAphspRun object into IAEAphsp files
0117
0118 void IAEAphspRun::DumpToIAEAphspFiles(const G4IAEAphspWriterStack* phspStack)
0119 {
0120   // Get info from vectors in G4IAEAphspWriterStack
0121   if (phspStack) {
0122     auto localPdgMtrx   = phspStack->GetPDGMtrx();
0123     auto localPosMtrx   = phspStack->GetPosMtrx();
0124     auto localMomMtrx   = phspStack->GetMomMtrx();
0125     auto localEneMtrx   = phspStack->GetEneMtrx();
0126     auto localWtMtrx    = phspStack->GetWtMtrx();
0127     auto localNstatMtrx = phspStack->GetNstatMtrx();
0128
0129     const size_t nPhsp = phspStack->GetZphspVec()->size();
0130     // -- DEBUG!!
0131     // G4cout << "IAEAphspRun: This run has " << nPhsp << " phsp planes stored."
0132     //     << G4endl;
0133
0134     if (nPhsp != localPdgMtrx->size() || nPhsp != localPosMtrx->size() ||
0135     nPhsp != localMomMtrx->size() || nPhsp != localEneMtrx->size() ||
0136     nPhsp != localWtMtrx->size()  || nPhsp != localNstatMtrx->size() ) {
0137       G4ExceptionDescription msg;
0138       msg << "Number of zphsp stored != size of vectors storing phsp data."
0139       << " MERGING IGNORED!" << G4endl;
0140       G4Exception("IAEAphspRun::DumpToIAEAphspFiles()",
0141           "IAEAphspRun001", JustWarning, msg);
0142     }
0143     else {
0144       size_t jj = 0; // phsp plane counter
0145       for (const auto& phspPdgVec : *localPdgMtrx) {
0146     size_t nPart = phspPdgVec->size(); // Get number of particles
0147     // -- DEBUG!!
0148     // G4cout << "\tPhsp #" << jj << " stores " << nPart << " particles"
0149     //        << G4endl;
0150
0151     if (nPart != (*localNstatMtrx)[jj]->size() ||
0152         nPart != (*localPosMtrx)[jj]->size() ||
0153         nPart != (*localMomMtrx)[jj]->size() ||
0154         nPart != (*localEneMtrx)[jj]->size() ||
0155         nPart != (*localWtMtrx)[jj]->size() ) {
0156       G4ExceptionDescription msg;
0157       msg << "Number of stored particles does not match in this "
0158           << "thread-local run for phps plane #" << jj
0159           << ". Merging ignored!" << G4endl;
0160       G4Exception("IAEAphspRun::DumpToIAEAphspFiles()",
0161               "IAEAphspRun002", JustWarning, msg);
0162     }
0163     else {
0164       // Everything OK to dump particles into the IAEAphsp output files
0165
0166       // 1. If the G4IAEAphspWriter object was not created yet,
0167       // create it, take data from G4IAEAphspWriterStack and open files
0168       if (!fIAEAphspWriter) {
0169         const G4String namePrefix = phspStack->GetFileName();
0170         fIAEAphspWriter = new G4IAEAphspWriter(namePrefix);
0171         fIAEAphspWriter->SetDataFromIAEAStack(phspStack);
0172         fIAEAphspWriter->OpenIAEAphspOutFiles(this);
0173       }
0174
0175       // 2. Loop over all vectors of this phsp to get the dynamic info
0176       // and write it into the corresponding IAEAphsp file
0177
0178       size_t ii = 0;
0179       for (const auto& pdg : *phspPdgVec) {
0180         G4int nStat = (*(*localNstatMtrx)[jj])[ii];
0181         G4double kinE = (*(*localEneMtrx)[jj])[ii];
0182         G4double wt = (*(*localWtMtrx)[jj])[ii];
0183         G4ThreeVector pos = (*(*localPosMtrx)[jj])[ii];
0184         G4ThreeVector momDir = (*(*localMomMtrx)[jj])[ii];
0185         fIAEAphspWriter->WriteIAEAParticle(jj, nStat, pdg, kinE, wt,
0186                            pos, momDir);
0187         ii++;
0188
0189         // -- DEBUG!!
0190         // G4cout << "\tPDG = " << pdg
0191         //     << "  n_stat = " << (*(*localNstatMtrx)[jj])[ii]
0192         //     << "  pos/cm = " << ((*(*localPosMtrx)[jj])[ii]) /cm
0193         //     << "  momDir = " << (*(*localMomMtrx)[jj])[ii]
0194         //     << "  kinE/MeV = " << ((*(*localEneMtrx)[jj])[ii]) /MeV
0195         //     << "  wt = " << (*(*localWtMtrx)[jj])[ii] << G4endl;
0196       }
0197     }
0198     jj++;
0199       }
0200     }
0201   }
0202   else {
0203     G4ExceptionDescription msg;
0204     msg << "This function is not meant to be called if no "
0205     << "G4IAEAphspWriterStack object has been defined."
0206     << G4endl;
0207     G4Exception("IAEAphspRun::DumpToIAEAphspFiles()",
0208         "IAEAphspRun003", FatalException, msg);
0209   }
0210 }