| 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 /// \file SAXSDetectorConstruction.cc
27 /// \brief Implementation of the SAXSDetectorConstruction class
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30
31 #include "SAXSDetectorConstruction.hh"
32
33 #include "SAXSSensitiveDetector.hh"
34
35 #include "G4AssemblyVolume.hh"
36 #include "G4Box.hh"
37 #include "G4Colour.hh"
38 #include "G4Cons.hh"
39 #include "G4Element.hh"
40 #include "G4ElementTable.hh"
41 #include "G4ExtendedMaterial.hh"
42 #include "G4GeometryManager.hh"
43 #include "G4IntersectionSolid.hh"
44 #include "G4LogicalVolume.hh"
45 #include "G4LogicalVolumeStore.hh"
46 #include "G4MIData.hh"
47 #include "G4Material.hh"
48 #include "G4MaterialTable.hh"
49 #include "G4MultiFunctionalDetector.hh"
50 #include "G4NistManager.hh"
51 #include "G4PSDoseDeposit.hh"
52 #include "G4PSEnergyDeposit.hh"
53 #include "G4PVPlacement.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4PhysicalVolumeStore.hh"
56 #include "G4Polyhedra.hh"
57 #include "G4RotationMatrix.hh"
58 #include "G4RunManager.hh"
59 #include "G4SDManager.hh"
60 #include "G4SDParticleFilter.hh"
61 #include "G4SolidStore.hh"
62 #include "G4Sphere.hh"
63 #include "G4SubtractionSolid.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4ThreeVector.hh"
66 #include "G4Trd.hh"
67 #include "G4Tubs.hh"
68 #include "G4UniformMagField.hh"
69 #include "G4UnionSolid.hh"
70 #include "G4VSensitiveDetector.hh"
71 #include "G4VisAttributes.hh"
72 #include "globals.hh"
73
74 #include <cmath>
75 #include <sstream>
76 #include <vector>
77
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79
80 SAXSDetectorConstruction::SAXSDetectorConstruction() : G4VUserDetectorConstruction(), fWorldLogic(0)
81 {
82 G4cout << "### DetectorConstruction Instantiated ###" << G4endl;
83
84 // instantiate the messenger (set methods will be called after construct)
85 fMessenger = new SAXSDetectorConstructionMessenger(this);
86
87 // set geometrical variables
88 SetGeometricalVariables();
89
90 // Initialization
91 fPhantomMaterialIndex = 1;
92
93 fComp0 = 0.0; // components of the "Medical Material (MedMat)"
94 fComp1 = 1.0; // Through macro I can set one Medical Material only
95 fComp2 = 0.0;
96 fComp3 = 0.0;
97
98 fCustomMatDensity = 1.00; // g/mol
99 fCustomMatHmassfract = 0.1119;
100 fCustomMatCmassfract = 0.;
101 fCustomMatNmassfract = 0.;
102 fCustomMatOmassfract = 0.8881;
103 fCustomMatNamassfract = 0.;
104 fCustomMatPmassfract = 0.;
105 fCustomMatSmassfract = 0.;
106 fCustomMatClmassfract = 0.;
107 fCustomMatKmassfract = 0.;
108 fCustomMatCamassfract = 0.;
109
110 fCustomMatFF = ""; // MIFF filename for a custom (extended) material
111
112 fSensitiveVolume = 0;
113
114 fIWantSlits = false;
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
118
119 SAXSDetectorConstruction::~SAXSDetectorConstruction() {}
120
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122
123 void SAXSDetectorConstruction::DefineMaterials()
124 {
125 // Define the NIST manager
126 G4NistManager* NistMan = G4NistManager::Instance();
127
128 // Define the required elements for compounds
129 G4Element* elH = NistMan->FindOrBuildElement("H");
130 G4Element* elC = NistMan->FindOrBuildElement("C");
131 G4Element* elN = NistMan->FindOrBuildElement("N");
132 G4Element* elO = NistMan->FindOrBuildElement("O");
133 G4Element* elNa = NistMan->FindOrBuildElement("Na");
134 G4Element* elP = NistMan->FindOrBuildElement("P");
135 G4Element* elS = NistMan->FindOrBuildElement("S");
136 G4Element* elCl = NistMan->FindOrBuildElement("Cl");
137 G4Element* elK = NistMan->FindOrBuildElement("K");
138 G4Element* elCa = NistMan->FindOrBuildElement("Ca");
139
140 // variable definition
141 G4double d; // density
142 G4int nel; // number of elements
143 G4String matname;
144
145 // Air
146 d = 1.29 * mg / cm3;
147 nel = 2;
148 G4double tAir = 293.15 * CLHEP::kelvin; // 20° Celsius
149 G4double pAir = 1. * CLHEP::atmosphere; // 1 atm
150 fAir = new G4Material("Air", d, nel, kStateGas, tAir, pAir);
151 fAir->AddElement(elN, 0.7);
152 fAir->AddElement(elO, 0.3);
153
154 // Fat (Tartari2002) (FF from Tartari2002)
155 G4double d_Fat = 0.923 * g / cm3;
156 nel = 3;
157 matname = "Fat_MI";
158 fFat = new G4Material(matname, d_Fat, nel);
159 fFat->AddElement(elH, 0.119);
160 fFat->AddElement(elC, 0.772);
161 fFat->AddElement(elO, 0.109);
162
163 // Water (FF from Tartari2002)
164 G4double d_Water = 1. * g / cm3;
165 nel = 2;
166 matname = "Water_MI";
167 fWater = new G4Material(matname, d_Water, nel);
168 fWater->AddElement(elH, 2);
169 fWater->AddElement(elO, 1);
170
171 // BoneMatrix (Collagen) (FF from Tartari2002)
172 G4double d_BoneMatrix = 1.263 * g / cm3;
173 nel = 4;
174 matname = "BoneMatrix_MI";
175 fBoneMatrix = new G4Material(matname, d_BoneMatrix, nel);
176 fBoneMatrix->AddElement(elH, 0.0344);
177 fBoneMatrix->AddElement(elC, 0.7140);
178 fBoneMatrix->AddElement(elN, 0.1827);
179 fBoneMatrix->AddElement(elO, 0.0689);
180
181 // Mineral (Hydroxyapatite) (Tartari2002) (FF from Tartari2002)
182 G4double d_Mineral = 2.74 * g / cm3;
183 nel = 4;
184 matname = "Mineral_MI";
185 fMineral = new G4Material(matname, d_Mineral, nel);
186 fMineral->AddElement(elH, 0.002);
187 fMineral->AddElement(elO, 0.414);
188 fMineral->AddElement(elP, 0.185);
189 fMineral->AddElement(elCa, 0.399);
190
191 // Medical Material (compostion of Water, Fat, BoneMatrix and Mineral)
192 G4double comp[] = {fComp0, fComp1, fComp2, fComp3};
193 G4double d_MedMat =
194 1 / (comp[0] / d_Fat + comp[1] / d_Water + comp[2] / d_BoneMatrix + comp[3] / d_Mineral);
195 G4int n_MedMat = 0;
196 for (size_t i = 0; i < 4; i++) {
197 if (comp[i] > 0) n_MedMat++;
198 if (comp[i] < 0 || comp[i] > 1) {
199 G4String excep = "Error in Medical Material composition: comp[i]<0 or comp[i]>1";
200 G4Exception("DetectorConstuction::DefineMaterials()", "dc0001", FatalException, excep);
201 return;
202 }
203 }
204 std::stringstream ss0, ss1, ss2, ss3;
205 ss0 << comp[0];
206 ss1 << comp[1];
207 ss2 << comp[2];
208 ss3 << comp[3];
209 if (comp[0] == 0 || comp[0] == 1) ss0 << ".00";
210 if (comp[1] == 0 || comp[1] == 1) ss1 << ".00";
211 if (comp[2] == 0 || comp[2] == 1) ss2 << ".00";
212 if (comp[3] == 0 || comp[3] == 1) ss3 << ".00";
213 if (ss0.str().size() < 4) ss0 << "0";
214 if (ss1.str().size() < 4) ss1 << "0";
215 if (ss2.str().size() < 4) ss2 << "0";
216 if (ss3.str().size() < 4) ss3 << "0";
217 if (ss0.str().size() != 4 || ss1.str().size() != 4 || ss2.str().size() != 4
218 || ss3.str().size() != 4)
219 {
220 G4String excep = "Error in MedMaterial composition: check the digits of the elements of comp";
221 G4Exception("DetectorConstuction::DefineMaterials()", "dc0002", FatalException, excep);
222 return;
223 }
224 matname = "MedMat_" + ss0.str() + "_" + ss1.str() + "_" + ss2.str() + "_" + ss3.str();
225 fMedMat = new G4Material(matname, d_MedMat, n_MedMat);
226 if (comp[0]) fMedMat->AddMaterial(fFat, comp[0]);
227 if (comp[1]) fMedMat->AddMaterial(fWater, comp[1]);
228 if (comp[2]) fMedMat->AddMaterial(fBoneMatrix, comp[2]);
229 if (comp[3]) fMedMat->AddMaterial(fMineral, comp[3]);
230 if (comp[0] + comp[1] + comp[2] + comp[3] != 1) {
231 G4String excep = "Error in Medical Material composition: sum(comp) != 1";
232 G4Exception("DetectorConstuction::DefineMaterials()", "dc0003", FatalException, excep);
233 return;
234 }
235 // If the user wants to use more than one MedMat, he has to create the mix
236 // and label the material properly, such as "MedMat_0.55_0.25_0.05_0.15".
237 // Such a name enables the automatic form factor calculation.
238
239 // PMMA (FF from Tartari2002)
240 d = 1.18 * g / cm3;
241 nel = 3;
242 matname = "PMMA_MI";
243 fPMMA = new G4Material(matname, d, nel);
244 fPMMA->AddElement(elH, 8);
245 fPMMA->AddElement(elC, 5);
246 fPMMA->AddElement(elO, 2);
247
248 // Adipose (Poletti2002) (FF from Poletti2002)
249 d = 0.92 * g / cm3;
250 nel = 4;
251 matname = "adipose_MI";
252 fAdipose = new G4Material(matname, d, nel);
253 fAdipose->AddElement(elH, 0.124);
254 fAdipose->AddElement(elC, 0.765);
255 fAdipose->AddElement(elN, 0.004);
256 fAdipose->AddElement(elO, 0.107);
257
258 // Glandular (Poletti2002) (FF from Poletti2002)
259 d = 1.04 * g / cm3;
260 nel = 4;
261 matname = "glandular_MI";
262 fGlandular = new G4Material(matname, d, nel);
263 fGlandular->AddElement(elH, 0.093);
264 fGlandular->AddElement(elC, 0.184);
265 fGlandular->AddElement(elN, 0.044);
266 fGlandular->AddElement(elO, 0.679);
267
268 // human breast 50/50 (ICRU44) (FF from Peplow1998)
269 d = 0.96 * g / cm3;
270 nel = 3;
271 matname = "breast5050_MI";
272 fBreast5050 = new G4Material(matname, d, nel);
273 fBreast5050->AddElement(elH, 0.115);
274 fBreast5050->AddElement(elC, 0.387);
275 fBreast5050->AddElement(elO, 0.498);
276
277 // Liver (ICRU46) (FF_pork_liver_Peplow1998)
278 d = 1.06 * g / cm3;
279 nel = 9;
280 matname = "liver_MI";
281 fliver = new G4Material(matname, d, nel);
282 fliver->AddElement(elH, 0.102);
283 fliver->AddElement(elC, 0.139);
284 fliver->AddElement(elN, 0.030);
285 fliver->AddElement(elO, 0.716);
286 fliver->AddElement(elNa, 0.002);
287 fliver->AddElement(elP, 0.003);
288 fliver->AddElement(elS, 0.003);
289 fliver->AddElement(elCl, 0.002);
290 fliver->AddElement(elK, 0.003);
291
292 // Kidney (ICRU46) (FF_pork_kidney_Peplow1998)
293 d = 1.05 * g / cm3;
294 nel = 10;
295 matname = "kidney_MI";
296 fkidney = new G4Material(matname, d, nel);
297 fkidney->AddElement(elH, 0.103);
298 fkidney->AddElement(elC, 0.132);
299 fkidney->AddElement(elN, 0.030);
300 fkidney->AddElement(elO, 0.724);
301 fkidney->AddElement(elNa, 0.002);
302 fkidney->AddElement(elP, 0.002);
303 fkidney->AddElement(elS, 0.002);
304 fkidney->AddElement(elCl, 0.002);
305 fkidney->AddElement(elK, 0.002);
306 fkidney->AddElement(elCa, 0.001);
307
308 // Lexan (Polycarbonate) (FF from Peplow1998)
309 d = 1.221 * g / cm3;
310 nel = 3;
311 matname = "Lexan_MI";
312 fLexan = new G4Material(matname, d, nel);
313 fLexan->AddElement(elH, 14);
314 fLexan->AddElement(elC, 16);
315 fLexan->AddElement(elO, 3);
316
317 // Kapton (FF from Peplow1998)
318 d = 1.42 * g / cm3;
319 nel = 4;
320 matname = "Kapton_MI";
321 fKapton = new G4Material(matname, d, nel);
322 fKapton->AddElement(elH, 28);
323 fKapton->AddElement(elC, 35);
324 fKapton->AddElement(elN, 2);
325 fKapton->AddElement(elO, 7);
326
327 // Carcinoma (muscle ICRU44) (FF from Kidane1999)
328 d = 1.05 * g / cm3; // check the density
329 nel = 9;
330 matname = "carcinoma_MI";
331 fcarcinoma = new G4Material(matname, d, nel);
332 fcarcinoma->AddElement(elH, 0.102);
333 fcarcinoma->AddElement(elC, 0.143);
334 fcarcinoma->AddElement(elN, 0.034);
335 fcarcinoma->AddElement(elO, 0.710);
336 fcarcinoma->AddElement(elNa, 0.001);
337 fcarcinoma->AddElement(elP, 0.002);
338 fcarcinoma->AddElement(elS, 0.003);
339 fcarcinoma->AddElement(elCl, 0.001);
340 fcarcinoma->AddElement(elK, 0.004);
341
342 // Nylon (FF from Kosanetzky1987)
343 d = 1.15 * g / cm3;
344 nel = 4;
345 matname = "Nylon_MI";
346 fNylon = new G4Material(matname, d, nel);
347 fNylon->AddElement(elH, 11);
348 fNylon->AddElement(elC, 6);
349 fNylon->AddElement(elN, 1);
350 fNylon->AddElement(elO, 1);
351
352 // Polyethylene (FF from Kosanetzky1987)
353 d = 0.94 * g / cm3; // MDPE => 0.92*g/cm3, HDPE => 0.94*g/cm3
354 nel = 2;
355 matname = "Polyethylene_MI";
356 fPolyethylene = new G4Material(matname, d, nel);
357 fPolyethylene->AddElement(elH, 4);
358 fPolyethylene->AddElement(elC, 2);
359
360 // Polystyrene (FF from Kosanetzky1987)
361 d = 1.05 * g / cm3;
362 nel = 2;
363 matname = "Polystyrene_MI";
364 fPolystyrene = new G4Material(matname, d, nel);
365 fPolystyrene->AddElement(elH, 8);
366 fPolystyrene->AddElement(elC, 8);
367
368 // GrayMatter (DeFelici2008) (FF from DeFelici2008)
369 d = 0.991 * g / cm3;
370 nel = 3;
371 matname = "grayMatter_MI";
372 fGrayMatter = new G4Material(matname, d, nel);
373 fGrayMatter->AddElement(elH, 0.1127);
374 fGrayMatter->AddElement(elC, 0.0849);
375 fGrayMatter->AddElement(elO, 0.8024);
376
377 // WhiteMatter (DeFelici2008) (FF from DeFelici2008)
378 d = 0.983 * g / cm3;
379 nel = 3;
380 matname = "whiteMatter_MI";
381 fWhiteMatter = new G4Material(matname, d, nel);
382 fWhiteMatter->AddElement(elH, 0.1134);
383 fWhiteMatter->AddElement(elC, 0.1621);
384 fWhiteMatter->AddElement(elO, 0.7245);
385
386 // Blood (beef) (ICRU46) (FF from Peplow1998)
387 d = 1.06 * g / cm3;
388 nel = 9;
389 matname = "blood_MI";
390 fbeefBlood = new G4Material(matname, d, nel);
391 fbeefBlood->AddElement(elH, 0.102);
392 fbeefBlood->AddElement(elC, 0.11);
393 fbeefBlood->AddElement(elN, 0.033);
394 fbeefBlood->AddElement(elO, 0.746);
395 fbeefBlood->AddElement(elNa, 0.001);
396 fbeefBlood->AddElement(elP, 0.001);
397 fbeefBlood->AddElement(elS, 0.002);
398 fbeefBlood->AddElement(elCl, 0.003);
399 fbeefBlood->AddElement(elK, 0.002);
400
401 // Formaline (FF from Peplow1998)
402 d = 1.083 * g / cm3;
403 nel = 3;
404 matname = "Formaline_MI";
405 fFormaline = new G4Material(matname, d, nel);
406 fFormaline->AddElement(elH, 2);
407 fFormaline->AddElement(elC, 1);
408 fFormaline->AddElement(elO, 1);
409
410 // Acetone (FF from Cozzini2010)
411 d = 0.7845 * g / cm3;
412 nel = 3;
413 matname = "Acetone_MI";
414 fAcetone = new G4Material(matname, d, nel);
415 fAcetone->AddElement(elH, 6);
416 fAcetone->AddElement(elC, 3);
417 fAcetone->AddElement(elO, 1);
418
419 // Hperoxide (FF from Cozzini2010)
420 d = 1.11 * g / cm3;
421 nel = 2;
422 matname = "Hperoxide_MI";
423 fHperoxide = new G4Material(matname, d, nel);
424 fHperoxide->AddElement(elH, 2);
425 fHperoxide->AddElement(elO, 2);
426
427 // CIRS30-70 (Poletti2002) (FF from Poletti2002)
428 d = 0.97 * g / cm3;
429 nel = 5;
430 matname = "CIRS30-70_MI";
431 fCIRS3070 = new G4Material(matname, d, nel);
432 fCIRS3070->AddElement(elH, 0.1178);
433 fCIRS3070->AddElement(elC, 0.7512);
434 fCIRS3070->AddElement(elN, 0.0066);
435 fCIRS3070->AddElement(elO, 0.1214);
436 fCIRS3070->AddElement(elCa, 0.0030);
437
438 // CIRS50-50 (Poletti2002) (FF from Poletti2002)
439 d = 0.98 * g / cm3;
440 nel = 5;
441 matname = "CIRS50-50_MI";
442 fCIRS5050 = new G4Material(matname, d, nel);
443 fCIRS5050->AddElement(elH, 0.1110);
444 fCIRS5050->AddElement(elC, 0.7274);
445 fCIRS5050->AddElement(elN, 0.0104);
446 fCIRS5050->AddElement(elO, 0.1482);
447 fCIRS5050->AddElement(elCa, 0.0030);
448
449 // CIRS70-30 (Poletti2002) (FF from Poletti2002)
450 d = 1.01 * g / cm3;
451 nel = 5;
452 matname = "CIRS70-30_MI";
453 fCIRS7030 = new G4Material(matname, d, nel);
454 fCIRS7030->AddElement(elH, 0.1172);
455 fCIRS7030->AddElement(elC, 0.7378);
456 fCIRS7030->AddElement(elN, 0.0130);
457 fCIRS7030->AddElement(elO, 0.1244);
458 fCIRS7030->AddElement(elCa, 0.0076);
459
460 // RMI454 (Poletti2002) (FF from Poletti2002)
461 d = 0.98 * g / cm3;
462 nel = 4;
463 matname = "RMI454_MI";
464 fRMI454 = new G4Material(matname, d, nel);
465 fRMI454->AddElement(elH, 0.0924);
466 fRMI454->AddElement(elC, 0.6935);
467 fRMI454->AddElement(elN, 0.0198);
468 fRMI454->AddElement(elO, 0.1943);
469
470 // Bone (King2011 decomposition) (FF from King2011)
471 d = 1.344 * g / cm3;
472 nel = 6;
473 matname = "bone_MI";
474 fBone = new G4Material(matname, d, nel);
475 fBone->AddElement(elH, 0.0582);
476 fBone->AddElement(elC, 0.3055);
477 fBone->AddElement(elN, 0.0347);
478 fBone->AddElement(elO, 0.3856);
479 fBone->AddElement(elP, 0.0684);
480 fBone->AddElement(elCa, 0.1476);
481
482 // FatLowX (Tartari2002) (FF_fat_Tartari2002_joint_lowXdata_ESRF2003)
483 nel = 3;
484 matname = "FatLowX_MI";
485 ffatLowX = new G4Material(matname, d_Fat, nel);
486 ffatLowX->AddElement(elH, 0.119);
487 ffatLowX->AddElement(elC, 0.772);
488 ffatLowX->AddElement(elO, 0.109);
489
490 // BonematrixLowX (Collagen)
491 //(Tartari2002) (FF_bonematrix_Tartari2002_joint_lowXdata)
492 nel = 4;
493 matname = "BoneMatrixLowX_MI";
494 fbonematrixLowX = new G4Material(matname, d_BoneMatrix, nel);
495 fbonematrixLowX->AddElement(elH, 0.0344);
496 fbonematrixLowX->AddElement(elC, 0.7140);
497 fbonematrixLowX->AddElement(elN, 0.1827);
498 fbonematrixLowX->AddElement(elO, 0.0689);
499
500 // dryBoneLowX
501 //(Tartari2002) (FF_dryBone_Tartari2002_joint_lowXdata_ESRF2003)
502 d = 2.06 * g / cm3;
503 nel = 6;
504 matname = "dryBoneLowX_MI";
505 fdryBoneLowX = new G4Material(matname, d, nel);
506 fdryBoneLowX->AddElement(elH, 0.0112);
507 fdryBoneLowX->AddElement(elC, 0.2013);
508 fdryBoneLowX->AddElement(elN, 0.0515);
509 fdryBoneLowX->AddElement(elO, 0.3148);
510 fdryBoneLowX->AddElement(elP, 0.1327);
511 fdryBoneLowX->AddElement(elCa, 0.2885);
512
513 // CustomMat (FF read from file)
514 nel = 0;
515 G4double sumMF = 0.;
516 G4cout << "CustomMat composition: " << G4endl;
517 if (fCustomMatHmassfract) {
518 G4cout << "CustomMatHmassfract: " << fCustomMatHmassfract << G4endl;
519 nel++;
520 sumMF += fCustomMatHmassfract;
521 }
522 if (fCustomMatCmassfract) {
523 G4cout << "CustomMatCmassfract: " << fCustomMatCmassfract << G4endl;
524 nel++;
525 sumMF += fCustomMatCmassfract;
526 }
527 if (fCustomMatNmassfract) {
528 G4cout << "CustomMatNmassfract: " << fCustomMatNmassfract << G4endl;
529 nel++;
530 sumMF += fCustomMatNmassfract;
531 }
532 if (fCustomMatOmassfract) {
533 G4cout << "CustomMatOmassfract: " << fCustomMatOmassfract << G4endl;
534 nel++;
535 sumMF += fCustomMatOmassfract;
536 }
537 if (fCustomMatNamassfract) {
538 G4cout << "CustomMatNamassfract: " << fCustomMatNamassfract << G4endl;
539 nel++;
540 sumMF += fCustomMatNamassfract;
541 }
542 if (fCustomMatPmassfract) {
543 G4cout << "CustomMatPmassfract: " << fCustomMatPmassfract << G4endl;
544 nel++;
545 sumMF += fCustomMatPmassfract;
546 }
547 if (fCustomMatSmassfract) {
548 G4cout << "CustomMatSmassfract: " << fCustomMatSmassfract << G4endl;
549 nel++;
550 sumMF += fCustomMatSmassfract;
551 }
552 if (fCustomMatClmassfract) {
553 G4cout << "CustomMatClmassfract: " << fCustomMatClmassfract << G4endl;
554 nel++;
555 sumMF += fCustomMatClmassfract;
556 }
557 if (fCustomMatKmassfract) {
558 G4cout << "CustomMatKmassfract: " << fCustomMatKmassfract << G4endl;
559 nel++;
560 sumMF += fCustomMatKmassfract;
561 }
562 if (fCustomMatCamassfract) {
563 G4cout << "CustomMatCamassfract: " << fCustomMatCamassfract << G4endl;
564 nel++;
565 sumMF += fCustomMatCamassfract;
566 }
567 if (sumMF == 0.) {
568 // set a default material (water),
569 // otherwiswe an error appears in the interactive mode
570 fCustomMatDensity = 1.00; // g/cm3
571 fCustomMatHmassfract = 0.1119;
572 G4cout << "CustomMat set, but not used!" << G4endl;
573 G4cout << "CustomMatHmassfract: " << fCustomMatHmassfract << G4endl;
574 fCustomMatOmassfract = 0.8881;
575 G4cout << "CustomMatOmassfract: " << fCustomMatOmassfract << G4endl;
576 nel = 2;
577 sumMF = 1;
578 }
579 if (sumMF != 1.) {
580 G4String excep = "Error in Custom Material composition: check elemental mass fractions";
581 G4Exception("DetectorConstuction::DefineMaterials()", "dc0004", FatalException, excep);
582 return;
583 }
584
585 d = fCustomMatDensity * g / cm3;
586 matname = "CustomMat";
587 // Notice: this is an extended material
588 fCustomMat = new G4ExtendedMaterial(matname, d, nel);
589 if (fCustomMatHmassfract) fCustomMat->AddElement(elH, fCustomMatHmassfract);
590 if (fCustomMatCmassfract) fCustomMat->AddElement(elC, fCustomMatCmassfract);
591 if (fCustomMatNmassfract) fCustomMat->AddElement(elN, fCustomMatNmassfract);
592 if (fCustomMatOmassfract) fCustomMat->AddElement(elO, fCustomMatOmassfract);
593 if (fCustomMatNamassfract) fCustomMat->AddElement(elNa, fCustomMatNamassfract);
594 if (fCustomMatPmassfract) fCustomMat->AddElement(elP, fCustomMatPmassfract);
595 if (fCustomMatSmassfract) fCustomMat->AddElement(elS, fCustomMatSmassfract);
596 if (fCustomMatClmassfract) fCustomMat->AddElement(elCl, fCustomMatClmassfract);
597 if (fCustomMatKmassfract) fCustomMat->AddElement(elK, fCustomMatKmassfract);
598 if (fCustomMatCamassfract) fCustomMat->AddElement(elCa, fCustomMatCamassfract);
599 // Register MI extension
600 fCustomMat->RegisterExtension(std::unique_ptr<G4MIData>(new G4MIData("MI")));
601 G4MIData* dataMICustomMat = (G4MIData*)fCustomMat->RetrieveExtension("MI");
602 dataMICustomMat->SetFilenameFF(fCustomMatFF);
603
604 // Nist Materials
605 fLead = NistMan->FindOrBuildMaterial("G4_Pb");
606 fTungsten = NistMan->FindOrBuildMaterial("G4_W");
607 fGe = NistMan->FindOrBuildMaterial("G4_Ge");
608 }
609
610 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
611
612 void SAXSDetectorConstruction::SetGeometricalVariables()
613 {
614 // World
615 fWorldSize = 10000. * mm;
616
617 // Phantom
618 fPhantomDiameter = 50. * mm;
619 fPhantomHeight = 15. * mm;
620 fPhantomZ = 0. * mm;
621
622 // setup angle (rad)
623 fthetaSetup = 0.;
624
625 // Slits
626 fSlitSize = 50. * mm;
627 fSlit1Thickness = 5. * mm;
628 fSlit2Thickness = 5. * mm;
629 fSlit3Thickness = 5. * mm;
630 fSlit4Thickness = 5. * mm;
631 fSlit1SampleDistance = 350. * mm;
632 fSlit2SampleDistance = 50. * mm;
633 fSlit3SampleDistance = 50. * mm;
634 fSlit4SampleDistance = 450. * mm;
635 fSlit1xAperture = 4. * mm;
636 fSlit2xAperture = 4. * mm;
637 fSlit3xAperture = 4. * mm;
638 fSlit4xAperture = 4. * mm;
639 fSlit1yAperture = 4. * mm;
640 fSlit2yAperture = 4. * mm;
641 fSlit3yAperture = 4. * mm;
642 fSlit4yAperture = 4. * mm;
643
644 // Detector
645 fDetectorSize = 20. * mm;
646 fDetectorThickness = 20. * mm;
647 fDetectorSampleDistance = 500. * mm;
648
649 // Shielding
650 fShieldingThickness = 4. * mm;
651 }
652
653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
654
655 G4VPhysicalVolume* SAXSDetectorConstruction::Construct()
656 {
657 // define custom materials
658 DefineMaterials();
659
660 // World
661 fWorldSolid = new G4Box("WorldSolid", fWorldSize * 0.5, fWorldSize * 0.5, fWorldSize * 0.5);
662
663 fWorldLogic = new G4LogicalVolume(fWorldSolid, fAir, "WorldLogic");
664
665 fWorldLogic->SetVisAttributes(G4VisAttributes::GetInvisible());
666
667 fWorldPhysical =
668 new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), fWorldLogic, "WorldPhysical", 0, false, 0);
669
670 // choose the phantom material
671 switch (fPhantomMaterialIndex) {
672 case (1):
673 fPhantomMaterial = fWater;
674 break;
675 case (2):
676 fPhantomMaterial = fMedMat;
677 break;
678 case (3):
679 fPhantomMaterial = fPMMA;
680 break;
681 case (4):
682 fPhantomMaterial = fAdipose;
683 break;
684 case (5):
685 fPhantomMaterial = fGlandular;
686 break;
687 case (6):
688 fPhantomMaterial = fBreast5050;
689 break;
690 case (7):
691 fPhantomMaterial = fcarcinoma;
692 break;
693 case (8):
694 fPhantomMaterial = fkidney;
695 break;
696 case (9):
697 fPhantomMaterial = fliver;
698 break;
699 case (10):
700 fPhantomMaterial = fFat;
701 break;
702 case (11):
703 fPhantomMaterial = fBoneMatrix;
704 break;
705 case (12):
706 fPhantomMaterial = fMineral;
707 break;
708 case (13):
709 fPhantomMaterial = fBone;
710 break;
711 case (14):
712 fPhantomMaterial = ffatLowX;
713 break;
714 case (15):
715 fPhantomMaterial = fbonematrixLowX;
716 break;
717 case (16):
718 fPhantomMaterial = fdryBoneLowX;
719 break;
720 case (17):
721 fPhantomMaterial = fLexan;
722 break;
723 case (18):
724 fPhantomMaterial = fKapton;
725 break;
726 case (19):
727 fPhantomMaterial = fNylon;
728 break;
729 case (20):
730 fPhantomMaterial = fPolyethylene;
731 break;
732 case (21):
733 fPhantomMaterial = fPolystyrene;
734 break;
735 case (22):
736 fPhantomMaterial = fFormaline;
737 break;
738 case (23):
739 fPhantomMaterial = fAcetone;
740 break;
741 case (24):
742 fPhantomMaterial = fHperoxide;
743 break;
744 case (25):
745 fPhantomMaterial = fCIRS3070;
746 break;
747 case (26):
748 fPhantomMaterial = fCIRS5050;
749 break;
750 case (27):
751 fPhantomMaterial = fCIRS7030;
752 break;
753 case (28):
754 fPhantomMaterial = fRMI454;
755 break;
756 case (29):
757 fPhantomMaterial = fAir;
758 break;
759 case (30):
760 fPhantomMaterial = fCustomMat;
761 break;
762 }
763
764 // Phantom (cylinder with axis orthogonal to the X-ray beam axis)
765 G4Tubs* PhantomSolid = new G4Tubs("PhantomSolid", 0., fPhantomDiameter * 0.5,
766 fPhantomHeight * 0.5, 0. * deg, 360. * deg);
767
768 fPhantomLogic = new G4LogicalVolume(PhantomSolid, fPhantomMaterial, "PhantomLogic");
769
770 G4VisAttributes* PhantomVisAttribute = new G4VisAttributes(G4Colour(1., 1., 1.));
771 PhantomVisAttribute->SetForceSolid(true);
772 fPhantomLogic->SetVisAttributes(PhantomVisAttribute);
773
774 G4cout << "Phantom material: " << fPhantomMaterial->GetName() << G4endl;
775 G4cout << "Phantom density: " << fPhantomMaterial->GetDensity() / (g / cm3) << " g/cm3" << G4endl;
776 G4cout << "Phantom mass: " << fPhantomLogic->GetMass() / g << " g" << G4endl;
777
778 G4double rotAngle = 90. * CLHEP::deg;
779 G4RotationMatrix* PhantomRotationMatrix = new G4RotationMatrix(0., 0., 0.);
780 PhantomRotationMatrix->rotateX(rotAngle);
781
782 fPhantomPhysical = new G4PVPlacement(PhantomRotationMatrix, G4ThreeVector(0., 0., fPhantomZ),
783 fPhantomLogic, "PhantomPhysical", fWorldLogic, false, 0);
784
785 // setup rotation matrix (downstream of the phantom/sample)
786 G4RotationMatrix* SetupRotationMatrix = new G4RotationMatrix();
787 SetupRotationMatrix->rotateY(-fthetaSetup);
788
789 // Slits
790 G4Box* Slit1OutSolid =
791 new G4Box("Slit1OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit1Thickness * 0.5);
792
793 G4Box* Slit2OutSolid =
794 new G4Box("Slit2OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit2Thickness * 0.5);
795
796 G4Box* Slit3OutSolid =
797 new G4Box("Slit3OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit3Thickness * 0.5);
798
799 G4Box* Slit4OutSolid =
800 new G4Box("Slit4OutSolid", fSlitSize * 0.5, fSlitSize * 0.5, fSlit4Thickness * 0.5);
801
802 G4Box* Hole1Solid =
803 new G4Box("Hole1Solid", fSlit1xAperture * 0.5, fSlit1yAperture * 0.5, fSlit1Thickness * 0.51);
804
805 G4Box* Hole2Solid =
806 new G4Box("Hole2Solid", fSlit2xAperture * 0.5, fSlit2yAperture * 0.5, fSlit2Thickness * 0.51);
807
808 G4Box* Hole3Solid =
809 new G4Box("Hole3Solid", fSlit3xAperture * 0.5, fSlit3yAperture * 0.5, fSlit3Thickness * 0.51);
810
811 G4Box* Hole4Solid =
812 new G4Box("Hole4Solid", fSlit4xAperture * 0.5, fSlit4yAperture * 0.5, fSlit4Thickness * 0.51);
813
814 G4SubtractionSolid* Slit1Solid = new G4SubtractionSolid("Slit1Solid", Slit1OutSolid, Hole1Solid);
815
816 G4SubtractionSolid* Slit2Solid = new G4SubtractionSolid("Slit1Solid", Slit2OutSolid, Hole2Solid);
817
818 G4SubtractionSolid* Slit3Solid = new G4SubtractionSolid("Slit3Solid", Slit3OutSolid, Hole3Solid);
819
820 G4SubtractionSolid* Slit4Solid = new G4SubtractionSolid("Slit4Solid", Slit4OutSolid, Hole4Solid);
821
822 fSlit1Logic = new G4LogicalVolume(Slit1Solid, fTungsten, "Slit1Logic");
823 fSlit2Logic = new G4LogicalVolume(Slit2Solid, fTungsten, "Slit2Logic");
824 fSlit3Logic = new G4LogicalVolume(Slit3Solid, fTungsten, "Slit3Logic");
825 fSlit4Logic = new G4LogicalVolume(Slit4Solid, fTungsten, "Slit4Logic");
826
827 if (fIWantSlits) {
828 G4cout << "Slit material: Tungsten" << G4endl;
829 G4cout << "Slit1 thickness: " << fSlit1Thickness / mm << " mm" << G4endl;
830 G4cout << "Slit2 thickness: " << fSlit2Thickness / mm << " mm" << G4endl;
831 G4cout << "Slit3 thickness: " << fSlit3Thickness / mm << " mm" << G4endl;
832 G4cout << "Slit4 thickness: " << fSlit4Thickness / mm << " mm" << G4endl;
833 G4cout << "Slit1 aperture: " << fSlit1xAperture / mm << " x " << fSlit1yAperture / mm << " mm2"
834 << G4endl;
835 G4cout << "Slit2 aperture: " << fSlit2xAperture / mm << " x " << fSlit2yAperture / mm << " mm2"
836 << G4endl;
837 G4cout << "Slit3 aperture: " << fSlit3xAperture / mm << " x " << fSlit3yAperture / mm << " mm2"
838 << G4endl;
839 G4cout << "Slit4 aperture: " << fSlit4xAperture / mm << " x " << fSlit4yAperture / mm << " mm2"
840 << G4endl;
841 }
842
843 G4VisAttributes* SlitlVisAttribute = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));
844 SlitlVisAttribute->SetForceSolid(true);
845 fSlit1Logic->SetVisAttributes(SlitlVisAttribute);
846 fSlit2Logic->SetVisAttributes(SlitlVisAttribute);
847 fSlit3Logic->SetVisAttributes(SlitlVisAttribute);
848 fSlit4Logic->SetVisAttributes(SlitlVisAttribute);
849
850 G4double Slit1z = fPhantomZ - fSlit1SampleDistance;
851 G4ThreeVector Slit1PositionVector = G4ThreeVector(0., 0., Slit1z);
852
853 G4double Slit2z = fPhantomZ - fSlit2SampleDistance;
854 G4ThreeVector Slit2PositionVector = G4ThreeVector(0., 0., Slit2z);
855
856 G4double Slit3x = fSlit3SampleDistance * std::sin(fthetaSetup);
857 G4double Slit3z = fPhantomZ + fSlit3SampleDistance * std::cos(fthetaSetup);
858 G4ThreeVector Slit3PositionVector = G4ThreeVector(Slit3x, 0., Slit3z);
859
860 G4double Slit4x = fSlit4SampleDistance * std::sin(fthetaSetup);
861 G4double Slit4z = fPhantomZ + fSlit4SampleDistance * std::cos(fthetaSetup);
862 G4ThreeVector Slit4PositionVector = G4ThreeVector(Slit4x, 0., Slit4z);
863
864 if (fIWantSlits) {
865 fSlit1Physical = new G4PVPlacement(0, Slit1PositionVector, fSlit1Logic, "Slit1Physical",
866 fWorldLogic, false, 0);
867
868 fSlit2Physical = new G4PVPlacement(0, Slit2PositionVector, fSlit2Logic, "Slit2Physical",
869 fWorldLogic, false, 0);
870
871 fSlit3Physical = new G4PVPlacement(SetupRotationMatrix, Slit3PositionVector, fSlit3Logic,
872 "Slit3Physical", fWorldLogic, false, 0);
873
874 fSlit4Physical = new G4PVPlacement(SetupRotationMatrix, Slit4PositionVector, fSlit4Logic,
875 "Slit4Physical", fWorldLogic, false, 0);
876 }
877
878 // Detector (with shielding)
879 G4Tubs* DetectorSolid = new G4Tubs("DetectorSolid", 0., fDetectorSize * 0.5,
880 fDetectorThickness * 0.5, 0. * deg, 360. * deg);
881
882 fDetectorLogic = new G4LogicalVolume(DetectorSolid, fGe, "DetectorLogic");
883
884 G4VisAttributes* DetectorVisAttribute = new G4VisAttributes(G4Colour(0., 0.5, 0.));
885 DetectorVisAttribute->SetForceSolid(true);
886 fDetectorLogic->SetVisAttributes(DetectorVisAttribute);
887
888 G4double Detx = fDetectorSampleDistance * std::sin(fthetaSetup);
889 G4double Detz = fPhantomZ + fDetectorSampleDistance * std::cos(fthetaSetup);
890 G4ThreeVector DetectorPositionVector = G4ThreeVector(Detx, 0., Detz);
891
892 fDetectorPhysical = new G4PVPlacement(SetupRotationMatrix, DetectorPositionVector, fDetectorLogic,
893 "DetectorPhysical", fWorldLogic, false, 0);
894
895 // Shielding
896 G4double ShieldingSize = fDetectorSize + 2 * fShieldingThickness;
897
898 G4double margin = 2.;
899 G4double ShieldingLength = fDetectorThickness + fShieldingThickness * margin;
900
901 G4double ShieldingSampleDistance = fDetectorSampleDistance + fDetectorThickness * 0.5
902 + fShieldingThickness - ShieldingLength * 0.5;
903
904 G4Tubs* ShieldingSolid = new G4Tubs("ShieldingSolid", fDetectorSize * 0.5, ShieldingSize * 0.5,
905 ShieldingLength * 0.5, 0. * deg, 360. * deg);
906
907 fShieldingLogic = new G4LogicalVolume(ShieldingSolid, fLead, "ShieldingLogic");
908
909 G4cout << "Shielding material: Lead" << G4endl;
910 G4cout << "Shielding thickness: " << fShieldingThickness / mm << " mm" << G4endl;
911
912 G4VisAttributes* ShieldingVisAttribute = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
913 ShieldingVisAttribute->SetForceSolid(true);
914 fShieldingLogic->SetVisAttributes(ShieldingVisAttribute);
915
916 G4double Shieldx = ShieldingSampleDistance * std::sin(fthetaSetup);
917 G4double Shieldz = fPhantomZ + ShieldingSampleDistance * std::cos(fthetaSetup);
918 G4ThreeVector ShieldingPositionVector = G4ThreeVector(Shieldx, 0., Shieldz);
919
920 G4double ShieldingBackSampleDistance =
921 fDetectorSampleDistance + fDetectorThickness * 0.5 + fShieldingThickness * 0.5;
922
923 G4Tubs* ShieldingBackSolid = new G4Tubs("ShieldingBackSolid", 0., fDetectorSize * 0.5,
924 fShieldingThickness * 0.5, 0. * deg, 360. * deg);
925
926 fShieldingBackLogic = new G4LogicalVolume(ShieldingBackSolid, fLead, "ShieldingBackLogic");
927
928 fShieldingBackLogic->SetVisAttributes(ShieldingVisAttribute);
929
930 G4double ShieldBackx = ShieldingBackSampleDistance * std::sin(fthetaSetup);
931 G4double ShieldBackz = fPhantomZ + ShieldingBackSampleDistance * std::cos(fthetaSetup);
932 G4ThreeVector ShieldingBackPositionVector = G4ThreeVector(ShieldBackx, 0., ShieldBackz);
933
934 fShieldingPhysical =
935 new G4PVPlacement(SetupRotationMatrix, ShieldingPositionVector, fShieldingLogic,
936 "ShieldingPhysical", fWorldLogic, false, 0);
937
938 fShieldingBackPhysical =
939 new G4PVPlacement(SetupRotationMatrix, ShieldingBackPositionVector, fShieldingBackLogic,
940 "ShieldingBackPhysical", fWorldLogic, false, 0);
941
942 return fWorldPhysical;
943 }
944
945 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
946
947 void SAXSDetectorConstruction::ConstructSDandField()
948 {
949 // Sensitive Volume
950 G4VSensitiveDetector* vDetector = new SAXSSensitiveDetector("det");
951 G4SDManager::GetSDMpointer()->AddNewDetector(vDetector);
952 fSensitiveVolume = fDetectorLogic;
953 fSensitiveVolume->SetSensitiveDetector(vDetector);
954 }
955
956 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
957