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