STCyclotron/src/STCyclotronDetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Author: F. Poignant, floriane.poignant@gmail.com
0027 //
0028 //
0029 //
0030 //
0031 //    ******************************************
0032 //    *                                        *
0033 //    *    STCyclotronDetectorConstruction.cc  *
0034 //    *                                        *
0035 //    ******************************************
0036 //
0037 //
0038 #include "STCyclotronDetectorConstruction.hh"
0039
0040 #include "G4RunManager.hh"
0041 #include "G4SDManager.hh"
0042 #include "G4NistManager.hh"
0043
0044 #include "G4Element.hh"
0045 #include "G4Material.hh"
0046
0047 #include "G4Box.hh"
0048 #include "G4Tubs.hh"
0049 #include "G4Cons.hh"
0050 #include "G4Polyhedra.hh"
0051 #include "G4LogicalVolume.hh"
0052 #include "G4PVPlacement.hh"
0053 #include "G4Transform3D.hh"
0054 #include "G4RotationMatrix.hh"
0055 #include "G4UnionSolid.hh"
0056 #include "G4SubtractionSolid.hh"
0057 #include "G4Region.hh"
0058 #include "G4Isotope.hh"
0059 #include "G4Element.hh"
0060
0061 #include "G4PhysicalConstants.hh"
0062 #include "G4SystemOfUnits.hh"
0063 #include "G4UnitsTable.hh"
0064
0065 #include "STCyclotronRun.hh"
0066 #include "STCyclotronSensitiveTarget.hh"
0067 #include "STCyclotronSensitiveFoil.hh"
0068 #include "STCyclotronDetectorMessenger.hh"
0069 #include "STCyclotronRunAction.hh"
0070
0071
0072 STCyclotronDetectorConstruction::STCyclotronDetectorConstruction()
0073  :fTarget_diameter(0),fIsotopeName(0),fIsotopeZ(0),fIsotopeN(0),fIsotopeA(0),
0074   fElementName(0),fElementSymbole(0),fElementNComponents(0),fElementAbundance(0),
0075   fNaturalElementName(0),fNaturalMaterialFractionMass(0), fDensity_target(0),
0076   fTarget_NComponents(0), fMaterialFractionMass(0),fIsotopeNameFoil(0),
0077   fIsotopeZFoil(0),fIsotopeNFoil(0),fIsotopeAFoil(0), fElementNameFoil(0),
0078   fElementSymboleFoil(0),fElementNComponentsFoil(0),fElementAbundanceFoil(0),
0079   fNaturalElementNameFoil(0),fNaturalMaterialFractionMassFoil(0),
0080   fDensity_foil(0), fFoil_NComponents(0), fMaterialFractionMassFoil(0),
0081   fTarget_thickness(0), fFoil_thickness(0), fTarget_Material(0), fFoil_Material(0),
0082   fZ_foil_position(0), fSolidFoil(nullptr),fLogicFoil(nullptr), fPhysFoil(nullptr),
0083   fLogicWorld(nullptr),
0084   fLayer_z_position_PART3(0), fPhysLayer_PART3(nullptr), fPhysTube_PART3(nullptr),
0085   fTube_outerRadius_PART4(0), fTube_length_PART4(0),fLayer_z_position_PART4(0),
0086   fPhysTube_PART4(nullptr), fPhysLayer_PART4(nullptr),
0087   fLayer1_z_position_PART4(0),fPhysLayer1_PART4(nullptr),
0088   fLogicTarget(nullptr), fTarget_z_position(0),fSolidTarget(nullptr),
0089   fPhysTarget(nullptr),
0090   fLayer1_z_position_PART5(0), fPhysLayer1_PART5(nullptr),
0091   fLayer2_z_position_PART5(0), fPhysLayer2_PART5(nullptr),
0092   fLayer3_z_position_PART5(0),fPhysLayer3_PART5(nullptr),
0093   fRegionTarget(nullptr), fRegionFoil(nullptr), fTargetVolume(0), fFoilVolume(0)
0094 {
0095   fDetectorMessenger = new STCyclotronDetectorMessenger(this);
0096 }
0097
0098 STCyclotronDetectorConstruction::~STCyclotronDetectorConstruction()
0099 {
0100   delete fDetectorMessenger;
0101 }
0102
0103 G4VPhysicalVolume* STCyclotronDetectorConstruction::Construct()
0104 {
0105   //Initialization of messenger parameters
0106   fTarget_diameter = 7.*mm;
0107   fDensity_target = 8.9*g/cm3;
0108   fTarget_thickness = 0.35*mm;
0109   fFoil_thickness = 0.000001*mm;
0110   fDensity_foil = 2.7*g/cm3;
0111
0112
0113   //Get nist material manager
0114   G4NistManager* nist = G4NistManager::Instance();
0115
0116   // Option to switch on/off checking of volumes overlaps
0117   G4bool checkOverlaps = true;
0118
0119
0120   //Create the world
0121   G4double world_hx = 1.*m;
0122   G4double world_hy = 1.*m;
0123   G4double world_hz = 1.*m;
0124   G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");
0125
0126   G4Box* solidWorld
0127     = new G4Box("World",
0128         world_hx,
0129         world_hy,
0130         world_hz);
0131
0132   fLogicWorld
0133     = new G4LogicalVolume(solidWorld,
0134               world_mat,
0135               "World");
0136
0137   G4VPhysicalVolume* physWorld
0138     = new G4PVPlacement(0,                     //no rotation
0139             G4ThreeVector(),       //at (0,0,0)
0140             fLogicWorld,            //its logical volume
0141             "World",               //its name
0142             0,                     //its mother  volume
0143             false,                 //no boolean operation
0144             0);                    //copy number
0145
0146
0147
0148   //////////////////////////////////////////////////////////////////////////////
0149   ///////////////////////////Create the detector////////////////////////////////
0150   //////////////////////////////////////////////////////////////////////////////
0151
0152
0153   //Overall parameters
0154   G4double startAngle     = 0.*deg;
0155   G4double spanningAngle  = 360.*deg;
0156
0157   ////////////////////////////////////
0158   /////////Define materials///////////
0159   ////////////////////////////////////
0160
0161
0162   //ALUMINIUM//
0163
0164   G4Material* al = nist->FindOrBuildMaterial("G4_Al");
0165
0166
0167
0168   //Create vacuum around the beam//
0169   G4double vacuum_atomic_number, vacuum_mass_of_mole, vacuum_density,vacuum_pressure,vacuum_temperature;
0170   vacuum_atomic_number = 1.;
0171   vacuum_mass_of_mole = 1.008*g/mole;
0172   vacuum_density      = 1.e-25*g/cm3;
0173   vacuum_pressure     = 1.e-8*bar;
0174   vacuum_temperature  = 293.*kelvin;     //from PhysicalConstants.h
0175   G4Material* vacuum_beam = new G4Material("vacuumBeam", vacuum_atomic_number,
0176                        vacuum_mass_of_mole, vacuum_density,
0177                        kStateGas,vacuum_temperature,
0178                        vacuum_pressure);
0179
0180
0181
0182   //HELIUM
0183   G4double helium_Z, helium_A, helium_density, helium_pressure, helium_temperature;
0184   helium_Z = 2.;
0185   helium_A = 4*g/mole;
0186   helium_density = 0.1785e-3*g/cm3; // with T=0°, 1 atm. To modify !
0187   helium_pressure = 2.*bar;
0188   helium_temperature = 293.*kelvin; //15 to 20°
0189   G4Material* helium = new G4Material("helium", helium_Z, helium_A, helium_density, kStateGas, helium_temperature, helium_pressure);
0190
0191
0192   //PLATINIUM
0193   G4Material* pt =  nist->FindOrBuildMaterial("G4_Pt");
0194
0195
0196   //TARGET MATERIAL INITIALIZATION
0197   /*G4String name;
0198   G4String symbole;
0199   G4int ncomponents;
0200   G4int n;
0201   G4double z_isotope;
0202   G4double abundance;
0203   G4double fractionmass;
0204   G4double a;*/
0205
0206
0207   /*
0208   //Pure Ni64
0209
0210   G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
0211
0212   G4Element* pureNi64 = new G4Element(name="pureNi64",symbole="64Ni", ncomponents=1);
0213   pureNi64->AddIsotope(Ni64, abundance = 100.*perCent);
0214
0215
0216   fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
0217   fTarget_Material->AddElement(pureNi64, fractionmass=1.);
0218  */
0219   /*
0220   //Ni64 94% enriched - Sz
0221
0222   G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
0223   G4Isotope* Ni58 = new G4Isotope(name="Zi58", z_isotope=28., n=58, a=58.*g/mole);
0224   G4Isotope* Ni60 = new G4Isotope(name="Zi60", z_isotope=28., n=60, a=60.*g/mole);
0225   G4Isotope* Ni61 = new G4Isotope(name="Zi61", z_isotope=28., n=61, a=61.*g/mole);
0226   G4Isotope* Ni62 = new G4Isotope(name="Zi62", z_isotope=28., n=62, a=62.*g/mole);
0227
0228
0229   G4Element* Ni64enriched95 = new G4Element(name="Ni64enriched95",symbole="64Ni_95", ncomponents=5);
0230   Ni64enriched95->AddIsotope(Ni64, abundance = 95.*perCent);
0231   Ni64enriched95->AddIsotope(Ni58, abundance = 2.6*perCent);
0232   Ni64enriched95->AddIsotope(Ni60, abundance = 1.72*perCent);
0233   Ni64enriched95->AddIsotope(Ni61, abundance = 0.15*perCent);
0234   Ni64enriched95->AddIsotope(Ni62, abundance = 0.53*perCent);
0235
0236
0237   fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
0238   fTarget_Material->AddElement(Ni64enriched95, fractionmass=1.);
0239
0240   //fTarget_Material =  nist->FindOrBuildMaterial("G4_Y");
0241   */
0242
0243   /*
0244
0245 //Ni64 95% enriched - Obata
0246
0247   G4Isotope* Ni64 = new G4Isotope(name="Zi64", z_isotope=28., n=64, a=64.*g/mole);
0248   G4Isotope* Ni58 = new G4Isotope(name="Zi58", z_isotope=28., n=58, a=58.*g/mole);
0249   G4Isotope* Ni60 = new G4Isotope(name="Zi60", z_isotope=28., n=60, a=60.*g/mole);
0250   G4Isotope* Ni61 = new G4Isotope(name="Zi61", z_isotope=28., n=61, a=61.*g/mole);
0251   G4Isotope* Ni62 = new G4Isotope(name="Zi62", z_isotope=28., n=62, a=62.*g/mole);
0252
0253
0254   G4Element* Ni64enriched95 = new G4Element(name="Ni64enriched95",symbole="64Ni_95", ncomponents=5);
0255   Ni64enriched95->AddIsotope(Ni64, abundance = 94.8*perCent);
0256   Ni64enriched95->AddIsotope(Ni58, abundance = 2.67*perCent);
0257   Ni64enriched95->AddIsotope(Ni60, abundance = 1.75*perCent);
0258   Ni64enriched95->AddIsotope(Ni61, abundance = 0.11*perCent);
0259   Ni64enriched95->AddIsotope(Ni62, abundance = 0.67*perCent);
0260
0261
0262   fTarget_Material = new G4Material("FTarget_Material",fDensity_target,ncomponents=1);
0263   fTarget_Material->AddElement(Ni64enriched95, fractionmass=1.);
0264   */
0265
0266   fTarget_Material =  nist->FindOrBuildMaterial("G4_Ni");
0267
0268   //FOIL MATERIAL
0269   fFoil_Material =  nist->FindOrBuildMaterial("G4_Al");
0270
0271
0272   ////////////////////////////////////////////////////////////////////
0273   /////////////////////////////PART1//////////////////////////////////
0274   ////////////////////////////////////////////////////////////////////
0275
0276   ////////////////////////////////////////////////////////////////////
0277   //////////////////////////LAYER PART 1//////////////////////////////
0278   ////////////////////////////////////////////////////////////////////
0279
0280
0281   //Create the (external) layer around the beam vacuum tube
0282
0283   G4double layer1_length_PART1         = 98.9*mm;
0284   G4double layer1_innerRadius_PART1    = 11.5*mm;
0285   G4double layer1_outerRadius_PART1    = 16.*mm;
0286   G4double layer1_hz_PART1             = 0.5*layer1_length_PART1;
0287
0288   G4Tubs* solidLayer1_PART1
0289     = new G4Tubs("Layer1_PART1",
0290          layer1_innerRadius_PART1,
0291          layer1_outerRadius_PART1,
0292          layer1_hz_PART1,
0293          startAngle,
0294          spanningAngle);
0295
0296
0297   G4double layer2_length_PART1         = 124.6*mm;
0298   G4double layer2_innerRadius_PART1    = 7.5*mm;
0299   G4double layer2_outerRadius_PART1    = 11.5*mm;
0300   G4double layer2_hz_PART1             = 0.5*layer2_length_PART1;
0301
0302   G4Tubs* solidLayer2_PART1
0303     = new G4Tubs("Layer2_PART1",
0304          layer2_innerRadius_PART1,
0305          layer2_outerRadius_PART1,
0306          layer2_hz_PART1,
0307          startAngle,
0308          spanningAngle);
0309
0310
0311
0312   G4RotationMatrix rot_layer_PART1;
0313   G4double z_layer_translation_PART1 = layer2_length_PART1-layer1_length_PART1;
0314   G4ThreeVector placement_layer_PART1 = G4ThreeVector(0.*mm,0.*mm,0.5*z_layer_translation_PART1);
0315   G4Transform3D transform_layer_PART1(rot_layer_PART1,placement_layer_PART1);
0316
0317   G4UnionSolid* solidLayer_PART1 = new G4UnionSolid("Layer_PART1", solidLayer2_PART1, solidLayer1_PART1, transform_layer_PART1);
0318
0319
0320   G4LogicalVolume* logicLayer_PART1
0321     = new G4LogicalVolume(solidLayer_PART1,
0322               al,
0323               "Layer_PART1");
0324
0325   /*  G4VPhysicalVolume* physLayer1_PART1 = */
0326   new G4PVPlacement(0,
0327             G4ThreeVector(0.*mm,0.*mm,0.5*layer2_length_PART1),
0328             logicLayer_PART1,
0329             "Layer_PART1",
0330             fLogicWorld,
0331             false,
0332             0,
0333             checkOverlaps);
0334
0335
0336   //////////////////////////////////////////////////////////////////
0337   ////////////////Create the beam vacuum tube///////////////////////
0338   //////////////////////////////////////////////////////////////////
0339
0340
0341   G4double tube_length_PART1 = layer2_length_PART1;
0342
0343   G4double tube_innerRadius_PART1    = 0.*mm;
0344   G4double tube_outerRadius_PART1    = 7.5*mm;
0345   G4double tube_hz_PART1             = 0.5*tube_length_PART1;
0346
0347   G4Tubs* solidTube_PART1
0348     = new G4Tubs("Tube_PART1",
0349          tube_innerRadius_PART1,
0350          tube_outerRadius_PART1,
0351          tube_hz_PART1,
0352          startAngle,
0353          spanningAngle);
0354
0355   G4LogicalVolume* logicTube_PART1
0356     = new G4LogicalVolume(solidTube_PART1,
0357               vacuum_beam,
0358               "Tube_PART1");
0359
0360   /*  G4VPhysicalVolume* physTube_PART1 = */
0361   new G4PVPlacement(0,
0362             G4ThreeVector(0.*mm,0.*mm,0.5*tube_length_PART1),
0363             logicTube_PART1,
0364             "Tube_PART1",
0365             fLogicWorld,
0366             false,
0367             0,
0368             checkOverlaps);
0369
0370
0371
0372
0373
0374
0375   ////////////////////////////////////////////////////////////////////////////////
0376   /////////////////////////// PART 2 = FOIL PART /////////////////////////////////
0377   ////////////////////////////////////////////////////////////////////////////////
0378
0379   ////////////////////////////////////////////////////////////////////
0380   //////////////////////////LAYER PART 2//////////////////////////////
0381   ////////////////////////////////////////////////////////////////////
0382
0383   //Create the (external) layer around the foil part
0384
0385   G4double layer_length_PART2          = 12.*mm;
0386   G4double layer_innerRadius_PART2     = 7.5*mm;
0387   G4double layer_outerRadius_PART2     = 15.*mm;
0388   G4double layer_hz_PART2              = 0.5*layer_length_PART2;
0389
0390   G4Tubs* solidLayer_PART2
0391     = new G4Tubs("Layer_PART2",
0392          layer_innerRadius_PART2,
0393          layer_outerRadius_PART2,
0394          layer_hz_PART2,
0395          startAngle,
0396          spanningAngle);
0397
0398
0399    G4LogicalVolume* logicLayer_PART2
0400     = new G4LogicalVolume(solidLayer_PART2,
0401               al,
0402               "Layer_PART2");
0403
0404    G4double layer_z_position_PART2 = tube_length_PART1 + 0.5*layer_length_PART2;
0405
0406    /*   G4VPhysicalVolume* physLayer1_PART2 = */
0407    new G4PVPlacement(0,
0408              G4ThreeVector(0.*mm,0.*mm,layer_z_position_PART2),
0409              logicLayer_PART2,
0410              "Layer_PART2",
0411              fLogicWorld,
0412              false,
0413              0,
0414              checkOverlaps);
0415
0416
0417
0418   ////////////////////////////////////////////////////////////////////
0419   //////////////////////////TUBE PART 2///////////////////////////////
0420   ////////////////////////////////////////////////////////////////////
0421
0422   //Create the tube before the foil
0423
0424   G4double tube_length_PART2 = layer_length_PART2;
0425
0426   G4double tube_innerRadius_PART2     = 0.*mm;
0427   G4double tube_outerRadius_PART2     = 7.5*mm;
0428   G4double tube_hz_PART2              = 0.5*tube_length_PART2;
0429
0430   G4Tubs* solidTube_PART2
0431     = new G4Tubs("Tube_PART2",
0432          tube_innerRadius_PART2,
0433          tube_outerRadius_PART2,
0434          tube_hz_PART2,
0435          startAngle,
0436          spanningAngle);
0437
0438   G4LogicalVolume* logicTube_PART2
0439    = new G4LogicalVolume(solidTube_PART2,
0440              vacuum_beam,
0441              "Tube_PART2");
0442
0443
0444   /*  G4VPhysicalVolume* physTube_PART2 = */
0445   new G4PVPlacement(0,
0446             G4ThreeVector(0.*mm,0.*mm, layer_z_position_PART2),
0447             logicTube_PART2,
0448             "Tube_PART2",
0449             fLogicWorld,
0450             false,
0451             0,
0452             checkOverlaps);
0453
0454
0455   /////////////////////////////////////////////////////////
0456   //////////////////////// GRID PART //////////////////////
0457   /////////////////////////////////////////////////////////
0458
0459
0460   //HEXAGONE
0461   G4double hexagone_length = layer_length_PART2 ;
0462   G4int numSide = 6;
0463   G4int numZPlanes = 2;
0464   const G4double zPlane[]={-0.5*hexagone_length,0.5*hexagone_length};
0465   const G4double rInner[]={3.*mm,3.*mm};
0466   const G4double rOuter[]={3.3*mm,3.3*mm};
0467
0468   G4Polyhedra* solidHexagone_0
0469     = new G4Polyhedra("Grid",
0470               0.*deg,
0471               360.*deg,
0472               numSide,
0473               numZPlanes,
0474               zPlane,
0475               rInner,
0476               rOuter);
0477
0478
0479   G4LogicalVolume* logicHexagone
0480     = new G4LogicalVolume(solidHexagone_0,
0481               al,
0482               "Grid");
0483
0484   /*  G4VPhysicalVolume* physHexagone = */
0485   new G4PVPlacement(0,
0486            G4ThreeVector(0.*mm,0.*mm,0.*mm),
0487            logicHexagone,
0488            "Grid",
0489            logicTube_PART2,
0490            false,
0491            0,
0492            checkOverlaps);
0493
0494
0495
0496
0497
0498   /////////////////////////////////////////////////////////////////////////
0499   //////////////////////////CREATE THE FOIL ///////////////////////////////
0500   /////////////////////////////////////////////////////////////////////////
0501
0502
0503   fZ_foil_position = 0.5*fFoil_thickness + tube_length_PART1 + layer_length_PART2;
0504
0505
0506
0507   G4double foil_innerRadius    = 0.*mm;
0508   G4double foil_outerRadius    = 16.*mm;
0509
0510   fSolidFoil
0511     = new G4Tubs("Foil",
0512          foil_innerRadius,
0513          foil_outerRadius,
0514          0.5*fFoil_thickness,
0515          startAngle,
0516          spanningAngle);
0517
0518   fFoilVolume = fSolidFoil->GetCubicVolume();
0519
0520   fLogicFoil
0521     = new G4LogicalVolume(fSolidFoil,
0522               fFoil_Material,
0523               "Foil");
0524
0525
0526
0527   fPhysFoil
0528     = new G4PVPlacement(0,
0529             G4ThreeVector(0.*mm,0.*mm,fZ_foil_position),
0530             fLogicFoil,
0531             "Foil",
0532             fLogicWorld,
0533             false,
0534             0,
0535             checkOverlaps);
0536
0537
0538
0539
0540
0541   /////////////////////////////////////////////////////////////////////////////////////
0542   /////////////////////////// PART 3 = AFTER THE FOIL /////////////////////////////////
0543   /////////////////////////////////////////////////////////////////////////////////////
0544
0545
0546   ////////////////////////////////////////////////////////////////////
0547   //////////////////////////LAYER PART 3//////////////////////////////
0548   ////////////////////////////////////////////////////////////////////
0549
0550   //Create the (external) layer around the beam
0551
0552
0553   G4double layer_length_PART3 = 38.32*mm;
0554
0555   G4double layer_innerRadius_PART3     = 7.5*mm;
0556   G4double layer_outerRadius_PART3     = 16.*mm;
0557   G4double layer_hz_PART3              = 0.5*layer_length_PART3;
0558
0559   G4Tubs* solidLayer_PART3
0560     = new G4Tubs("Layer_PART3",
0561          layer_innerRadius_PART3,
0562          layer_outerRadius_PART3,
0563          layer_hz_PART3,
0564          startAngle,
0565          spanningAngle);
0566
0567   G4LogicalVolume* logicLayer_PART3
0568     = new G4LogicalVolume(solidLayer_PART3,
0569               al,
0570               "Layer_PART3");
0571
0572   fLayer_z_position_PART3 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + 0.5*layer_length_PART3;
0573
0574   fPhysLayer_PART3
0575     = new G4PVPlacement(0,
0576             G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3),
0577             logicLayer_PART3,
0578             "Layer_PART3",
0579             fLogicWorld,
0580             false,
0581             0,
0582             checkOverlaps);
0583
0584
0585
0586   ////////////////////////////////////////////////////////////////////
0587   //////////////////////////TUBE PART 3///////////////////////////////
0588   ////////////////////////////////////////////////////////////////////
0589
0590   //Create the tube after the foil, MATERIAL = HELIUM
0591
0592   G4double tube_length_PART3 = layer_length_PART3;
0593
0594   G4double tube_innerRadius_PART3     = 0.*mm;
0595   G4double tube_outerRadius_PART3     = 7.5*mm;
0596   G4double tube_hz_PART3              = 0.5*tube_length_PART3;
0597
0598   G4Tubs* solidTube_PART3
0599     = new G4Tubs("Tube_PART3",
0600          tube_innerRadius_PART3,
0601          tube_outerRadius_PART3,
0602          tube_hz_PART3,
0603          startAngle,
0604          spanningAngle);
0605
0606   G4LogicalVolume* logicTube_PART3
0607     = new G4LogicalVolume(solidTube_PART3,
0608               helium,
0609               "Tube_PART3");
0610
0611
0612   fPhysTube_PART3
0613     = new G4PVPlacement(0,
0614             G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART3),
0615             logicTube_PART3,
0616             "Tube_PART3",
0617             fLogicWorld,
0618             false,
0619             0,
0620             checkOverlaps);
0621
0622
0623   ////////////////////////////////////////////////////////////////////////////////////////
0624   /////////////////////////// PART 4 = BEFORE THE TARGET/////////////////////////////////
0625   ////////////////////////////////////////////////////////////////////////////////////////
0626
0627   ////////////////////////////////////////////////////////////////////
0628   //////////////////////////LAYER PART4///////////////////////////////
0629   ////////////////////////////////////////////////////////////////////
0630
0631
0632   G4double layer1_length_PART4    = 11.5*mm; //Length of the gold part
0633
0634   G4double layer2_length_PART4    = 4.6*mm;
0635   G4double layer2_Rmin1_PART4     = 8.5*mm;
0636   G4double layer2_Rmax1_PART4     = 9.*mm;
0637   G4double layer2_Rmin2_PART4     = 8.5*mm;
0638   G4double layer2_Rmax2_PART4     = 14.*mm;
0639   G4double layer2_hz_PART4        = 0.5*layer2_length_PART4;
0640
0641   G4Cons* solidLayer2_PART4
0642     = new G4Cons("Layer2_PART4",
0643          layer2_Rmin1_PART4,
0644          layer2_Rmax1_PART4,
0645          layer2_Rmin2_PART4,
0646          layer2_Rmax2_PART4,
0647          layer2_hz_PART4,
0648          startAngle,
0649          spanningAngle);
0650
0651   //Create the (external) aluminium layer around the beam before the target
0652
0653   G4double layer3_length_PART4 = layer1_length_PART4-layer2_length_PART4;
0654   G4double layer3_innerRadius_PART4     = 8.5*mm;
0655   G4double layer3_outerRadius_PART4     = 14.*mm;
0656   G4double layer3_hz_PART4              = 0.5*layer3_length_PART4;
0657
0658   G4Tubs* solidLayer3_PART4
0659     = new G4Tubs("Layer3_PART4",
0660          layer3_innerRadius_PART4,
0661          layer3_outerRadius_PART4,
0662          layer3_hz_PART4,
0663          startAngle,
0664          spanningAngle);
0665
0666
0667
0668   G4RotationMatrix rot_layer_PART4;
0669   G4double z_layer_translation_PART4 = 0.5*(layer2_length_PART4+layer3_length_PART4);
0670   G4ThreeVector placement_layer_PART4 = G4ThreeVector(0.*mm,0.*mm,z_layer_translation_PART4);
0671   G4Transform3D transform_layer_PART4(rot_layer_PART4,placement_layer_PART4);
0672
0673   G4UnionSolid* solidLayer_PART4 = new G4UnionSolid("Layer_PART1", solidLayer2_PART4, solidLayer3_PART4, transform_layer_PART4);
0674
0675
0676
0677   G4LogicalVolume* logicLayer_PART4
0678     = new G4LogicalVolume(solidLayer_PART4,
0679               al,
0680               "Layer_PART4");
0681
0682   fLayer_z_position_PART4 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +0.5*layer2_length_PART4;
0683
0684   fPhysLayer_PART4
0685     = new G4PVPlacement(0,
0686             G4ThreeVector(0.*mm,0.*mm,fLayer_z_position_PART4),
0687             logicLayer_PART4,
0688             "Layer_PART4",
0689             fLogicWorld,
0690             false,
0691             0,
0692             checkOverlaps);
0693
0694
0695
0696
0697
0698    //Create the (internal) gold layer around the beam
0699
0700   G4double layer1_innerRadius_PART4     = 8.*mm;
0701   G4double layer1_outerRadius_PART4     = 8.5*mm;
0702   G4double layer1_hz_PART4              = 0.5*layer1_length_PART4;
0703
0704   G4Tubs* solidLayer1_PART4
0705     = new G4Tubs("Layer1_PART4",
0706          layer1_innerRadius_PART4,
0707          layer1_outerRadius_PART4,
0708          layer1_hz_PART4,
0709          startAngle,
0710          spanningAngle);
0711
0712
0713   G4LogicalVolume* logicLayer1_PART4
0714     = new G4LogicalVolume(solidLayer1_PART4,
0715               pt,
0716               "Layer1_PART4");
0717
0718   fLayer1_z_position_PART4 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +0.5*layer1_length_PART4;
0719
0720
0721   fPhysLayer1_PART4
0722     = new G4PVPlacement(0,
0723             G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4),
0724             logicLayer1_PART4,
0725             "Layer1_PART4",
0726             fLogicWorld,
0727             false,
0728             0,
0729             checkOverlaps);
0730
0731   ////////////////////////////////////////////////////////////////////
0732   //////////////////////////TUBE PART 4///////////////////////////////
0733   ////////////////////////////////////////////////////////////////////
0734
0735
0736   //Create the tube before the target
0737
0738
0739   fTube_length_PART4 = layer1_length_PART4;
0740
0741   G4double tube_innerRadius_PART4     = 0.*mm;
0742   fTube_outerRadius_PART4     = 8.*mm;
0743   G4double tube_hz_PART4              = 0.5*fTube_length_PART4;
0744
0745   G4Tubs* solidTube_PART4
0746     = new G4Tubs("Tube_PART4",
0747          tube_innerRadius_PART4,
0748          fTube_outerRadius_PART4,
0749          tube_hz_PART4,
0750          startAngle,
0751          spanningAngle);
0752
0753   G4LogicalVolume* logicTube_PART4
0754     = new G4LogicalVolume(solidTube_PART4,
0755               helium,
0756               "Tube_PART4");
0757
0758
0759   fPhysTube_PART4
0760     = new G4PVPlacement(0,
0761             G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART4),
0762             logicTube_PART4,
0763             "Tube_PART4",
0764             fLogicWorld,
0765             false,
0766             0,
0767             checkOverlaps);
0768
0769
0770
0771
0772   ////////////////////////////////////////////////////////////////////
0773   //////////////////////////// TARGET ////////////////////////////////
0774   ////////////////////////////////////////////////////////////////////
0775
0776
0777
0778
0779   //Design the target inside the tube
0780
0781   G4double target_innerRadius     = 0.*mm;
0782   G4double target_outerRadius     = 0.5*fTarget_diameter;
0783   G4double target_hz              = 0.5*fTarget_thickness;
0784
0785   fSolidTarget
0786     = new G4Tubs("Target",
0787          target_innerRadius,
0788          target_outerRadius,
0789          target_hz,
0790          startAngle,
0791          spanningAngle);
0792
0793   fTargetVolume = fSolidTarget->GetCubicVolume();
0794
0795   fLogicTarget
0796     = new G4LogicalVolume(fSolidTarget,
0797               fTarget_Material,
0798               "Target");
0799
0800   fTarget_z_position = 0.5*fTube_length_PART4-0.5*fTarget_thickness;  //inside the tube!
0801
0802
0803   fPhysTarget
0804     = new G4PVPlacement(0,
0805             G4ThreeVector(0.*mm,0.*mm, fTarget_z_position),
0806             fLogicTarget,
0807             "Target",
0808             logicTube_PART4,
0809             false,
0810             0,
0811             checkOverlaps);
0812
0813
0814
0815   ///////////////////////////////////////////////////////////////////////////////////////
0816   /////////////////////////// PART 5 = AFTER THE TARGET /////////////////////////////////
0817   ///////////////////////////////////////////////////////////////////////////////////////
0818
0819   ////////////////////////////////////////////////////////////////////
0820   //////////////////////////LAYER PART 5//////////////////////////////
0821   ////////////////////////////////////////////////////////////////////
0822
0823   //Create the (internal) platinium layer
0824
0825   G4double layer1_length_PART5 = 0.5*mm;
0826
0827   G4double layer1_innerRadius_PART5     = 0.*mm;
0828   G4double layer1_outerRadius_PART5     = 8.5*mm;
0829   G4double layer1_hz_PART5              = 0.5*layer1_length_PART5;
0830
0831   G4Tubs* solidLayer1_PART5
0832     = new G4Tubs("Layer1_PART5",
0833          layer1_innerRadius_PART5,
0834          layer1_outerRadius_PART5,
0835          layer1_hz_PART5,
0836          startAngle,
0837          spanningAngle);
0838
0839   G4LogicalVolume* logicLayer1_PART5
0840     = new G4LogicalVolume(solidLayer1_PART5,
0841               pt,
0842               "Layer1_PART5");
0843
0844   fLayer1_z_position_PART5 = tube_length_PART1 + layer_length_PART2 + fFoil_thickness + layer_length_PART3 +layer1_length_PART4 + 0.5*layer1_length_PART5;
0845
0846   fPhysLayer1_PART5
0847     = new G4PVPlacement(0,
0848             G4ThreeVector(0.*mm,0.*mm,fLayer1_z_position_PART5),
0849             logicLayer1_PART5,
0850             "Layer1_PART5",
0851             fLogicWorld,
0852             false,
0853             0,
0854             checkOverlaps);
0855
0856
0857
0858   //Create the (external) aluminium layer after the target
0859
0860   G4double layer2_length_PART5 = 0.5*mm;
0861
0862   G4double layer2_innerRadius_PART5     = 8.5*mm;
0863   G4double layer2_outerRadius_PART5     = 14.*mm;
0864   G4double layer2_hz_PART5              = 0.5*layer2_length_PART5;
0865
0866   G4Tubs* solidLayer2_PART5
0867     = new G4Tubs("Layer2_PART5",
0868          layer2_innerRadius_PART5,
0869          layer2_outerRadius_PART5,
0870          layer2_hz_PART5,
0871          startAngle,
0872          spanningAngle);
0873
0874   G4LogicalVolume* logicLayer2_PART5
0875     = new G4LogicalVolume(solidLayer2_PART5,
0876               al,
0877               "Layer2_PART5");
0878
0879   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;
0880
0881   fPhysLayer2_PART5
0882     = new G4PVPlacement(0,
0883             G4ThreeVector(0.*mm,0.*mm,fLayer2_z_position_PART5),
0884             logicLayer2_PART5,
0885             "Layer2_PART5",
0886             fLogicWorld,
0887             false,
0888             0,
0889             checkOverlaps);
0890
0891
0892
0893   //Create the end of the beam target machine
0894
0895   G4double layer3_length_PART5 = 3.*mm;
0896
0897   G4double layer3_innerRadius_PART5     = 0.*mm;
0898   G4double layer3_outerRadius_PART5     = 14.*mm;
0899   G4double layer3_hz_PART5              = 0.5*layer3_length_PART5;
0900
0901   G4Tubs* solidLayer3_PART5
0902     = new G4Tubs("Layer3_PART5",
0903          layer3_innerRadius_PART5,
0904          layer3_outerRadius_PART5,
0905          layer3_hz_PART5,
0906          startAngle,
0907          spanningAngle);
0908
0909   G4LogicalVolume* logicLayer3_PART5
0910     = new G4LogicalVolume(solidLayer3_PART5,
0911               al,
0912               "Layer3_PART5");
0913
0914   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;
0915
0916   fPhysLayer3_PART5
0917     = new G4PVPlacement(0,
0918             G4ThreeVector(0.*mm,0.*mm,fLayer3_z_position_PART5),
0919             logicLayer3_PART5,
0920             "Layer3_PART5",
0921             fLogicWorld,
0922             false,
0923             0,
0924             checkOverlaps);
0925
0926
0927
0928
0929   //////////////////////////////////////////////
0930   /////////// Set sensitive region /////////////
0931   //////////////////////////////////////////////
0932
0933
0934   if(!fRegionTarget)
0935     {
0936       fRegionTarget = new G4Region("Target");
0937       fLogicTarget -> SetRegion(fRegionTarget);
0938       fRegionTarget -> AddRootLogicalVolume(fLogicTarget);
0939     }
0940
0941   if(!fRegionFoil&&fPhysFoil!=nullptr)
0942     {
0943       fRegionFoil = new G4Region("Foil");
0944       fLogicFoil -> SetRegion(fRegionFoil);
0945       fRegionFoil -> AddRootLogicalVolume(fLogicFoil);
0946     }
0947
0948
0949
0950
0951   //Always return physical world//
0952
0953   return physWorld;
0954
0955
0956 }
0957
0958 void STCyclotronDetectorConstruction::ConstructSDandField()
0959 {
0960
0961   if(fLogicTarget != nullptr){
0962     STCyclotronSensitiveTarget* TargetSD = new STCyclotronSensitiveTarget("TargetSD", this);
0963     SetSensitiveDetector("Target",TargetSD);
0964   }
0965
0966
0967   if(fLogicFoil != nullptr){
0968     STCyclotronSensitiveFoil* FoilSD = new STCyclotronSensitiveFoil("FoilSD", this);
0969     SetSensitiveDetector("Foil",FoilSD);
0970   }
0971 }
0972
0973 void STCyclotronDetectorConstruction::SetTargetDiameter(G4double targetDiameter)
0974 {
0975
0976   if(fTarget_diameter!=targetDiameter){
0977     if(targetDiameter/2.>fTube_outerRadius_PART4){ G4cout << "Error : the diameter is bigger than the tube" << G4endl;}
0978     else{
0979       fTarget_diameter = targetDiameter;
0980       if(fSolidTarget) fSolidTarget->SetOuterRadius(0.5*fTarget_diameter);
0981       G4cout << "The new diameter of the target is " << fTarget_diameter << " mm." << G4endl;
0982     }
0983   }
0984
0985   G4RunManager::GetRunManager() -> GeometryHasBeenModified();
0986   G4cout << "... Geometry is notified .... " << G4endl;
0987
0988 }
0989
0990 //SET MATERIAL METHODS//
0991
0992 void STCyclotronDetectorConstruction::SetTargetIsotopeName(G4String name){
0993   fIsotopeName.push_back(name);
0994   G4RunManager::GetRunManager() -> GeometryHasBeenModified();
0995   G4cout << "... Geometry is notified .... " << G4endl;
0996 }
0997
0998 void STCyclotronDetectorConstruction::SetTargetIsotopeZ(G4double z){
0999   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