Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Please cite the following paper if you use 26 // Please cite the following paper if you use this software 27 // Nucl.Instrum.Meth.B260:20-27, 2007 27 // Nucl.Instrum.Meth.B260:20-27, 2007 28 // 28 // 29 // Based on purging magnet advanced example. 29 // Based on purging magnet advanced example. 30 // 30 // 31 31 32 #include "DetectorConstruction.hh" 32 #include "DetectorConstruction.hh" 33 33 34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh" 35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "G4NistManager.hh" 36 #include "G4NistManager.hh" 37 #include "G4RunManager.hh" 37 #include "G4RunManager.hh" 38 38 39 // Field 39 // Field 40 #include "G4Mag_UsualEqRhs.hh" 40 #include "G4Mag_UsualEqRhs.hh" 41 #include "G4TransportationManager.hh" 41 #include "G4TransportationManager.hh" 42 #include "G4ClassicalRK4.hh" 42 #include "G4ClassicalRK4.hh" 43 #include "G4PropagatorInField.hh" 43 #include "G4PropagatorInField.hh" 44 44 45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 46 47 G4ThreadLocal TabulatedField3D* DetectorConstr 47 G4ThreadLocal TabulatedField3D* DetectorConstruction::fField = 0; 48 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 50 51 DetectorConstruction::DetectorConstruction() 51 DetectorConstruction::DetectorConstruction() 52 { 52 { 53 fDetectorMessenger = new DetectorMessenger(th 53 fDetectorMessenger = new DetectorMessenger(this); 54 54 55 // Default values (square field, coef calcula 55 // Default values (square field, coef calculation, profile) 56 56 57 fModel=1; 57 fModel=1; 58 fG1=-11.964623; 58 fG1=-11.964623; 59 fG2=16.494652; 59 fG2=16.494652; 60 fG3=9.866770; 60 fG3=9.866770; 61 fG4=-6.244493; 61 fG4=-6.244493; 62 fCoef=0; 62 fCoef=0; 63 fProfile=1; 63 fProfile=1; 64 fGrid=0; 64 fGrid=0; 65 65 66 } 66 } 67 67 68 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 69 70 DetectorConstruction::~DetectorConstruction() 70 DetectorConstruction::~DetectorConstruction() 71 { delete fDetectorMessenger;} 71 { delete fDetectorMessenger;} 72 72 73 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 74 74 75 G4VPhysicalVolume* DetectorConstruction::Const 75 G4VPhysicalVolume* DetectorConstruction::Construct() 76 76 77 { 77 { 78 DefineMaterials(); 78 DefineMaterials(); 79 return ConstructVolumes(); 79 return ConstructVolumes(); 80 } 80 } 81 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 83 84 void DetectorConstruction::DefineMaterials() 84 void DetectorConstruction::DefineMaterials() 85 { 85 { 86 G4String name, symbol; 86 G4String name, symbol; 87 G4double density; 87 G4double density; 88 88 89 G4double z, a; 89 G4double z, a; 90 90 91 // Vacuum standard definition... 91 // Vacuum standard definition... 92 density = universe_mean_density; 92 density = universe_mean_density; 93 G4Material* vacuum = new G4Material(name="Va 93 G4Material* vacuum = new G4Material(name="Vacuum", z=1., a=1.01*g/mole, 94 density); 94 density); 95 95 96 // NIST 96 // NIST 97 G4NistManager *man=G4NistManager::Instance() 97 G4NistManager *man=G4NistManager::Instance(); 98 man->SetVerbose(1); 98 man->SetVerbose(1); 99 99 100 // 100 // 101 101 102 G4cout << G4endl << *(G4Material::GetMateria 102 G4cout << G4endl << *(G4Material::GetMaterialTable()) << G4endl; 103 103 104 // Default materials in setup. 104 // Default materials in setup. 105 fDefaultMaterial = vacuum; 105 fDefaultMaterial = vacuum; 106 fGridMaterial = man->FindOrBuildMaterial("G4 106 fGridMaterial = man->FindOrBuildMaterial("G4_Ni"); 107 } 107 } 108 108 109 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 110 110 111 G4VPhysicalVolume* DetectorConstruction::Const 111 G4VPhysicalVolume* DetectorConstruction::ConstructVolumes() 112 { 112 { 113 113 114 fSolidWorld = new G4Box("World", //it 114 fSolidWorld = new G4Box("World", //its name 115 12*m/2,12*m/2,22*m/2); //its size 115 12*m/2,12*m/2,22*m/2); //its size 116 116 117 117 118 fLogicWorld = new G4LogicalVolume(fSolidWorl 118 fLogicWorld = new G4LogicalVolume(fSolidWorld, //its solid 119 fDefaultMaterial, //its material 119 fDefaultMaterial, //its material 120 "World"); //its name 120 "World"); //its name 121 121 122 fPhysiWorld = new G4PVPlacement(0, //no 122 fPhysiWorld = new G4PVPlacement(0, //no rotation 123 G4ThreeVector(), //at (0,0,0) 123 G4ThreeVector(), //at (0,0,0) 124 "World", // 124 "World", //its name 125 fLogicWorld, 125 fLogicWorld, //its logical volume 126 NULL, // 126 NULL, //its mother volume 127 false, // 127 false, //no boolean operation 128 0); //co 128 0); //copy number 129 129 130 130 131 // MAGNET VOLUME 131 // MAGNET VOLUME 132 132 133 fSolidVol = new G4Box("Vol", //its na 133 fSolidVol = new G4Box("Vol", //its name 134 10*m/2,10*m/2,9.120*m/2); //its si 134 10*m/2,10*m/2,9.120*m/2); //its size 135 135 136 136 137 fLogicVol = new G4LogicalVolume(fSolidVol, 137 fLogicVol = new G4LogicalVolume(fSolidVol, //its solid 138 fDefaultMaterial, //its material 138 fDefaultMaterial, //its material 139 "Vol"); //its name 139 "Vol"); //its name 140 140 141 fPhysiVol = new G4PVPlacement(0, //no r 141 fPhysiVol = new G4PVPlacement(0, //no rotation 142 G4ThreeVector(0,0,-4310*mm), //at ( 142 G4ThreeVector(0,0,-4310*mm), //at (0,0,0) 143 "Vol", // 143 "Vol", //its name 144 fLogicVol, 144 fLogicVol, //its logical volume 145 fPhysiWorld, 145 fPhysiWorld, //its mother volume 146 false, // 146 false, //no boolean operation 147 0); //co 147 0); //copy number 148 148 149 // GRID 149 // GRID 150 150 151 if (fGrid==1) 151 if (fGrid==1) 152 { 152 { 153 153 154 G4cout << G4endl; 154 G4cout << G4endl; 155 155 156 G4cout << " ********************** " << G4en 156 G4cout << " ********************** " << G4endl; 157 G4cout << " **** GRID IN PLACE *** " << G4en 157 G4cout << " **** GRID IN PLACE *** " << G4endl; 158 G4cout << " ********************** " << G4en 158 G4cout << " ********************** " << G4endl; 159 159 160 G4double x_grid=5.0*mm; 160 G4double x_grid=5.0*mm; 161 G4double y_grid=5.0*mm; 161 G4double y_grid=5.0*mm; 162 G4double grid_Zpos=(250+200)*mm; // 250 162 G4double grid_Zpos=(250+200)*mm; // 250+10 mm for object size of 50µm diam 163 163 164 //G4double thickness_grid=10*micrometer; 164 //G4double thickness_grid=10*micrometer; 165 G4double thickness_grid=100*micrometer; 165 G4double thickness_grid=100*micrometer; 166 166 167 G4double z_grid=thickness_grid/2.0; 167 G4double z_grid=thickness_grid/2.0; 168 168 169 fSolidGridVol= new G4Box("GridVolume",x_grid 169 fSolidGridVol= new G4Box("GridVolume",x_grid,y_grid,z_grid); //its size 170 170 171 fLogicGridVol = new G4LogicalVolume(fSolidGr 171 fLogicGridVol = new G4LogicalVolume(fSolidGridVol, //its solid 172 fGridMaterial, //its 172 fGridMaterial, //its material 173 "GridVolume"); //its name 173 "GridVolume"); //its name 174 174 175 fPhysiGridVol = new G4PVPlacement(0, 175 fPhysiGridVol = new G4PVPlacement(0, //no rotation 176 G4ThreeVector(0,0,grid_Zpos), // o 176 G4ThreeVector(0,0,grid_Zpos), // origin 177 fLogicGridVol 177 fLogicGridVol, //its logical volume 178 "GridVolume", 178 "GridVolume", //its name 179 fLogicWorld, 179 fLogicWorld, //its mother volume 180 false, 180 false, //no boolean operation 181 0); 181 0); 182 182 183 // Holes in grid 183 // Holes in grid 184 184 185 G4double holeSize= 9e-3*mm; 185 G4double holeSize= 9e-3*mm; 186 G4double pix_grid=1.3e-2*mm; 186 G4double pix_grid=1.3e-2*mm; 187 G4int num_half_grid=100; 187 G4int num_half_grid=100; 188 188 189 fSolidGridVol_Hole= new G4Box("GridHole",hol 189 fSolidGridVol_Hole= new G4Box("GridHole",holeSize/2,holeSize/2,z_grid); //its size 190 190 191 fLogicGridVol_Hole = new G4LogicalVolume(fSo 191 fLogicGridVol_Hole = new G4LogicalVolume(fSolidGridVol_Hole, //its solid 192 fDefaultMaterial, 192 fDefaultMaterial, //its material 193 "GridHole"); 193 "GridHole"); //its name 194 194 195 195 196 for(int i=-num_half_grid;i<num_half_grid;i++ 196 for(int i=-num_half_grid;i<num_half_grid;i++) 197 { 197 { 198 for (int j=-num_half_grid;j<num_half_gri 198 for (int j=-num_half_grid;j<num_half_grid;j++) 199 { 199 { 200 200 201 G4double x0_grid,y0_grid,z0_grid; 201 G4double x0_grid,y0_grid,z0_grid; 202 G4int number_index_grid; 202 G4int number_index_grid; 203 203 204 x0_grid=pix_grid*i; 204 x0_grid=pix_grid*i; 205 y0_grid=pix_grid*j; 205 y0_grid=pix_grid*j; 206 z0_grid=0.0*mm; 206 z0_grid=0.0*mm; 207 207 208 number_index_grid=(i+num_half_grid)*1000+( 208 number_index_grid=(i+num_half_grid)*1000+(j+num_half_grid); 209 209 210 fPhysiGridVol_Hole = new G4PVPlacement( 210 fPhysiGridVol_Hole = new G4PVPlacement(0, //no rotation 211 G4ThreeVector(x0_grid,y0_grid,z0_gr 211 G4ThreeVector(x0_grid,y0_grid,z0_grid),//origin 212 fLogicGridVol 212 fLogicGridVol_Hole, //its logical volume 213 "GridHole", //its name 213 "GridHole", //its name 214 fLogicGridVol 214 fLogicGridVol, //its mother volume 215 false, 215 false, //no boolean operation 216 number_index_ 216 number_index_grid); 217 } 217 } 218 } 218 } 219 219 220 // Grid imaging plane 220 // Grid imaging plane 221 221 222 G4double ContVolSizeXY = 1*m; 222 G4double ContVolSizeXY = 1*m; 223 G4double ImPlaneWidth = 0.001*mm; 223 G4double ImPlaneWidth = 0.001*mm; 224 224 225 fSolidControlVol_GridShadow = 225 fSolidControlVol_GridShadow = 226 new G4Box 226 new G4Box 227 ("ControlVol_GridShadow", ContVolSizeXY/2, 227 ("ControlVol_GridShadow", ContVolSizeXY/2, ContVolSizeXY/2 , ImPlaneWidth/2); 228 228 229 fLogicControlVol_GridShadow = 229 fLogicControlVol_GridShadow = 230 new G4LogicalVolume 230 new G4LogicalVolume 231 (fSolidControlVol_GridShadow, fDefaultMate 231 (fSolidControlVol_GridShadow, fDefaultMaterial, "ControlVol_GridShadow"); 232 232 233 fPhysiControlVol_GridShadow = 233 fPhysiControlVol_GridShadow = 234 new G4PVPlacement 234 new G4PVPlacement 235 ( 0, G4ThreeVector(0,0,(250+300)*mm), fLog 235 ( 0, G4ThreeVector(0,0,(250+300)*mm), fLogicControlVol_GridShadow, "ControlVol_GridShadow", 236 fLogicWorld, false, 0); 236 fLogicWorld, false, 0); 237 237 238 238 239 } // end GRID 239 } // end GRID 240 240 241 // STEP MINIMUM SIZE 241 // STEP MINIMUM SIZE 242 fLogicVol->SetUserLimits(new G4UserLimits(1* 242 fLogicVol->SetUserLimits(new G4UserLimits(1*mm)); 243 243 244 return fPhysiWorld; 244 return fPhysiWorld; 245 } 245 } 246 246 247 //....oooOO0OOooo........oooOO0OOooo........oo 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 248 248 249 void DetectorConstruction::SetG1(G4float value 249 void DetectorConstruction::SetG1(G4float value) 250 { 250 { 251 fG1 = value; 251 fG1 = value; 252 G4RunManager::GetRunManager()->ReinitializeG 252 G4RunManager::GetRunManager()->ReinitializeGeometry(); 253 } 253 } 254 254 255 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 256 256 257 void DetectorConstruction::SetG2(G4float value 257 void DetectorConstruction::SetG2(G4float value) 258 { 258 { 259 fG2 = value; 259 fG2 = value; 260 G4RunManager::GetRunManager()->ReinitializeG 260 G4RunManager::GetRunManager()->ReinitializeGeometry(); 261 } 261 } 262 262 263 //....oooOO0OOooo........oooOO0OOooo........oo 263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 264 264 265 void DetectorConstruction::SetG3(G4float value 265 void DetectorConstruction::SetG3(G4float value) 266 { 266 { 267 fG3 = value; 267 fG3 = value; 268 G4RunManager::GetRunManager()->ReinitializeG 268 G4RunManager::GetRunManager()->ReinitializeGeometry(); 269 } 269 } 270 270 271 //....oooOO0OOooo........oooOO0OOooo........oo 271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 272 272 273 void DetectorConstruction::SetG4(G4float value 273 void DetectorConstruction::SetG4(G4float value) 274 { 274 { 275 fG4 = value; 275 fG4 = value; 276 G4RunManager::GetRunManager()->ReinitializeG 276 G4RunManager::GetRunManager()->ReinitializeGeometry(); 277 } 277 } 278 278 279 //....oooOO0OOooo........oooOO0OOooo........oo 279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 280 280 281 void DetectorConstruction::SetModel(G4int mode 281 void DetectorConstruction::SetModel(G4int modelChoice) 282 { 282 { 283 if (modelChoice==1) fModel=1; 283 if (modelChoice==1) fModel=1; 284 if (modelChoice==2) fModel=2; 284 if (modelChoice==2) fModel=2; 285 if (modelChoice==3) fModel=3; 285 if (modelChoice==3) fModel=3; 286 G4RunManager::GetRunManager()->ReinitializeG 286 G4RunManager::GetRunManager()->ReinitializeGeometry(); 287 } 287 } 288 288 289 //....oooOO0OOooo........oooOO0OOooo........oo 289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 290 290 291 void DetectorConstruction::SetCoef(G4int val) 291 void DetectorConstruction::SetCoef(G4int val) 292 { 292 { 293 fCoef=val; 293 fCoef=val; 294 G4RunManager::GetRunManager()->ReinitializeG 294 G4RunManager::GetRunManager()->ReinitializeGeometry(); 295 } 295 } 296 296 297 //....oooOO0OOooo........oooOO0OOooo........oo 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 298 298 299 G4int DetectorConstruction::GetCoef() 299 G4int DetectorConstruction::GetCoef() 300 { 300 { 301 return fCoef; 301 return fCoef; 302 } 302 } 303 303 304 //....oooOO0OOooo........oooOO0OOooo........oo 304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 305 305 306 void DetectorConstruction::SetProfile(G4int my 306 void DetectorConstruction::SetProfile(G4int myProfile) 307 { 307 { 308 fProfile=myProfile; 308 fProfile=myProfile; 309 G4RunManager::GetRunManager()->ReinitializeG 309 G4RunManager::GetRunManager()->ReinitializeGeometry(); 310 } 310 } 311 311 312 //....oooOO0OOooo........oooOO0OOooo........oo 312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 313 313 314 void DetectorConstruction::SetGrid(G4int myGri 314 void DetectorConstruction::SetGrid(G4int myGrid) 315 { 315 { 316 fGrid=myGrid; 316 fGrid=myGrid; 317 G4RunManager::GetRunManager()->ReinitializeG 317 G4RunManager::GetRunManager()->ReinitializeGeometry(); 318 } 318 } 319 319 320 //....oooOO0OOooo........oooOO0OOooo........oo 320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 321 321 322 void DetectorConstruction::ConstructSDandField 322 void DetectorConstruction::ConstructSDandField() 323 { 323 { 324 fField = new TabulatedField3D(fG1, fG2, 324 fField = new TabulatedField3D(fG1, fG2, fG3, fG4, fModel); 325 325 326 //This is thread-local 326 //This is thread-local 327 G4FieldManager* fFieldMgr = 327 G4FieldManager* fFieldMgr = 328 G4TransportationManager::GetTransportationMa 328 G4TransportationManager::GetTransportationManager()->GetFieldManager(); 329 329 330 G4Mag_UsualEqRhs* fEquation = new G4Mag_ 330 G4Mag_UsualEqRhs* fEquation = new G4Mag_UsualEqRhs (fField); 331 331 332 G4ClassicalRK4* fStepper = new G4Classic 332 G4ClassicalRK4* fStepper = new G4ClassicalRK4 (fEquation); 333 333 334 G4ChordFinder* fChordFinder = new G4Chor 334 G4ChordFinder* fChordFinder = new G4ChordFinder(fField,1e-9*m,fStepper); 335 335 336 fFieldMgr->SetChordFinder(fChordFinder); 336 fFieldMgr->SetChordFinder(fChordFinder); 337 fFieldMgr->SetDetectorField(fField); 337 fFieldMgr->SetDetectorField(fField); 338 338 339 // SI: 01-07-2018 : following settings w 339 // SI: 01-07-2018 : following settings were initially set to 1e-9*m 340 // instead of 1e-9*m, but they now indu 340 // instead of 1e-9*m, but they now induce warnings as 341 // *** G4Exception : GeomNav1002 341 // *** G4Exception : GeomNav1002 342 // issued by : G4PropagatorInField::Com 342 // issued by : G4PropagatorInField::ComputeStep 343 343 344 fFieldMgr->GetChordFinder()->SetDeltaCho 344 fFieldMgr->GetChordFinder()->SetDeltaChord(1e-7*m); 345 fFieldMgr->SetDeltaIntersection(1e-7*m); 345 fFieldMgr->SetDeltaIntersection(1e-7*m); 346 fFieldMgr->SetDeltaOneStep(1e-7*m); 346 fFieldMgr->SetDeltaOneStep(1e-7*m); 347 347 348 // 348 // 349 349 350 // To avoid G4MagIntegratorDriver::OneGo 350 // To avoid G4MagIntegratorDriver::OneGoodStep:Stepsize underflows in Stepper 351 351 352 if (fCoef==1) 352 if (fCoef==1) 353 { 353 { 354 G4PropagatorInField* fPropInField = 354 G4PropagatorInField* fPropInField = 355 G4TransportationManager::GetTranspor 355 G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); 356 fPropInField->SetMinimumEpsilonStep(1e 356 fPropInField->SetMinimumEpsilonStep(1e-11); 357 fPropInField->SetMaximumEpsilonStep(1e 357 fPropInField->SetMaximumEpsilonStep(1e-10); 358 358 359 } 359 } 360 else 360 else 361 { 361 { 362 G4PropagatorInField* fPropInField = 362 G4PropagatorInField* fPropInField = 363 G4TransportationManager::GetTranspor 363 G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); 364 fPropInField->SetMinimumEpsilonStep(1e 364 fPropInField->SetMinimumEpsilonStep(1e-9); 365 fPropInField->SetMaximumEpsilonStep(1e 365 fPropInField->SetMaximumEpsilonStep(1e-8); 366 } 366 } 367 367 368 } 368 } 369 369 370 370