| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Author: F. Poignant, floriane.poignant@gmail.com
27 //
28 //
29 //
30 //
31 // ******************************************
32 // * *
33 // * STCyclotronDetectorConstruction.cc *
34 // * *
35 // ******************************************
36 //
37 //
38 #include "STCyclotronDetectorConstruction.hh"
39
40 #include "G4RunManager.hh"
41 #include "G4SDManager.hh"
42 #include "G4NistManager.hh"
43
44 #include "G4Element.hh"
45 #include "G4Material.hh"
46
47 #include "G4Box.hh"
48 #include "G4Tubs.hh"
49 #include "G4Cons.hh"
50 #include "G4Polyhedra.hh"
51 #include "G4LogicalVolume.hh"
52 #include "G4PVPlacement.hh"
53 #include "G4Transform3D.hh"
54 #include "G4RotationMatrix.hh"
55 #include "G4UnionSolid.hh"
56 #include "G4SubtractionSolid.hh"
57 #include "G4Region.hh"
58 #include "G4Isotope.hh"
59 #include "G4Element.hh"
60
61 #include "G4PhysicalConstants.hh"
62 #include "G4SystemOfUnits.hh"
63 #include "G4UnitsTable.hh"
64
65 #include "STCyclotronRun.hh"
66 #include "STCyclotronSensitiveTarget.hh"
67 #include "STCyclotronSensitiveFoil.hh"
68 #include "STCyclotronDetectorMessenger.hh"
69 #include "STCyclotronRunAction.hh"
70
71
72 STCyclotronDetectorConstruction::STCyclotronDetectorConstruction()
73 :fTarget_diameter(0),fIsotopeName(0),fIsotopeZ(0),fIsotopeN(0),fIsotopeA(0),
74 fElementName(0),fElementSymbole(0),fElementNComponents(0),fElementAbundance(0),
75 fNaturalElementName(0),fNaturalMaterialFractionMass(0), fDensity_target(0),
76 fTarget_NComponents(0), fMaterialFractionMass(0),fIsotopeNameFoil(0),
77 fIsotopeZFoil(0),fIsotopeNFoil(0),fIsotopeAFoil(0), fElementNameFoil(0),
78 fElementSymboleFoil(0),fElementNComponentsFoil(0),fElementAbundanceFoil(0),
79 fNaturalElementNameFoil(0),fNaturalMaterialFractionMassFoil(0),
80 fDensity_foil(0), fFoil_NComponents(0), fMaterialFractionMassFoil(0),
81 fTarget_thickness(0), fFoil_thickness(0), fTarget_Material(0), fFoil_Material(0),
82 fZ_foil_position(0), fSolidFoil(nullptr),fLogicFoil(nullptr), fPhysFoil(nullptr),
83 fLogicWorld(nullptr),
84 fLayer_z_position_PART3(0), fPhysLayer_PART3(nullptr), fPhysTube_PART3(nullptr),
85 fTube_outerRadius_PART4(0), fTube_length_PART4(0),fLayer_z_position_PART4(0),
86 fPhysTube_PART4(nullptr), fPhysLayer_PART4(nullptr),
87 fLayer1_z_position_PART4(0),fPhysLayer1_PART4(nullptr),
88 fLogicTarget(nullptr), fTarget_z_position(0),fSolidTarget(nullptr),
89 fPhysTarget(nullptr),
90 fLayer1_z_position_PART5(0), fPhysLayer1_PART5(nullptr),
91 fLayer2_z_position_PART5(0), fPhysLayer2_PART5(nullptr),
92 fLayer3_z_position_PART5(0),fPhysLayer3_PART5(nullptr),
93 fRegionTarget(nullptr), fRegionFoil(nullptr), fTargetVolume(0), fFoilVolume(0)
94 {
95 fDetectorMessenger = new STCyclotronDetectorMessenger(this);
96 }
97
98 STCyclotronDetectorConstruction::~STCyclotronDetectorConstruction()
99 {
100 delete fDetectorMessenger;
101 }
102
103 G4VPhysicalVolume* STCyclotronDetectorConstruction::Construct()
104 {
105 //Initialization of messenger parameters
106 fTarget_diameter = 7.*mm;
107 fDensity_target = 8.9*g/cm3;
108 fTarget_thickness = 0.35*mm;
109 fFoil_thickness = 0.000001*mm;
110 fDensity_foil = 2.7*g/cm3;
111
112
113 //Get nist material manager
114 G4NistManager* nist = G4NistManager::Instance();
115
116 // Option to switch on/off checking of volumes overlaps
117 G4bool checkOverlaps = true;
118
119
120 //Create the world
121 G4double world_hx = 1.*m;
122 G4double world_hy = 1.*m;
123 G4double world_hz = 1.*m;
124 G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");
125
126 G4Box* solidWorld
127 = new G4Box("World",
128 world_hx,
129 world_hy,
130 world_hz);
131
132 fLogicWorld
133 = new G4LogicalVolume(solidWorld,
134 world_mat,
135 "World");
136
137 G4VPhysicalVolume* physWorld
138 = new G4PVPlacement(0, //no rotation
139 G4ThreeVector(), //at (0,0,0)
140 fLogicWorld, //its logical volume
141 "World", //its name
142 0, //its mother volume
143 false, //no boolean operation
144 0); //copy number
145
146
147
148 //////////////////////////////////////////////////////////////////////////////
149 ///////////////////////////Create the detector////////////////////////////////
150 //////////////////////////////////////////////////////////////////////////////
151
152
153 //Overall parameters
154 G4double startAngle = 0.*deg;
155 G4double spanningAngle = 360.*deg;
156
157 ////////////////////////////////////
158 /////////Define materials///////////
159 ////////////////////////////////////
160
161
162 //ALUMINIUM//
163
164 G4Material* al = nist->FindOrBuildMaterial("G4_Al");
165
166
167
168 //Create vacuum around the beam//
169 G4double vacuum_atomic_number, vacuum_mass_of_mole, vacuum_density,vacuum_pressure,vacuum_temperature;
170 vacuum_atomic_number = 1.;
171 vacuum_mass_of_mole = 1.008*g/mole;
172 vacuum_density = 1.e-25*g/cm3;
173 vacuum_pressure = 1.e-8*bar;
174 vacuum_temperature = 293.*kelvin; //from PhysicalConstants.h
175 G4Material* vacuum_beam = new G4Material("vacuumBeam", vacuum_atomic_number,
176 vacuum_mass_of_mole, vacuum_density,
177 kStateGas,vacuum_temperature,
178 vacuum_pressure);
179
180
181
182 //HELIUM
183 G4double helium_Z, helium_A, helium_density, helium_pressure, helium_temperature;
184 helium_Z = 2.;
185 helium_A = 4*g/mole;
186 helium_density = 0.1785e-3*g/cm3; // with T=0°, 1 atm. To modify !
187 helium_pressure = 2.*bar;
188 helium_temperature = 293.*kelvin; //15 to 20°
189 G4Material* helium = new G4Material("helium", helium_Z, helium_A, helium_density, kStateGas, helium_temperature, helium_pressure);
190
191
192 //PLATINIUM
193 G4Material* pt = nist->FindOrBuildMaterial("G4_Pt");
194
195
196 //TARGET MATERIAL INITIALIZATION
197 /*G4String name;
198 G4String symbole;
199 G4int ncomponents;
200 G4int n;
201 G4double z_isotope;
202 G4double abundance;
203 G4double fractionmass;
204 G4double a;*/
205
206
207 /*
208 //Pure Ni64
209
210 G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
211
212 G4Element* pureNi64 = new G4Element(name="pureNi64",symbole="64Ni", ncomponents=1);
213 pureNi64->AddIsotope(Ni64, abundance = 100.*perCent);
214
215
216 fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
217 fTarget_Material->AddElement(pureNi64, fractionmass=1.);
218 */
219 /*
220 //Ni64 94% enriched - Sz
221
222 G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
223 G4Isotope* Ni58 = new G4Isotope(name="Zi58", z_isotope=28., n=58, a=58.*g/mole);
224 G4Isotope* Ni60 = new G4Isotope(name="Zi60", z_isotope=28., n=60, a=60.*g/mole);
225 G4Isotope* Ni61 = new G4Isotope(name="Zi61", z_isotope=28., n=61, a=61.*g/mole);
226 G4Isotope* Ni62 = new G4Isotope(name="Zi62", z_isotope=28., n=62, a=62.*g/mole);
227
228
229 G4Element* Ni64enriched95 = new G4Element(name="Ni64enriched95",symbole="64Ni_95", ncomponents=5);
230 Ni64enriched95->AddIsotope(Ni64, abundance = 95.*perCent);
231 Ni64enriched95->AddIsotope(Ni58, abundance = 2.6*perCent);
232 Ni64enriched95->AddIsotope(Ni60, abundance = 1.72*perCent);
233 Ni64enriched95->AddIsotope(Ni61, abundance = 0.15*perCent);
234 Ni64enriched95->AddIsotope(Ni62, abundance = 0.53*perCent);
235
236
237 fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
238 fTarget_Material->AddElement(Ni64enriched95, fractionmass=1.);
239
240 //fTarget_Material = nist->FindOrBuildMaterial("G4_Y");
241 */
242
243 /*
244
245 //Ni64 95% enriched - Obata
246
247 G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
248 G4Isotope* Ni58 = new G4Isotope(name="Zi58", z_isotope=28., n=58, a=58.*g/mole);
249 G4Isotope* Ni60 = new G4Isotope(name="Zi60", z_isotope=28., n=60, a=60.*g/mole);
250 G4Isotope* Ni61 = new G4Isotope(name="Zi61", z_isotope=28., n=61, a=61.*g/mole);
251 G4Isotope* Ni62 = new G4Isotope(name="Zi62", z_isotope=28., n=62, a=62.*g/mole);
252
253
254 G4Element* Ni64enriched95 = new G4Element(name="Ni64enriched95",symbole="64Ni_95", ncomponents=5);
255 Ni64enriched95->AddIsotope(Ni64, abundance = 94.8*perCent);
256 Ni64enriched95->AddIsotope(Ni58, abundance = 2.67*perCent);
257 Ni64enriched95->AddIsotope(Ni60, abundance = 1.75*perCent);
258 Ni64enriched95->AddIsotope(Ni61, abundance = 0.11*perCent);
259 Ni64enriched95->AddIsotope(Ni62, abundance = 0.67*perCent);
260
261
262 fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
263 fTarget_Material->AddElement(Ni64enriched95, fractionmass=1.);
264 */
265
266 fTarget_Material = nist->FindOrBuildMaterial("G4_Ni");
267
268 //FOIL MATERIAL
269 fFoil_Material = nist->FindOrBuildMaterial("G4_Al");
270
271
272 ////////////////////////////////////////////////////////////////////
273 /////////////////////////////PART1//////////////////////////////////
274 ////////////////////////////////////////////////////////////////////
275
276 ////////////////////////////////////////////////////////////////////
277 //////////////////////////LAYER PART 1//////////////////////////////
278 ////////////////////////////////////////////////////////////////////
279
280
281 //Create the (external) layer around the beam vacuum tube
282
283 G4double layer1_length_PART1 = 98.9*mm;
284 G4double layer1_innerRadius_PART1 = 11.5*mm;
285 G4double layer1_outerRadius_PART1 = 16.*mm;
286 G4double layer1_hz_PART1 = 0.5*layer1_length_PART1;
287
288 G4Tubs* solidLayer1_PART1
289 = new G4Tubs("Layer1_PART1",
290 layer1_innerRadius_PART1,
291 layer1_outerRadius_PART1,
292 layer1_hz_PART1,
293 startAngle,
294 spanningAngle);
295
296
297 G4double layer2_length_PART1 = 124.6*mm;
298 G4double layer2_innerRadius_PART1 = 7.5*mm;
299 G4double layer2_outerRadius_PART1 = 11.5*mm;
300 G4double layer2_hz_PART1 = 0.5*layer2_length_PART1;
301
302 G4Tubs* solidLayer2_PART1
303 = new G4Tubs("Layer2_PART1",
304 layer2_innerRadius_PART1,
305 layer2_outerRadius_PART1,
306 layer2_hz_PART1,
307 startAngle,
308 spanningAngle);
309
310
311
312 G4RotationMatrix rot_layer_PART1;
313 G4double z_layer_translation_PART1 = layer2_length_PART1-layer1_length_PART1;
314 G4ThreeVector placement_layer_PART1 = G4ThreeVector(0.*mm,0.*mm,0.5*z_layer_translation_PART1);
315 G4Transform3D transform_layer_PART1(rot_layer_PART1,placement_layer_PART1);
316
317 G4UnionSolid* solidLayer_PART1 = new G4UnionSolid("Layer_PART1", solidLayer2_PART1, solidLayer1_PART1, transform_layer_PART1);
318
319
320 G4LogicalVolume* logicLayer_PART1
321 = new G4LogicalVolume(solidLayer_PART1,
322 al,
323 "Layer_PART1");
324
325 /* G4VPhysicalVolume* physLayer1_PART1 = */
326 new G4PVPlacement(0,
327 G4ThreeVector(0.*mm,0.*mm,0.5*layer2_length_PART1),
328 logicLayer_PART1,
329 "Layer_PART1",
330 fLogicWorld,
331 false,
332 0,
333 checkOverlaps);
334
335
336 //////////////////////////////////////////////////////////////////
337 ////////////////Create the beam vacuum tube///////////////////////
338 //////////////////////////////////////////////////////////////////
339
340
341 G4double tube_length_PART1 = layer2_length_PART1;
342
343 G4double tube_innerRadius_PART1 = 0.*mm;
344 G4double tube_outerRadius_PART1 = 7.5*mm;
345 G4double tube_hz_PART1 = 0.5*tube_length_PART1;
346
347 G4Tubs* solidTube_PART1
348 = new G4Tubs("Tube_PART1",
349 tube_innerRadius_PART1,
350 tube_outerRadius_PART1,
351 tube_hz_PART1,
352 startAngle,
353 spanningAngle);
354
355 G4LogicalVolume* logicTube_PART1
356 = new G4LogicalVolume(solidTube_PART1,
357 vacuum_beam,
358 "Tube_PART1");
359
360 /* G4VPhysicalVolume* physTube_PART1 = */
361 new G4PVPlacement(0,
362 G4ThreeVector(0.*mm,0.*mm,0.5*tube_length_PART1),
363 logicTube_PART1,
364 "Tube_PART1",
365 fLogicWorld,
366 false,
367 0,
368 checkOverlaps);
369
370
371
372
373
374
375 ////////////////////////////////////////////////////////////////////////////////
376 /////////////////////////// PART 2 = FOIL PART /////////////////////////////////
377 ////////////////////////////////////////////////////////////////////////////////
378
379 ////////////////////////////////////////////////////////////////////
380 //////////////////////////LAYER PART 2//////////////////////////////
381 ////////////////////////////////////////////////////////////////////
382
383 //Create the (external) layer around the foil part
384
385 G4double layer_length_PART2 = 12.*mm;
386 G4double layer_innerRadius_PART2 = 7.5*mm;
387 G4double layer_outerRadius_PART2 = 15.*mm;
388 G4double layer_hz_PART2 = 0.5*layer_length_PART2;
389
390 G4Tubs* solidLayer_PART2
391 = new G4Tubs("Layer_PART2",
392 layer_innerRadius_PART2,
393 layer_outerRadius_PART2,
394 layer_hz_PART2,
395 startAngle,
396 spanningAngle);
397
398
399 G4LogicalVolume* logicLayer_PART2
400 = new G4LogicalVolume(solidLayer_PART2,
401 al,
402 "Layer_PART2");
403
404 G4double layer_z_position_PART2 = tube_length_PART1 + 0.5*layer_length_PART2;
405
406 /* G4VPhysicalVolume* physLayer1_PART2 = */
407 new G4PVPlacement(0,
408 G4ThreeVector(0.*mm,0.*mm,layer_z_position_PART2),
409 logicLayer_PART2,
410 "Layer_PART2",
411 fLogicWorld,
412 false,
413 0,
414 checkOverlaps);
415
416
417
418 ////////////////////////////////////////////////////////////////////
419 //////////////////////////TUBE PART 2///////////////////////////////
420 ////////////////////////////////////////////////////////////////////
421
422 //Create the tube before the foil
423
424 G4double tube_length_PART2 = layer_length_PART2;
425
426 G4double tube_innerRadius_PART2 = 0.*mm;
427 G4double tube_outerRadius_PART2 = 7.5*mm;
428 G4double tube_hz_PART2 = 0.5*tube_length_PART2;
429
430 G4Tubs* solidTube_PART2
431 = new G4Tubs("Tube_PART2",
432 tube_innerRadius_PART2,
433 tube_outerRadius_PART2,
434 tube_hz_PART2,
435 startAngle,
436 spanningAngle);
437
438 G4LogicalVolume* logicTube_PART2
439 = new G4LogicalVolume(solidTube_PART2,
440 vacuum_beam,
441 "Tube_PART2");
442
443
444 /* G4VPhysicalVolume* physTube_PART2 = */
445 new G4PVPlacement(0,
446 G4ThreeVector(0.*mm,0.*mm, layer_z_position_PART2),
447 logicTube_PART2,
448 "Tube_PART2",
449 fLogicWorld,
450 false,
451 0,
452 checkOverlaps);
453
454
455 /////////////////////////////////////////////////////////
456 //////////////////////// GRID PART //////////////////////
457 /////////////////////////////////////////////////////////
458
459
460 //HEXAGONE
461 G4double hexagone_length = layer_length_PART2 ;
462 G4int numSide = 6;
463 G4int numZPlanes = 2;
464 const G4double zPlane[]={-0.5*hexagone_length,0.5*hexagone_length};
465 const G4double rInner[]={3.*mm,3.*mm};
466 const G4double rOuter[]={3.3*mm,3.3*mm};
467
468 G4Polyhedra* solidHexagone_0
469 = new G4Polyhedra("Grid",
470 0.*deg,
471 360.*deg,
472 numSide,
473 numZPlanes,
474 zPlane,
475 rInner,
476 rOuter);
477
478
479 G4LogicalVolume* logicHexagone
480 = new G4LogicalVolume(solidHexagone_0,
481 al,
482 "Grid");
483
484 /* G4VPhysicalVolume* physHexagone = */
485 new G4PVPlacement(0,
486 G4ThreeVector(0.*mm,0.*mm,0.*mm),
487 logicHexagone,
488 "Grid",
489 logicTube_PART2,
490 false,
491 0,
492 checkOverlaps);
493
494
495
496
497
498 /////////////////////////////////////////////////////////////////////////
499 //////////////////////////CREATE THE FOIL ///////////////////////////////
500 /////////////////////////////////////////////////////////////////////////
501
502
503 fZ_foil_position = 0.5*fFoil_thickness + tube_length_PART1 + layer_length_PART2;
504
505
506
507 G4double foil_innerRadius = 0.*mm;
508 G4double foil_outerRadius = 16.*mm;
509
510 fSolidFoil
511 = new G4Tubs("Foil",
512 foil_innerRadius,
513 foil_outerRadius,
514 0.5*fFoil_thickness,
515 startAngle,
516 spanningAngle);
517
518 fFoilVolume = fSolidFoil->GetCubicVolume();
519
520 fLogicFoil
521 = new G4LogicalVolume(fSolidFoil,
522 fFoil_Material,
523 "Foil");
524
525
526
527 fPhysFoil
528 = new G4PVPlacement(0,
529 G4ThreeVector(0.*mm,0.*mm,fZ_foil_position),
530 fLogicFoil,
531 "Foil",
532 fLogicWorld,
533 false,
534 0,
535 checkOverlaps);
536
537
538
539
540
541 /////////////////////////////////////////////////////////////////////////////////////
542 /////////////////////////// PART 3 = AFTER THE FOIL /////////////////////////////////
543 /////////////////////////////////////////////////////////////////////////////////////
544
545
546 ////////////////////////////////////////////////////////////////////
547 //////////////////////////LAYER PART 3//////////////////////////////
548 ////////////////////////////////////////////////////////////////////
549
550 //Create the (external) layer around the beam
551
552
553 G4double layer_length_PART3 = 38.32*mm;
554
555 G4double layer_innerRadius_PART3 = 7.5*mm;
556 G4double layer_outerRadius_PART3 = 16.*mm;
557 G4double layer_hz_PART3 = 0.5*layer_length_PART3;
558
559 G4Tubs* solidLayer_PART3
560 = new G4Tubs("Layer_PART3",
561 layer_innerRadius_PART3,
562 layer_outerRadius_PART3,
563 layer_hz_PART3,
564 startAngle,
565 spanningAngle);
566
567 G4LogicalVolume* logicLayer_PART3
568 = new G4LogicalVolume(solidLayer_PART3,
569 al,
570 "Layer_PART3");
571
572 fLayer_z_position_PART3 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + 0.5*layer_length_PART3;
573
574 fPhysLayer_PART3
575 = new G4PVPlacement(0,
576 G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3),
577 logicLayer_PART3,
578 "Layer_PART3",
579 fLogicWorld,
580 false,
581 0,
582 checkOverlaps);
583
584
585
586 ////////////////////////////////////////////////////////////////////
587 //////////////////////////TUBE PART 3///////////////////////////////
588 ////////////////////////////////////////////////////////////////////
589
590 //Create the tube after the foil, MATERIAL = HELIUM
591
592 G4double tube_length_PART3 = layer_length_PART3;
593
594 G4double tube_innerRadius_PART3 = 0.*mm;
595 G4double tube_outerRadius_PART3 = 7.5*mm;
596 G4double tube_hz_PART3 = 0.5*tube_length_PART3;
597
598 G4Tubs* solidTube_PART3
599 = new G4Tubs("Tube_PART3",
600 tube_innerRadius_PART3,
601 tube_outerRadius_PART3,
602 tube_hz_PART3,
603 startAngle,
604 spanningAngle);
605
606 G4LogicalVolume* logicTube_PART3
607 = new G4LogicalVolume(solidTube_PART3,
608 helium,
609 "Tube_PART3");
610
611
612 fPhysTube_PART3
613 = new G4PVPlacement(0,
614 G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3),
615 logicTube_PART3,
616 "Tube_PART3",
617 fLogicWorld,
618 false,
619 0,
620 checkOverlaps);
621
622
623 ////////////////////////////////////////////////////////////////////////////////////////
624 /////////////////////////// PART 4 = BEFORE THE TARGET/////////////////////////////////
625 ////////////////////////////////////////////////////////////////////////////////////////
626
627 ////////////////////////////////////////////////////////////////////
628 //////////////////////////LAYER PART4///////////////////////////////
629 ////////////////////////////////////////////////////////////////////
630
631
632 G4double layer1_length_PART4 = 11.5*mm; //Length of the gold part
633
634 G4double layer2_length_PART4 = 4.6*mm;
635 G4double layer2_Rmin1_PART4 = 8.5*mm;
636 G4double layer2_Rmax1_PART4 = 9.*mm;
637 G4double layer2_Rmin2_PART4 = 8.5*mm;
638 G4double layer2_Rmax2_PART4 = 14.*mm;
639 G4double layer2_hz_PART4 = 0.5*layer2_length_PART4;
640
641 G4Cons* solidLayer2_PART4
642 = new G4Cons("Layer2_PART4",
643 layer2_Rmin1_PART4,
644 layer2_Rmax1_PART4,
645 layer2_Rmin2_PART4,
646 layer2_Rmax2_PART4,
647 layer2_hz_PART4,
648 startAngle,
649 spanningAngle);
650
651 //Create the (external) aluminium layer around the beam before the target
652
653 G4double layer3_length_PART4 = layer1_length_PART4-layer2_length_PART4;
654 G4double layer3_innerRadius_PART4 = 8.5*mm;
655 G4double layer3_outerRadius_PART4 = 14.*mm;
656 G4double layer3_hz_PART4 = 0.5*layer3_length_PART4;
657
658 G4Tubs* solidLayer3_PART4
659 = new G4Tubs("Layer3_PART4",
660 layer3_innerRadius_PART4,
661 layer3_outerRadius_PART4,
662 layer3_hz_PART4,
663 startAngle,
664 spanningAngle);
665
666
667
668 G4RotationMatrix rot_layer_PART4;
669 G4double z_layer_translation_PART4 = 0.5*(layer2_length_PART4+layer3_length_PART4);
670 G4ThreeVector placement_layer_PART4 = G4ThreeVector(0.*mm,0.*mm,z_layer_translation_PART4);
671 G4Transform3D transform_layer_PART4(rot_layer_PART4,placement_layer_PART4);
672
673 G4UnionSolid* solidLayer_PART4 = new G4UnionSolid("Layer_PART1", solidLayer2_PART4, solidLayer3_PART4, transform_layer_PART4);
674
675
676
677 G4LogicalVolume* logicLayer_PART4
678 = new G4LogicalVolume(solidLayer_PART4,
679 al,
680 "Layer_PART4");
681
682 fLayer_z_position_PART4 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +0.5*layer2_length_PART4;
683
684 fPhysLayer_PART4
685 = new G4PVPlacement(0,
686 G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART4),
687 logicLayer_PART4,
688 "Layer_PART4",
689 fLogicWorld,
690 false,
691 0,
692 checkOverlaps);
693
694
695
696
697
698 //Create the (internal) gold layer around the beam
699
700 G4double layer1_innerRadius_PART4 = 8.*mm;
701 G4double layer1_outerRadius_PART4 = 8.5*mm;
702 G4double layer1_hz_PART4 = 0.5*layer1_length_PART4;
703
704 G4Tubs* solidLayer1_PART4
705 = new G4Tubs("Layer1_PART4",
706 layer1_innerRadius_PART4,
707 layer1_outerRadius_PART4,
708 layer1_hz_PART4,
709 startAngle,
710 spanningAngle);
711
712
713 G4LogicalVolume* logicLayer1_PART4
714 = new G4LogicalVolume(solidLayer1_PART4,
715 pt,
716 "Layer1_PART4");
717
718 fLayer1_z_position_PART4 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +0.5*layer1_length_PART4;
719
720
721 fPhysLayer1_PART4
722 = new G4PVPlacement(0,
723 G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4),
724 logicLayer1_PART4,
725 "Layer1_PART4",
726 fLogicWorld,
727 false,
728 0,
729 checkOverlaps);
730
731 ////////////////////////////////////////////////////////////////////
732 //////////////////////////TUBE PART 4///////////////////////////////
733 ////////////////////////////////////////////////////////////////////
734
735
736 //Create the tube before the target
737
738
739 fTube_length_PART4 = layer1_length_PART4;
740
741 G4double tube_innerRadius_PART4 = 0.*mm;
742 fTube_outerRadius_PART4 = 8.*mm;
743 G4double tube_hz_PART4 = 0.5*fTube_length_PART4;
744
745 G4Tubs* solidTube_PART4
746 = new G4Tubs("Tube_PART4",
747 tube_innerRadius_PART4,
748 fTube_outerRadius_PART4,
749 tube_hz_PART4,
750 startAngle,
751 spanningAngle);
752
753 G4LogicalVolume* logicTube_PART4
754 = new G4LogicalVolume(solidTube_PART4,
755 helium,
756 "Tube_PART4");
757
758
759 fPhysTube_PART4
760 = new G4PVPlacement(0,
761 G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4),
762 logicTube_PART4,
763 "Tube_PART4",
764 fLogicWorld,
765 false,
766 0,
767 checkOverlaps);
768
769
770
771
772 ////////////////////////////////////////////////////////////////////
773 //////////////////////////// TARGET ////////////////////////////////
774 ////////////////////////////////////////////////////////////////////
775
776
777
778
779 //Design the target inside the tube
780
781 G4double target_innerRadius = 0.*mm;
782 G4double target_outerRadius = 0.5*fTarget_diameter;
783 G4double target_hz = 0.5*fTarget_thickness;
784
785 fSolidTarget
786 = new G4Tubs("Target",
787 target_innerRadius,
788 target_outerRadius,
789 target_hz,
790 startAngle,
791 spanningAngle);
792
793 fTargetVolume = fSolidTarget->GetCubicVolume();
794
795 fLogicTarget
796 = new G4LogicalVolume(fSolidTarget,
797 fTarget_Material,
798 "Target");
799
800 fTarget_z_position = 0.5*fTube_length_PART4-0.5*fTarget_thickness; //inside the tube!
801
802
803 fPhysTarget
804 = new G4PVPlacement(0,
805 G4ThreeVector(0.*mm,0.*mm, fTarget_z_position),
806 fLogicTarget,
807 "Target",
808 logicTube_PART4,
809 false,
810 0,
811 checkOverlaps);
812
813
814
815 ///////////////////////////////////////////////////////////////////////////////////////
816 /////////////////////////// PART 5 = AFTER THE TARGET /////////////////////////////////
817 ///////////////////////////////////////////////////////////////////////////////////////
818
819 ////////////////////////////////////////////////////////////////////
820 //////////////////////////LAYER PART 5//////////////////////////////
821 ////////////////////////////////////////////////////////////////////
822
823 //Create the (internal) platinium layer
824
825 G4double layer1_length_PART5 = 0.5*mm;
826
827 G4double layer1_innerRadius_PART5 = 0.*mm;
828 G4double layer1_outerRadius_PART5 = 8.5*mm;
829 G4double layer1_hz_PART5 = 0.5*layer1_length_PART5;
830
831 G4Tubs* solidLayer1_PART5
832 = new G4Tubs("Layer1_PART5",
833 layer1_innerRadius_PART5,
834 layer1_outerRadius_PART5,
835 layer1_hz_PART5,
836 startAngle,
837 spanningAngle);
838
839 G4LogicalVolume* logicLayer1_PART5
840 = new G4LogicalVolume(solidLayer1_PART5,
841 pt,
842 "Layer1_PART5");
843
844 fLayer1_z_position_PART5 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +layer1_length_PART4 + 0.5*layer1_length_PART5;
845
846 fPhysLayer1_PART5
847 = new G4PVPlacement(0,
848 G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART5),
849 logicLayer1_PART5,
850 "Layer1_PART5",
851 fLogicWorld,
852 false,
853 0,
854 checkOverlaps);
855
856
857
858 //Create the (external) aluminium layer after the target
859
860 G4double layer2_length_PART5 = 0.5*mm;
861
862 G4double layer2_innerRadius_PART5 = 8.5*mm;
863 G4double layer2_outerRadius_PART5 = 14.*mm;
864 G4double layer2_hz_PART5 = 0.5*layer2_length_PART5;
865
866 G4Tubs* solidLayer2_PART5
867 = new G4Tubs("Layer2_PART5",
868 layer2_innerRadius_PART5,
869 layer2_outerRadius_PART5,
870 layer2_hz_PART5,
871 startAngle,
872 spanningAngle);
873
874 G4LogicalVolume* logicLayer2_PART5
875 = new G4LogicalVolume(solidLayer2_PART5,
876 al,
877 "Layer2_PART5");
878
879 fLayer2_z_position_PART5 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +layer2_length_PART4+layer3_length_PART4 + 0.5*layer2_length_PART5;
880
881 fPhysLayer2_PART5
882 = new G4PVPlacement(0,
883 G4ThreeVector(0.*mm,0.*mm,fLayer2_z_position_PART5),
884 logicLayer2_PART5,
885 "Layer2_PART5",
886 fLogicWorld,
887 false,
888 0,
889 checkOverlaps);
890
891
892
893 //Create the end of the beam target machine
894
895 G4double layer3_length_PART5 = 3.*mm;
896
897 G4double layer3_innerRadius_PART5 = 0.*mm;
898 G4double layer3_outerRadius_PART5 = 14.*mm;
899 G4double layer3_hz_PART5 = 0.5*layer3_length_PART5;
900
901 G4Tubs* solidLayer3_PART5
902 = new G4Tubs("Layer3_PART5",
903 layer3_innerRadius_PART5,
904 layer3_outerRadius_PART5,
905 layer3_hz_PART5,
906 startAngle,
907 spanningAngle);
908
909 G4LogicalVolume* logicLayer3_PART5
910 = new G4LogicalVolume(solidLayer3_PART5,
911 al,
912 "Layer3_PART5");
913
914 fLayer3_z_position_PART5 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +layer2_length_PART4+layer3_length_PART4 + layer2_length_PART5 + 0.5*layer3_length_PART5;
915
916 fPhysLayer3_PART5
917 = new G4PVPlacement(0,
918 G4ThreeVector(0.*mm,0.*mm,fLayer3_z_position_PART5),
919 logicLayer3_PART5,
920 "Layer3_PART5",
921 fLogicWorld,
922 false,
923 0,
924 checkOverlaps);
925
926
927
928
929 //////////////////////////////////////////////
930 /////////// Set sensitive region /////////////
931 //////////////////////////////////////////////
932
933
934 if(!fRegionTarget)
935 {
936 fRegionTarget = new G4Region("Target");
937 fLogicTarget -> SetRegion(fRegionTarget);
938 fRegionTarget -> AddRootLogicalVolume(fLogicTarget);
939 }
940
941 if(!fRegionFoil&&fPhysFoil!=nullptr)
942 {
943 fRegionFoil = new G4Region("Foil");
944 fLogicFoil -> SetRegion(fRegionFoil);
945 fRegionFoil -> AddRootLogicalVolume(fLogicFoil);
946 }
947
948
949
950
951 //Always return physical world//
952
953 return physWorld;
954
955
956 }
957
958 void STCyclotronDetectorConstruction::ConstructSDandField()
959 {
960
961 if(fLogicTarget != nullptr){
962 STCyclotronSensitiveTarget* TargetSD = new STCyclotronSensitiveTarget("TargetSD", this);
963 SetSensitiveDetector("Target",TargetSD);
964 }
965
966
967 if(fLogicFoil != nullptr){
968 STCyclotronSensitiveFoil* FoilSD = new STCyclotronSensitiveFoil("FoilSD", this);
969 SetSensitiveDetector("Foil",FoilSD);
970 }
971 }
972
973 void STCyclotronDetectorConstruction::SetTargetDiameter(G4double targetDiameter)
974 {
975
976 if(fTarget_diameter!=targetDiameter){
977 if(targetDiameter/2.>fTube_outerRadius_PART4){ G4cout << "Error : the diameter is bigger than the tube" << G4endl;}
978 else{
979 fTarget_diameter = targetDiameter;
980 if(fSolidTarget) fSolidTarget->SetOuterRadius(0.5*fTarget_diameter);
981 G4cout << "The new diameter of the target is " << fTarget_diameter << " mm." << G4endl;
982 }
983 }
984
985 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
986 G4cout << "... Geometry is notified .... " << G4endl;
987
988 }
989
990 //SET MATERIAL METHODS//
991
992 void STCyclotronDetectorConstruction::SetTargetIsotopeName(G4String name){
993 fIsotopeName.push_back(name);
994 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
995 G4cout << "... Geometry is notified .... " << G4endl;
996 }
997
998 void STCyclotronDetectorConstruction::SetTargetIsotopeZ(G4double z){
999 fIsotopeZ.push_back(z);
1000 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1001 G4cout << "... Geometry is notified .... " << G4endl;
1002 }
1003
1004 void STCyclotronDetectorConstruction::SetTargetIsotopeN(G4int n){
1005 fIsotopeN.push_back(n);
1006 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1007 G4cout << "... Geometry is notified .... " << G4endl;
1008 }
1009
1010 void STCyclotronDetectorConstruction::SetTargetIsotopeA(G4double a){
1011 fIsotopeA.push_back(a);
1012 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1013 G4cout << "... Geometry is notified .... " << G4endl;
1014 }
1015
1016 void STCyclotronDetectorConstruction::SetTargetElementName(G4String name){
1017 fElementName.push_back(name);
1018 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1019 G4cout << "... Geometry is notified .... " << G4endl;
1020 }
1021
1022 void STCyclotronDetectorConstruction::SetTargetElementSymbole(G4String symbole){
1023 fElementSymbole.push_back(symbole);
1024 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1025 G4cout << "... Geometry is notified .... " << G4endl;
1026 }
1027
1028 void STCyclotronDetectorConstruction::SetTargetElementNComponents(G4int n){
1029 fElementNComponents.push_back(n);
1030 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1031 G4cout << "... Geometry is notified .... " << G4endl;
1032 }
1033
1034 void STCyclotronDetectorConstruction::SetTargetElementAbundance(G4double abundance){
1035 fElementAbundance.push_back(abundance);
1036 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1037 G4cout << "... Geometry is notified .... " << G4endl;
1038 }
1039
1040 void STCyclotronDetectorConstruction::SetTargetMaterialDensity(G4double materialDensity){
1041 fDensity_target = materialDensity;
1042 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1043 G4cout << "... Geometry is notified .... " << G4endl;
1044 }
1045
1046 void STCyclotronDetectorConstruction::SetTargetMaterialNComponents(G4int n){
1047 fTarget_NComponents = n;
1048 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1049 G4cout << "... Geometry is notified .... " << G4endl;
1050 }
1051
1052 void STCyclotronDetectorConstruction::SetTargetMaterialFractionMass(G4double fractionMass){
1053 fMaterialFractionMass.push_back(fractionMass);
1054 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1055 G4cout << "... Geometry is notified .... " << G4endl;
1056 }
1057
1058 void STCyclotronDetectorConstruction::SetTargetNaturalElement(G4String name){
1059 fNaturalElementName.push_back(name);
1060 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1061 G4cout << "... Geometry is notified .... " << G4endl;
1062 }
1063
1064 void STCyclotronDetectorConstruction::SetTargetNaturalMaterialFractionMass(G4double fractionMass){
1065 fNaturalMaterialFractionMass.push_back(fractionMass);
1066 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1067 G4cout << "... Geometry is notified .... " << G4endl;
1068 }
1069
1070 //UPDATE THE MATERIAL TO APPLY THE MODIFICATIONS
1071
1072 G4bool STCyclotronDetectorConstruction::UpdateMaterial(){
1073
1074 G4int nElements = fTarget_NComponents;
1075
1076 G4double density = fDensity_target*g/cm3;
1077 G4Material* material = new G4Material("FTarget_Material", density, nElements);
1078
1079 G4double checkFractionMass = 0;
1080
1081
1082 //CHECK THAT NUMBER OF ELEMENTS IN THE MATERIAL IS EQUAL TO THE NUMBER OF ELEMENTS DEFINED USING ISOTOPES PLUS THE NUMBER OF ELEMENT DEFINED USING NIST
1083
1084 G4int counterElement = 0;
1085
1086 //std::vector<G4String>::size_type sz = fElementName.size();
1087 std::ptrdiff_t const sizeElementName = fElementName.size();
1088 std::ptrdiff_t const sizeNaturalElementName = fNaturalElementName.size();
1089
1090
1091 if(sizeElementName!=0){
1092
1093 //ELEMENT LOOP
1094 if(nElements != sizeElementName + sizeNaturalElementName){
1095 G4cout << "Error : the number of elements in the target was set up at " << nElements << " but you defined only " << fElementName.size()+fNaturalElementName.size() << "elements." << G4endl ;
1096 return false;
1097 }
1098
1099 for(int i = 0; i<sizeElementName ; i++){
1100
1101 checkFractionMass = checkFractionMass + fMaterialFractionMass[i];
1102 G4Element* elementi = new G4Element(fElementName[i],fElementSymbole[i],fElementNComponents[i]);
1103
1104 G4double checkAbundance = 0.;
1105
1106 //ISOTOPE LOOP FOR THE ELEMENT
1107 std::ptrdiff_t const sizeIsotopeName = fIsotopeName.size();
1108
1109 if(fElementNComponents[i] != sizeIsotopeName){
1110 G4cout << "Error : the number of isotopes defined in the target's element" << fElementName[i] << "was set up at " << fElementNComponents[i] << " but you defined only " << fIsotopeName.size() << "elements." << G4endl;
1111 return false;
1112 }
1113 for(G4int j=counterElement;j<counterElement+fElementNComponents[i];++j){
1114 G4double A = fIsotopeA[j]*g/mole;
1115 G4Isotope* isotopeij = new G4Isotope(fIsotopeName[j], G4int(fIsotopeZ[j]), fIsotopeN[j], A);
1116 checkAbundance = checkAbundance + fElementAbundance[j];
1117 elementi->AddIsotope(isotopeij,fElementAbundance[j]);
1118 }
1119 if((1.-checkAbundance)>1E-5){
1120 G4cout << "Error : the total abundance of isotopes in your target's element " << fElementName[i] << " is equal to " << checkAbundance << ". It must be equal to 1." << G4endl;
1121 return false;
1122 }
1123
1124 counterElement = counterElement + fElementNComponents[i];
1125 material->AddElement(elementi, fMaterialFractionMass[i]);
1126
1127 }
1128 }
1129
1130
1131 if(sizeNaturalElementName!=0){
1132 for(int i=0;i<sizeNaturalElementName;i++){
1133 checkFractionMass = checkFractionMass + fNaturalMaterialFractionMass[i];
1134 G4NistManager* man = G4NistManager::Instance();
1135 G4Element* element = man->FindOrBuildElement(fNaturalElementName[i]);
1136 material->AddElement(element, fNaturalMaterialFractionMass[i]);
1137 }
1138 }
1139
1140 if((1.-checkFractionMass)>1E-5){
1141 G4cout << "Error : the sum of the fraction mass of each element in the target is equal to " << checkFractionMass << ". It must be equal to 1." << G4endl;
1142 return false;
1143 }
1144
1145 fTarget_Material = material;
1146 G4cout << "You succesfully changed the material of the target." << G4endl;
1147 G4cout << "Here is the new material : " << G4endl;
1148 G4cout << "Number of elements : " << material->GetNumberOfElements() << G4endl;
1149
1150 std::ptrdiff_t const sizeMaterial = material->GetNumberOfElements();
1151
1152
1153 for(int i = 0; i<sizeMaterial;i++){
1154 const G4Element* element = material->GetElement(i);
1155 G4cout << element->GetName() << " having the following isotopes : " << G4endl;
1156
1157 std::ptrdiff_t const sizeIsotope = element->GetNumberOfIsotopes();
1158
1159
1160 for(int j=0; j<sizeIsotope; j++){
1161 const G4Isotope* isotope = element->GetIsotope(j);
1162 G4cout << " - " << (element->GetRelativeAbundanceVector())[j]*100 << "% of " << isotope->GetName() << " Z = " << isotope->GetZ() << " N = " << isotope->GetN() << "." << G4endl;
1163 }
1164 }
1165 fLogicTarget->SetMaterial(fTarget_Material);
1166
1167 //
1168 fIsotopeName.clear();
1169 fIsotopeZ.clear();
1170 fIsotopeN.clear();
1171 fIsotopeA.clear();
1172 fElementName.clear();
1173 fElementSymbole.clear();
1174 fElementNComponents.clear();
1175 fElementAbundance.clear();
1176 fNaturalElementName.clear();
1177 fNaturalMaterialFractionMass.clear();
1178 fMaterialFractionMass.clear();
1179 //
1180
1181 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1182 G4cout << "... Geometry is notified .... " << G4endl;
1183
1184
1185 return true;
1186
1187
1188
1189 }
1190
1191 //SET TARGET MATERIAL WITH PHYSICS NIST LIST//
1192
1193 void STCyclotronDetectorConstruction::SetTargetMaterial(G4String materialName){
1194
1195 G4NistManager* nistManager = G4NistManager::Instance();
1196 G4Material* targetMaterial = nistManager->FindOrBuildMaterial(materialName);
1197 if(fTarget_Material!=targetMaterial){
1198
1199 if(targetMaterial){
1200
1201 fTarget_Material = targetMaterial;
1202 fLogicTarget->SetMaterial(fTarget_Material);
1203
1204 G4cout << "The new material of the target is : " << fTarget_Material << G4endl;
1205
1206 }
1207 else{ G4cout << "This material wasn't found in the NIST list." << G4endl; }
1208
1209 }
1210
1211 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1212 G4cout << "... Geometry is notified .... " << G4endl;
1213
1214 }
1215
1216
1217 void STCyclotronDetectorConstruction::SetFoilIsotopeName(G4String name){
1218 fIsotopeNameFoil.push_back(name);
1219 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1220 G4cout << "... Geometry is notified .... " << G4endl;
1221 }
1222
1223 void STCyclotronDetectorConstruction::SetFoilIsotopeZ(G4double z){
1224 fIsotopeZFoil.push_back(z);
1225 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1226 G4cout << "... Geometry is notified .... " << G4endl;
1227 }
1228
1229 void STCyclotronDetectorConstruction::SetFoilIsotopeN(G4int n){
1230 fIsotopeNFoil.push_back(n);
1231 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1232 G4cout << "... Geometry is notified .... " << G4endl;
1233 }
1234
1235 void STCyclotronDetectorConstruction::SetFoilIsotopeA(G4double a){
1236 fIsotopeAFoil.push_back(a);
1237 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1238 G4cout << "... Geometry is notified .... " << G4endl;
1239 }
1240
1241 void STCyclotronDetectorConstruction::SetFoilElementName(G4String name){
1242 fElementNameFoil.push_back(name);
1243 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1244 G4cout << "... Geometry is notified .... " << G4endl;
1245 }
1246
1247 void STCyclotronDetectorConstruction::SetFoilElementSymbole(G4String symbole){
1248 fElementSymboleFoil.push_back(symbole);
1249 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1250 G4cout << "... Geometry is notified .... " << G4endl;
1251 }
1252
1253 void STCyclotronDetectorConstruction::SetFoilElementNComponents(G4int n){
1254 fElementNComponentsFoil.push_back(n);
1255 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1256 G4cout << "... Geometry is notified .... " << G4endl;
1257 }
1258
1259 void STCyclotronDetectorConstruction::SetFoilElementAbundance(G4double abundance){
1260 fElementAbundanceFoil.push_back(abundance);
1261 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1262 G4cout << "... Geometry is notified .... " << G4endl;
1263 }
1264
1265 void STCyclotronDetectorConstruction::SetFoilMaterialDensity(G4double materialDensity){
1266 fDensity_foil = materialDensity;
1267 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1268 G4cout << "... Geometry is notified .... " << G4endl;
1269 }
1270
1271 void STCyclotronDetectorConstruction::SetFoilMaterialNComponents(G4int n){
1272 fFoil_NComponents = n;
1273 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1274 G4cout << "... Geometry is notified .... " << G4endl;
1275 }
1276
1277 void STCyclotronDetectorConstruction::SetFoilMaterialFractionMass(G4double fractionMass){
1278 fMaterialFractionMassFoil.push_back(fractionMass);
1279 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1280 G4cout << "... Geometry is notified .... " << G4endl;
1281 }
1282
1283 void STCyclotronDetectorConstruction::SetFoilNaturalElement(G4String name){
1284 fNaturalElementNameFoil.push_back(name);
1285 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1286 G4cout << "... Geometry is notified .... " << G4endl;
1287 }
1288
1289 void STCyclotronDetectorConstruction::SetFoilNaturalMaterialFractionMass(G4double fractionMass){
1290 fNaturalMaterialFractionMassFoil.push_back(fractionMass);
1291 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1292 G4cout << "... Geometry is notified .... " << G4endl;
1293 }
1294
1295 G4bool STCyclotronDetectorConstruction::UpdateFoilMaterial(){
1296
1297 G4int nElements = fFoil_NComponents;
1298 G4double density = fDensity_foil*g/cm3;
1299
1300 G4Material* material = new G4Material("FFoil_Material", density,nElements);
1301
1302 G4double checkFractionMass = 0;
1303
1304
1305 //CHECK THAT NUMBER OF ELEMENTS IN THE MATERIAL IS EQUAL TO THE NUMBER OF ELEMENTS DEFINED USING ISOTOPES PLUS THE NUMBER OF ELEMENT DEFINED USING NIST
1306
1307 G4int counterElement = 0;
1308
1309 std::ptrdiff_t const sizeElementName = fElementNameFoil.size();
1310 std::ptrdiff_t const sizeNaturalElementName = fNaturalElementNameFoil.size();
1311
1312
1313
1314 if(sizeElementName!=0){
1315 //ELEMENT LOOP
1316 if(nElements != sizeElementName + sizeNaturalElementName){
1317 G4cout << "Error : the number of elements was set up at " << nElements << " but you defined only " << fElementNameFoil.size()+fNaturalElementNameFoil.size() << "elements." << G4endl ;
1318 return false;
1319 }
1320 for(int i = 0; i<sizeElementName ; i++){
1321
1322 checkFractionMass = checkFractionMass + fMaterialFractionMassFoil[i];
1323 G4Element* elementi = new G4Element(fElementNameFoil[i],fElementSymboleFoil[i],fElementNComponentsFoil[i]);
1324
1325 G4double checkAbundance = 0.;
1326 //ISOTOPE LOOP FOR THE ELEMENT
1327
1328 std::ptrdiff_t const sizeIsotopeName = fIsotopeNameFoil.size();
1329
1330 if(fElementNComponentsFoil[i] != sizeIsotopeName){
1331 G4cout << "Error : the number of isotopes in element" << fElementNameFoil[i] << " of the foil was set up at " << fElementNComponentsFoil[i] << " but you defined only " << fIsotopeNameFoil.size() << "elements." << G4endl;
1332 return false;
1333 }
1334 for(G4int j=counterElement;j<counterElement+fElementNComponentsFoil[i];++j){
1335 G4double A = fIsotopeAFoil[j]*g/cm3;
1336 G4Isotope* isotopeij = new G4Isotope(fIsotopeNameFoil[j], G4int(fIsotopeZFoil[j]), fIsotopeNFoil[j],A);
1337 checkAbundance = checkAbundance + fElementAbundanceFoil[j];
1338 elementi->AddIsotope(isotopeij,fElementAbundanceFoil[j]);
1339 }
1340 if((1.-checkAbundance)>1E-5){
1341 G4cout << "Error : the total abundance of isotopes in your foil's element " << fElementNameFoil[i] << " is equal to " << checkAbundance << ". It must be equal to 1." << G4endl;
1342 return false;
1343 }
1344
1345 counterElement = counterElement + fElementNComponentsFoil[i];
1346 material->AddElement(elementi, fMaterialFractionMassFoil[i]);
1347
1348 }
1349 }
1350
1351 if(sizeNaturalElementName!=0){
1352 for(int i=0;i<sizeNaturalElementName;i++){
1353 checkFractionMass = checkFractionMass + fNaturalMaterialFractionMassFoil[i];
1354 G4NistManager* man = G4NistManager::Instance();
1355 G4Element* element = man->FindOrBuildElement(fNaturalElementNameFoil[i]);
1356 material->AddElement(element, fNaturalMaterialFractionMassFoil[i]);
1357 }
1358 }
1359
1360 if((1.-checkFractionMass)>1E-5){
1361 G4cout << "Error : the sum of the fraction mass of each element in the foil is equal to " << checkFractionMass << ". It must be equal to 1." << G4endl;
1362 return false;
1363 }
1364
1365 fFoil_Material = material;
1366 G4cout << "You succesfully changed the material of the foil." << G4endl;
1367 G4cout << "Here is the new material : " << G4endl;
1368 G4cout << "Number of elements : " << material->GetNumberOfElements() << G4endl;
1369 std::ptrdiff_t const sizeMaterial = material->GetNumberOfElements();
1370
1371
1372 for(int i = 0; i<sizeMaterial;i++){
1373 const G4Element* element = material->GetElement(i);
1374 G4cout << element->GetName() << " having the following isotopes : " << G4endl;
1375
1376 std::ptrdiff_t const sizeIsotope = element->GetNumberOfIsotopes();
1377
1378 for(int j=0; j<sizeIsotope; j++){
1379 const G4Isotope* isotope = element->GetIsotope(j);
1380 G4cout << " - " << (element->GetRelativeAbundanceVector())[j]*100 << "% of " << isotope->GetName() << " Z = " << isotope->GetZ() << " N = " << isotope->GetN() << "." << G4endl;
1381 }
1382 }
1383 fLogicFoil->SetMaterial(fFoil_Material);
1384
1385 //
1386 fIsotopeNameFoil.clear();
1387 fIsotopeZFoil.clear();
1388 fIsotopeNFoil.clear();
1389 fIsotopeAFoil.clear();
1390 fElementNameFoil.clear();
1391 fElementSymboleFoil.clear();
1392 fElementNComponentsFoil.clear();
1393 fElementAbundanceFoil.clear();
1394 fNaturalElementNameFoil.clear();
1395 fNaturalMaterialFractionMassFoil.clear();
1396 fMaterialFractionMassFoil.clear();
1397 //
1398
1399 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1400 G4cout << "... Geometry is notified .... " << G4endl;
1401
1402 return true;
1403 }
1404
1405
1406 void STCyclotronDetectorConstruction::SetFoilMaterial(G4String materialName){
1407
1408 G4NistManager* nistManager = G4NistManager::Instance();
1409 G4Material* foilMaterial = nistManager->FindOrBuildMaterial(materialName);
1410 if(fFoil_Material!=foilMaterial){
1411
1412 if(foilMaterial){
1413
1414 fFoil_Material = foilMaterial;
1415 fLogicFoil->SetMaterial(fFoil_Material);
1416
1417 G4cout << "The new material of the foil is : " << fFoil_Material << G4endl;
1418
1419 }
1420 else{ G4cout << "This material wasn't found in the NIST list." << G4endl; }
1421
1422 }
1423 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1424 G4cout << "... Geometry is notified .... " << G4endl;
1425
1426 }
1427
1428 //SET PARAMETERS FOR THE TARGET
1429
1430 void STCyclotronDetectorConstruction::SetTargetThickness(G4double targetThickness)
1431 {
1432
1433 if(fTarget_thickness!=targetThickness){
1434 if(targetThickness > fTube_length_PART4){ G4cout << "Error : the target thickness is longer than the tube length" << G4endl; }
1435 else {
1436 //Modify the position of the target
1437 fTarget_z_position = fTarget_z_position + 0.5*fTarget_thickness;
1438 fTarget_thickness = targetThickness;
1439 fTarget_z_position = fTarget_z_position - 0.5*fTarget_thickness;
1440 fPhysTarget->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fTarget_z_position*mm));
1441 //Modify the thickness of the target
1442 fSolidTarget->SetZHalfLength(0.5*targetThickness);
1443 G4cout << "The new thickness of the target is " << fTarget_thickness << " mm." << G4endl;
1444 G4cout << "Position 1 = " << GetTargetPosition1() << " -- Position 2 = " << GetTargetPosition2() << G4endl;
1445
1446 }
1447 }
1448
1449 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1450 G4cout << "... Geometry is notified .... " << G4endl;
1451 }
1452
1453
1454 void STCyclotronDetectorConstruction::SetFoilThickness(G4double foilThickness)
1455 {
1456
1457 if(fFoil_thickness != foilThickness){
1458
1459
1460 //Change de position of the detector parts after the foil
1461 //Change the position to set it so there is no foil
1462 fLayer_z_position_PART3 = fLayer_z_position_PART3 - fFoil_thickness;
1463 fLayer_z_position_PART4 = fLayer_z_position_PART4 - fFoil_thickness;
1464 fLayer1_z_position_PART4 = fLayer1_z_position_PART4 - fFoil_thickness;
1465 fLayer1_z_position_PART5 = fLayer1_z_position_PART5 - fFoil_thickness;
1466 fLayer2_z_position_PART5 = fLayer2_z_position_PART5 - fFoil_thickness;
1467 fLayer3_z_position_PART5 = fLayer3_z_position_PART5 - fFoil_thickness;
1468 fZ_foil_position = fZ_foil_position - 0.5*fFoil_thickness;
1469
1470 fFoil_thickness = foilThickness;
1471
1472 //Change the position using the new foil thickness
1473 fLayer_z_position_PART3 = fLayer_z_position_PART3 + fFoil_thickness;
1474 fLayer_z_position_PART4 = fLayer_z_position_PART4 + fFoil_thickness;
1475 fLayer1_z_position_PART4 = fLayer1_z_position_PART4 + fFoil_thickness;
1476 fLayer1_z_position_PART5 = fLayer1_z_position_PART5 + fFoil_thickness;
1477 fLayer2_z_position_PART5 = fLayer2_z_position_PART5 + fFoil_thickness;
1478 fLayer3_z_position_PART5 = fLayer3_z_position_PART5 + fFoil_thickness;
1479 fZ_foil_position = fZ_foil_position + 0.5*fFoil_thickness;
1480
1481 fPhysLayer_PART3->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3*mm));
1482 fPhysTube_PART3->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3*mm));
1483 fPhysLayer_PART4->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART4*mm));
1484 fPhysLayer1_PART4->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4*mm));
1485 fPhysTube_PART4->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4*mm));
1486 fPhysLayer1_PART5->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART5*mm));
1487 fPhysLayer2_PART5->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer2_z_position_PART5*mm));
1488 fPhysLayer3_PART5->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fLayer3_z_position_PART5*mm));
1489
1490
1491 if(fPhysFoil) fPhysFoil->SetTranslation(G4ThreeVector(0.*mm,0.*mm,fZ_foil_position*mm));
1492 //Change the thickness
1493 if(fSolidFoil) fSolidFoil->SetZHalfLength(0.5*foilThickness*mm);
1494 G4cout << "The new thickness of the foil is " << fFoil_thickness << " mm." << G4endl;
1495
1496 }
1497
1498 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
1499 G4cout << "... Geometry is notified .... " << G4endl;
1500
1501 }
1502
1503
1504