| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file PrimaryGeneratorSourceGRASCSV.cc
28 /// \brief Primary Generator source class implementation for GRAS CSV phase space
29 /// for Molecular DNA simulation
30
31 #include "PrimaryGeneratorSourceGRASCSV.hh"
32
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
34
35 PrimaryGeneratorSourceGRASCSV::PrimaryGeneratorSourceGRASCSV(const G4String& filename)
36 {
37 // set a default buffer size
38 fEndOfFile = false;
39
40 // open file
41 fInputFile.open(filename);
42 if (fInputFile.is_open()) {
43 G4cout << "*** Opening particle source file " << filename << " ***" << G4endl;
44 // find the starting particle source box
45 G4String line = "";
46 G4int totalevents = 0;
47 while (std::getline(fInputFile, line)) {
48 if (line.find("TOTAL_EVENTS") != std::string::npos) {
49 // Create a stringstream of the current line
50 std::stringstream ss(line);
51 G4String tmp = "";
52 ss >> tmp >> tmp >> totalevents;
53 }
54 if (line.find
55 ("TWO-STAGE PARTICLE PHASE SPACE OUTPUT FILE: PARTICLE PHASE-SPACE INFORMATION")
56 != std::string::npos)
57 {
58 for (G4int i = 0; i < 18; i++) {
59 std::getline(fInputFile, line);
60 // Print the header information
61 G4cout << line << G4endl;
62 }
63 break;
64 }
65 }
66 // set the number of particles
67 fnParticles = totalevents;
68 }
69 else {
70 G4cout << "*** Warning: can't open particle source file in reading mode ***" << G4endl;
71 }
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75
76 PrimaryGeneratorSourceGRASCSV::~PrimaryGeneratorSourceGRASCSV()
77 {
78 // close file and clear structures
79 if (fInputFile.is_open()) {
80 fInputFile.close();
81 }
82 fPrimaryList.clear();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86
87 G4double PrimaryGeneratorSourceGRASCSV::RecomputeNParticles()
88 {
89 // to be implemented
90 // the variable in CSV has to be recomputed reading all the file
91 return fnParticles;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95
96 Primary* PrimaryGeneratorSourceGRASCSV::GetPrimary()
97 {
98 if (fInputFile.peek() == EOF || fEndOfFile == true) return nullptr;
99 // Read in primaries if the primary list is empty
100 if (fPrimaryList.size() == 0) {
101 // Read 100 primaries at a time
102 G4int num = fBufferSize;
103 if (num > fnParticles) {
104 num = fnParticles;
105 }
106
107 for (G4int i = 0; i < num; i++) {
108 if (fInputFile.peek() == EOF) {
109 // End of file
110 fEndOfFile = true;
111 break;
112 }
113 G4String line = "";
114 std::getline(fInputFile, line);
115 if (line.find
116 ("TWO-STAGE PARTICLE PHASE SPACE OUTPUT FILE: PARTICLE PHASE-SPACE INFORMATION")
117 != std::string::npos)
118 {
119 const G4int lines = 18;
120 for (G4int j = 0; j < lines; j++) {
121 std::getline(fInputFile, line);
122 }
123 i--;
124 continue;
125 }
126 if (line.find("Block") != std::string::npos || line.find("File") != std::string::npos
127 || line.find("End") != std::string::npos || line.find("*") != std::string::npos)
128 {
129 // End of Block, skip
130 i--;
131 continue;
132 }
133 std::stringstream ss(line);
134 std::vector<G4String> tempVect;
135 G4String tempString = "";
136 while (ss >> tempString)
137 tempVect.push_back(tempString);
138
139 try {
140 // Define position and momentum direction vectors
141 G4ThreeVector pos = G4ThreeVector(std::stod(tempVect.at(10)),
142 std::stod(tempVect.at(11)),
143 std::stod(tempVect.at(12)));
144 G4ThreeVector momDir = G4ThreeVector(std::stod(tempVect.at(16)),
145 std::stod(tempVect.at(17)),
146 std::stod(tempVect.at(18)));
147 G4int PDGEncoding = std::stoi(tempVect.at(2));
148 G4double KE = 0;
149 KE = std::stod(tempVect.at(8));
150
151 // Generate primary particle
152 auto* primary = new Primary();
153 primary->SetName(PDGEncoding);
154 primary->SetPosition(pos);
155 primary->SetMomentumDirection(momDir);
156 primary->SetEnergy(KE);
157
158 // Populate the primaries list
159 fPrimaryList.push_back(primary);
160 }
161 catch (const std::invalid_argument&) {
162 G4cerr << "Phase Space reading error: invalid_argument" << G4endl;
163 }
164 catch (const std::out_of_range&) {
165 G4cerr << "Phase Space reading error: out_of_range" << G4endl;
166 }
167 }
168 }
169
170 // Get first element and delete it so it's not reused
171 Primary* primary = nullptr;
172 if (fPrimaryList.size() > 0) {
173 primary = fPrimaryList.front();
174 fPrimaryList.pop_front();
175 }
176 return primary;
177 }
178