| 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 // MONTE CARLO SIMULATION OF REALISTIC GEOMETRY FROM MICROSCOPES IMAGES
28 //
29 // Authors and contributors:
30 // P. Barberet, S. Incerti, N. H. Tran, L. Morelli
31 //
32 // University of Bordeaux, CNRS, LP2i, UMR5797, Gradignan, France
33 //
34 // If you use this code, please cite the following publication:
35 // P. Barberet et al.,
36 // "Monte-Carlo dosimetry on a realistic cell monolayer
37 // geometry exposed to alpha particles."
38 // Ph. Barberet et al 2012 Phys. Med. Biol. 57 2189
39 // doi: 110.1088/0031-9155/57/8/2189
40 // --------------------------------------------------------------------------------
41
42 #include "CellParameterisation.hh"
43
44 #include "G4Material.hh"
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 CellParameterisation *CellParameterisation::gInstance = nullptr;
49
50 CellParameterisation::CellParameterisation
51 (G4String fileName,
52 G4Material *RedMat, G4Material *GreenMat, G4Material *BlueMat,
53 G4double shiftX, G4double shiftY, G4double shiftZ
54 )
55 :fRedMaterial(RedMat), fGreenMaterial(GreenMat), fBlueMaterial(BlueMat),
56 fShiftX(shiftX), fShiftY(shiftY), fShiftZ(shiftZ)
57 {
58 Initialize(fileName);
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63 void CellParameterisation::Initialize(const G4String &fileName)
64 {
65 G4int ncols, l, mat;
66 G4int pixelX, pixelY, pixelZ;
67 G4double x, y, z, den1, den2, den3;
68
69 ncols = 0;
70 l = 0;
71
72 // Read phantom
73
74 FILE *fMap;
75 fMap = fopen(fileName, "r");
76
77 fRedMass = 0;
78 fGreenMass = 0;
79 fBlueMass = 0;
80
81 ncols = fscanf(fMap, "%d %d %d %d", &fPhantomTotalPixels, &fRedTotalPixels, &fGreenTotalPixels,
82 &fBlueTotalPixels);
83 ncols = fscanf(fMap, "%lf %lf %lf %s", &fSizeRealX, &fSizeRealY, &fSizeRealZ, &fRealUnit);
84 ncols = fscanf(fMap, "%lf %lf %lf %s", &fDimCellBoxX, &fDimCellBoxY, &fDimCellBoxZ, &fRealUnit);
85
86 fMapCell = new G4ThreeVector[fPhantomTotalPixels]; //geant4 coordinates space
87 fMapCellPxl = new G4ThreeVector[fPhantomTotalPixels]; //voxel space
88 fMapCellOriginal = new G4ThreeVector[fPhantomTotalPixels]; //original coordinates space
89 fMaterial = new G4int[fPhantomTotalPixels];
90
91 fDimCellBoxX = fDimCellBoxX * um;
92 fDimCellBoxY = fDimCellBoxY * um;
93 fDimCellBoxZ = fDimCellBoxZ * um;
94
95 den1 = fRedMaterial->GetDensity();
96 den2 = fGreenMaterial->GetDensity();
97 den3 = fBlueMaterial->GetDensity();
98
99 fOffsetX = -fSizeRealX / 2 *um;
100 fOffsetY = -fSizeRealY / 2 *um;
101 fOffsetZ = -fSizeRealZ / 2 *um;
102
103 G4cout << G4endl;
104 G4cout << " #########################################################################" << G4endl;
105 G4cout << " Phantom placement and density " << G4endl;
106 G4cout << " #########################################################################" << G4endl;
107 G4cout << G4endl;
108 G4cout << " ==========> Phantom origin - X (um) = " << (fOffsetX + fShiftX)/um << G4endl;
109 G4cout << " ==========> Phantom origin - Y (um) = " << (fOffsetY + fShiftY)/um << G4endl;
110 G4cout << " ==========> Phantom origin - Z (um) = " << (fOffsetZ + fShiftZ)/um << G4endl;
111 G4cout << G4endl;
112 G4cout << " ==========> Red density (g/cm3) = " << den1/(g/cm3) << G4endl;
113 G4cout << " ==========> Green density (g/cm3) = " << den2/(g/cm3) << G4endl;
114 G4cout << " ==========> Blue density (g/cm3) = " << den3/(g/cm3) << G4endl;
115 G4cout << G4endl;
116 G4cout << " #########################################################################" << G4endl;
117 G4cout << G4endl;
118
119 while (1)
120 {
121 ncols = fscanf(fMap, "%lf %lf %lf %d", &x, &y, &z, &mat);
122 if (ncols < 0) break;
123
124 G4ThreeVector v( x*um + fOffsetX + fShiftX, // phantom shift
125 -(y*um + fOffsetY + fShiftY),
126 z*um + fOffsetZ + fShiftZ );
127
128 // Pixel coordinates
129 pixelX = (x*um)/fDimCellBoxX;
130 pixelY = (y*um)/fDimCellBoxY;
131 pixelZ = (z*um)/fDimCellBoxZ;
132
133 G4ThreeVector w(pixelX, pixelY, pixelZ);
134
135 G4ThreeVector v_original(x*um, y*um, z*um);
136
137 fMapCell[l] = v;
138 fMapCellPxl[l] = w;
139 fMapCellOriginal[l] = v_original;
140
141 fMaterial[l] = mat;
142
143 if (mat == 1){
144 fRedMass += den1 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
145 }
146 else if (mat == 2){
147 fGreenMass += den2 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
148 }
149 else if (mat == 3){
150 fBlueMass += den3 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
151 }
152 l++;
153 }
154
155 fclose(fMap);
156
157 fRedAttributes = new G4VisAttributes;
158 fRedAttributes->SetColour(G4Colour(1, 0, 0));
159 fRedAttributes->SetForceSolid(false);
160
161 fGreenAttributes = new G4VisAttributes;
162 fGreenAttributes->SetColour(G4Colour(0, 1, 0));
163 fGreenAttributes->SetForceSolid(false);
164
165 fBlueAttributes = new G4VisAttributes;
166 fBlueAttributes->SetColour(G4Colour(0, 0, 1));
167 fBlueAttributes->SetForceSolid(false);
168
169 gInstance = this;
170 }
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172
173 CellParameterisation::~CellParameterisation()
174 {
175 delete[] fMapCell;
176 delete[] fMapCellPxl;
177 delete[] fMaterial;
178 }
179
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181
182 void CellParameterisation::ComputeTransformation
183 (const G4int copyNo, G4VPhysicalVolume *physVol) const
184 {
185 if(fMapCell == nullptr)
186 {
187 G4ExceptionDescription ex;
188 ex<< "fMapCell == nullptr ";
189 G4Exception("CellParameterisation::ComputeTransformation",
190 "CellParameterisation001",
191 FatalException,
192 ex);
193 }
194 else
195 {
196 G4ThreeVector
197 origin(fMapCell[copyNo].x(), fMapCell[copyNo].y(), fMapCell[copyNo].z());
198
199 physVol->SetTranslation(origin);
200 }
201 }
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203
204 G4Material *
205 CellParameterisation::ComputeMaterial(const G4int copyNo,
206 G4VPhysicalVolume *physVol,
207 const G4VTouchable *)
208 {
209 if (fMaterial[copyNo] == 3) // fMaterial 3 is blue
210 {
211 physVol->SetName("physicalMat3");
212 physVol->GetLogicalVolume()->SetVisAttributes(fBlueAttributes);
213 return fBlueMaterial;
214 }
215 else if (fMaterial[copyNo] == 2) // fMaterial 2 is green
216 {
217 physVol->SetName("physicalMat2");
218 physVol->GetLogicalVolume()->SetVisAttributes(fGreenAttributes);
219 return fGreenMaterial;
220 }
221 else if (fMaterial[copyNo] == 1) // fMaterial 1 is red
222 {
223 physVol->SetName("physicalMat1");
224 physVol->GetLogicalVolume()->SetVisAttributes(fRedAttributes);
225 return fRedMaterial;
226 }
227
228 return physVol->GetLogicalVolume()->GetMaterial();
229 }
230