ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PurgMagDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PurgMagDetectorConstruction.cc
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 }