| 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 // Code developed by:
28 // S.Guatelli, M. Large and A. Malaroda, University of Wollongong
29 //
30 //Original code from geant4/examples/extended/runAndEvent/RE03
31 //
32 #include <vector>
33 #include <map>
34 #include "ICRP110UserScoreWriter.hh"
35 #include "ICRP110ScoreWriterMessenger.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4SDParticleFilter.hh"
38 #include "G4VPrimitiveScorer.hh"
39 #include "G4VScoringMesh.hh"
40
41 ICRP110UserScoreWriter::ICRP110UserScoreWriter():
42 G4VScoreWriter()
43 {
44 fMessenger = new ICRP110ScoreWriterMessenger(this);
45 fSex = "female"; //Default phantom sex is female
46 fSection = "head"; // Default phantom section is head
47 }
48
49 ICRP110UserScoreWriter::~ICRP110UserScoreWriter()
50 {
51 delete fMessenger;
52 }
53
54 void ICRP110UserScoreWriter::DumpQuantityToFile(const G4String & psName, const G4String & fileName, const G4String & option)
55 {
56 using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
57
58 if(verboseLevel > 0)
59 {
60 G4cout << "ICRP110UserScorer-defined DumpQuantityToFile() method is invoked." << G4endl;
61 }
62
63 // change the option string into lowercase to the case-insensitive.
64 G4String opt = option;
65 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
66
67 // confirm the option
68 if(opt.size() == 0) opt = "csv";
69
70 //--------------------------------------------------------------------//
71 //----------------Create Scoring Mesh Output Text File----------------//
72 //--------------------------------------------------------------------//
73 // First we create use the scoring mesh to create a default output text
74 // file containing 4 columns: voxel number along x, y, z, and edep deposited
75 // in that voxel (in J). This file is to be called "PhantomMesh_Edep.txt".
76
77 std::ofstream ofile(fileName);
78
79 if(!ofile)
80 {
81 G4cerr << "ERROR : DumpToFile : File open error -> " << fileName << G4endl;
82 return;
83 }
84 ofile << "# mesh name: " << fScoringMesh -> GetWorldName() << G4endl;
85
86 // retrieve the map
87 MeshScoreMap fSMap = fScoringMesh -> GetScoreMap();
88
89 auto msMapItr = fSMap.find(psName);
90
91 if(msMapItr == fSMap.end())
92 {
93 G4cerr << "ERROR : DumpToFile : Unknown quantity, \""<< psName
94 << "\"." << G4endl;
95 return;
96 }
97
98 std::map<G4int, G4StatDouble*> * score = msMapItr -> second-> GetMap();
99
100 ofile << "# primitive scorer name: " << msMapItr -> first << G4endl;
101
102 // declare dose array and initialize to zero.
103 std::vector<double> ScoringMeshEdep;
104 for(G4int y = 0; y < fNMeshSegments[0]*fNMeshSegments[1]*fNMeshSegments[2]; y++) ScoringMeshEdep.push_back(0.);
105
106 ofile << std::setprecision(16); // for double value with 8 bytes
107
108 for(G4int x = 0; x < fNMeshSegments[0]; x++) {
109 for(G4int y = 0; y < fNMeshSegments[1]; y++) {
110 for(G4int z = 0; z < fNMeshSegments[2]; z++){
111 // Retrieve dose in each scoring mesh bin/voxel
112 G4int idx = GetIndex(x, y, z);
113 std::map<G4int, G4StatDouble*>::iterator value = score -> find(idx);
114 if (value != score -> end()) ScoringMeshEdep[idx] += (value->second->sum_wx())/(joule);
115 }
116 }
117 }
118
119 ofile << std::setprecision(6);
120
121 ofile << std::setprecision(16); // for double value with 8 bytes
122
123 for(G4int x = 0; x < fNMeshSegments[0]; x++) {
124 for(G4int y = 0; y < fNMeshSegments[1]; y++) {
125 for(G4int z = 0; z < fNMeshSegments[2]; z++){
126
127 G4int idx = GetIndex(x, y, z);
128
129 if (ScoringMeshEdep[idx] != 0){
130 ofile << x << '\t' << y << '\t' << z << '\t' << ScoringMeshEdep[idx] << G4endl;
131 }
132 //Store x,y,z and dose for each voxel in output text file.
133
134 }
135 }
136 }
137
138 // Close the output ASCII file
139 ofile.close();
140
141 //----------------------------------------------------------------------------------------//
142 //-----Read Data.dat File to determine the name and number of slice files to open---------//
143 //----------------------------------------------------------------------------------------//
144 // Using the macro commands from the .in files, the UserScoreWriter identifies which
145 // phantom sex and section has been constructed. We then store the names of the individual z-slices
146 // which have been called upon in the detector construction when creating the phantom.
147
148 G4int NSlices = 0;
149 G4int NXVoxels = 0;
150 G4int NYVoxels = 0;
151 std::ifstream DataFile;
152
153 G4cout << "Phantom Sex: " << fSex << G4endl;
154 G4cout << "Phantom Section: " << fSection << G4endl;
155
156 G4String male = "male";
157 G4String female = "female";
158 //G4String
159
160 //Determine Phantom Sex and Section which was Simulated
161 if (fSex == "male")
162 {
163 if (fSection == "head")
164 {
165 DataFile.open("ICRPdata/MaleHead.dat");
166 G4cout << "Selecting data file ICRPdata/MaleHead.dat..." << G4endl;
167 }
168 if (fSection == "trunk")
169 {
170 DataFile.open("ICRPdata/MaleTrunk.dat");
171 G4cout << "Selecting data file ICRPdata/MaleTrunk.dat..." << G4endl;
172 }
173 if (fSection == "full")
174 {
175 DataFile.open("ICRPdata/MaleData.dat");
176 G4cout << "Selecting data file ICRPdata/MaleData.dat..." << G4endl;
177 }
178 }
179 else if(fSex == "female")
180 {
181 if (fSection == "head")
182 {
183 DataFile.open("ICRPdata/FemaleHead.dat");
184 G4cout << "Selecting data file ICRPdata/FemaleHead.dat..." << G4endl;
185 }
186 if (fSection == "trunk")
187 {
188 DataFile.open("ICRPdata/FemaleTrunk.dat");
189 G4cout << "Selecting data file ICRPdata/FemaleTrunk.dat..." << G4endl;
190 }
191 if (fSection == "full")
192 {
193 DataFile.open("ICRPdata/FemaleData.dat");
194 G4cout << "Selecting data file ICRPdata/FemaleData.dat..." << G4endl;
195 }
196 }
197 else
198 {
199 G4cout << "Phantom Sex or section not correctly specified to ICRP110UserScoreWriter" << G4endl;
200 }
201
202 //Check if file opens
203 if(DataFile.good() != 1 )
204 {
205 G4cout << "Problem Reading Data File" << G4endl;
206 }
207 else
208 {
209 G4cout << "Opening Data.dat File..." << G4endl;
210 }
211
212 DataFile >> NSlices;
213 G4cout << "Number of Phantom Slices Simulated = " << NSlices << G4endl;
214
215 DataFile >> NXVoxels >> NYVoxels;
216 G4cout << "Number of X Voxels per slice = " << NXVoxels << G4endl;
217 G4cout << "Number of Y Voxels per slice = " << NYVoxels << G4endl;
218
219
220 G4int VoxelsPerSlice = 0;
221 VoxelsPerSlice = NXVoxels * NYVoxels;
222
223 //Skip lines 4-60 of Data.dat as they hold no useful information for this code
224 for (G4int i = 0; i < 58; i++){
225 DataFile.ignore(256, '\n');
226 }
227
228 // Read file names to open from Data.dat (i.e. those that were used in simulation)
229
230 std::vector<G4String> SliceName; //char SliceName[NSlices][20]; //Stores name and number of phantom slices used in simulation
231
232 for(G4int i = 0; i < NSlices; i++){
233 SliceName.push_back("empty");
234 }
235
236 for (G4int i = 0; i < NSlices; i++){
237 DataFile >> SliceName[i];
238 }
239
240
241 //--------------------------------------------------------------------//
242 //----------------------Read Phantom Slice Files----------------------//
243 //--------------------------------------------------------------------//
244 // Reads each of the phantom z-slice files identified by the above code
245 // and stores the organIDs within each voxel sequentially.
246 // Later, we will compare the dose in each voxel to the organ ID of each
247 // voxel to calculate total dose in each organ.
248
249 G4int ARRAY_SIZE = VoxelsPerSlice * NSlices;
250
251 std::vector<G4int> OrganIDs;
252
253 std::ifstream PhantomFile;
254
255 for (G4int ii = 0; ii < NSlices; ii++){
256
257 G4String sliceVar = SliceName[ii];
258 G4String slice;
259
260 if (fSex == "male")
261 {
262 slice = "ICRPdata/ICRP110_g4dat/AM/"+sliceVar;
263 }
264 else if(fSex == "female")
265 {
266 slice = "ICRPdata/ICRP110_g4dat/AF/"+sliceVar;
267 }
268
269 PhantomFile.open(slice.c_str());
270
271 //Check if file opens
272 if(PhantomFile.good() != 1 )
273 {
274 G4cout << "Problem Reading Phantom Slice File:" << SliceName[ii] << G4endl;
275 }
276 else
277 {
278 G4cout << "Opening Phantom Slice File: " << SliceName[ii] << G4endl;
279 }
280
281 G4int FNVoxelX;
282 G4int FNVoxelY;
283 G4int FNVoxelZ;
284
285 //Input first 3 numbers of file - They are not organ IDs
286 PhantomFile >> FNVoxelX >> FNVoxelY >> FNVoxelZ;
287
288 G4int aa = 0;
289
290 for (G4int i = 0; i < VoxelsPerSlice; i++)
291 {
292 PhantomFile >> aa;
293 OrganIDs.push_back(aa);
294 }
295
296 PhantomFile.close();
297 }
298
299 //---------------------------------------------------------------//
300 //--------------Read Data from Scoring Mesh Text File------------//
301 //---------------------------------------------------------------//
302 // Opens and reads initial/default scoring mesh output file which
303 // was created at the beginning of the code. We now store the edep in each
304 // voxel to cross reference against the organ ID of each voxel for
305 // calculations of total dose in each organ.
306
307 std::ifstream MeshFile(fileName);
308
309 //Check if file opens
310 if(MeshFile.good() != 1 )
311 {
312 G4cout << "Problem Reading Data File: " << fileName << G4endl;
313 }
314 else {
315 G4cout << "Opening File: " << fileName << G4endl;
316 }
317
318
319 //-----Reads Phantom Mesh Text File and Stores Data in 6 different Vectors-------//
320
321 G4int col = 4;
322 G4int lines = VoxelsPerSlice*NSlices;
323
324 //Ignore first 2 lines of PhantomMesh.txt as they are text headers
325 MeshFile.ignore(256, '\n');
326 MeshFile.ignore(256, '\n');
327
328 std::vector<G4int> X_MeshID; //Stores X-position of all scoring mesh voxels
329 std::vector<G4int> Y_MeshID; //Stores Y-position
330 std::vector<G4int> Z_MeshID; //Stores Z-position
331 std::vector<G4double> Edep; //Stores edep in voxels
332
333 G4int nX = 0; //Number along X of scoring mesh voxel
334 G4int nY = 0; //Number along Y
335 G4int nZ = 0; //Number along Z
336 G4double EDep = 0.0; //edep deposited in individual voxels
337
338 for (G4int i=0; i< lines; i++){
339 for (G4int j=0; j < col;){
340
341 MeshFile >> nX; //Reads number along X of current scoring mesh voxel
342 X_MeshID.push_back(nX); //Stores it sequentially in vector X_MeshID
343 j++;
344
345 MeshFile >> nY; // Reads number along Y
346 Y_MeshID.push_back(nY); // Stores in vector
347 j++;
348
349 MeshFile >> nZ; // Reads number along Z
350 Z_MeshID.push_back(nZ); // Stores in vector
351 j++;
352
353 MeshFile >> EDep; // Reads edep in each voxel
354 Edep.push_back(EDep); // Stores in vector
355 j++;
356
357 }
358 }
359
360 MeshFile.close();
361
362 //--------------------------------------------------------------------//
363 //---------Reads AM/AF_organs.dat file and stores info about----------//
364 //------------------------the phantom organs--------------------------//
365 //--------------------------------------------------------------------//
366
367 std::ifstream PhantomOrganNames;
368
369 if (strcmp(fSex.c_str(), male.c_str()) == 0)
370 {
371 PhantomOrganNames.open ("ICRPdata/ICRP110_g4dat/P110_data_V1.2/AM/AM_organs.dat");
372 //Check if file opens
373 if(PhantomOrganNames.good() != 1 )
374 {
375 G4cout << "Problem reading AM_organs.dat" << G4endl;
376 }
377 else
378 {
379 G4cout << "Reading AM_organs.dat" << G4endl;
380 }
381 }
382 else if(strcmp(fSex.c_str(), female.c_str()) == 0)
383 {
384 PhantomOrganNames.open ("ICRPdata/ICRP110_g4dat/P110_data_V1.2/AF/AF_organs.dat");
385 //Check if file opens
386 if(PhantomOrganNames.good() != 1 )
387 {
388 G4cout << "Problem reading AF_organs.dat" << G4endl;
389 }
390 else
391 {
392 G4cout << "Reading AF_organs.dat" << G4endl;
393 }
394 }
395
396
397 PhantomOrganNames.ignore(256, '\n');
398 PhantomOrganNames.ignore(256, '\n');
399 PhantomOrganNames.ignore(256, '\n');
400 PhantomOrganNames.ignore(256, '\n');
401
402 G4String str;
403 std::vector<G4String> OrganNames;
404
405 OrganNames.push_back("0 Air"); // Register air surrounding phantom
406 //Fill with organ IDs of the phantom from ICRP data files (located in /ICRPdata/ICRP110_g4dat/P110_data_V1.2/)
407 while(getline(PhantomOrganNames, str))
408 {
409 OrganNames.push_back(str);
410 }
411 OrganNames.push_back("141 Phantom Top/Bottom Skin Layer"); //Registers top and bottom slices of phantom made entirely of
412 // skin. The skin in these layers has organ ID 141 to differentiate it from other skin, and is given its
413 // own organ ID so that the user can choose whether to include it or not.
414
415 //-------------------------------------------------------------------//
416 //---------Reads OrganMasses.dat file and stores info about----------//
417 //--------------------the phantom organ massess----------------------//
418 //-------------------------------------------------------------------//
419
420 G4int NOrganIDs = OrganNames.size();
421
422 std::ifstream OrganMasses;
423
424 OrganMasses.open ("ICRPdata/OrganMasses.dat");
425 //Check if file opens
426 if(OrganMasses.good() != 1 )
427 {
428 G4cout << "Problem reading OrganMasses.dat" << G4endl;
429 }
430 else
431 {
432 G4cout << "Reading OrganMasses.dat" << G4endl;
433 }
434
435
436 OrganMasses.ignore(256, '\n'); //Igonore first line as it is a header
437
438 std::vector<G4int> iteratorID;
439 std::vector<G4double> MaleOrganMasses;
440 std::vector<G4double> FemaleOrganMasses;
441
442 G4int itID = 0;
443 G4double massM = 0.0;
444 G4double massF = 0.0;
445
446 for (G4int i = 0; i < NOrganIDs; i++)
447 {
448 OrganMasses >> itID;
449 iteratorID.push_back(itID);
450
451 OrganMasses >> massM;
452 MaleOrganMasses.push_back(massM);
453
454 OrganMasses >> massF;
455 FemaleOrganMasses.push_back(massF);
456 }
457
458 //---------------------------------------------------------------------------//
459 //----------------------Writes Outputs of code to File-----------------------//
460 //-------------------------------OrganDeps.out------------------------------//
461 //---------------------------------------------------------------------------//
462 // As the final step, we compare the edep in each voxel with the voxels organID
463 // and sum the edep in voxels with identical organIDs to obtain total edep in
464 // each organ. We then divide total edep in each organ by their respective organ
465 // mass to give total dose received in each organ (in Gy).
466 // All this information is then output to the file "ICRP110.out".
467
468 std::ofstream OutputFile2;
469
470 G4int VoxelNumber = 0;
471 G4int OrganIndex = 0;
472 G4cout << "NOrganIDs: " << NOrganIDs << G4endl;
473 std::vector <G4double> OrganDep;
474 std::vector <G4double> OrganDose;
475
476 G4double a = 0.0;
477 G4double b = 0.0;
478
479 for (G4int i = 0; i < NOrganIDs; i++)
480 {
481 OrganDep.push_back(a);
482 OrganDose.push_back(b);
483 }
484
485
486 for (G4int i = 0; i < ARRAY_SIZE; i++){
487 VoxelNumber = X_MeshID[i] + NXVoxels * Y_MeshID[i] + VoxelsPerSlice * Z_MeshID[i];
488 OrganIndex = OrganIDs[VoxelNumber];
489 OrganDep[OrganIndex] += Edep[i];
490 }
491
492 //Calculate dose in each organ by dividing edep in each organ by organ masses (in kg)
493 if (strcmp(fSex.c_str(), male.c_str()) == 0)
494 {
495 for (G4int i = 0; i < NOrganIDs; i++)
496 {
497 OrganDose[i] = (MaleOrganMasses[i] == 0 ) ? 0 : OrganDep[i]/(MaleOrganMasses[i] * 1e-3);
498 }
499 }
500 else if(strcmp(fSex.c_str(), female.c_str()) == 0)
501 {
502 for (G4int i = 0; i < NOrganIDs; i++)
503 {
504 OrganDose[i] = (FemaleOrganMasses[i] == 0 ) ? 0 : OrganDep[i]/(FemaleOrganMasses[i] * 1e-3);
505 }
506 }
507
508 OutputFile2.open ("ICRP110.out");
509
510 //Check if file opens
511 if(OutputFile2.good() != 1 )
512 {
513 G4cout << "Problem writing output to ICRP110.out" << G4endl;
514 }
515 else {
516 G4cout << "Writing output to ICRP110.out" << G4endl;
517 }
518
519
520 G4double TotalDep = 0.0;
521 G4double TotalDose = 0.0;
522
523 OutputFile2 << G4endl;
524 OutputFile2 << '\t' << "-------------------------------- " << G4endl;
525 OutputFile2 << '\t' << "OrganID" << '\t' << "Edep (J)" << '\t' << "Dose (Gy)" << G4endl;
526 OutputFile2 << '\t' << "-------------------------------- " << G4endl;
527
528
529 for (G4int i = 1; i < NOrganIDs; i++)
530 {
531 if (OrganDep[i] != 0)
532 {
533 if (i != 140) //Skip dose deposited in air inside body
534 {
535 OutputFile2 << '\t' << i << " |" << '\t' << '\t' << OrganDep[i] << '\t' << OrganDose[i] << G4endl;
536 }
537 }
538 }
539
540 OutputFile2 << "----------------------------------------------------------------------------" << G4endl;
541 OutputFile2 << "-------------------------------ORGAN INFO-----------------------------------" << G4endl;
542 OutputFile2 << "-----------------(of organs where edep/dose was recorded)-------------------" << G4endl;
543 OutputFile2 << "----------------------------------------------------------------------------" << G4endl;
544 OutputFile2 << "ID" << '\t' << '\t' << "Organ Name" << '\t' << " " << '\t' << " " << '\t' << " " << '\t' << "Material ID" << '\t' << '\t' << "Density (g/cm^3) " << G4endl;
545
546 for(G4int i = 1; i < NOrganIDs; i++)
547 {
548 if (OrganDep[i] != 0)
549 {
550 if (i != 140) //Skip dose deposited in air inside body
551 {
552 OutputFile2 << OrganNames[i] << G4endl;
553 }
554 }
555 }
556
557 //Sum total dose over all organs
558 for (G4int i = 1; i < NOrganIDs; i++)
559 {
560 if (i != 140) //Skip dose deposited in air inside body
561 {
562 TotalDep += OrganDep[i];
563 TotalDose += OrganDose[i];
564 }
565 }
566
567 OutputFile2 << G4endl;
568 OutputFile2 << "Total Edep over all organs = " << TotalDep << " J" << G4endl;
569 OutputFile2 << "Total dose absorbed over all organs = " << TotalDose << " Gy" << G4endl;
570
571
572 OutputFile2 << G4endl;
573 OutputFile2 << "----------------------------------------------------------------------------" << G4endl;
574 OutputFile2 << "----------------ORGAN ENERGY DEPOSITIONS AND ABSORBED DOSE------------------" << G4endl;
575 OutputFile2 << "-----------(for all organs [includes air - OrganIDs = 0, 140])--------------" << G4endl;
576 OutputFile2 << "--------------([and top/bottom skin layer - OrganIDs = 141])----------------" << G4endl;
577 OutputFile2 << "----------------------------------------------------------------------------" << G4endl;
578 OutputFile2 << "OrganID" << '\t' << "Edep (J) " << '\t' << "Dose (Gy) " << G4endl;
579 OutputFile2 << "-------------------------------" << G4endl;
580
581 for (G4int i = 0; i < NOrganIDs; i++)
582 {
583 OutputFile2 << i << " | " << '\t' << '\t' << OrganDep[i] << '\t' << OrganDose[i] << G4endl;
584 }
585
586 OutputFile2 << "Total energy depositied over all organs = " << TotalDep << " J" << G4endl;
587 OutputFile2 << "Total absorbed dose over all organs = " << TotalDose << " Gy " << G4endl;
588 G4cout << "Total energy deposited over all Organs within the Phantom is " << TotalDep << " J" << G4endl;
589 G4cout << "Total absorbed dose over all phantom organs is " << TotalDose << " Gy " << G4endl;
590
591 OutputFile2.close();
592 }
593
594
595 // Sets the sex of the phantom as defined through the messenger class
596 void ICRP110UserScoreWriter::SetPhantomSex(G4String newSex)
597 {
598 fSex = newSex;
599
600 if (fSex == "male")
601 {
602 G4cout << ">> Male Phantom identified by UserScoreWriter." << G4endl;
603 }
604 if (fSex == "female")
605 {
606 G4cout << ">> Female Phantom identified by UserScoreWriter." << G4endl;
607 }
608 if ((fSex != "female") && (fSex != "male"))
609 G4cout << fSex << " can not be defined!" << G4endl;
610 }
611
612 // Sets the section of the phantom as defined through the messenger class
613 void ICRP110UserScoreWriter::SetPhantomSection(G4String newSection)
614 {
615 fSection = newSection;
616
617 if (fSection == "head")
618 {
619 G4cout << ">> Partial Head Phantom identified by UserScoreWriter." << G4endl;
620 }
621 if (fSection == "trunk")
622 {
623 G4cout << ">> Partial Trunk Phantom identified by UserScoreWriter." << G4endl;
624 }
625 if (fSection == "full")
626 {
627 G4cout << ">> Custom/Full Phantom identified by UserScoreWriter." << G4endl;
628 }
629 if ((fSection != "head") && (fSection != "trunk") && (fSection != "full"))
630 G4cout << fSection << " can not be defined!" << G4endl;
631 }
632