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 // ------------------------------------------------------------------- 27 // Nucl.Instrum.Meth.B260:20-27, 2007 << 27 // $Id$ 28 // << 28 // ------------------------------------------------------------------- 29 // Based on purging magnet advanced example. << 30 // << 31 29 32 #include "DetectorConstruction.hh" 30 #include "DetectorConstruction.hh" 33 << 34 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh" 35 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" 36 #include "G4NistManager.hh" 33 #include "G4NistManager.hh" 37 #include "G4RunManager.hh" << 38 << 39 // Field << 40 #include "G4Mag_UsualEqRhs.hh" << 41 #include "G4TransportationManager.hh" << 42 #include "G4ClassicalRK4.hh" << 43 #include "G4PropagatorInField.hh" << 44 << 45 //....oooOO0OOooo........oooOO0OOooo........oo << 46 << 47 G4ThreadLocal TabulatedField3D* DetectorConstr << 48 34 49 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 36 51 DetectorConstruction::DetectorConstruction() 37 DetectorConstruction::DetectorConstruction() 52 { 38 { 53 fDetectorMessenger = new DetectorMessenger(th << 39 detectorMessenger = new DetectorMessenger(this); 54 << 40 gradientsInitialized=false; 55 // Default values (square field, coef calcula << 41 G1=0; G2=0; G3=0; G4=0; coef=0; profile=0; grid=0; 56 << 57 fModel=1; << 58 fG1=-11.964623; << 59 fG2=16.494652; << 60 fG3=9.866770; << 61 fG4=-6.244493; << 62 fCoef=0; << 63 fProfile=1; << 64 fGrid=0; << 65 << 66 } 42 } 67 43 68 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 45 70 DetectorConstruction::~DetectorConstruction() 46 DetectorConstruction::~DetectorConstruction() 71 { delete fDetectorMessenger;} << 47 { delete detectorMessenger;} 72 48 73 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 74 50 75 G4VPhysicalVolume* DetectorConstruction::Const 51 G4VPhysicalVolume* DetectorConstruction::Construct() 76 52 77 { 53 { 78 DefineMaterials(); 54 DefineMaterials(); 79 return ConstructVolumes(); 55 return ConstructVolumes(); 80 } 56 } 81 57 82 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 59 84 void DetectorConstruction::DefineMaterials() 60 void DetectorConstruction::DefineMaterials() 85 { 61 { 86 G4String name, symbol; 62 G4String name, symbol; 87 G4double density; 63 G4double density; 88 64 89 G4double z, a; 65 G4double z, a; 90 66 91 // Vacuum standard definition... 67 // Vacuum standard definition... 92 density = universe_mean_density; 68 density = universe_mean_density; 93 G4Material* vacuum = new G4Material(name="Va 69 G4Material* vacuum = new G4Material(name="Vacuum", z=1., a=1.01*g/mole, 94 density); 70 density); 95 71 96 // NIST 72 // NIST 97 G4NistManager *man=G4NistManager::Instance() 73 G4NistManager *man=G4NistManager::Instance(); 98 man->SetVerbose(1); 74 man->SetVerbose(1); 99 75 100 // << 101 << 102 G4cout << G4endl << *(G4Material::GetMateria 76 G4cout << G4endl << *(G4Material::GetMaterialTable()) << G4endl; 103 77 104 // Default materials in setup. 78 // Default materials in setup. 105 fDefaultMaterial = vacuum; << 79 defaultMaterial = vacuum; 106 fGridMaterial = man->FindOrBuildMaterial("G4 << 80 gridMaterial = man->FindOrBuildMaterial("G4_Ni"); 107 } 81 } 108 82 109 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 110 84 111 G4VPhysicalVolume* DetectorConstruction::Const 85 G4VPhysicalVolume* DetectorConstruction::ConstructVolumes() 112 { 86 { 113 87 114 fSolidWorld = new G4Box("World", //it << 88 static G4bool fieldIsInitialized = false; >> 89 if(!fieldIsInitialized && gradientsInitialized) >> 90 { >> 91 G4FieldManager* pFieldMgr; >> 92 G4MagIntegratorStepper* pStepper; >> 93 G4Mag_UsualEqRhs* pEquation; >> 94 >> 95 G4MagneticField* Field= new TabulatedField3D(G1, G2, G3, G4, model); >> 96 >> 97 pEquation = new G4Mag_UsualEqRhs (Field); >> 98 pStepper = new G4ClassicalRK4 (pEquation); >> 99 pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); >> 100 >> 101 G4ChordFinder *pChordFinder = new G4ChordFinder(Field,1e-9*m,pStepper); >> 102 pFieldMgr->SetChordFinder( pChordFinder ); >> 103 >> 104 pFieldMgr->SetDetectorField(Field); >> 105 >> 106 fieldIsInitialized = true; >> 107 >> 108 // tuned parameters >> 109 pFieldMgr->GetChordFinder()->SetDeltaChord(1.e-9*m); >> 110 pFieldMgr->SetDeltaIntersection(1.e-9*m); >> 111 pFieldMgr->SetDeltaOneStep(1.e-9*m); >> 112 >> 113 G4PropagatorInField *propInField; >> 114 propInField = >> 115 G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); >> 116 propInField->SetMinimumEpsilonStep(1e-11); >> 117 propInField->SetMaximumEpsilonStep(1e-10); >> 118 >> 119 } >> 120 >> 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 122 >> 123 solidWorld = new G4Box("World", //its name 115 12*m/2,12*m/2,22*m/2); //its size 124 12*m/2,12*m/2,22*m/2); //its size 116 125 117 126 118 fLogicWorld = new G4LogicalVolume(fSolidWorl << 127 logicWorld = new G4LogicalVolume(solidWorld, //its solid 119 fDefaultMaterial, //its material << 128 defaultMaterial, //its material 120 "World"); //its name << 129 "World"); //its name 121 130 122 fPhysiWorld = new G4PVPlacement(0, //no << 131 physiWorld = new G4PVPlacement(0, //no rotation 123 G4ThreeVector(), //at (0,0,0) 132 G4ThreeVector(), //at (0,0,0) 124 "World", // 133 "World", //its name 125 fLogicWorld, << 134 logicWorld, //its logical volume 126 NULL, // 135 NULL, //its mother volume 127 false, // 136 false, //no boolean operation 128 0); //co 137 0); //copy number 129 138 130 139 131 // MAGNET VOLUME 140 // MAGNET VOLUME 132 141 133 fSolidVol = new G4Box("Vol", //its na << 142 solidVol = new G4Box("Vol", //its name 134 10*m/2,10*m/2,9.120*m/2); //its si 143 10*m/2,10*m/2,9.120*m/2); //its size 135 144 136 145 137 fLogicVol = new G4LogicalVolume(fSolidVol, << 146 logicVol = new G4LogicalVolume(solidVol, //its solid 138 fDefaultMaterial, //its material << 147 defaultMaterial, //its material 139 "Vol"); //its name << 148 "Vol"); //its name 140 149 141 fPhysiVol = new G4PVPlacement(0, //no r << 150 physiVol = new G4PVPlacement(0, //no rotation 142 G4ThreeVector(0,0,-4310*mm), //at ( 151 G4ThreeVector(0,0,-4310*mm), //at (0,0,0) 143 "Vol", // 152 "Vol", //its name 144 fLogicVol, << 153 logicVol, //its logical volume 145 fPhysiWorld, << 154 physiWorld, //its mother volume 146 false, // 155 false, //no boolean operation 147 0); //co 156 0); //copy number 148 157 149 // GRID 158 // GRID 150 159 151 if (fGrid==1) << 160 if (grid==1) 152 { 161 { 153 162 154 G4cout << G4endl; 163 G4cout << G4endl; 155 164 156 G4cout << " ********************** " << G4en 165 G4cout << " ********************** " << G4endl; 157 G4cout << " **** GRID IN PLACE *** " << G4en 166 G4cout << " **** GRID IN PLACE *** " << G4endl; 158 G4cout << " ********************** " << G4en 167 G4cout << " ********************** " << G4endl; 159 168 160 G4double x_grid=5.0*mm; 169 G4double x_grid=5.0*mm; 161 G4double y_grid=5.0*mm; 170 G4double y_grid=5.0*mm; 162 G4double grid_Zpos=(250+200)*mm; // 250 << 171 G4double grid_Zpos=(250+200)*mm; // 250+10 mm for object size of 50µm diam 163 172 164 //G4double thickness_grid=10*micrometer; 173 //G4double thickness_grid=10*micrometer; 165 G4double thickness_grid=100*micrometer; 174 G4double thickness_grid=100*micrometer; 166 175 167 G4double z_grid=thickness_grid/2.0; 176 G4double z_grid=thickness_grid/2.0; 168 177 169 fSolidGridVol= new G4Box("GridVolume",x_grid << 178 solidGridVol= new G4Box("GridVolume",x_grid,y_grid,z_grid); //its size 170 179 171 fLogicGridVol = new G4LogicalVolume(fSolidGr << 180 logicGridVol = new G4LogicalVolume(solidGridVol, //its solid 172 fGridMaterial, //its << 181 gridMaterial, //its material 173 "GridVolume"); //its name << 182 "GridVolume"); //its name 174 183 175 fPhysiGridVol = new G4PVPlacement(0, << 184 physiGridVol = new G4PVPlacement(0, //no rotation 176 G4ThreeVector(0,0,grid_Zpos), // o 185 G4ThreeVector(0,0,grid_Zpos), // origin 177 fLogicGridVol << 186 logicGridVol, //its logical volume 178 "GridVolume", 187 "GridVolume", //its name 179 fLogicWorld, << 188 logicWorld, //its mother volume 180 false, 189 false, //no boolean operation 181 0); 190 0); 182 191 183 // Holes in grid 192 // Holes in grid 184 193 185 G4double holeSize= 9e-3*mm; 194 G4double holeSize= 9e-3*mm; 186 G4double pix_grid=1.3e-2*mm; 195 G4double pix_grid=1.3e-2*mm; 187 G4int num_half_grid=100; 196 G4int num_half_grid=100; 188 197 189 fSolidGridVol_Hole= new G4Box("GridHole",hol << 198 solidGridVol_Hole= new G4Box("GridHole",holeSize/2,holeSize/2,z_grid); //its size 190 199 191 fLogicGridVol_Hole = new G4LogicalVolume(fSo << 200 logicGridVol_Hole = new G4LogicalVolume(solidGridVol_Hole, //its solid 192 fDefaultMaterial, << 201 defaultMaterial, //its material 193 "GridHole"); << 202 "GridHole"); //its name 194 203 195 204 196 for(int i=-num_half_grid;i<num_half_grid;i++ 205 for(int i=-num_half_grid;i<num_half_grid;i++) 197 { 206 { 198 for (int j=-num_half_grid;j<num_half_gri 207 for (int j=-num_half_grid;j<num_half_grid;j++) 199 { 208 { 200 209 201 G4double x0_grid,y0_grid,z0_grid; 210 G4double x0_grid,y0_grid,z0_grid; 202 G4int number_index_grid; 211 G4int number_index_grid; 203 212 204 x0_grid=pix_grid*i; 213 x0_grid=pix_grid*i; 205 y0_grid=pix_grid*j; 214 y0_grid=pix_grid*j; 206 z0_grid=0.0*mm; 215 z0_grid=0.0*mm; 207 216 208 number_index_grid=(i+num_half_grid)*1000+( 217 number_index_grid=(i+num_half_grid)*1000+(j+num_half_grid); 209 218 210 fPhysiGridVol_Hole = new G4PVPlacement( << 219 physiGridVol_Hole = new G4PVPlacement(0, //no rotation 211 G4ThreeVector(x0_grid,y0_grid,z0_gr 220 G4ThreeVector(x0_grid,y0_grid,z0_grid),//origin 212 fLogicGridVol << 221 logicGridVol_Hole, //its logical volume 213 "GridHole", //its name 222 "GridHole", //its name 214 fLogicGridVol << 223 logicGridVol, //its mother volume 215 false, 224 false, //no boolean operation 216 number_index_ 225 number_index_grid); 217 } 226 } 218 } 227 } 219 228 220 // Grid imaging plane 229 // Grid imaging plane 221 230 222 G4double ContVolSizeXY = 1*m; 231 G4double ContVolSizeXY = 1*m; 223 G4double ImPlaneWidth = 0.001*mm; 232 G4double ImPlaneWidth = 0.001*mm; 224 233 225 fSolidControlVol_GridShadow = << 234 solidControlVol_GridShadow = 226 new G4Box 235 new G4Box 227 ("ControlVol_GridShadow", ContVolSizeXY/2, 236 ("ControlVol_GridShadow", ContVolSizeXY/2, ContVolSizeXY/2 , ImPlaneWidth/2); 228 237 229 fLogicControlVol_GridShadow = << 238 logicControlVol_GridShadow = 230 new G4LogicalVolume 239 new G4LogicalVolume 231 (fSolidControlVol_GridShadow, fDefaultMate << 240 (solidControlVol_GridShadow, defaultMaterial, "ControlVol_GridShadow"); 232 241 233 fPhysiControlVol_GridShadow = << 242 physiControlVol_GridShadow = 234 new G4PVPlacement 243 new G4PVPlacement 235 ( 0, G4ThreeVector(0,0,(250+300)*mm), fLog << 244 ( 0, G4ThreeVector(0,0,(250+300)*mm), logicControlVol_GridShadow, "ControlVol_GridShadow",logicWorld, false, 0); 236 fLogicWorld, false, 0); << 237 245 238 246 239 } // end GRID 247 } // end GRID 240 248 241 // STEP MINIMUM SIZE 249 // STEP MINIMUM SIZE 242 fLogicVol->SetUserLimits(new G4UserLimits(1* << 250 logicVol->SetUserLimits(new G4UserLimits(1*mm)); 243 251 244 return fPhysiWorld; << 252 return physiWorld; 245 } 253 } 246 254 247 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 248 256 249 void DetectorConstruction::SetG1(G4float value 257 void DetectorConstruction::SetG1(G4float value) 250 { 258 { 251 fG1 = value; << 259 G1 = value; 252 G4RunManager::GetRunManager()->ReinitializeG << 253 } 260 } 254 261 255 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 256 263 257 void DetectorConstruction::SetG2(G4float value 264 void DetectorConstruction::SetG2(G4float value) 258 { 265 { 259 fG2 = value; << 266 G2 = value; 260 G4RunManager::GetRunManager()->ReinitializeG << 261 } 267 } 262 268 263 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 264 270 265 void DetectorConstruction::SetG3(G4float value 271 void DetectorConstruction::SetG3(G4float value) 266 { 272 { 267 fG3 = value; << 273 G3 = value; 268 G4RunManager::GetRunManager()->ReinitializeG << 269 } 274 } 270 275 271 //....oooOO0OOooo........oooOO0OOooo........oo 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 272 277 273 void DetectorConstruction::SetG4(G4float value 278 void DetectorConstruction::SetG4(G4float value) 274 { 279 { 275 fG4 = value; << 280 G4 = value; 276 G4RunManager::GetRunManager()->ReinitializeG << 277 } 281 } 278 282 279 //....oooOO0OOooo........oooOO0OOooo........oo 283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 280 284 281 void DetectorConstruction::SetModel(G4int mode 285 void DetectorConstruction::SetModel(G4int modelChoice) 282 { 286 { 283 if (modelChoice==1) fModel=1; << 287 if (modelChoice==1) model=1; 284 if (modelChoice==2) fModel=2; << 288 if (modelChoice==2) model=2; 285 if (modelChoice==3) fModel=3; << 289 if (modelChoice==3) model=3; 286 G4RunManager::GetRunManager()->ReinitializeG << 287 } 290 } 288 291 289 //....oooOO0OOooo........oooOO0OOooo........oo 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 290 293 291 void DetectorConstruction::SetCoef(G4int val) << 294 #include "G4RunManager.hh" >> 295 >> 296 void DetectorConstruction::UpdateGeometry() 292 { 297 { 293 fCoef=val; << 298 gradientsInitialized=true; 294 G4RunManager::GetRunManager()->ReinitializeG << 299 G4RunManager::GetRunManager()->DefineWorldVolume(ConstructVolumes()); 295 } 300 } 296 301 >> 302 297 //....oooOO0OOooo........oooOO0OOooo........oo 303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 298 304 299 G4int DetectorConstruction::GetCoef() << 305 void DetectorConstruction::SetCoef() 300 { 306 { 301 return fCoef; << 307 coef=1; 302 } 308 } 303 309 304 //....oooOO0OOooo........oooOO0OOooo........oo 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 305 311 306 void DetectorConstruction::SetProfile(G4int my << 312 G4int DetectorConstruction::GetCoef() 307 { 313 { 308 fProfile=myProfile; << 314 return coef; 309 G4RunManager::GetRunManager()->ReinitializeG << 310 } 315 } 311 316 312 //....oooOO0OOooo........oooOO0OOooo........oo 317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 313 318 314 void DetectorConstruction::SetGrid(G4int myGri << 319 void DetectorConstruction::SetProfile(G4int myProfile) 315 { 320 { 316 fGrid=myGrid; << 321 profile=myProfile; 317 G4RunManager::GetRunManager()->ReinitializeG << 318 } 322 } 319 323 320 //....oooOO0OOooo........oooOO0OOooo........oo 324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 321 325 322 void DetectorConstruction::ConstructSDandField << 326 void DetectorConstruction::SetGrid(G4int myGrid) 323 { 327 { 324 fField = new TabulatedField3D(fG1, fG2, << 328 grid=myGrid; 325 << 326 //This is thread-local << 327 G4FieldManager* fFieldMgr = << 328 G4TransportationManager::GetTransportationMa << 329 << 330 G4Mag_UsualEqRhs* fEquation = new G4Mag_ << 331 << 332 G4ClassicalRK4* fStepper = new G4Classic << 333 << 334 G4ChordFinder* fChordFinder = new G4Chor << 335 << 336 fFieldMgr->SetChordFinder(fChordFinder); << 337 fFieldMgr->SetDetectorField(fField); << 338 << 339 // SI: 01-07-2018 : following settings w << 340 // instead of 1e-9*m, but they now indu << 341 // *** G4Exception : GeomNav1002 << 342 // issued by : G4PropagatorInField::Com << 343 << 344 fFieldMgr->GetChordFinder()->SetDeltaCho << 345 fFieldMgr->SetDeltaIntersection(1e-7*m); << 346 fFieldMgr->SetDeltaOneStep(1e-7*m); << 347 << 348 // << 349 << 350 // To avoid G4MagIntegratorDriver::OneGo << 351 << 352 if (fCoef==1) << 353 { << 354 G4PropagatorInField* fPropInField = << 355 G4TransportationManager::GetTranspor << 356 fPropInField->SetMinimumEpsilonStep(1e << 357 fPropInField->SetMaximumEpsilonStep(1e << 358 << 359 } << 360 else << 361 { << 362 G4PropagatorInField* fPropInField = << 363 G4TransportationManager::GetTranspor << 364 fPropInField->SetMinimumEpsilonStep(1e << 365 fPropInField->SetMaximumEpsilonStep(1e << 366 } << 367 << 368 } 329 } 369 330 370 331