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