purging_magnet/src/PurgMagDetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Code developed by:
0027 //  S.Larsson
0028 //
0029 //    *****************************************
0030 //    *                                       *
0031 //    *    PurgMagDetectorConstruction.cc     *
0032 //    *                                       *
0033 //    *****************************************
0034 //
0035 //
0036 #include "PurgMagDetectorConstruction.hh"
0037 #include "PurgMagTabulatedField3D.hh"
0038 #include "globals.hh"
0039 #include "G4PhysicalConstants.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4ThreeVector.hh"
0042 #include "G4Material.hh"
0043 #include "G4Box.hh"
0044 #include "G4Trd.hh"
0045 #include "G4Tubs.hh"
0046 #include "G4LogicalVolume.hh"
0047 #include "G4PVPlacement.hh"
0048 #include "G4PVReplica.hh"
0049 #include "G4PVParameterised.hh"
0050 #include "G4Mag_UsualEqRhs.hh"
0051 #include "G4FieldManager.hh"
0052 #include "G4TransportationManager.hh"
0053 #include "G4EqMagElectricField.hh"
0054
0055 #include "G4ChordFinder.hh"
0056 #include "G4UniformMagField.hh"
0057 #include "G4ExplicitEuler.hh"
0058 #include "G4ImplicitEuler.hh"
0059 #include "G4SimpleRunge.hh"
0060 #include "G4SimpleHeum.hh"
0061 #include "G4ClassicalRK4.hh"
0062 #include "G4HelixExplicitEuler.hh"
0063 #include "G4HelixImplicitEuler.hh"
0064 #include "G4HelixSimpleRunge.hh"
0065 #include "G4CashKarpRKF45.hh"
0066 #include "G4RKG3_Stepper.hh"
0067
0068 #include "G4VisAttributes.hh"
0069 #include "G4Colour.hh"
0070 #include "G4UnitsTable.hh"
0071 #include "G4ios.hh"
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0074 // Possibility to turn off (0) magnetic field and measurement volume.
0075 #define GAP 1          // Magnet geometric volume
0076 #define MAG 1          // Magnetic field grid
0077 #define MEASUREVOL 1   // Volume for measurement
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0080
0081 PurgMagDetectorConstruction::PurgMagDetectorConstruction()
0082
0083   :physiWorld(NULL), logicWorld(NULL), solidWorld(NULL),
0084    physiGap1(NULL), logicGap1(NULL), solidGap1(NULL),
0085    physiGap2(NULL), logicGap2(NULL), solidGap2(NULL),
0086    physiMeasureVolume(NULL), logicMeasureVolume(NULL),
0087    solidMeasureVolume(NULL),
0088    WorldMaterial(NULL),
0089    GapMaterial(NULL)
0090
0091 {
0092   fField.Put(0);
0093   WorldSizeXY=WorldSizeZ=0;
0094   GapSizeX1=GapSizeX2=GapSizeY1=GapSizeY2=GapSizeZ=0;
0095   MeasureVolumeSizeXY=MeasureVolumeSizeZ=0;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0099
0100 PurgMagDetectorConstruction::~PurgMagDetectorConstruction()
0101 {}
0102
0103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0104
0105 G4VPhysicalVolume* PurgMagDetectorConstruction::Construct()
0106
0107 {
0108   DefineMaterials();
0109   return ConstructCalorimeter();
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0113
0114 void PurgMagDetectorConstruction::DefineMaterials()
0115 {
0116   //This function illustrates the possible ways to define materials.
0117   //Density and mass per mole taken from Physics Handbook for Science
0118   //and engineering, sixth edition. This is a general material list
0119   //with extra materials for other examples.
0120
0121   G4String name, symbol;
0122   G4double density;
0123
0124   G4int ncomponents, natoms;
0125   G4double fractionmass;
0126   G4double temperature, pressure;
0127
0128   // Define Elements
0129   // Example: G4Element* Notation  = new G4Element ("Element", "Notation", z, a);
0130   G4Element*   H  = new G4Element ("Hydrogen", "H", 1. ,  1.01*g/mole);
0131   G4Element*   N  = new G4Element ("Nitrogen", "N", 7., 14.01*g/mole);
0132   G4Element*   O  = new G4Element ("Oxygen"  , "O", 8. , 16.00*g/mole);
0133   G4Element*   Ar = new G4Element ("Argon" , "Ar", 18., 39.948*g/mole );
0134
0135
0136   // Define Material
0137   // Example: G4Material* Notation = new G4Material("Material", z, a, density);
0138   /* Not used in this setup, will be used in further development.
0139   G4Material* He = new G4Material("Helium", 2., 4.00*g/mole, 0.178*mg/cm3);
0140   G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3);
0141   G4Material* W  = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
0142   G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
0143   */
0144   G4Material* Fe = new G4Material("Iron", 26., 55.84*g/mole, 7.87*g/cm3);
0145
0146   // Define materials from elements.
0147
0148   // Case 1: chemical molecule
0149   // Water
0150   density = 1.000*g/cm3;
0151   G4Material* H2O = new G4Material(name="H2O"  , density, ncomponents=2);
0152   H2O->AddElement(H, natoms=2);
0153   H2O->AddElement(O, natoms=1);
0154
0155   // Case 2: mixture by fractional mass.
0156   // Air
0157   density = 1.290*mg/cm3;
0158   G4Material* Air = new G4Material(name="Air"  , density, ncomponents=2);
0159   Air->AddElement(N, fractionmass=0.7);
0160   Air->AddElement(O, fractionmass=0.3);
0161
0162   // Vacuum
0163   density     = 1.e-5*g/cm3;
0164   pressure    = 2.e-2*bar;
0165   temperature = STP_Temperature;         //from PhysicalConstants.h
0166   G4Material* vacuum = new G4Material(name="vacuum", density, ncomponents=1,
0167                                       kStateGas,temperature,pressure);
0168   vacuum->AddMaterial(Air, fractionmass=1.);
0169
0170
0171   // Laboratory vacuum: Dry air (average composition)
0172   density = 1.7836*mg/cm3 ;       // STP
0173   G4Material* Argon = new G4Material(name="Argon", density, ncomponents=1);
0174   Argon->AddElement(Ar, 1);
0175
0176   density = 1.25053*mg/cm3 ;       // STP
0177   G4Material* Nitrogen = new G4Material(name="N2", density, ncomponents=1);
0178   Nitrogen->AddElement(N, 2);
0179
0180   density = 1.4289*mg/cm3 ;       // STP
0181   G4Material* Oxygen = new G4Material(name="O2", density, ncomponents=1);
0182   Oxygen->AddElement(O, 2);
0183
0184
0185   density  = 1.2928*mg/cm3 ;       // STP
0186   density *= 1.0e-8 ;              // pumped vacuum
0187
0188   temperature = STP_Temperature;
0189   pressure = 1.0e-8*STP_Pressure;
0190
0191   G4Material* LaboratoryVacuum = new G4Material(name="LaboratoryVacuum",
0192                         density,ncomponents=3,
0193                         kStateGas,temperature,pressure);
0194   LaboratoryVacuum->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
0195   LaboratoryVacuum->AddMaterial( Oxygen,   fractionmass = 0.2315 ) ;
0196   LaboratoryVacuum->AddMaterial( Argon,    fractionmass = 0.0128 ) ;
0197
0198
0199   G4cout << G4endl << *(G4Material::GetMaterialTable()) << G4endl;
0200
0201
0202   // Default materials in setup.
0203   WorldMaterial = LaboratoryVacuum;
0204   GapMaterial = Fe;
0205
0206
0207   G4cout << "end material"<< G4endl;
0208 }
0209
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0211 G4VPhysicalVolume* PurgMagDetectorConstruction::ConstructCalorimeter()
0212 {
0213   // Complete the parameters definition
0214
0215   //The World
0216   WorldSizeXY  = 300.*cm;  // Cube
0217   WorldSizeZ   = 300.*cm;
0218
0219   //Measurement volume
0220   MeasureVolumeSizeXY = 280.*cm;  // Cubic slice
0221   MeasureVolumeSizeZ  = 1.*cm;
0222
0223   // Position of measurement volume.
0224   // SSD is Source to Surface Distance. Source in origo and measurements 50 cm
0225   // below in the z-direction (symbolizin a patient at SSD = 50 cm)
0226
0227   SSD = 50.*cm;
0228   MeasureVolumePosition = -(SSD + MeasureVolumeSizeZ/2);
0229
0230
0231   // Geometric definition of the gap of the purging magnet. Approximation of
0232   // the shape of the pole gap.
0233
0234   GapSizeY1 = 10.*cm;    // length along x at the surface positioned at -dz
0235   GapSizeY2 = 10.*cm;    // length along x at the surface positioned at +dz
0236   GapSizeX1 = 10.*cm;    // length along y at the surface positioned at -dz
0237   GapSizeX2 = 18.37*cm;  // length along y at the surface positioned at +dz
0238   GapSizeZ  = 11.5*cm;   // length along z axis
0239
0240   Gap1PosY = 0.*cm;
0241   Gap1PosX = -9.55*cm;
0242   Gap1PosZ = -6.89*cm;
0243
0244   Gap2PosY = 0.*cm;
0245   Gap2PosX = 9.55*cm;
0246   Gap2PosZ = -6.89*cm;
0247
0248
0249   // Coordinate correction for field grif.
0250   // Gap opening at z = -11.4 mm.
0251   // In field grid coordonates gap at z = -0.007m in field from z = 0.0m to
0252   // z = 0.087m.
0253   // -> zOffset = -11.4-(-7) = 4.4 mm
0254
0255   zOffset = 4.4*mm;
0256
0257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0258 //
0259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0260
0261 // Some out prints of the setup.
0262
0263   G4cout << "\n-----------------------------------------------------------"
0264      << "\n      Geometry and materials"
0265      << "\n-----------------------------------------------------------"
0266      << "\n ---> World:"
0267      << "\n ---> " << WorldMaterial->GetName() << " in World"
0268      << "\n ---> " << "WorldSizeXY: " << G4BestUnit(WorldSizeXY,"Length")
0269      << "\n ---> " << "WorldSizeZ: " << G4BestUnit(WorldSizeZ,"Length");
0270
0271 #if GAP
0272   G4cout << "\n-----------------------------------------------------------"
0273      << "\n ---> Purging Magnet:"
0274      << "\n ---> " << "Gap made of "<< GapMaterial->GetName()
0275      << "\n ---> " << "GapSizeY1: " << G4BestUnit(GapSizeY1,"Length")
0276      << "\n ---> " << "GapSizeY2: " << G4BestUnit(GapSizeY2,"Length")
0277      << "\n ---> " << "GapSizeX1: " << G4BestUnit(GapSizeX1,"Length")
0278      << "\n ---> " << "GapSizeX2: " << G4BestUnit(GapSizeX2,"Length");
0279 #endif
0280
0281 #if MEASUREVOL
0282   G4cout << "\n-----------------------------------------------------------"
0283      << "\n ---> Measurement Volume:"
0284      << "\n ---> " << WorldMaterial->GetName() << " in Measurement volume"
0285      << "\n ---> " << "MeasureVolumeXY: " << G4BestUnit(MeasureVolumeSizeXY,"Length")
0286      << "\n ---> " << "MeasureVolumeZ: " << G4BestUnit(MeasureVolumeSizeZ,"Length")
0287      << "\n ---> " << "At SSD =  " << G4BestUnit(MeasureVolumePosition,"Length");
0288 #endif
0289
0290   G4cout << "\n-----------------------------------------------------------\n";
0291
0292
0293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0294   //
0295   // World
0296   //
0297
0298
0299   solidWorld = new G4Box("World",                      //its name
0300                WorldSizeXY/2,WorldSizeXY/2,WorldSizeZ/2);  //its size
0301
0302
0303   logicWorld = new G4LogicalVolume(solidWorld,          //its solid
0304                    WorldMaterial,   //its material
0305                    "World");        //its name
0306
0307   physiWorld = new G4PVPlacement(0,         //no rotation
0308                  G4ThreeVector(),   //at (0,0,0)
0309                                  "World",       //its name
0310                                  logicWorld,        //its logical volume
0311                                  NULL,          //its mother  volume
0312                                  false,         //no boolean operation
0313                                  0);            //copy number
0314
0315   // Visualization attributes
0316   G4VisAttributes* simpleWorldVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0)); //White
0317   simpleWorldVisAtt->SetVisibility(true);
0318   logicWorld->SetVisAttributes(simpleWorldVisAtt);
0319
0320
0321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0322   //
0323   // Measurement Volume
0324   //
0325
0326 #if MEASUREVOL
0327
0328   solidMeasureVolume = new G4Box("MeasureVolume",                      //its name
0329                    MeasureVolumeSizeXY/2,MeasureVolumeSizeXY/2,MeasureVolumeSizeZ/2);  //its size
0330
0331   logicMeasureVolume = new G4LogicalVolume(solidMeasureVolume,  //its solid
0332                                    WorldMaterial,           //its material
0333                                    "MeasureVolume");        //its name
0334
0335   physiMeasureVolume = new G4PVPlacement(0,                      //no rotation
0336                  G4ThreeVector(0.,0.,MeasureVolumePosition), //at (0,0,0)
0337                                  "MeasureVolume",                    //its name
0338                                  logicMeasureVolume,                     //its logical volume
0339                                  physiWorld,                         //its mother  volume
0340                                  false,                              //no boolean operation
0341                                  0);                                 //copy number
0342
0343   // Visualization attributes
0344   G4VisAttributes* simpleMeasureVolumeVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0)); //White
0345   simpleMeasureVolumeVisAtt->SetVisibility(true);
0346   simpleMeasureVolumeVisAtt->SetForceSolid(true);
0347   logicMeasureVolume->SetVisAttributes(simpleMeasureVolumeVisAtt);
0348
0349 #endif
0350
0351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0352   //
0353   //Gap cone. Opening 20 deg. Two separate trapezoids. Iron.
0354   //
0355
0356 #if GAP
0357
0358   //Gap part 1, placed in negative x-direction.
0359
0360   solidGap1 = new G4Trd("Gap1",
0361             GapSizeX1/2,  // Half-length along x at the surface positioned at -dz
0362             GapSizeX2/2,  // Half-length along x at the surface positioned at +dz
0363             GapSizeY1/2,  // Half-length along y at the surface positioned at -dz
0364             GapSizeY2/2,  // Half-length along y at the surface positioned at +dz
0365             GapSizeZ/2 ); // Half-length along z axis
0366
0367   logicGap1 = new G4LogicalVolume(solidGap1,            //its solid
0368                   GapMaterial,          //its material
0369                   "Gap1");              //its name
0370
0371   physiGap1 = new G4PVPlacement(0,                              //90 deg rotation
0372                 G4ThreeVector(Gap1PosX,Gap1PosY,Gap1PosZ),  //position
0373                 "Gap1",                             //its name
0374                 logicGap1,                          //its logical volume
0375                 physiWorld,                         //its mother  volume
0376                 false,                              //no boolean operation
0377                 0);                             //copy number
0378
0379   //Gap part 2, placed in positive x-direction.
0380
0381   solidGap2 = new G4Trd("Gap2",
0382                     GapSizeX1/2,  // Half-length along x at the surface positioned at -dz
0383                 GapSizeX2/2,  // Half-length along x at the surface positioned at +dz
0384                     GapSizeY1/2,  // Half-length along y at the surface positioned at -dz
0385             GapSizeY2/2,  // Half-length along y at the surface positioned at +dz
0386                     GapSizeZ/2 ); // Half-length along z axis
0387
0388   logicGap2 = new G4LogicalVolume(solidGap2,            //its solid
0389                   GapMaterial,          //its material
0390                   "Gap2");              //its name
0391
0392   physiGap2 = new G4PVPlacement(0,                              //no rotation
0393                 G4ThreeVector(Gap2PosX,Gap2PosY,Gap2PosZ),  //position
0394                 "Gap2",                             //its name
0395                 logicGap2,                          //its logical volume
0396                 physiWorld,                         //its mother  volume
0397                 false,                              //no boolean operation
0398                 0);                             //copy number
0399
0400   // Visualization attributes
0401   G4VisAttributes* simpleGap1VisAtt= new G4VisAttributes(G4Colour(0.0,0.0,1.0)); //yellow
0402   simpleGap1VisAtt->SetVisibility(true);
0403   simpleGap1VisAtt->SetForceSolid(true);
0404   logicGap1->SetVisAttributes(simpleGap1VisAtt);
0405
0406   G4VisAttributes* simpleGap2VisAtt= new G4VisAttributes(G4Colour(0.0,0.0,1.0)); //yellow
0407   simpleGap2VisAtt->SetVisibility(true);
0408   simpleGap2VisAtt->SetForceSolid(true);
0409   logicGap2->SetVisAttributes(simpleGap2VisAtt);
0410
0411 #endif
0412
0413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0414
0415   return physiWorld;
0416 }
0417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0418
0419 void PurgMagDetectorConstruction::ConstructSDandField()
0420 {
0421 //  Magnetic Field - Purging magnet
0422 //
0423 #if MAG
0424
0425   if (fField.Get() == 0)
0426     {
0427       //Field grid in A9.TABLE. File must be in accessible from run urn directory.
0428       G4MagneticField* PurgMagField= new PurgMagTabulatedField3D("PurgMag3D.TABLE", zOffset);
0429       fField.Put(PurgMagField);
0430
0431       //This is thread-local
0432       G4FieldManager* pFieldMgr =
0433     G4TransportationManager::GetTransportationManager()->GetFieldManager();
0434
0435       G4cout<< "DeltaStep "<<pFieldMgr->GetDeltaOneStep()/mm <<"mm" <<G4endl;
0436       //G4ChordFinder *pChordFinder = new G4ChordFinder(PurgMagField);
0437
0438       pFieldMgr->SetDetectorField(fField.Get());
0439       pFieldMgr->CreateChordFinder(fField.Get());
0440
0441     }
0442 #endif
0443 }