source: trunk/examples/advanced/purging_magnet/src/PurgMagDetectorConstruction.cc@ 866

Last change on this file since 866 was 807, checked in by garnier, 17 years ago

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