| 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 // S.Larsson
28 //
29 // *****************************************
30 // * *
31 // * PurgMagDetectorConstruction.cc *
32 // * *
33 // *****************************************
34 //
35 //
36 #include "PurgMagDetectorConstruction.hh"
37 #include "PurgMagTabulatedField3D.hh"
38 #include "globals.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ThreeVector.hh"
42 #include "G4Material.hh"
43 #include "G4Box.hh"
44 #include "G4Trd.hh"
45 #include "G4Tubs.hh"
46 #include "G4LogicalVolume.hh"
47 #include "G4PVPlacement.hh"
48 #include "G4PVReplica.hh"
49 #include "G4PVParameterised.hh"
50 #include "G4Mag_UsualEqRhs.hh"
51 #include "G4FieldManager.hh"
52 #include "G4TransportationManager.hh"
53 #include "G4EqMagElectricField.hh"
54
55 #include "G4ChordFinder.hh"
56 #include "G4UniformMagField.hh"
57 #include "G4ExplicitEuler.hh"
58 #include "G4ImplicitEuler.hh"
59 #include "G4SimpleRunge.hh"
60 #include "G4SimpleHeum.hh"
61 #include "G4ClassicalRK4.hh"
62 #include "G4HelixExplicitEuler.hh"
63 #include "G4HelixImplicitEuler.hh"
64 #include "G4HelixSimpleRunge.hh"
65 #include "G4CashKarpRKF45.hh"
66 #include "G4RKG3_Stepper.hh"
67
68 #include "G4VisAttributes.hh"
69 #include "G4Colour.hh"
70 #include "G4UnitsTable.hh"
71 #include "G4ios.hh"
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 // Possibility to turn off (0) magnetic field and measurement volume.
75 #define GAP 1 // Magnet geometric volume
76 #define MAG 1 // Magnetic field grid
77 #define MEASUREVOL 1 // Volume for measurement
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81 PurgMagDetectorConstruction::PurgMagDetectorConstruction()
82
83 :physiWorld(NULL), logicWorld(NULL), solidWorld(NULL),
84 physiGap1(NULL), logicGap1(NULL), solidGap1(NULL),
85 physiGap2(NULL), logicGap2(NULL), solidGap2(NULL),
86 physiMeasureVolume(NULL), logicMeasureVolume(NULL),
87 solidMeasureVolume(NULL),
88 WorldMaterial(NULL),
89 GapMaterial(NULL)
90
91 {
92 fField.Put(0);
93 WorldSizeXY=WorldSizeZ=0;
94 GapSizeX1=GapSizeX2=GapSizeY1=GapSizeY2=GapSizeZ=0;
95 MeasureVolumeSizeXY=MeasureVolumeSizeZ=0;
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99
100 PurgMagDetectorConstruction::~PurgMagDetectorConstruction()
101 {}
102
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105 G4VPhysicalVolume* PurgMagDetectorConstruction::Construct()
106
107 {
108 DefineMaterials();
109 return ConstructCalorimeter();
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113
114 void PurgMagDetectorConstruction::DefineMaterials()
115 {
116 //This function illustrates the possible ways to define materials.
117 //Density and mass per mole taken from Physics Handbook for Science
118 //and engineering, sixth edition. This is a general material list
119 //with extra materials for other examples.
120
121 G4String name, symbol;
122 G4double density;
123
124 G4int ncomponents, natoms;
125 G4double fractionmass;
126 G4double temperature, pressure;
127
128 // Define Elements
129 // Example: G4Element* Notation = new G4Element ("Element", "Notation", z, a);
130 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
131 G4Element* N = new G4Element ("Nitrogen", "N", 7., 14.01*g/mole);
132 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
133 G4Element* Ar = new G4Element ("Argon" , "Ar", 18., 39.948*g/mole );
134
135
136 // Define Material
137 // Example: G4Material* Notation = new G4Material("Material", z, a, density);
138 /* Not used in this setup, will be used in further development.
139 G4Material* He = new G4Material("Helium", 2., 4.00*g/mole, 0.178*mg/cm3);
140 G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3);
141 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
142 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
143 */
144 G4Material* Fe = new G4Material("Iron", 26., 55.84*g/mole, 7.87*g/cm3);
145
146 // Define materials from elements.
147
148 // Case 1: chemical molecule
149 // Water
150 density = 1.000*g/cm3;
151 G4Material* H2O = new G4Material(name="H2O" , density, ncomponents=2);
152 H2O->AddElement(H, natoms=2);
153 H2O->AddElement(O, natoms=1);
154
155 // Case 2: mixture by fractional mass.
156 // Air
157 density = 1.290*mg/cm3;
158 G4Material* Air = new G4Material(name="Air" , density, ncomponents=2);
159 Air->AddElement(N, fractionmass=0.7);
160 Air->AddElement(O, fractionmass=0.3);
161
162 // Vacuum
163 density = 1.e-5*g/cm3;
164 pressure = 2.e-2*bar;
165 temperature = STP_Temperature; //from PhysicalConstants.h
166 G4Material* vacuum = new G4Material(name="vacuum", density, ncomponents=1,
167 kStateGas,temperature,pressure);
168 vacuum->AddMaterial(Air, fractionmass=1.);
169
170
171 // Laboratory vacuum: Dry air (average composition)
172 density = 1.7836*mg/cm3 ; // STP
173 G4Material* Argon = new G4Material(name="Argon", density, ncomponents=1);
174 Argon->AddElement(Ar, 1);
175
176 density = 1.25053*mg/cm3 ; // STP
177 G4Material* Nitrogen = new G4Material(name="N2", density, ncomponents=1);
178 Nitrogen->AddElement(N, 2);
179
180 density = 1.4289*mg/cm3 ; // STP
181 G4Material* Oxygen = new G4Material(name="O2", density, ncomponents=1);
182 Oxygen->AddElement(O, 2);
183
184
185 density = 1.2928*mg/cm3 ; // STP
186 density *= 1.0e-8 ; // pumped vacuum
187
188 temperature = STP_Temperature;
189 pressure = 1.0e-8*STP_Pressure;
190
191 G4Material* LaboratoryVacuum = new G4Material(name="LaboratoryVacuum",
192 density,ncomponents=3,
193 kStateGas,temperature,pressure);
194 LaboratoryVacuum->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
195 LaboratoryVacuum->AddMaterial( Oxygen, fractionmass = 0.2315 ) ;
196 LaboratoryVacuum->AddMaterial( Argon, fractionmass = 0.0128 ) ;
197
198
199 G4cout << G4endl << *(G4Material::GetMaterialTable()) << G4endl;
200
201
202 // Default materials in setup.
203 WorldMaterial = LaboratoryVacuum;
204 GapMaterial = Fe;
205
206
207 G4cout << "end material"<< G4endl;
208 }
209
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 G4VPhysicalVolume* PurgMagDetectorConstruction::ConstructCalorimeter()
212 {
213 // Complete the parameters definition
214
215 //The World
216 WorldSizeXY = 300.*cm; // Cube
217 WorldSizeZ = 300.*cm;
218
219 //Measurement volume
220 MeasureVolumeSizeXY = 280.*cm; // Cubic slice
221 MeasureVolumeSizeZ = 1.*cm;
222
223 // Position of measurement volume.
224 // SSD is Source to Surface Distance. Source in origo and measurements 50 cm
225 // below in the z-direction (symbolizin a patient at SSD = 50 cm)
226
227 SSD = 50.*cm;
228 MeasureVolumePosition = -(SSD + MeasureVolumeSizeZ/2);
229
230
231 // Geometric definition of the gap of the purging magnet. Approximation of
232 // the shape of the pole gap.
233
234 GapSizeY1 = 10.*cm; // length along x at the surface positioned at -dz
235 GapSizeY2 = 10.*cm; // length along x at the surface positioned at +dz
236 GapSizeX1 = 10.*cm; // length along y at the surface positioned at -dz
237 GapSizeX2 = 18.37*cm; // length along y at the surface positioned at +dz
238 GapSizeZ = 11.5*cm; // length along z axis
239
240 Gap1PosY = 0.*cm;
241 Gap1PosX = -9.55*cm;
242 Gap1PosZ = -6.89*cm;
243
244 Gap2PosY = 0.*cm;
245 Gap2PosX = 9.55*cm;
246 Gap2PosZ = -6.89*cm;
247
248
249 // Coordinate correction for field grif.
250 // Gap opening at z = -11.4 mm.
251 // In field grid coordonates gap at z = -0.007m in field from z = 0.0m to
252 // z = 0.087m.
253 // -> zOffset = -11.4-(-7) = 4.4 mm
254
255 zOffset = 4.4*mm;
256
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258 //
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261 // Some out prints of the setup.
262
263 G4cout << "\n-----------------------------------------------------------"
264 << "\n Geometry and materials"
265 << "\n-----------------------------------------------------------"
266 << "\n ---> World:"
267 << "\n ---> " << WorldMaterial->GetName() << " in World"
268 << "\n ---> " << "WorldSizeXY: " << G4BestUnit(WorldSizeXY,"Length")
269 << "\n ---> " << "WorldSizeZ: " << G4BestUnit(WorldSizeZ,"Length");
270
271 #if GAP
272 G4cout << "\n-----------------------------------------------------------"
273 << "\n ---> Purging Magnet:"
274 << "\n ---> " << "Gap made of "<< GapMaterial->GetName()
275 << "\n ---> " << "GapSizeY1: " << G4BestUnit(GapSizeY1,"Length")
276 << "\n ---> " << "GapSizeY2: " << G4BestUnit(GapSizeY2,"Length")
277 << "\n ---> " << "GapSizeX1: " << G4BestUnit(GapSizeX1,"Length")
278 << "\n ---> " << "GapSizeX2: " << G4BestUnit(GapSizeX2,"Length");
279 #endif
280
281 #if MEASUREVOL
282 G4cout << "\n-----------------------------------------------------------"
283 << "\n ---> Measurement Volume:"
284 << "\n ---> " << WorldMaterial->GetName() << " in Measurement volume"
285 << "\n ---> " << "MeasureVolumeXY: " << G4BestUnit(MeasureVolumeSizeXY,"Length")
286 << "\n ---> " << "MeasureVolumeZ: " << G4BestUnit(MeasureVolumeSizeZ,"Length")
287 << "\n ---> " << "At SSD = " << G4BestUnit(MeasureVolumePosition,"Length");
288 #endif
289
290 G4cout << "\n-----------------------------------------------------------\n";
291
292
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 //
295 // World
296 //
297
298
299 solidWorld = new G4Box("World", //its name
300 WorldSizeXY/2,WorldSizeXY/2,WorldSizeZ/2); //its size
301
302
303 logicWorld = new G4LogicalVolume(solidWorld, //its solid
304 WorldMaterial, //its material
305 "World"); //its name
306
307 physiWorld = new G4PVPlacement(0, //no rotation
308 G4ThreeVector(), //at (0,0,0)
309 "World", //its name
310 logicWorld, //its logical volume
311 NULL, //its mother volume
312 false, //no boolean operation
313 0); //copy number
314
315 // Visualization attributes
316 G4VisAttributes* simpleWorldVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0)); //White
317 simpleWorldVisAtt->SetVisibility(true);
318 logicWorld->SetVisAttributes(simpleWorldVisAtt);
319
320
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322 //
323 // Measurement Volume
324 //
325
326 #if MEASUREVOL
327
328 solidMeasureVolume = new G4Box("MeasureVolume", //its name
329 MeasureVolumeSizeXY/2,MeasureVolumeSizeXY/2,MeasureVolumeSizeZ/2); //its size
330
331 logicMeasureVolume = new G4LogicalVolume(solidMeasureVolume, //its solid
332 WorldMaterial, //its material
333 "MeasureVolume"); //its name
334
335 physiMeasureVolume = new G4PVPlacement(0, //no rotation
336 G4ThreeVector(0.,0.,MeasureVolumePosition), //at (0,0,0)
337 "MeasureVolume", //its name
338 logicMeasureVolume, //its logical volume
339 physiWorld, //its mother volume
340 false, //no boolean operation
341 0); //copy number
342
343 // Visualization attributes
344 G4VisAttributes* simpleMeasureVolumeVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0)); //White
345 simpleMeasureVolumeVisAtt->SetVisibility(true);
346 simpleMeasureVolumeVisAtt->SetForceSolid(true);
347 logicMeasureVolume->SetVisAttributes(simpleMeasureVolumeVisAtt);
348
349 #endif
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 //
353 //Gap cone. Opening 20 deg. Two separate trapezoids. Iron.
354 //
355
356 #if GAP
357
358 //Gap part 1, placed in negative x-direction.
359
360 solidGap1 = new G4Trd("Gap1",
361 GapSizeX1/2, // Half-length along x at the surface positioned at -dz
362 GapSizeX2/2, // Half-length along x at the surface positioned at +dz
363 GapSizeY1/2, // Half-length along y at the surface positioned at -dz
364 GapSizeY2/2, // Half-length along y at the surface positioned at +dz
365 GapSizeZ/2 ); // Half-length along z axis
366
367 logicGap1 = new G4LogicalVolume(solidGap1, //its solid
368 GapMaterial, //its material
369 "Gap1"); //its name
370
371 physiGap1 = new G4PVPlacement(0, //90 deg rotation
372 G4ThreeVector(Gap1PosX,Gap1PosY,Gap1PosZ), //position
373 "Gap1", //its name
374 logicGap1, //its logical volume
375 physiWorld, //its mother volume
376 false, //no boolean operation
377 0); //copy number
378
379 //Gap part 2, placed in positive x-direction.
380
381 solidGap2 = new G4Trd("Gap2",
382 GapSizeX1/2, // Half-length along x at the surface positioned at -dz
383 GapSizeX2/2, // Half-length along x at the surface positioned at +dz
384 GapSizeY1/2, // Half-length along y at the surface positioned at -dz
385 GapSizeY2/2, // Half-length along y at the surface positioned at +dz
386 GapSizeZ/2 ); // Half-length along z axis
387
388 logicGap2 = new G4LogicalVolume(solidGap2, //its solid
389 GapMaterial, //its material
390 "Gap2"); //its name
391
392 physiGap2 = new G4PVPlacement(0, //no rotation
393 G4ThreeVector(Gap2PosX,Gap2PosY,Gap2PosZ), //position
394 "Gap2", //its name
395 logicGap2, //its logical volume
396 physiWorld, //its mother volume
397 false, //no boolean operation
398 0); //copy number
399
400 // Visualization attributes
401 G4VisAttributes* simpleGap1VisAtt= new G4VisAttributes(G4Colour(0.0,0.0,1.0)); //yellow
402 simpleGap1VisAtt->SetVisibility(true);
403 simpleGap1VisAtt->SetForceSolid(true);
404 logicGap1->SetVisAttributes(simpleGap1VisAtt);
405
406 G4VisAttributes* simpleGap2VisAtt= new G4VisAttributes(G4Colour(0.0,0.0,1.0)); //yellow
407 simpleGap2VisAtt->SetVisibility(true);
408 simpleGap2VisAtt->SetForceSolid(true);
409 logicGap2->SetVisAttributes(simpleGap2VisAtt);
410
411 #endif
412
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414
415 return physiWorld;
416 }
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
418
419 void PurgMagDetectorConstruction::ConstructSDandField()
420 {
421 // Magnetic Field - Purging magnet
422 //
423 #if MAG
424
425 if (fField.Get() == 0)
426 {
427 //Field grid in A9.TABLE. File must be in accessible from run urn directory.
428 G4MagneticField* PurgMagField= new PurgMagTabulatedField3D("PurgMag3D.TABLE", zOffset);
429 fField.Put(PurgMagField);
430
431 //This is thread-local
432 G4FieldManager* pFieldMgr =
433 G4TransportationManager::GetTransportationManager()->GetFieldManager();
434
435 G4cout<< "DeltaStep "<<pFieldMgr->GetDeltaOneStep()/mm <<"mm" <<G4endl;
436 //G4ChordFinder *pChordFinder = new G4ChordFinder(PurgMagField);
437
438 pFieldMgr->SetDetectorField(fField.Get());
439 pFieldMgr->CreateChordFinder(fField.Get());
440
441 }
442 #endif
443 }
444