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