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

Last change on this file since 1321 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 16.5 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: geant4-09-03-cand-01 $
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.