moleculardna/src/PrimaryGeneratorSourceGRASCSV.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 /// \file PrimaryGeneratorSourceGRASCSV.cc
0028 /// \brief Primary Generator source class implementation for GRAS CSV phase space
0029 /// for Molecular DNA simulation
0030
0031 #include "PrimaryGeneratorSourceGRASCSV.hh"
0032
0033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0034
0035 PrimaryGeneratorSourceGRASCSV::PrimaryGeneratorSourceGRASCSV(const G4String& filename)
0036 {
0037   // set a default buffer size
0038   fEndOfFile = false;
0039
0040   // open file
0041   fInputFile.open(filename);
0042   if (fInputFile.is_open()) {
0043     G4cout << "*** Opening particle source file " << filename << " ***" << G4endl;
0044     // find the starting particle source box
0045     G4String line = "";
0046     G4int totalevents = 0;
0047     while (std::getline(fInputFile, line)) {
0048       if (line.find("TOTAL_EVENTS") != std::string::npos) {
0049         // Create a stringstream of the current line
0050         std::stringstream ss(line);
0051         G4String tmp = "";
0052         ss >> tmp >> tmp >> totalevents;
0053       }
0054       if (line.find
0055           ("TWO-STAGE PARTICLE PHASE SPACE OUTPUT FILE: PARTICLE PHASE-SPACE INFORMATION")
0056           != std::string::npos)
0057       {
0058         for (G4int i = 0; i < 18; i++) {
0059           std::getline(fInputFile, line);
0060           // Print the header information
0061           G4cout << line << G4endl;
0062         }
0063         break;
0064       }
0065     }
0066     // set the number of particles
0067     fnParticles = totalevents;
0068   }
0069   else {
0070     G4cout << "*** Warning: can't open particle source file in reading mode ***" << G4endl;
0071   }
0072 }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 PrimaryGeneratorSourceGRASCSV::~PrimaryGeneratorSourceGRASCSV()
0077 {
0078   // close file and clear structures
0079   if (fInputFile.is_open()) {
0080     fInputFile.close();
0081   }
0082   fPrimaryList.clear();
0083 }
0084
0085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0086
0087 G4double PrimaryGeneratorSourceGRASCSV::RecomputeNParticles()
0088 {
0089   // to be implemented
0090   // the variable in CSV has to be recomputed reading all the file
0091   return fnParticles;
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095
0096 Primary* PrimaryGeneratorSourceGRASCSV::GetPrimary()
0097 {
0098   if (fInputFile.peek() == EOF || fEndOfFile == true) return nullptr;
0099   // Read in primaries if the primary list is empty
0100   if (fPrimaryList.size() == 0) {
0101     // Read 100 primaries at a time
0102     G4int num = fBufferSize;
0103     if (num > fnParticles) {
0104       num = fnParticles;
0105     }
0106
0107     for (G4int i = 0; i < num; i++) {
0108       if (fInputFile.peek() == EOF) {
0109         // End of file
0110         fEndOfFile = true;
0111         break;
0112       }
0113       G4String line = "";
0114       std::getline(fInputFile, line);
0115       if (line.find
0116           ("TWO-STAGE PARTICLE PHASE SPACE OUTPUT FILE: PARTICLE PHASE-SPACE INFORMATION")
0117           != std::string::npos)
0118       {
0119         const G4int lines = 18;
0120         for (G4int j = 0; j < lines; j++) {
0121           std::getline(fInputFile, line);
0122         }
0123         i--;
0124         continue;
0125       }
0126       if (line.find("Block") != std::string::npos || line.find("File") != std::string::npos
0127           || line.find("End") != std::string::npos || line.find("*") != std::string::npos)
0128       {
0129         // End of Block, skip
0130         i--;
0131         continue;
0132       }
0133       std::stringstream ss(line);
0134       std::vector<G4String> tempVect;
0135       G4String tempString = "";
0136       while (ss >> tempString)
0137         tempVect.push_back(tempString);
0138
0139       try {
0140         // Define position and momentum direction vectors
0141         G4ThreeVector pos = G4ThreeVector(std::stod(tempVect.at(10)),
0142                                           std::stod(tempVect.at(11)),
0143                                           std::stod(tempVect.at(12)));
0144         G4ThreeVector momDir = G4ThreeVector(std::stod(tempVect.at(16)),
0145                                              std::stod(tempVect.at(17)),
0146                                              std::stod(tempVect.at(18)));
0147         G4int PDGEncoding = std::stoi(tempVect.at(2));
0148         G4double KE = 0;
0149         KE = std::stod(tempVect.at(8));
0150
0151         // Generate primary particle
0152         auto* primary = new Primary();
0153         primary->SetName(PDGEncoding);
0154         primary->SetPosition(pos);
0155         primary->SetMomentumDirection(momDir);
0156         primary->SetEnergy(KE);
0157
0158         // Populate the primaries list
0159         fPrimaryList.push_back(primary);
0160       }
0161       catch (const std::invalid_argument&) {
0162         G4cerr << "Phase Space reading error: invalid_argument" << G4endl;
0163       }
0164       catch (const std::out_of_range&) {
0165         G4cerr << "Phase Space reading error: out_of_range" << G4endl;
0166       }
0167     }
0168   }
0169
0170   // Get first element and delete it so it's not reused
0171   Primary* primary = nullptr;
0172   if (fPrimaryList.size() > 0) {
0173     primary = fPrimaryList.front();
0174     fPrimaryList.pop_front();
0175   }
0176   return primary;
0177 }