| 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 ChemGeoImport.cc
28 /// \brief Implementation of the ChemGeoImport class
29
30 #include "ChemGeoImport.hh"
31 #include "G4Filesystem.hh"
32 #include "G4DNAMolecule.hh"
33
34 ChemGeoImport::ChemGeoImport()
35 {
36 GetVoxelDefFilePathList();
37 fpGun = new UserMoleculeGun();
38 }
39
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41
42 ChemGeoImport::~ChemGeoImport()
43 {
44 if(fpGun)
45 delete fpGun;
46 }
47
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50 void ChemGeoImport::InsertMoleculeInWorld()
51 {
52 // The idea is to add all the molecules specified in the input files
53
54 if(fIsParsed)
55 {
56 // Create the molecules
57
58 // Loop on all the parsed molecules
59 for(G4int i=0, ie=fMolecules.size(); i<ie; i++)
60 {
61 // Retrieve general molecule informations
62 //
63 G4String name = fMolecules[i].fName;
64
65 G4ThreeVector moleculePosition = fMolecules[i].fPosition;
66 G4int copyNum = fMolecules[i].fCopyNumber;
67
68 G4int strand = fMolecules[i].fStrand;
69 ChemMolecule Bmolecule(name, copyNum, moleculePosition, strand, -1, -1, -1);
70
71 if(name=="phosphate1" || name=="phosphate2")
72 {
73 name=G4Phosphate::Definition()->GetName();
74 Bmolecule.fName=name;
75 }
76 else if(name=="deoxyribose1" || name=="deoxyribose2")
77 {
78 name=G4Deoxyribose::Definition()->GetName();
79 Bmolecule.fName=name;
80 }
81 else if(name=="base_adenine")
82 {
83 name=G4Adenine::Definition()->GetName();
84 Bmolecule.fName=name;
85 }
86 else if(name=="base_guanine")
87 {
88 name=G4Guanine::Definition()->GetName();
89 Bmolecule.fName=name;
90 }
91 else if(name=="base_thymine")
92 {
93 name=G4Thymine::Definition()->GetName();
94 Bmolecule.fName=name;
95 }
96 else if(name=="base_cytosine")
97 {
98 name=G4Cytosine::Definition()->GetName();
99 Bmolecule.fName=name;
100 }
101 else if(name=="histone")
102 {
103 name=G4Histone::Definition()->GetName();
104 Bmolecule.fName=name;
105 }
106 else if(name=="solvatedElectron")
107 {
108 name=G4Electron_aq::Definition()->GetName();
109 }
110 else if(name=="water")
111 {
112 name=G4H2O::Definition()->GetName();
113 }
114 else
115 {
116 G4String msg =
117 "The name "+ name+" is not specified in the listed chemical molecules";
118 G4Exception("ChemGeoImport::BuildGeometry", "", FatalException, msg);
119 }
120
121 // Check if the molecule is on the "remove list"
122 G4bool toBeRemoved = IsMoleculeInTheRemoveTable(Bmolecule);
123 if(!toBeRemoved)
124 {
125 // Molecule is not in the "remove list" and we can add it to the simulation
126 // Check the molecule to be added is not a water molecule (special case)
127 if(name != G4H2O::Definition()->GetName() )
128 fpGun->AddMolecule(name, moleculePosition, 1.e-12*s, copyNum, strand);
129 else // Water molecule case
130 fpGun->AddWaterMolecule(moleculePosition, fMolecules.at(i).fTrackId,
131 ElectronicModification(fMolecules.at(i).fState),
132 fMolecules.at(i).fElectronicLevel);
133 }
134
135 }
136 G4DNAChemistryManager::Instance()->SetGun(fpGun);
137 }
138 else
139 {
140 G4String msg =
141 "ChemGeoImport::InsertMoleculeInWorld: The parse method needs to be called first.";
142 G4Exception("ChemGeoImport::ChemGeoImport::InsertMoleculeInWorld", "",FatalException, msg);
143 }
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147
148 void ChemGeoImport::Reset()
149 {
150 // Clear the containers
151 if(fpGun){
152 delete fpGun;
153 fpGun = new UserMoleculeGun();
154 }
155 fMolecules.clear();
156 fMolecules.shrink_to_fit();
157 fToBeRemovedMol.clear();
158 fIsParsed = false;
159 fFactor = 1.;
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163
164 void ChemGeoImport::ParseFiles(const G4String& chemInputFile)
165 {
166 G4fs::path aP{std::string(chemInputFile)};
167 if (G4fs::exists(aP)) {
168 Reset();
169 ParseChemInputFile(chemInputFile);
170 auto geoPathFileName = GetVoxelDefFilePath(fGeoNameFromChemInput);
171
172 ParseGeoFile(geoPathFileName);
173 fIsParsed = true;
174 }
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178
179 void ChemGeoImport::ParseChemInputFile(const G4String& fileName)
180 {
181 // Setup the input stream
182 std::ifstream file;
183 file.open(fileName.c_str() );
184
185 if(!file.good() )
186 {
187 // Geant4 exception
188 G4String msg = fileName+" could not be opened";
189 G4Exception("ChemGeoImport::ParseChemInputFile", "", FatalException, msg);
190 }
191
192 // Define the line string variable
193 G4String line;
194
195 // Read the file line per line
196 while(std::getline(file, line) )
197 {
198 // Check the line to determine if it is empty
199 if(line.empty() )
200 continue; // skip the line if it is empty
201
202 // Data string stream
203 std::istringstream issLine(line);
204
205 // String to determine the first letter/word
206 G4String firstItem;
207
208 // Put the first letter/word within the string
209 issLine >> firstItem;
210
211 // Check first letter to determine if the line is data or comment
212 if(firstItem=="#")
213 continue; // skip the line if it is comment
214
215 else if(firstItem=="_input")
216 {
217 G4int type(-1), state(-1), electronicLevel(-1), parentTrackId(-1);
218 G4double x, y, z;
219 issLine >> type >> state >> electronicLevel;
220 issLine >> x >> y >> z;
221 issLine >> parentTrackId;
222
223 x *= fFactor*nm;
224 y *= fFactor*nm;
225 z *= fFactor*nm;
226
227 G4String name;
228 if(type==1)
229 name="water";
230 else if(type==2)
231 name="solvatedElectron";
232 else
233 {
234 G4ExceptionDescription description;
235 description << "The type " << type <<" is not recognized";
236 G4Exception("ChemGeoImport::ParseFile", "Fatal", FatalException, description, "");
237 }
238
239 ChemMolecule molecule(name, -1, G4ThreeVector(x,y,z), -1,
240 state, electronicLevel, parentTrackId);
241
242 fMolecules.push_back(molecule);
243 }
244
245 else if(firstItem=="_remove")
246 {
247 G4String name;
248 issLine >> name;
249
250 G4int copyNumber;
251 issLine >> copyNumber;
252
253 G4int strand;
254 issLine >> strand;
255
256 fToBeRemovedMol.push_back(ChemMolecule(name,copyNumber,G4ThreeVector(),strand,-1,-1,-1));
257 }
258
259 else if(firstItem=="_eventNum")
260 {
261 // Nothing
262 }
263
264 else if(firstItem=="_voxelType")
265 {
266 issLine >> fGeoNameFromChemInput;
267 }
268
269 else if(firstItem=="_voxelCopyNumber")
270 {
271 // Nothing
272 }
273
274 else if(firstItem=="_Version")
275 {
276 // Nothing
277 }
278
279 else
280 {
281 // Geant4 exception
282 G4String msg =
283 firstItem+" is not defined in the parser. Check the input file: "+fileName+".";
284 G4Exception("ChemGeoImport::ParseChemInputFile", "Geo_WrongParse",FatalException, msg);
285 }
286 }
287 file.close();
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291
292 void ChemGeoImport::ParseGeoFile(const G4String& fileName)
293 {
294 // Setup the input stream
295 std::ifstream file(fileName.c_str());
296
297 // Check if the file was correctly opened
298 if(!file.is_open() )
299 {
300 // Geant4 exception
301 G4String msg = fileName+" could not be opened";
302 G4Exception("ChemGeoImport::ParseGeoFile", "", FatalException, msg);
303 }
304
305 // Define the line string variable
306 G4String line;
307
308 // Read the file line per line
309 while(std::getline(file, line) )
310 {
311 // Check the line to determine if it is empty
312 if(line.empty() )
313 continue; // skip the line if it is empty
314
315 // Data string stream
316 std::istringstream issLine(line);
317
318 // String to determine the first letter/word
319 G4String firstItem;
320
321 // Put the first letter/word within the string
322 issLine >> firstItem;
323
324 // Check first letter to determine if the line is data or comment
325 if(firstItem=="#")
326 continue; // skip the line if it is comment
327
328 // Use the file
329 else if(firstItem=="_Name")
330 {
331 G4String name;
332 issLine >> name;
333 }
334 else if(firstItem=="_Size")
335 {
336 G4double size;
337 issLine >> size;
338 size *= fFactor*nm;
339
340 fSize = size;
341 }
342 else if(firstItem=="_Number")
343 {
344 // Nothing
345 }
346 else if(firstItem=="_Radius")
347 {
348 // Nothing
349 }
350 else if(firstItem=="_Version")
351 {
352 // Nothing
353 }
354 else if(firstItem=="_pl")
355 {
356 G4String name;
357 issLine >> name;
358
359 G4String material;
360 issLine >> material;
361
362 G4int strand;
363 issLine >> strand;
364
365 G4int copyNumber;
366 issLine >> copyNumber;
367
368 G4double x;
369 issLine >> x;
370 x *= fFactor*nm;
371
372 G4double y;
373 issLine >> y;
374 y *= fFactor*nm;
375
376 G4double z;
377 issLine >> z;
378 z *= fFactor*nm;
379
380 ChemMolecule molecule(name, copyNumber, G4ThreeVector(x, y, z), strand, -1, -1, -1);
381
382 fMolecules.push_back(molecule);
383 }
384
385 else
386 {
387 // Geant4 exception
388 G4String msg =
389 firstItem+" is not defined in the parser. Check the input file: "+fileName+".";
390 G4Exception("ChemGeoImport::ParseGeoFile", "Geo_WrongParse", FatalException, msg);
391 }
392 }
393 file.close();
394 }
395
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
397
398 G4bool ChemGeoImport::IsMoleculeInTheRemoveTable(const ChemMolecule& molecule)
399 {
400 if(std::find(fToBeRemovedMol.begin(),fToBeRemovedMol.end(),molecule) != fToBeRemovedMol.end())
401 return true;
402 else
403 return false;
404 }
405
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407
408 G4String ChemGeoImport::GetVoxelDefFilePath(G4String bareName)
409 {
410 G4String strRes = "";
411 for (auto const &entry : fVoxelDefFilesList) {
412 G4fs::path voxelP{std::string(entry)};
413 if (voxelP.stem().string() == bareName) {
414 strRes = entry;
415 }
416 }
417 return strRes;
418 }
419
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
421
422 void ChemGeoImport::GetVoxelDefFilePathList()
423 {
424 G4fs::path thisP = G4fs::current_path();
425 G4bool doesWantedFileExist = false;
426 for (const auto &entry : G4fs::directory_iterator(thisP)){
427 if (entry.path().filename() == "imp.info") {
428 std::ifstream file(entry.path().c_str());
429 if(!file.good() ){
430 G4String msg =
431 "File imp.info is broken. Check its content or try to rerun the PhysicalStage?";
432 G4Exception("ChemGeoImport::GetVoxelDefFilePathList()", "", FatalException, msg);
433 }
434 doesWantedFileExist = true;
435 G4String line;
436 while(std::getline(file, line) ){
437 std::istringstream iss(line);
438 G4String flag;
439 G4String voxelDefFile;
440 iss >> flag;
441 if ( flag == "_geovolxelpath") {
442 iss >> voxelDefFile;
443 fVoxelDefFilesList.insert(voxelDefFile);
444 }
445 }
446 file.close();
447 }
448 }
449
450 if (!doesWantedFileExist) {
451 G4String msg = "File imp.info does not exist. Did you run the Physical Stage?";
452 G4Exception("ChemGeoImport::GetVoxelDefFilePathList()", "", FatalException, msg);
453 }
454 }
455
456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......