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