ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
STCyclotronDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file STCyclotronDetectorConstruction.cc
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 
151 
152 
153  //Overall parameters
154  G4double startAngle = 0.*deg;
155  G4double spanningAngle = 360.*deg;
156 
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-30*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 
275 
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 
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 
378 
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 
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 
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 
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 
544 
545 
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 
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 
626 
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 
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 
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 
818 
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 
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