| 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 // Code developed by
27 // Silvia Pozzi (1), silvia.pozzi@iss.it
28 // Barbara Caccia (1), barbara.caccia@iss.it
29 // Carlo Mancini Terracciano (2), carlo.mancini.terracciano@roma1.infn.it
30 // (1) Istituto Superiore di Sanita' and INFN Roma, Italy
31 // (2) Univ. La Sapienza and INFN Roma, Italy
32
33 #include "DetectorConstruction.hh"
34 #include "DetectorMessenger.hh"
35
36 #include <G4LogicalVolume.hh>
37 #include <G4PVPlacement.hh>
38 #include <G4PVReplica.hh>
39 #include <G4ProductionCuts.hh>
40 #include <G4NistManager.hh>
41 #include <G4SystemOfUnits.hh>
42 #include <G4VisAttributes.hh>
43 #include <G4Box.hh>
44 #include <G4VSolid.hh>
45 #include <G4Cons.hh>
46 #include <G4Tubs.hh>
47 #include <G4BooleanSolid.hh>
48 #include <G4SubtractionSolid.hh>
49 #include <G4SDManager.hh>
50 #include <G4GlobalMagFieldMessenger.hh>
51 #include <G4UnionSolid.hh>
52 #include <G4SubtractionSolid.hh>
53 #include <G4RunManager.hh>
54 #include <G4GeometryManager.hh>
55 #include <G4PhysicalVolumeStore.hh>
56 #include <G4LogicalVolumeStore.hh>
57 #include <G4SolidStore.hh>
58 #include <G4Transform3D.hh>
59 #include <CLHEP/Vector/Rotation.h>
60
61
62 G4VPhysicalVolume* DetectorConstruction::Construct()
63 {
64 detectorMessenger = new DetectorMessenger(this);
65
66 //// Saturn 43
67 G4Material* air = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
68
69 G4double worldSizeX = 3.5 * m;
70 G4double worldSizeY = 3.5 * m;
71 G4double worldSizeZ = 3.5 * m;
72
73 G4VSolid* worldBox = new G4Box("world", worldSizeX/2., worldSizeY/2., worldSizeZ/2.);
74
75 G4LogicalVolume* worldLog = new G4LogicalVolume(worldBox, air, "world");
76
77 G4VisAttributes* visAttr = new G4VisAttributes();
78 visAttr->SetVisibility(false);
79 worldLog->SetVisAttributes(visAttr);
80
81 world_phys = new G4PVPlacement(nullptr, {}, worldLog, "world", nullptr, false, 0);
82
83 ConstructMaterials();
84 ConstructAccelerator();
85 ConstructPhantom();
86
87 return world_phys;
88 }
89
90 void DetectorConstruction::ConstructPhantom()
91 {
92 G4Material* water = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
93
94 G4ThreeVector phantomHalfDim = {phantomSideDim/2., phantomSideDim/2., phantomSideDim/2.};
95
96 G4VSolid* Phantom = new G4Box("Phantom", phantomHalfDim.getX(), phantomHalfDim.getY(), phantomHalfDim.getZ());
97 G4LogicalVolume* Phantom_log = new G4LogicalVolume(Phantom, water, "Phantom_log", 0, 0, 0);
98
99 G4VSolid* spessore = new G4Box("Spessore", phantomHalfDim.getX(), phantomHalfDim.getY(), 1.25*mm);
100 G4LogicalVolume* spessore_log = new G4LogicalVolume(spessore, water, "spessore_log", 0, 0, 0);
101 new G4PVPlacement(nullptr, {0.,0., 1.25*mm}, "spessore_phys", spessore_log, world_phys, false, 1);
102
103 // axes origin on the face of the phantom
104 G4ThreeVector phantomPos = {0., 0., phantomHalfDim.getZ()+voxelDepthDim/2};
105
106 phantom_phys = new G4PVPlacement(nullptr, phantomPos,
107 "phantom_phys", Phantom_log, world_phys, false, 1);
108
109 PhantomSegmentation(Phantom_log);
110
111 // Region for cuts
112 G4Region *regVol= new G4Region("fullWaterPhantomR");
113 G4ProductionCuts* cuts = new G4ProductionCuts;
114 cuts->SetProductionCut(0.1*mm);
115 regVol->SetProductionCuts(cuts);
116
117 Phantom_log -> SetRegion(regVol);
118 regVol->AddRootLogicalVolume(Phantom_log);
119
120 //--------- Visualization attributes -------------------------------
121 /// Phantom
122 G4VisAttributes* simpleH2OVisAtt= new G4VisAttributes(G4Colour::Cyan());
123 simpleH2OVisAtt->SetVisibility(true);
124 simpleH2OVisAtt->SetForceSolid(false);
125 Phantom_log->SetVisAttributes(simpleH2OVisAtt);
126
127 }
128
129 G4Material* DetectorConstruction::GetMaterial(const G4String materialName)
130 {
131 G4Material* material = nullptr;
132
133 if (materialName == "kapton")
134 {
135 material = mat_Kapton;
136 }
137 else if (materialName == "steel")
138 {
139 material = mat_Ssteel;
140 }
141 else if (materialName == "xc10")
142 {
143 material = mat_XC10;
144 }
145 else if (materialName == "wnicu")
146 {
147 material = mat_WNICU;
148 }
149 else
150 {
151 G4cerr << "* DetectorConstruction::GetMaterial: material " <<
152 materialName << " not found" << G4endl;
153 }
154
155 return material;
156 }
157
158 void DetectorConstruction::ConstructMaterials()
159 {
160 G4NistManager* nist = G4NistManager::Instance();
161
162 G4Material* el_H = nist -> FindOrBuildMaterial("G4_H");
163 G4Material* el_C = nist -> FindOrBuildMaterial("G4_C");
164 G4Material* el_N = nist -> FindOrBuildMaterial("G4_N");
165 G4Material* el_O = nist -> FindOrBuildMaterial("G4_O");
166 G4Material* el_Si = nist -> FindOrBuildMaterial("G4_Si");
167 G4Material* el_Cr = nist -> FindOrBuildMaterial("G4_Cr");
168 G4Material* el_Mn = nist -> FindOrBuildMaterial("G4_Mn");
169 G4Material* el_Fe = nist -> FindOrBuildMaterial("G4_Fe");
170 G4Material* el_Ni = nist -> FindOrBuildMaterial("G4_Ni");
171 G4Material* el_Cu = nist -> FindOrBuildMaterial("G4_Cu");
172 G4Material* el_W = nist -> FindOrBuildMaterial("G4_W");
173
174 mat_Ssteel = new G4Material("StainlessSteel", 7.8 *g/cm3, 6);
175 mat_Ssteel -> AddMaterial(el_Fe, 0.6898);
176 mat_Ssteel -> AddMaterial(el_Cr, 0.18);
177 mat_Ssteel -> AddMaterial(el_Ni, 0.10);
178 mat_Ssteel -> AddMaterial(el_Mn, 0.02);
179 mat_Ssteel -> AddMaterial(el_Si, 0.01);
180 mat_Ssteel -> AddMaterial(el_C, 0.0002);
181
182 mat_XC10 = new G4Material("CARBON_STEEL", 7.8 *g/cm3, 3);
183 mat_XC10 -> AddMaterial(el_Fe, 0.993);
184 mat_XC10 -> AddMaterial(el_Mn, 0.006);
185 mat_XC10 -> AddMaterial(el_C, 0.001);
186
187 mat_WNICU = new G4Material("Denal(WNICU)", 16.8*g/cm3, 3);
188 mat_WNICU -> AddMaterial(el_W, 0.905);
189 mat_WNICU -> AddMaterial(el_Ni, 0.07);
190 mat_WNICU -> AddMaterial(el_Cu, 0.025);
191
192 mat_Kapton = new G4Material("Kapton", 1.42*g/cm3, 4);
193 mat_Kapton -> AddMaterial(el_C, 69.1133 *perCent);
194 mat_Kapton -> AddMaterial(el_O, 20.9235 *perCent);
195 mat_Kapton -> AddMaterial(el_N, 7.3270 *perCent);
196 mat_Kapton -> AddMaterial(el_H, 2.6362 *perCent);
197
198 }
199
200 void DetectorConstruction::ConstructAccelerator()
201 {
202 G4bool checkOverlap = true;
203
204 G4Material* Vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_Galactic");
205
206 accHalfSize = {600.*mm, 600.*mm, 600.*mm};
207 G4ThreeVector initialCentre = {0.*mm, 0.*mm, -std::abs(sourceToSkinDistance)};
208
209 G4Box* accWorldBox = new G4Box("accWorld", accHalfSize.getX(), accHalfSize.getY(), accHalfSize.getZ());
210 G4LogicalVolume* accWorldLV = new G4LogicalVolume(accWorldBox, Vacuum, "accWorldL", 0, 0, 0);
211 accWorld_phys = new G4PVPlacement(nullptr, initialCentre, "acceleratorBox", accWorldLV, world_phys, false, checkOverlap);
212
213 G4VisAttributes* simpleAlSVisAtt = new G4VisAttributes(G4Colour::White());
214 simpleAlSVisAtt -> SetVisibility(false);
215 accWorldLV -> SetVisAttributes(simpleAlSVisAtt);
216
217 ConstructTarget();
218 ConstructPrimaryCollimator();
219 ConstructVacuumWindow();
220 ConstructFlatteningFilter();
221 ConstructIonizationChamber();
222 ConstructJawsX();
223 ConstructJawsY();
224 }
225
226 void DetectorConstruction::ConstructTarget()
227 {
228 G4bool checkOverlap = true;
229
230 G4ThreeVector targetCentre, boxHalfSide, tubeCentre;
231 G4double tubeRadius, tubeHalfLengthZ;
232 targetCentre.set(0, 0, -3.5*mm);
233 boxHalfSide.set(5.*mm, 5.*mm, 7.5*mm);
234
235 tubeRadius = 3.*mm;
236 tubeHalfLengthZ = 5.5*mm; // 11 mm tot
237
238 G4Material* diskMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_Ti");
239 G4Material* targetMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_W");
240
241 // Physical volumes
242 G4Box* box = new G4Box("targetBox", boxHalfSide.getX(), boxHalfSide.getY(), boxHalfSide.getZ());
243 G4Tubs* tube = new G4Tubs("targetTube", 0.*mm, tubeRadius, tubeHalfLengthZ, 0.*deg, 360.*deg);
244
245 G4SubtractionSolid* BoxMinusTube = new G4SubtractionSolid("TargetSolid", box, tube, nullptr, {0.,0.,-2*mm});
246 G4LogicalVolume* BoxMinusTubeLV = new G4LogicalVolume(BoxMinusTube, targetMaterial, "BoxMinusTubeLV",0,0,0);
247 new G4PVPlacement(0, targetCentre, "TargetPV", BoxMinusTubeLV, accWorld_phys, false, checkOverlap);
248
249 G4Tubs* disk = new G4Tubs("targetDisk", 0.*mm, 4.*mm, 0.025*mm, 0.*deg, 360.*deg);
250 G4LogicalVolume* diskLV = new G4LogicalVolume(disk, diskMaterial, "targetDiskLV",0,0,0);
251 new G4PVPlacement(0, {0.,0.,-15.*mm}, "diskPV", diskLV, accWorld_phys, false, checkOverlap);
252
253 // Visualization
254 G4VisAttributes* simpleAlSVisAtt;
255 simpleAlSVisAtt = new G4VisAttributes(G4Colour::Grey());
256 simpleAlSVisAtt -> SetVisibility(true);
257 simpleAlSVisAtt -> SetForceSolid(true);
258 BoxMinusTubeLV -> SetVisAttributes(simpleAlSVisAtt);
259 diskLV -> SetVisAttributes(simpleAlSVisAtt);
260
261 // Region for cuts
262 G4Region *regVol;
263 regVol = new G4Region("TargetR");
264 G4ProductionCuts* cuts = new G4ProductionCuts;
265 cuts -> SetProductionCut(1.*mm);
266 regVol -> SetProductionCuts(cuts);
267 BoxMinusTubeLV -> SetRegion(regVol);
268 regVol -> AddRootLogicalVolume(BoxMinusTubeLV);
269 diskLV -> SetRegion(regVol);
270 regVol -> AddRootLogicalVolume(diskLV);
271 }
272
273 void DetectorConstruction::ConstructPrimaryCollimator()
274 {
275 G4bool checkOverlap = true;
276
277 G4double bigTube1HalfLengthZ, bigTube1Radius, bigTube1Z;
278 G4double bigTube2HalfLengthZ, bigTube2Radius, bigTube2Z;
279 G4double bigTube3HalfLengthZ, bigTube3Radius, bigTube3Z;
280
281 G4double tube1Radius, tube1HalfLengthZ, tube2Radius, tube2HalfLengthZ;
282
283 G4double cone1RadiusMin, cone1RadiusMax, cone1HalfLengthZ;
284 G4double cone2RadiusMin, cone2RadiusMax, cone2HalfLengthZ;
285 G4double cone3RadiusMin, cone3RadiusMax, cone3HalfLengthZ;
286
287 // upper principal tube = bigTube1
288 bigTube1Radius = 141./2.*mm;
289 bigTube1HalfLengthZ = 79./2. * mm;
290 bigTube1Z = 5. * mm + bigTube1HalfLengthZ;
291 tube1Radius = 34./2. * mm;
292 tube1HalfLengthZ = 15.86/2. * mm;
293 cone1RadiusMin = 34./2.*mm; // upper cone
294 cone1RadiusMax = 54./2.*mm;
295 cone1HalfLengthZ = (bigTube1HalfLengthZ - tube1HalfLengthZ)*mm; //(79.-15.86)/2.*mm;
296
297 // middle principal tube = bigTube2
298 bigTube2Radius = 141./2. * mm;
299 bigTube2HalfLengthZ = 57.5/2. * mm;
300 bigTube2Z = bigTube1Z + bigTube1HalfLengthZ + bigTube2HalfLengthZ;
301 tube2Radius = 53./2. * mm;
302 tube2HalfLengthZ = 1.35/2. * mm;
303 cone2RadiusMin = 53./2.*mm; // middle cone
304 cone2RadiusMax = 81./2.*mm;
305 cone2HalfLengthZ = bigTube2HalfLengthZ - tube2HalfLengthZ; // (57.50 - 1.35)/2.*mm;
306
307 // lower principal tube = bigTube3
308 bigTube3Radius = 241./2.*mm;
309 bigTube3HalfLengthZ = 40.50/2.*mm;
310 bigTube3Z = bigTube2Z + bigTube2HalfLengthZ + 19.5 + bigTube3HalfLengthZ;
311 cone3RadiusMin = 88.06/2.*mm; // lower cone
312 cone3RadiusMax = 109.76/2.*mm;
313 cone3HalfLengthZ = 40.50/2.*mm;
314
315 G4Material* upperTubeMaterial = GetMaterial("xc10");
316 G4Material* middleTubeMaterial = GetMaterial("wnicu");
317 G4Material* lowerTubeMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
318
319 // Physical volumes
320 G4Tubs* bigTube1 = new G4Tubs("PriCollBigTube1", 0.*mm, bigTube1Radius, bigTube1HalfLengthZ, 0.*deg, 360.*deg);
321 G4Tubs* tube1 = new G4Tubs("PriCollTube1", 0.*mm, tube1Radius, tube1HalfLengthZ, 0.*deg, 360.*deg);
322 G4Cons* cone1 = new G4Cons("PriCollCone1", 0., cone1RadiusMin, 0., cone1RadiusMax, cone1HalfLengthZ, 0, 360.*deg);
323
324 G4Tubs* bigTube2 = new G4Tubs("PriCollBigTube2", 0.*mm, bigTube2Radius, bigTube2HalfLengthZ, 0.*deg, 360.*deg);
325 G4Tubs* tube2 = new G4Tubs("PriCollTube2", 0.*mm, tube2Radius, tube2HalfLengthZ, 0.*deg, 360.*deg);
326 G4Cons* cone2 = new G4Cons("PriCollCone2", 0., cone2RadiusMin, 0., cone2RadiusMax, cone2HalfLengthZ, 0, 360.*deg);
327
328 G4Tubs* bigTube3 = new G4Tubs("PriCollBigTube3", 0.*mm, bigTube3Radius, bigTube3HalfLengthZ, 0.*deg, 360.*deg);
329 G4Cons* cone3 = new G4Cons("PriCollCone3", 0., cone3RadiusMin, 0., cone3RadiusMax, cone3HalfLengthZ, 0, 360.*deg);
330
331 G4UnionSolid* tube1AndCone1 = new G4UnionSolid("PriCollTube1AndCone1",
332 tube1, cone1, nullptr, {0., 0., (tube1HalfLengthZ + cone1HalfLengthZ)});
333 G4SubtractionSolid* bigTube1NotTube1AndCone1 = new G4SubtractionSolid(
334 "PriColltube1NotTube1AndCone1", bigTube1, tube1AndCone1, 0,
335 {0., 0., (-bigTube1HalfLengthZ + tube1HalfLengthZ)});
336 G4LogicalVolume* bigTube1NotTube1AndCone1LV = new G4LogicalVolume(
337 bigTube1NotTube1AndCone1, upperTubeMaterial, "PriCollBigTube1NotTube1AndCone1LV", 0, 0, 0);
338
339 G4UnionSolid* tube2AndCone2 = new G4UnionSolid("PriCollTube2AndCone2",
340 tube2, cone2, nullptr, {0., 0., (tube2HalfLengthZ + cone2HalfLengthZ)});
341 G4SubtractionSolid* bigTube2NotTube2AndCone2 = new G4SubtractionSolid(
342 "PriCollBigTube2NotTube2Cone2", bigTube2, tube2AndCone2, 0,
343 {0., 0., (-bigTube2HalfLengthZ + tube2HalfLengthZ)});
344 G4LogicalVolume* bigTube2NotTube2AndCone2LV = new G4LogicalVolume(
345 bigTube2NotTube2AndCone2, middleTubeMaterial, "PriCollTube2NotTube2Cone2LV", 0, 0, 0);
346
347 G4SubtractionSolid* bigTube3NotCone3 = new G4SubtractionSolid("PriCollTube4NotCone3", bigTube3, cone3);
348 G4LogicalVolume* bigTube3NotCone3LV = new G4LogicalVolume(bigTube3NotCone3,
349 lowerTubeMaterial, "PriCollBigTube3NotCone3LV",0,0,0);
350
351 new G4PVPlacement(nullptr, {0, 0, bigTube1Z}, "PriCollUpperPV",
352 bigTube1NotTube1AndCone1LV, accWorld_phys, false, checkOverlap);
353 new G4PVPlacement(nullptr, {0, 0, bigTube2Z}, "PriCollMiddlePV",
354 bigTube2NotTube2AndCone2LV, accWorld_phys, false, checkOverlap);
355 new G4PVPlacement(nullptr, {0, 0, bigTube3Z}, "PriCollLowerPV",
356 bigTube3NotCone3LV, accWorld_phys, false, checkOverlap);
357
358 // Visualization
359 G4VisAttributes* simpleAlSVisAtt1 = new G4VisAttributes(G4Colour::Green());
360 G4VisAttributes* simpleAlSVisAtt2 = new G4VisAttributes(G4Colour::Red());
361 G4VisAttributes* simpleAlSVisAtt3 = new G4VisAttributes(G4Colour::Blue());
362 simpleAlSVisAtt1 -> SetVisibility(true);
363 simpleAlSVisAtt1 -> SetForceSolid(false);
364 simpleAlSVisAtt1 -> SetForceWireframe(true);
365 bigTube1NotTube1AndCone1LV -> SetVisAttributes(simpleAlSVisAtt1); //primary collimator
366 simpleAlSVisAtt2 -> SetVisibility(true);
367 simpleAlSVisAtt2 -> SetForceSolid(false);
368 simpleAlSVisAtt2 -> SetForceWireframe(true);
369 bigTube2NotTube2AndCone2LV -> SetVisAttributes(simpleAlSVisAtt2); //primary collimator
370 simpleAlSVisAtt3 -> SetVisibility(true);
371 simpleAlSVisAtt3 -> SetForceSolid(false);
372 simpleAlSVisAtt3 -> SetForceWireframe(true);
373 bigTube3NotCone3LV -> SetVisAttributes(simpleAlSVisAtt3); //secondary collimator
374
375 // Region for cuts
376 G4Region* regVol = new G4Region("primaryCollimatorR");
377 G4ProductionCuts* cuts = new G4ProductionCuts;
378 cuts->SetProductionCut(1.*cm);
379 regVol->SetProductionCuts(cuts);
380 bigTube1NotTube1AndCone1LV -> SetRegion(regVol);
381 regVol -> AddRootLogicalVolume(bigTube1NotTube1AndCone1LV);
382 bigTube2NotTube2AndCone2LV -> SetRegion(regVol);
383 regVol -> AddRootLogicalVolume(bigTube2NotTube2AndCone2LV);
384 bigTube3NotCone3LV -> SetRegion(regVol);
385 regVol -> AddRootLogicalVolume(bigTube3NotCone3LV);
386 }
387
388 void DetectorConstruction::ConstructFlatteningFilter()
389 {
390 G4double tube1HalfLengthZ; //tube1Radius,
391 G4double cone1RadiusMin, cone1RadiusMax, cone1HalfLengthZ;
392 G4double cone2RadiusMin, cone2RadiusMax, cone2HalfLengthZ;
393 G4double cone3RadiusMin, cone3RadiusMax, cone3HalfLengthZ, ffZ;
394
395 tubeFFRadius = 108./2.*mm; // top principal tube
396 tube1HalfLengthZ = 7.5/2.*mm;
397
398 cone1RadiusMin = 54./2.*mm; // top cone
399 cone1RadiusMax = 76./2.*mm;
400 cone1HalfLengthZ = 13.7/2.*mm;
401
402 cone2RadiusMin = 8./2.*mm; // mid cone
403 cone2RadiusMax = 54./2.*mm;
404 cone2HalfLengthZ = (44.3-13.7)/2.*mm;
405
406 cone3RadiusMin = 0.000001*mm; // bottom cone
407 cone3RadiusMax = 8./2.*mm;
408 cone3HalfLengthZ = (46.8-44.3)/2.*mm;
409 tubeFFFirstFaceZ = 149.50*mm;
410
411 ffZ = tubeFFFirstFaceZ + tube1HalfLengthZ; // the half point is located ad the centre of the tube1 because of the solids unions and translations
412
413 G4Material* ffMaterial = GetMaterial("steel");
414
415 // Physical volumes
416 G4Tubs* tube1 = new G4Tubs("ffTube1", 0.*mm, tubeFFRadius, tube1HalfLengthZ, 0.*deg, 360.*deg);
417 G4Cons* cone1 = new G4Cons("ffCone1", 0., cone1RadiusMax, 0., cone1RadiusMin, cone1HalfLengthZ, 0, 360.*deg);
418 G4Cons* cone2 = new G4Cons("ffCone2", 0., cone2RadiusMax, 0., cone2RadiusMin, cone2HalfLengthZ, 0, 360.*deg);
419 G4Cons* cone3 = new G4Cons("ffCone3", 0., cone3RadiusMax, 0., cone3RadiusMin, cone3HalfLengthZ, 0, 360.*deg);
420
421 G4double pos = (cone1HalfLengthZ - tube1HalfLengthZ - 1.)*mm; // == (13.7/2. - 7.5/2. -1.)*mm
422 G4UnionSolid *tube1AndCone1 = new G4UnionSolid("ffTube1AndCone1", tube1,
423 cone1, 0, {0., 0., pos});
424
425 pos += cone1HalfLengthZ + cone2HalfLengthZ;
426 G4UnionSolid* tubeCone1AndCone2 = new G4UnionSolid("fftubeCone1AndCone2",
427 tube1AndCone1, cone2, 0, {0., 0., pos});
428
429 pos += cone2HalfLengthZ + cone3HalfLengthZ;
430 G4UnionSolid* tubeCone12AndCone3 = new G4UnionSolid("fftubeCone12AndCone3",
431 tubeCone1AndCone2, cone3, 0, {0., 0., pos});
432
433 G4LogicalVolume* ffLV = new G4LogicalVolume(tubeCone12AndCone3, ffMaterial, "ffLV",0,0,0);
434
435 new G4PVPlacement(0, {0,0,ffZ}, "ffPV", ffLV, accWorld_phys, false, 0);
436
437 // Visualization
438 G4VisAttributes* simpleAlSVisAtt= new G4VisAttributes(G4Colour::Magenta());
439 simpleAlSVisAtt->SetVisibility(true);
440 simpleAlSVisAtt->SetForceSolid(true);
441 ffLV->SetVisAttributes(simpleAlSVisAtt);
442
443 // Region for cuts
444 G4Region* regVol;
445 regVol= new G4Region("ffR");
446 G4ProductionCuts* cuts = new G4ProductionCuts;
447 cuts->SetProductionCut(0.1*cm);
448 regVol->SetProductionCuts(cuts);
449 ffLV->SetRegion(regVol);
450 regVol->AddRootLogicalVolume(ffLV);
451 }
452
453 void DetectorConstruction::ConstructIonizationChamber()
454 {
455 G4bool checkOverlap = true;
456
457 G4double tubeRadius, tubeA1Z, tubeA2Z, tubeA3Z, tubeA4Z, tubeA5Z, tubeA6Z;
458 G4double kaptonThickness, AlThickness1, AlThickness8;
459
460 kaptonThickness = 0.025/2. *mm;
461 AlThickness1 = 1.6e-5/2. *cm;
462 AlThickness8 = 8e-6/2. *cm;
463
464 tubeRadius = 110./2.*mm; // upper principal tube
465 tubeA1Z = (202.5 + AlThickness8)*mm;
466 tubeA2Z = (204.5 + AlThickness1)*mm;
467 tubeA3Z = (206.5 + AlThickness8)*mm;
468 tubeA4Z = (207.5 + AlThickness8)*mm;
469 tubeA5Z = (209.5 + AlThickness1)*mm;
470 tubeA6Z = (211.5 + AlThickness8)*mm;
471
472 G4double pos1 = (AlThickness1 + kaptonThickness)*mm;
473 G4double pos8 = (AlThickness8 + kaptonThickness)*mm;
474
475 G4Material* el_Al = G4NistManager::Instance()->FindOrBuildMaterial("G4_Al");
476 G4Material* Kapton = GetMaterial("kapton");
477
478 // Physical volumes
479 G4Tubs* tubeAl8 = new G4Tubs("ICTube1A", 0.*mm, tubeRadius, AlThickness8, 0.*deg, 360.*deg);
480 G4Tubs* tubeK = new G4Tubs("ICTube1B", 0.*mm, tubeRadius, kaptonThickness, 0.*deg, 360.*deg);
481 G4Tubs* tubeAl1 = new G4Tubs("ICTube2A", 0.*mm, tubeRadius, AlThickness1, 0.*deg, 360.*deg);
482
483 G4LogicalVolume* IC_Al8_LV = new G4LogicalVolume(tubeAl8, el_Al, "IC_Al8_LV", 0, 0, 0);
484 G4LogicalVolume* IC_Al1_LV = new G4LogicalVolume(tubeAl1, el_Al, "IC_Al1_LV", 0, 0, 0);
485 G4LogicalVolume* IC_K_LV = new G4LogicalVolume(tubeK, Kapton, "IC_Al_LV", 0, 0, 0);
486
487 new G4PVPlacement(0, {0,0,tubeA1Z}, "IC_AL8_PV1", IC_Al8_LV, accWorld_phys, false, checkOverlap);
488 new G4PVPlacement(0, {0,0,tubeA1Z + pos8}, "IC_K_PV1", IC_K_LV, accWorld_phys, false, checkOverlap);
489
490 new G4PVPlacement(0, {0,0,tubeA2Z}, "IC_AL1_PV2", IC_Al1_LV, accWorld_phys, false, checkOverlap);
491 new G4PVPlacement(0, {0,0,tubeA2Z + pos1}, "IC_K_PV2", IC_K_LV, accWorld_phys, false, checkOverlap);
492
493 new G4PVPlacement(0, {0,0,tubeA3Z}, "IC_AL8_PV3", IC_Al8_LV, accWorld_phys, false, checkOverlap);
494 new G4PVPlacement(0, {0,0,tubeA3Z + pos8}, "IC_K_PV3", IC_K_LV, accWorld_phys, false, checkOverlap);
495
496 new G4PVPlacement(0, {0,0,tubeA4Z}, "IC_AL8_PV4", IC_Al8_LV, accWorld_phys, false, checkOverlap);
497 new G4PVPlacement(0, {0,0,tubeA4Z + pos8}, "IC_K_PV4", IC_K_LV, accWorld_phys, false, checkOverlap);
498
499 new G4PVPlacement(0, {0,0,tubeA5Z}, "IC_AL1_PV5", IC_Al1_LV, accWorld_phys, false, checkOverlap);
500 new G4PVPlacement(0, {0,0,tubeA5Z + pos1}, "IC_K_PV5", IC_K_LV, accWorld_phys, false, checkOverlap);
501
502 new G4PVPlacement(0, {0,0,tubeA6Z}, "IC_AL8_PV6", IC_Al8_LV, accWorld_phys, false, checkOverlap);
503 new G4PVPlacement(0, {0,0,tubeA6Z + pos8}, "IC_K_PV6", IC_K_LV, accWorld_phys, false, checkOverlap);
504
505 // Visualization
506 G4VisAttributes* simpleAlSVisAttAl1 = new G4VisAttributes(G4Colour::Red());
507 simpleAlSVisAttAl1 -> SetVisibility(true);
508 simpleAlSVisAttAl1 -> SetForceSolid(true);
509 IC_Al1_LV -> SetVisAttributes(simpleAlSVisAttAl1);
510
511 G4VisAttributes* simpleAlSVisAttAl8 = new G4VisAttributes(G4Colour::Magenta());
512 simpleAlSVisAttAl8 -> SetVisibility(true);
513 simpleAlSVisAttAl8 -> SetForceSolid(true);
514 IC_Al8_LV -> SetVisAttributes(simpleAlSVisAttAl8);
515
516 G4VisAttributes* simpleAlSVisAttK = new G4VisAttributes(G4Colour::Blue());
517 simpleAlSVisAttK -> SetVisibility(true);
518 simpleAlSVisAttK -> SetForceSolid(true);
519 IC_K_LV -> SetVisAttributes(simpleAlSVisAttK);
520
521 // Region for cuts
522 G4Region * regVol = new G4Region("IonChamberR");
523 G4ProductionCuts* cuts = new G4ProductionCuts;
524 cuts -> SetProductionCut(0.1*mm);
525 regVol -> SetProductionCuts(cuts);
526
527 IC_Al1_LV -> SetRegion(regVol);
528 regVol -> AddRootLogicalVolume(IC_Al1_LV);
529 IC_Al8_LV -> SetRegion(regVol);
530 regVol -> AddRootLogicalVolume(IC_Al8_LV);
531 IC_K_LV -> SetRegion(regVol);
532 regVol -> AddRootLogicalVolume(IC_K_LV);
533 }
534
535 void DetectorConstruction::ConstructVacuumWindow()
536 {
537 G4double tubeRadius, tubeHalfLengthZ, tubeZ;
538
539 tubeRadius = 116.53/2.*mm; // upper principal tube
540 tubeHalfLengthZ = 2./2.*mm;
541 tubeZ = 215.75*mm + tubeHalfLengthZ;
542
543 G4Material* el_Al= G4NistManager::Instance()->FindOrBuildMaterial("G4_Al");
544
545 // Physical volumes
546 G4Tubs* tube = new G4Tubs("VacWindowTube", 0.*mm, tubeRadius, tubeHalfLengthZ, 0.*deg, 360.*deg);
547 G4LogicalVolume* tubeLV = new G4LogicalVolume(tube, el_Al, "VacWindowTubeLV",0,0,0);
548
549 new G4PVPlacement(0, {0,0,tubeZ}, "VacWindowTubePV", tubeLV, accWorld_phys, false, 0);
550
551 // Visualization
552 G4VisAttributes* simpleAlSVisAtt = new G4VisAttributes(G4Colour::Green());
553 simpleAlSVisAtt -> SetVisibility(true);
554 simpleAlSVisAtt -> SetForceSolid(true);
555 tubeLV -> SetVisAttributes(simpleAlSVisAtt);
556
557 // Region for cuts
558 G4Region* regVol= new G4Region("VacWindowR");
559 G4ProductionCuts* cuts = new G4ProductionCuts;
560 cuts->SetProductionCut(0.1*cm);
561 regVol->SetProductionCuts(cuts);
562 tubeLV -> SetRegion(regVol);
563 regVol -> AddRootLogicalVolume(tubeLV);
564 }
565 void DetectorConstruction::ConstructJawsX()
566 {
567 G4bool checkOverlap = false;
568
569 G4double boxSide = 101.*mm;
570 G4double box1HalfLengthZ = 3.*mm;
571 G4double box2HalfLengthZ = 35.*mm;
572 G4double box3HalfLengthZ = 35.*mm;
573 G4double box4HalfLengthZ = 27.*mm;
574 G4double box5HalfLengthZ = 10.*mm;
575 G4double box6HalfLengthZ = 5.*mm;
576
577 G4Material* el_Pb = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
578 G4Material* XC10 = GetMaterial("xc10");
579 G4Material* WNICU = GetMaterial("wnicu");
580
581 G4ThreeVector boxJawXSide = {boxSide, boxSide, 212.*mm};
582
583 G4Box* boxJaw1X = new G4Box("Jaw1Xbox", boxJawXSide.getX()/2., boxJawXSide.getY()/2., boxJawXSide.getZ()/2.);
584 G4LogicalVolume* boxJaw1XLV = new G4LogicalVolume(boxJaw1X,
585 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw1XboxLV", 0, 0, 0);
586 G4Box* boxJaw2X = new G4Box("Jaw2Xbox", boxJawXSide.getX()/2., boxJawXSide.getY()/2., boxJawXSide.getZ()/2.);
587 G4LogicalVolume* boxJaw2XLV = new G4LogicalVolume(boxJaw2X,
588 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw2XboxLV", 0, 0, 0);
589
590 jaw1XInitialPos.setX(boxJawXSide.getX()/2.);
591 jaw1XInitialPos.setY(0.);
592 jaw1XInitialPos.setZ(275.5*mm + boxJawXSide.getZ()/2.);
593
594 G4RotationMatrix* cRotationJaw1X = new G4RotationMatrix();
595 cRotationJaw1X -> rotateY(std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
596
597 G4ThreeVector position = *cRotationJaw1X * jaw1XInitialPos;
598 *cRotationJaw1X = cRotationJaw1X -> inverse();
599
600 boxJaw1X_phys =
601 new G4PVPlacement (cRotationJaw1X, position,
602 "Jaw1XPV", boxJaw1XLV, accWorld_phys, false, 0, checkOverlap);
603
604 jaw2XInitialPos.setX(-jaw1XInitialPos.getX());
605 jaw2XInitialPos.setY(jaw1XInitialPos.getY());
606 jaw2XInitialPos.setZ(jaw1XInitialPos.getZ());
607
608 G4RotationMatrix* cRotationJaw2X = new G4RotationMatrix();
609 cRotationJaw2X -> rotateY(-std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
610
611 position = *cRotationJaw2X * jaw2XInitialPos;
612 *cRotationJaw2X = cRotationJaw2X -> inverse();
613
614 boxJaw2X_phys =
615 new G4PVPlacement (cRotationJaw2X, position,
616 "Jaw2XPV", boxJaw2XLV, accWorld_phys, false, 0, 0);
617
618 // Physical volumes
619 G4Box* box1 = new G4Box("Jaws1XBox1", boxSide/2., boxSide/2., box1HalfLengthZ/2.);
620 G4Box* box2 = new G4Box("Jaws1XBox2", boxSide/2., boxSide/2., box2HalfLengthZ/2.);
621 G4Box* box3 = new G4Box("Jaws1XBox3", boxSide/2., boxSide/2., box3HalfLengthZ/2.);
622 G4Box* box4 = new G4Box("Jaws1XBox4", boxSide/2., boxSide/2., box4HalfLengthZ/2.);
623 G4Box* box5 = new G4Box("Jaws1XBox5", boxSide/2., boxSide/2., box5HalfLengthZ/2.);
624 G4Box* box6 = new G4Box("Jaws1XBox6", boxSide/2., boxSide/2., box6HalfLengthZ/2.);
625 G4LogicalVolume* box1LV1 = new G4LogicalVolume(box1, XC10, "Jaws1XLV1", 0, 0, 0);
626 G4LogicalVolume* box2LV1 = new G4LogicalVolume(box2, el_Pb, "Jaws1XLV2", 0, 0, 0);
627 G4LogicalVolume* box3LV1 = new G4LogicalVolume(box3, WNICU, "Jaws1XLV3", 0, 0, 0);
628 G4LogicalVolume* box4LV1 = new G4LogicalVolume(box4, el_Pb, "Jaws1XLV4", 0, 0, 0);
629 G4LogicalVolume* box5LV1 = new G4LogicalVolume(box5, el_Pb, "Jaws1XLV5", 0, 0, 0);
630 G4LogicalVolume* box6LV1 = new G4LogicalVolume(box6, WNICU, "Jaws1XLV6", 0, 0, 0);
631
632 G4LogicalVolume* box1LV2 = new G4LogicalVolume(box1, XC10, "Jaws2XLV1", 0, 0, 0);
633 G4LogicalVolume* box2LV2 = new G4LogicalVolume(box2, el_Pb, "Jaws2XLV2", 0, 0, 0);
634 G4LogicalVolume* box3LV2 = new G4LogicalVolume(box3, WNICU, "Jaws2XLV3", 0, 0, 0);
635 G4LogicalVolume* box4LV2 = new G4LogicalVolume(box4, el_Pb, "Jaws2XLV4", 0, 0, 0);
636 G4LogicalVolume* box5LV2 = new G4LogicalVolume(box5, el_Pb, "Jaws2XLV5", 0, 0, 0);
637 G4LogicalVolume* box6LV2 = new G4LogicalVolume(box6, WNICU, "Jaws2XLV6", 0, 0, 0);
638
639 G4ThreeVector centre = {0., 0., 0.};
640 G4double zCentreCurrentBox = -boxJawXSide.getZ()/2.+box1HalfLengthZ/2.*mm;
641
642 centre.setZ(zCentreCurrentBox);
643 new G4PVPlacement(nullptr, centre, "Jaws1XPV1", box1LV1, boxJaw1X_phys, false, 0, checkOverlap);
644 new G4PVPlacement(nullptr, centre, "Jaws2XPV1", box1LV2, boxJaw2X_phys, false, 0, checkOverlap);
645
646 zCentreCurrentBox += (box1HalfLengthZ+box2HalfLengthZ)/2.*mm;
647 centre.setZ(zCentreCurrentBox);
648 new G4PVPlacement(nullptr, centre, "Jaws1XPV2", box2LV1, boxJaw1X_phys, false, 0, checkOverlap);
649 new G4PVPlacement(nullptr, centre, "Jaws2XPV2", box2LV2, boxJaw2X_phys, false, 0, checkOverlap);
650
651 zCentreCurrentBox += (box2HalfLengthZ+box3HalfLengthZ)/2.*mm;
652 centre.setZ(zCentreCurrentBox);
653 new G4PVPlacement(nullptr, centre, "Jaws1XPV3", box3LV1, boxJaw1X_phys, false, 0, checkOverlap);
654 new G4PVPlacement(nullptr, centre, "Jaws2XPV3", box3LV2, boxJaw2X_phys, false, 0, checkOverlap);
655
656 zCentreCurrentBox += (box3HalfLengthZ+box4HalfLengthZ)/2.*mm;
657 centre.setZ(zCentreCurrentBox);
658 new G4PVPlacement(nullptr, centre, "Jaws1XPV4", box4LV1, boxJaw1X_phys, false, 0, checkOverlap);
659 new G4PVPlacement(nullptr, centre, "Jaws2XPV4", box4LV2, boxJaw2X_phys, false, 0, checkOverlap);
660
661 zCentreCurrentBox += (box4HalfLengthZ+box5HalfLengthZ)/2.+97.*mm;
662 centre.setZ(zCentreCurrentBox);
663 new G4PVPlacement(nullptr, centre, "Jaws1XPV5", box5LV1, boxJaw1X_phys, false, 0, checkOverlap);
664 new G4PVPlacement(nullptr, centre, "Jaws2XPV5", box5LV2, boxJaw2X_phys, false, 0, checkOverlap);
665
666 zCentreCurrentBox += (box5HalfLengthZ+box6HalfLengthZ)/2.;
667 centre.setZ(zCentreCurrentBox);
668 new G4PVPlacement(nullptr, centre, "Jaws1XPV6", box6LV1, boxJaw1X_phys, false, 0, checkOverlap);
669 new G4PVPlacement(nullptr, centre, "Jaws2XPV6", box6LV2, boxJaw2X_phys, false, 0, checkOverlap);
670
671 // Visualization
672 G4VisAttributes* visAtt = new G4VisAttributes(G4Colour::Yellow());
673 visAtt -> SetVisibility(false);
674 visAtt -> SetForceSolid(false);
675
676 G4VisAttributes* simpleAlSVisAttPb = new G4VisAttributes(G4Colour::Blue());
677 simpleAlSVisAttPb -> SetVisibility(true);
678 simpleAlSVisAttPb -> SetForceSolid(true);
679
680 G4VisAttributes* simpleAlSVisAttXC10 = new G4VisAttributes(G4Colour::Green());
681 simpleAlSVisAttXC10 -> SetVisibility(true);
682 simpleAlSVisAttXC10 ->SetForceSolid(true);
683
684 G4VisAttributes* simpleAlSVisAttWNICU = new G4VisAttributes(G4Colour::Red());
685 simpleAlSVisAttWNICU ->SetVisibility(true);
686 simpleAlSVisAttWNICU ->SetForceSolid(true);
687
688 box1LV1 -> SetVisAttributes(simpleAlSVisAttXC10);
689 box2LV1 -> SetVisAttributes(simpleAlSVisAttPb);
690 box3LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
691 box4LV1 -> SetVisAttributes(simpleAlSVisAttPb);
692 box5LV1 -> SetVisAttributes(simpleAlSVisAttPb);
693 box6LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
694
695 box1LV2 -> SetVisAttributes(simpleAlSVisAttXC10);
696 box2LV2 -> SetVisAttributes(simpleAlSVisAttPb);
697 box3LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
698 box4LV2 -> SetVisAttributes(simpleAlSVisAttPb);
699 box5LV2 -> SetVisAttributes(simpleAlSVisAttPb);
700 box6LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
701
702 boxJaw1XLV -> SetVisAttributes(visAtt);
703 boxJaw2XLV -> SetVisAttributes(visAtt);
704
705 // Region for cuts
706 G4Region* regVol = new G4Region("JawsXR");
707 G4ProductionCuts* cuts = new G4ProductionCuts;
708 cuts -> SetProductionCut(2.*cm);
709 regVol -> SetProductionCuts(cuts);
710 box1LV1 -> SetRegion(regVol);
711 regVol -> AddRootLogicalVolume(box1LV1);
712 box2LV1 -> SetRegion(regVol);
713 regVol -> AddRootLogicalVolume(box2LV1);
714 box3LV1 -> SetRegion(regVol);
715 regVol -> AddRootLogicalVolume(box3LV1);
716 box4LV1 -> SetRegion(regVol);
717 regVol -> AddRootLogicalVolume(box4LV1);
718 box5LV1 -> SetRegion(regVol);
719 regVol -> AddRootLogicalVolume(box5LV1);
720 box6LV1 -> SetRegion(regVol);
721 regVol -> AddRootLogicalVolume(box6LV1);
722
723 box1LV2 -> SetRegion(regVol);
724 regVol -> AddRootLogicalVolume(box1LV2);
725 box2LV2 -> SetRegion(regVol);
726 regVol -> AddRootLogicalVolume(box2LV2);
727 box3LV2 -> SetRegion(regVol);
728 regVol -> AddRootLogicalVolume(box3LV2);
729 box4LV2 -> SetRegion(regVol);
730 regVol -> AddRootLogicalVolume(box4LV2);
731 box5LV2 -> SetRegion(regVol);
732 regVol -> AddRootLogicalVolume(box5LV2);
733 box6LV2 -> SetRegion(regVol);
734 regVol -> AddRootLogicalVolume(box6LV2);
735 }
736
737 void DetectorConstruction::ConstructJawsY()
738 {
739 G4bool checkOverlap = false;
740
741 G4double boxSide = 101.*mm;
742 G4double box1HalfLengthZ = 15.*mm;
743 G4double box2HalfLengthZ = 31.*mm;
744 G4double box3HalfLengthZ = 35.*mm;
745 G4double box4HalfLengthZ = 21.*mm;
746
747 G4Material* el_Pb = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
748 G4Material* XC10 = GetMaterial("xc10");
749 G4Material* WNICU = GetMaterial("wnicu");
750
751 G4ThreeVector boxJawYSide = {boxSide, boxSide, 219.5*mm}; // 219, perché 219.5??
752
753 G4Box* boxJaw1Y = new G4Box("Jaw1Ybox", boxJawYSide.getX()/2., boxJawYSide.getY()/2., boxJawYSide.getZ()/2.);
754 G4LogicalVolume* boxJaw1YLV = new G4LogicalVolume(boxJaw1Y,
755 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw1YboxLV", 0, 0, 0);
756 G4Box* boxJaw2Y = new G4Box("Jaw2Ybox", boxJawYSide.getX()/2., boxJawYSide.getY()/2., boxJawYSide.getZ()/2.);
757 G4LogicalVolume* boxJaw2YLV = new G4LogicalVolume(boxJaw2Y,
758 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw2YboxLV", 0, 0, 0);
759
760 jaw1YInitialPos.setX(0.);
761 jaw1YInitialPos.setY(boxJawYSide.getX()/2.);
762 jaw1YInitialPos.setZ(248.5*mm + boxJawYSide.getZ()/2.);
763
764 G4RotationMatrix* cRotationJaw1Y = new G4RotationMatrix();
765 cRotationJaw1Y -> rotateX(-std::fabs(std::atan(jawAperture/2/-(-(sourceToSkinDistance+10*cm)))));
766
767 G4ThreeVector position = *cRotationJaw1Y * jaw1YInitialPos;
768 *cRotationJaw1Y = cRotationJaw1Y -> inverse();
769
770 boxJaw1Y_phys =
771 new G4PVPlacement (cRotationJaw1Y, position,
772 "Jaw1YPV", boxJaw1YLV, accWorld_phys, false, 0, checkOverlap);
773
774 jaw2YInitialPos.setX(jaw1YInitialPos.getX());
775 jaw2YInitialPos.setY(-jaw1YInitialPos.getY());
776 jaw2YInitialPos.setZ(jaw1YInitialPos.getZ());
777
778 G4RotationMatrix* cRotationJaw2Y = new G4RotationMatrix();
779 cRotationJaw2Y -> rotateX(std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
780
781 position = *cRotationJaw2Y * jaw2YInitialPos;
782 *cRotationJaw2Y = cRotationJaw2Y -> inverse();
783
784 boxJaw2Y_phys =
785 new G4PVPlacement (cRotationJaw2Y, position,
786 "Jaw2YPV", boxJaw2YLV, accWorld_phys, false, 0, checkOverlap);
787
788 // Physical volumes
789 G4Box *box1 = new G4Box("Jaws1YBox1", boxSide/2., boxSide/2., box1HalfLengthZ/2.);
790 G4Box *box2 = new G4Box("Jaws1YBox2", boxSide/2., boxSide/2., box2HalfLengthZ/2.);
791 G4Box *box3 = new G4Box("Jaws1YBox3", boxSide/2., boxSide/2., box3HalfLengthZ/2.);
792 G4Box *box4 = new G4Box("Jaws1YBox4", boxSide/2., boxSide/2., box4HalfLengthZ/2.);
793 G4LogicalVolume* box1LV1 = new G4LogicalVolume(box1, XC10, "Jaws1YLV1", 0,0,0);
794 G4LogicalVolume* box2LV1 = new G4LogicalVolume(box2, el_Pb, "Jaws1YLV2", 0,0,0);
795 G4LogicalVolume* box3LV1 = new G4LogicalVolume(box3, WNICU, "Jaws1YLV3", 0,0,0);
796 G4LogicalVolume* box4LV1 = new G4LogicalVolume(box4, el_Pb, "Jaws1YLV4", 0,0,0);
797
798 G4LogicalVolume* box1LV2 = new G4LogicalVolume(box1, XC10, "Jaws2YLV1", 0,0,0);
799 G4LogicalVolume* box2LV2 = new G4LogicalVolume(box2, el_Pb, "Jaws2YLV2", 0,0,0);
800 G4LogicalVolume* box3LV2 = new G4LogicalVolume(box3, WNICU, "Jaws2YLV3", 0,0,0);
801 G4LogicalVolume* box4LV2 = new G4LogicalVolume(box4, el_Pb, "Jaws2YLV4", 0,0,0);
802
803 G4ThreeVector centre = {0., 0., 0.};
804 G4double zCentreCurrentBox = -boxJawYSide.getZ()/2.+box1HalfLengthZ/2.*mm;
805
806 centre.setZ(zCentreCurrentBox);
807 new G4PVPlacement(nullptr, centre, "Jaws1YPV1", box1LV1, boxJaw1Y_phys, false, 0, checkOverlap);
808 new G4PVPlacement(nullptr, centre, "Jaws2YPV1", box1LV2, boxJaw2Y_phys, false, 0, checkOverlap);
809
810 zCentreCurrentBox += (box1HalfLengthZ+box2HalfLengthZ)/2.*mm + 117.5*mm;
811 centre.setZ(zCentreCurrentBox);
812 new G4PVPlacement(nullptr, centre, "Jaws1YPV2", box2LV1, boxJaw1Y_phys, false, 0, checkOverlap);
813 new G4PVPlacement(nullptr, centre, "Jaws2YPV2", box2LV2, boxJaw2Y_phys, false, 0, checkOverlap);
814
815 zCentreCurrentBox += (box2HalfLengthZ+box3HalfLengthZ)/2.*mm;
816 centre.setZ(zCentreCurrentBox);
817 new G4PVPlacement(nullptr, centre, "Jaws1YPV3", box3LV1, boxJaw1Y_phys, false, 0, checkOverlap);
818 new G4PVPlacement(nullptr, centre, "Jaws2YPV3", box3LV2, boxJaw2Y_phys, false, 0, checkOverlap);
819
820 zCentreCurrentBox += (box3HalfLengthZ+box4HalfLengthZ)/2.*mm;
821 centre.setZ(zCentreCurrentBox);
822 new G4PVPlacement(nullptr, centre, "Jaws1YPV4", box4LV1, boxJaw1Y_phys, false, 0, checkOverlap);
823 new G4PVPlacement(nullptr, centre, "Jaws2YPV4", box4LV2, boxJaw2Y_phys, false, 0, checkOverlap);
824
825 // Visualization
826 G4VisAttributes* visAtt = new G4VisAttributes(G4Colour::Yellow());
827 visAtt -> SetVisibility(false);
828 visAtt -> SetForceSolid(false);
829
830 G4VisAttributes* simpleAlSVisAttPb = new G4VisAttributes(G4Colour::Blue());
831 simpleAlSVisAttPb -> SetVisibility(true);
832 simpleAlSVisAttPb -> SetForceSolid(true);
833
834 G4VisAttributes* simpleAlSVisAttXC10 = new G4VisAttributes(G4Colour::Green());
835 simpleAlSVisAttXC10 -> SetVisibility(true);
836 simpleAlSVisAttXC10 -> SetForceSolid(true);
837
838 G4VisAttributes* simpleAlSVisAttWNICU = new G4VisAttributes(G4Colour::Red());
839 simpleAlSVisAttWNICU -> SetVisibility(true);
840 simpleAlSVisAttWNICU -> SetForceSolid(true);
841
842 box1LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
843 box2LV1 -> SetVisAttributes(simpleAlSVisAttPb);
844 box3LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
845 box4LV1 -> SetVisAttributes(simpleAlSVisAttPb);
846
847 box1LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
848 box2LV2 -> SetVisAttributes(simpleAlSVisAttPb);
849 box3LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
850 box4LV2 -> SetVisAttributes(simpleAlSVisAttPb);
851
852 boxJaw1YLV -> SetVisAttributes(visAtt);
853 boxJaw2YLV -> SetVisAttributes(visAtt);
854
855 // Region for cuts
856 G4Region* regVol = new G4Region("JawsYR");
857 G4ProductionCuts* cuts = new G4ProductionCuts;
858 cuts -> SetProductionCut(2.*cm);
859 regVol -> SetProductionCuts(cuts);
860
861 box1LV1 -> SetRegion(regVol);
862 regVol -> AddRootLogicalVolume(box1LV1);
863 box2LV1 -> SetRegion(regVol);
864 regVol -> AddRootLogicalVolume(box2LV1);
865 box3LV1 -> SetRegion(regVol);
866 regVol -> AddRootLogicalVolume(box3LV1);
867 box4LV1 -> SetRegion(regVol);
868 regVol -> AddRootLogicalVolume(box4LV1);
869
870 box1LV2 -> SetRegion(regVol);
871 regVol -> AddRootLogicalVolume(box1LV2);
872 box2LV2 -> SetRegion(regVol);
873 regVol -> AddRootLogicalVolume(box2LV2);
874 box3LV2 -> SetRegion(regVol);
875 regVol -> AddRootLogicalVolume(box3LV2);
876 box4LV2 -> SetRegion(regVol);
877 regVol -> AddRootLogicalVolume(box4LV2);
878 }
879
880 void DetectorConstruction::PhantomSegmentation(G4LogicalVolume* Phantom_log)
881 {
882
883 G4ThreeVector phantomHalfDim = {phantomSideDim/2., phantomSideDim/2., phantomSideDim/2.};
884
885 nSideCells = CheckPhantomSegmentation(phantomSideDim/voxelSideDim);
886 nDepthCells = CheckPhantomSegmentation(phantomSideDim/voxelDepthDim);
887
888 G4Material* water = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
889
890 G4String yRepName("RepY");
891 G4VSolid* solYRep = new G4Box(yRepName, phantomHalfDim.getX(),
892 phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ());
893 G4LogicalVolume* logYRep = new G4LogicalVolume(solYRep, water, yRepName);
894 new G4PVReplica(yRepName, logYRep, Phantom_log, kYAxis, nSideCells,
895 2.*phantomHalfDim.getY()/nSideCells);
896
897 G4String xRepName("RepX");
898 G4VSolid* solXRep = new G4Box(xRepName, phantomHalfDim.getX()/nSideCells,
899 phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ());
900 G4LogicalVolume* logXRep = new G4LogicalVolume(solXRep, water, xRepName);
901 new G4PVReplica(xRepName, logXRep, logYRep, kXAxis, nSideCells,
902 2.*phantomHalfDim.getX()/nSideCells);
903
904 G4String zVoxName("phantomSens");
905 G4VSolid* solVoxel = new G4Box(zVoxName, phantomHalfDim.getX()/nSideCells,
906 phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ()/nDepthCells);
907
908 LVPhantomSens = new G4LogicalVolume(solVoxel, water, zVoxName); // This is the Sensitive Volume
909 new G4PVReplica(zVoxName, LVPhantomSens, logXRep, kZAxis, nDepthCells,
910 2.*phantomHalfDim.getZ()/nDepthCells);
911
912 }
913
914 G4int DetectorConstruction::CheckPhantomSegmentation(G4int nCells)
915 {
916 // I want an odd number of voxels
917 // check whether nCells is integer
918 if ( nCells != (int)nCells )
919 {
920 nCells = std::ceil(nCells);
921 }
922 // check whether nCells is even
923 if ( std::fmod(nCells,2) == 0)
924 {
925 nCells += 1;
926 }
927
928 return nCells;
929 }
930
931 G4double DetectorConstruction::GetAccOriginPosition()
932 {
933 return sourceToSkinDistance;
934 }
935
936 G4double DetectorConstruction::GetFFilterRadius()
937 {
938 return tubeFFRadius;
939 }
940
941 G4double DetectorConstruction::GetFFilterZ()
942 {
943 return tubeFFFirstFaceZ;
944 }
945
946 G4double DetectorConstruction::GetVoxelDepthDim()
947 {
948 return voxelDepthDim;
949 }
950
951 G4int DetectorConstruction::GetNumberSideCells() const
952 {
953 return nSideCells;
954 }
955
956 G4int DetectorConstruction::GetNumberDepthCells() const
957 {
958 return nDepthCells;
959 }
960
961 G4int DetectorConstruction::GetPhantomDepth() const
962 {
963 return phantomSideDim;
964 }
965
966 void DetectorConstruction::SetJaws(G4double value)
967 {
968 jawAperture = value;
969 }
970
971 void DetectorConstruction::SetTargetPosition(G4double value)
972 {
973 sourceToSkinDistance = value;
974 }
975
976 void DetectorConstruction::SetPhantomSide(G4double value)
977 {
978 phantomSideDim = value;
979 }
980
981 void DetectorConstruction::SetVoxelSide(G4double value)
982 {
983 voxelSideDim = value;
984 }
985
986 void DetectorConstruction::SetVoxelDepth(G4double value)
987 {
988 voxelDepthDim = value;
989 }
990
991 void DetectorConstruction::UpdateGeometry(G4String string, G4double value)
992 {
993
994 if ( string == "fieldSide" )
995 {
996 //
997 G4double halfFieldSide = value/2;
998 G4ThreeVector position;
999 G4RotationMatrix* rMatrix1X = new G4RotationMatrix();
1000
1001 rMatrix1X -> rotateY(std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1002 position = *rMatrix1X * jaw1XInitialPos;
1003 boxJaw1X_phys -> SetTranslation(position);
1004 *rMatrix1X = rMatrix1X -> inverse();
1005 boxJaw1X_phys -> SetRotation(rMatrix1X);
1006
1007 G4RotationMatrix* rMatrix2X = new G4RotationMatrix();
1008 rMatrix2X -> rotateY(-std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1009
1010 position = *rMatrix2X * jaw2XInitialPos;
1011 boxJaw2X_phys -> SetTranslation(position);
1012 *rMatrix2X = rMatrix2X -> inverse();
1013 boxJaw2X_phys -> SetRotation(rMatrix2X);
1014
1015 G4RotationMatrix* rMatrix1Y = new G4RotationMatrix();
1016 rMatrix1Y -> rotateX(-std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1017 position = *rMatrix1Y * jaw1YInitialPos;
1018 boxJaw1Y_phys -> SetTranslation(position);
1019 *rMatrix1Y = rMatrix1Y -> inverse();
1020 boxJaw1Y_phys -> SetRotation(rMatrix1Y);
1021
1022 G4RotationMatrix* rMatrix2Y = new G4RotationMatrix();
1023 rMatrix2Y -> rotateX(std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1024
1025 position = *rMatrix2Y * jaw2YInitialPos;
1026 boxJaw2Y_phys -> SetTranslation(position);
1027 *rMatrix2Y = rMatrix2Y -> inverse();
1028 boxJaw2Y_phys -> SetRotation(rMatrix2Y);
1029 }
1030 else if ( string == "sourceToSkinDistance")
1031 {
1032 accWorld_phys -> SetTranslation({0., 0., value});
1033 }
1034 else if (string == "phantomSide")
1035 {
1036 G4Box* box = (G4Box *) phantom_phys -> GetLogicalVolume() -> GetSolid();
1037
1038 // change phantom side
1039 box -> SetXHalfLength(phantomSideDim/2.);
1040 box -> SetYHalfLength(phantomSideDim/2.);
1041 box -> SetZHalfLength(phantomSideDim/2.);
1042 phantom_phys -> SetTranslation({0., 0., phantomSideDim/2.});
1043
1044 // clear daughters and rebuild them
1045 phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1046 PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1047 }
1048 else if (string == "voxelSide")
1049 {
1050 phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1051 PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1052 }
1053 else if (string == "voxelDepth")
1054 {
1055 phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1056 PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1057 }
1058 else
1059 G4cerr << "*** DetectorMessenger::UpdateGeometry: Command not found" << G4endl;
1060
1061 G4RunManager::GetRunManager()->GeometryHasBeenModified();
1062 }
1063
1064
1065