Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:20:21 

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 // Code developed by
0027 // Silvia Pozzi (1), silvia.pozzi@iss.it
0028 // Barbara Caccia (1), barbara.caccia@iss.it
0029 // Carlo Mancini Terracciano (2), carlo.mancini.terracciano@roma1.infn.it
0030 // (1) Istituto Superiore di Sanita' and INFN Roma, Italy
0031 // (2) Univ. La Sapienza and INFN Roma, Italy
0032 
0033 #include "DetectorConstruction.hh"
0034 #include "DetectorMessenger.hh"
0035 
0036 #include <G4LogicalVolume.hh>
0037 #include <G4PVPlacement.hh>
0038 #include <G4PVReplica.hh>
0039 #include <G4ProductionCuts.hh>
0040 #include <G4NistManager.hh>
0041 #include <G4SystemOfUnits.hh>
0042 #include <G4VisAttributes.hh>
0043 #include <G4Box.hh>
0044 #include <G4VSolid.hh>
0045 #include <G4Cons.hh>
0046 #include <G4Tubs.hh>
0047 #include <G4BooleanSolid.hh>
0048 #include <G4SubtractionSolid.hh>
0049 #include <G4SDManager.hh>
0050 #include <G4GlobalMagFieldMessenger.hh>
0051 #include <G4UnionSolid.hh>
0052 #include <G4SubtractionSolid.hh>
0053 #include <G4RunManager.hh>
0054 #include <G4GeometryManager.hh>
0055 #include <G4PhysicalVolumeStore.hh>
0056 #include <G4LogicalVolumeStore.hh>
0057 #include <G4SolidStore.hh>
0058 #include <G4Transform3D.hh>
0059 #include <CLHEP/Vector/Rotation.h>
0060 
0061 
0062 G4VPhysicalVolume* DetectorConstruction::Construct()
0063 {
0064     detectorMessenger = new DetectorMessenger(this);
0065 
0066     //// Saturn 43
0067     G4Material* air = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
0068 
0069     G4double worldSizeX = 3.5 * m;
0070     G4double worldSizeY = 3.5 * m;
0071     G4double worldSizeZ = 3.5 * m;
0072 
0073     G4VSolid* worldBox = new G4Box("world", worldSizeX/2., worldSizeY/2., worldSizeZ/2.);
0074 
0075     G4LogicalVolume* worldLog = new G4LogicalVolume(worldBox, air, "world");
0076 
0077     G4VisAttributes* visAttr = new G4VisAttributes();
0078     visAttr->SetVisibility(false);
0079     worldLog->SetVisAttributes(visAttr);
0080 
0081     world_phys = new G4PVPlacement(nullptr, {}, worldLog, "world", nullptr, false, 0);
0082 
0083     ConstructMaterials();
0084     ConstructAccelerator();
0085     ConstructPhantom();
0086 
0087     return world_phys;
0088 }
0089 
0090 void DetectorConstruction::ConstructPhantom()
0091 {
0092     G4Material* water = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
0093 
0094     G4ThreeVector phantomHalfDim = {phantomSideDim/2., phantomSideDim/2., phantomSideDim/2.};
0095 
0096     G4VSolid* Phantom = new G4Box("Phantom", phantomHalfDim.getX(), phantomHalfDim.getY(), phantomHalfDim.getZ());
0097     G4LogicalVolume* Phantom_log = new G4LogicalVolume(Phantom, water, "Phantom_log", 0, 0, 0);
0098 
0099     G4VSolid* spessore = new G4Box("Spessore", phantomHalfDim.getX(), phantomHalfDim.getY(), 1.25*mm);
0100     G4LogicalVolume* spessore_log = new G4LogicalVolume(spessore, water, "spessore_log", 0, 0, 0);
0101     new G4PVPlacement(nullptr, {0.,0., 1.25*mm}, "spessore_phys", spessore_log,  world_phys, false, 1);
0102 
0103     // axes origin on the face of the phantom
0104     G4ThreeVector phantomPos = {0., 0., phantomHalfDim.getZ()+voxelDepthDim/2};
0105 
0106     phantom_phys = new G4PVPlacement(nullptr, phantomPos,
0107                  "phantom_phys", Phantom_log, world_phys, false, 1);
0108 
0109     PhantomSegmentation(Phantom_log);
0110 
0111     // Region for cuts
0112     G4Region *regVol= new G4Region("fullWaterPhantomR");
0113     G4ProductionCuts* cuts = new G4ProductionCuts;
0114     cuts->SetProductionCut(0.1*mm);
0115     regVol->SetProductionCuts(cuts);
0116 
0117     Phantom_log -> SetRegion(regVol);
0118     regVol->AddRootLogicalVolume(Phantom_log);
0119 
0120     //--------- Visualization attributes -------------------------------
0121     /// Phantom
0122     G4VisAttributes* simpleH2OVisAtt= new G4VisAttributes(G4Colour::Cyan());
0123     simpleH2OVisAtt->SetVisibility(true);
0124     simpleH2OVisAtt->SetForceSolid(false);
0125     Phantom_log->SetVisAttributes(simpleH2OVisAtt);
0126 
0127 }
0128 
0129 G4Material* DetectorConstruction::GetMaterial(const G4String materialName)
0130 {
0131     G4Material* material = nullptr;
0132 
0133         if (materialName == "kapton")
0134         {
0135             material = mat_Kapton;
0136         }
0137         else if (materialName == "steel")
0138         {
0139             material = mat_Ssteel;
0140         }
0141         else if (materialName == "xc10")
0142         {
0143             material = mat_XC10;
0144         }
0145         else if (materialName == "wnicu")
0146         {
0147             material = mat_WNICU;
0148         }
0149         else
0150         {
0151             G4cerr << "* DetectorConstruction::GetMaterial: material " <<
0152                     materialName << " not found" << G4endl;
0153         }
0154 
0155         return material;
0156 }
0157 
0158 void DetectorConstruction::ConstructMaterials()
0159 {
0160     G4NistManager* nist = G4NistManager::Instance();
0161 
0162     G4Material* el_H  = nist -> FindOrBuildMaterial("G4_H");
0163     G4Material* el_C  = nist -> FindOrBuildMaterial("G4_C");
0164     G4Material* el_N  = nist -> FindOrBuildMaterial("G4_N");
0165     G4Material* el_O  = nist -> FindOrBuildMaterial("G4_O");
0166     G4Material* el_Si = nist -> FindOrBuildMaterial("G4_Si");
0167     G4Material* el_Cr = nist -> FindOrBuildMaterial("G4_Cr");
0168     G4Material* el_Mn = nist -> FindOrBuildMaterial("G4_Mn");
0169     G4Material* el_Fe = nist -> FindOrBuildMaterial("G4_Fe");
0170     G4Material* el_Ni = nist -> FindOrBuildMaterial("G4_Ni");
0171     G4Material* el_Cu = nist -> FindOrBuildMaterial("G4_Cu");
0172     G4Material* el_W  = nist -> FindOrBuildMaterial("G4_W");
0173 
0174     mat_Ssteel = new G4Material("StainlessSteel", 7.8 *g/cm3, 6);
0175     mat_Ssteel -> AddMaterial(el_Fe, 0.6898);
0176     mat_Ssteel -> AddMaterial(el_Cr, 0.18);
0177     mat_Ssteel -> AddMaterial(el_Ni, 0.10);
0178     mat_Ssteel -> AddMaterial(el_Mn, 0.02);
0179     mat_Ssteel -> AddMaterial(el_Si, 0.01);
0180     mat_Ssteel -> AddMaterial(el_C,  0.0002);
0181 
0182     mat_XC10 = new G4Material("CARBON_STEEL", 7.8 *g/cm3, 3);
0183     mat_XC10 -> AddMaterial(el_Fe, 0.993);
0184     mat_XC10 -> AddMaterial(el_Mn, 0.006);
0185     mat_XC10 -> AddMaterial(el_C,  0.001);
0186 
0187     mat_WNICU = new G4Material("Denal(WNICU)", 16.8*g/cm3, 3);
0188     mat_WNICU -> AddMaterial(el_W,  0.905);
0189     mat_WNICU -> AddMaterial(el_Ni, 0.07);
0190     mat_WNICU -> AddMaterial(el_Cu, 0.025);
0191 
0192     mat_Kapton = new G4Material("Kapton", 1.42*g/cm3, 4);
0193     mat_Kapton -> AddMaterial(el_C, 69.1133 *perCent);
0194     mat_Kapton -> AddMaterial(el_O, 20.9235 *perCent);
0195     mat_Kapton -> AddMaterial(el_N, 7.3270 *perCent);
0196     mat_Kapton -> AddMaterial(el_H, 2.6362 *perCent);
0197 
0198 }
0199 
0200 void DetectorConstruction::ConstructAccelerator()
0201 {
0202     G4bool checkOverlap = true;
0203 
0204     G4Material* Vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_Galactic");
0205 
0206     accHalfSize = {600.*mm, 600.*mm, 600.*mm};
0207     G4ThreeVector initialCentre = {0.*mm, 0.*mm, -std::abs(sourceToSkinDistance)};
0208 
0209     G4Box* accWorldBox = new G4Box("accWorld", accHalfSize.getX(), accHalfSize.getY(), accHalfSize.getZ());
0210     G4LogicalVolume* accWorldLV = new G4LogicalVolume(accWorldBox, Vacuum, "accWorldL", 0, 0, 0);
0211     accWorld_phys = new G4PVPlacement(nullptr, initialCentre, "acceleratorBox", accWorldLV, world_phys, false, checkOverlap);
0212 
0213     G4VisAttributes* simpleAlSVisAtt = new G4VisAttributes(G4Colour::White());
0214     simpleAlSVisAtt -> SetVisibility(false);
0215     accWorldLV -> SetVisAttributes(simpleAlSVisAtt);
0216 
0217     ConstructTarget();
0218     ConstructPrimaryCollimator();
0219     ConstructVacuumWindow();
0220     ConstructFlatteningFilter();
0221     ConstructIonizationChamber();
0222     ConstructJawsX();
0223     ConstructJawsY();
0224 }
0225 
0226 void DetectorConstruction::ConstructTarget()
0227 {
0228     G4bool checkOverlap = true;
0229 
0230     G4ThreeVector targetCentre, boxHalfSide, tubeCentre;
0231     G4double tubeRadius, tubeHalfLengthZ;
0232     targetCentre.set(0, 0, -3.5*mm);
0233     boxHalfSide.set(5.*mm, 5.*mm, 7.5*mm);
0234 
0235     tubeRadius      = 3.*mm;
0236     tubeHalfLengthZ = 5.5*mm; // 11 mm tot
0237 
0238     G4Material* diskMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_Ti");
0239     G4Material* targetMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_W");
0240 
0241     // Physical volumes
0242     G4Box* box = new G4Box("targetBox", boxHalfSide.getX(), boxHalfSide.getY(), boxHalfSide.getZ());
0243     G4Tubs* tube = new G4Tubs("targetTube", 0.*mm, tubeRadius, tubeHalfLengthZ, 0.*deg, 360.*deg);
0244 
0245     G4SubtractionSolid* BoxMinusTube = new G4SubtractionSolid("TargetSolid", box, tube, nullptr, {0.,0.,-2*mm});
0246     G4LogicalVolume* BoxMinusTubeLV = new G4LogicalVolume(BoxMinusTube, targetMaterial, "BoxMinusTubeLV",0,0,0);
0247     new G4PVPlacement(0, targetCentre, "TargetPV", BoxMinusTubeLV, accWorld_phys, false, checkOverlap);
0248 
0249     G4Tubs* disk = new G4Tubs("targetDisk", 0.*mm, 4.*mm, 0.025*mm, 0.*deg, 360.*deg);
0250     G4LogicalVolume* diskLV = new G4LogicalVolume(disk, diskMaterial, "targetDiskLV",0,0,0);
0251     new G4PVPlacement(0, {0.,0.,-15.*mm}, "diskPV", diskLV, accWorld_phys, false, checkOverlap);
0252 
0253     // Visualization
0254     G4VisAttributes* simpleAlSVisAtt;
0255     simpleAlSVisAtt = new G4VisAttributes(G4Colour::Grey());
0256     simpleAlSVisAtt -> SetVisibility(true);
0257     simpleAlSVisAtt -> SetForceSolid(true);
0258     BoxMinusTubeLV  -> SetVisAttributes(simpleAlSVisAtt);
0259     diskLV -> SetVisAttributes(simpleAlSVisAtt);
0260 
0261     // Region for cuts
0262     G4Region *regVol;
0263     regVol = new G4Region("TargetR");
0264     G4ProductionCuts* cuts = new G4ProductionCuts;
0265     cuts -> SetProductionCut(1.*mm);
0266     regVol -> SetProductionCuts(cuts);
0267     BoxMinusTubeLV -> SetRegion(regVol);
0268     regVol -> AddRootLogicalVolume(BoxMinusTubeLV);
0269     diskLV -> SetRegion(regVol);
0270     regVol -> AddRootLogicalVolume(diskLV);
0271 }
0272 
0273 void DetectorConstruction::ConstructPrimaryCollimator()
0274 {
0275     G4bool checkOverlap = true;
0276 
0277     G4double bigTube1HalfLengthZ, bigTube1Radius, bigTube1Z;
0278     G4double bigTube2HalfLengthZ, bigTube2Radius, bigTube2Z;
0279     G4double bigTube3HalfLengthZ, bigTube3Radius, bigTube3Z;
0280 
0281     G4double tube1Radius, tube1HalfLengthZ, tube2Radius, tube2HalfLengthZ;
0282 
0283     G4double cone1RadiusMin, cone1RadiusMax, cone1HalfLengthZ;
0284     G4double cone2RadiusMin, cone2RadiusMax, cone2HalfLengthZ;
0285     G4double cone3RadiusMin, cone3RadiusMax, cone3HalfLengthZ;
0286 
0287     // upper principal tube = bigTube1
0288     bigTube1Radius      = 141./2.*mm;
0289     bigTube1HalfLengthZ = 79./2. * mm;
0290     bigTube1Z           = 5. * mm + bigTube1HalfLengthZ;
0291     tube1Radius         = 34./2. * mm;
0292     tube1HalfLengthZ    = 15.86/2. * mm;
0293     cone1RadiusMin      = 34./2.*mm; // upper cone
0294     cone1RadiusMax      = 54./2.*mm;
0295     cone1HalfLengthZ    = (bigTube1HalfLengthZ - tube1HalfLengthZ)*mm; //(79.-15.86)/2.*mm;
0296 
0297     // middle principal tube = bigTube2
0298     bigTube2Radius      = 141./2. * mm;
0299     bigTube2HalfLengthZ = 57.5/2. * mm;
0300     bigTube2Z           = bigTube1Z + bigTube1HalfLengthZ + bigTube2HalfLengthZ;
0301     tube2Radius         = 53./2. * mm;
0302     tube2HalfLengthZ    = 1.35/2. * mm;
0303     cone2RadiusMin      = 53./2.*mm; // middle cone
0304     cone2RadiusMax      = 81./2.*mm;
0305     cone2HalfLengthZ    = bigTube2HalfLengthZ - tube2HalfLengthZ; // (57.50 - 1.35)/2.*mm;
0306 
0307     // lower principal tube = bigTube3
0308     bigTube3Radius      = 241./2.*mm;
0309     bigTube3HalfLengthZ = 40.50/2.*mm;
0310     bigTube3Z           = bigTube2Z + bigTube2HalfLengthZ + 19.5 + bigTube3HalfLengthZ;
0311     cone3RadiusMin      = 88.06/2.*mm; // lower cone
0312     cone3RadiusMax      = 109.76/2.*mm;
0313     cone3HalfLengthZ    = 40.50/2.*mm;
0314 
0315     G4Material* upperTubeMaterial = GetMaterial("xc10");
0316     G4Material* middleTubeMaterial = GetMaterial("wnicu");
0317     G4Material* lowerTubeMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
0318 
0319     // Physical volumes
0320     G4Tubs* bigTube1 = new G4Tubs("PriCollBigTube1", 0.*mm, bigTube1Radius, bigTube1HalfLengthZ, 0.*deg, 360.*deg);
0321     G4Tubs* tube1 = new G4Tubs("PriCollTube1", 0.*mm, tube1Radius, tube1HalfLengthZ, 0.*deg, 360.*deg);
0322     G4Cons* cone1 = new G4Cons("PriCollCone1", 0., cone1RadiusMin, 0., cone1RadiusMax, cone1HalfLengthZ, 0, 360.*deg);
0323 
0324     G4Tubs* bigTube2 = new G4Tubs("PriCollBigTube2", 0.*mm, bigTube2Radius, bigTube2HalfLengthZ, 0.*deg, 360.*deg);
0325     G4Tubs* tube2 = new G4Tubs("PriCollTube2", 0.*mm, tube2Radius, tube2HalfLengthZ, 0.*deg, 360.*deg);
0326     G4Cons* cone2 = new G4Cons("PriCollCone2", 0., cone2RadiusMin, 0., cone2RadiusMax, cone2HalfLengthZ, 0, 360.*deg);
0327 
0328     G4Tubs* bigTube3 = new G4Tubs("PriCollBigTube3", 0.*mm, bigTube3Radius, bigTube3HalfLengthZ, 0.*deg, 360.*deg);
0329     G4Cons* cone3 = new G4Cons("PriCollCone3", 0., cone3RadiusMin, 0., cone3RadiusMax, cone3HalfLengthZ, 0, 360.*deg);
0330 
0331     G4UnionSolid* tube1AndCone1 = new G4UnionSolid("PriCollTube1AndCone1",
0332                 tube1, cone1, nullptr, {0., 0., (tube1HalfLengthZ + cone1HalfLengthZ)});
0333     G4SubtractionSolid* bigTube1NotTube1AndCone1 = new G4SubtractionSolid(
0334             "PriColltube1NotTube1AndCone1", bigTube1, tube1AndCone1, 0,
0335             {0., 0., (-bigTube1HalfLengthZ + tube1HalfLengthZ)});
0336     G4LogicalVolume* bigTube1NotTube1AndCone1LV = new G4LogicalVolume(
0337             bigTube1NotTube1AndCone1, upperTubeMaterial, "PriCollBigTube1NotTube1AndCone1LV", 0, 0, 0);
0338 
0339     G4UnionSolid* tube2AndCone2 = new G4UnionSolid("PriCollTube2AndCone2",
0340                 tube2, cone2, nullptr, {0., 0., (tube2HalfLengthZ + cone2HalfLengthZ)});
0341     G4SubtractionSolid* bigTube2NotTube2AndCone2 = new G4SubtractionSolid(
0342             "PriCollBigTube2NotTube2Cone2", bigTube2, tube2AndCone2, 0,
0343             {0., 0., (-bigTube2HalfLengthZ + tube2HalfLengthZ)});
0344     G4LogicalVolume* bigTube2NotTube2AndCone2LV  = new G4LogicalVolume(
0345             bigTube2NotTube2AndCone2, middleTubeMaterial, "PriCollTube2NotTube2Cone2LV", 0, 0, 0);
0346 
0347     G4SubtractionSolid* bigTube3NotCone3 = new G4SubtractionSolid("PriCollTube4NotCone3", bigTube3, cone3);
0348     G4LogicalVolume* bigTube3NotCone3LV  = new G4LogicalVolume(bigTube3NotCone3,
0349             lowerTubeMaterial, "PriCollBigTube3NotCone3LV",0,0,0);
0350 
0351     new G4PVPlacement(nullptr, {0, 0, bigTube1Z}, "PriCollUpperPV",
0352             bigTube1NotTube1AndCone1LV, accWorld_phys, false, checkOverlap);
0353     new G4PVPlacement(nullptr, {0, 0, bigTube2Z}, "PriCollMiddlePV",
0354             bigTube2NotTube2AndCone2LV, accWorld_phys, false, checkOverlap);
0355     new G4PVPlacement(nullptr, {0, 0, bigTube3Z}, "PriCollLowerPV",
0356             bigTube3NotCone3LV, accWorld_phys, false, checkOverlap);
0357 
0358     // Visualization
0359     G4VisAttributes* simpleAlSVisAtt1 = new G4VisAttributes(G4Colour::Green());
0360     G4VisAttributes* simpleAlSVisAtt2 = new G4VisAttributes(G4Colour::Red());
0361     G4VisAttributes* simpleAlSVisAtt3 = new G4VisAttributes(G4Colour::Blue());
0362     simpleAlSVisAtt1 -> SetVisibility(true);
0363     simpleAlSVisAtt1 -> SetForceSolid(false);
0364     simpleAlSVisAtt1 -> SetForceWireframe(true);
0365     bigTube1NotTube1AndCone1LV -> SetVisAttributes(simpleAlSVisAtt1);  //primary collimator
0366     simpleAlSVisAtt2 -> SetVisibility(true);
0367     simpleAlSVisAtt2 -> SetForceSolid(false);
0368     simpleAlSVisAtt2 -> SetForceWireframe(true);
0369     bigTube2NotTube2AndCone2LV -> SetVisAttributes(simpleAlSVisAtt2); //primary collimator
0370     simpleAlSVisAtt3 -> SetVisibility(true);
0371     simpleAlSVisAtt3 -> SetForceSolid(false);
0372     simpleAlSVisAtt3 -> SetForceWireframe(true);
0373     bigTube3NotCone3LV -> SetVisAttributes(simpleAlSVisAtt3); //secondary collimator
0374 
0375     // Region for cuts
0376     G4Region* regVol = new G4Region("primaryCollimatorR");
0377     G4ProductionCuts* cuts = new G4ProductionCuts;
0378     cuts->SetProductionCut(1.*cm);
0379     regVol->SetProductionCuts(cuts);
0380     bigTube1NotTube1AndCone1LV -> SetRegion(regVol);
0381     regVol -> AddRootLogicalVolume(bigTube1NotTube1AndCone1LV);
0382     bigTube2NotTube2AndCone2LV -> SetRegion(regVol);
0383     regVol -> AddRootLogicalVolume(bigTube2NotTube2AndCone2LV);
0384     bigTube3NotCone3LV -> SetRegion(regVol);
0385     regVol -> AddRootLogicalVolume(bigTube3NotCone3LV);
0386 }
0387 
0388 void DetectorConstruction::ConstructFlatteningFilter()
0389 {
0390     G4double tube1HalfLengthZ; //tube1Radius,
0391     G4double cone1RadiusMin, cone1RadiusMax, cone1HalfLengthZ;
0392     G4double cone2RadiusMin, cone2RadiusMax, cone2HalfLengthZ;
0393     G4double cone3RadiusMin, cone3RadiusMax, cone3HalfLengthZ, ffZ;
0394 
0395     tubeFFRadius     = 108./2.*mm; // top principal tube
0396     tube1HalfLengthZ = 7.5/2.*mm;
0397 
0398     cone1RadiusMin   = 54./2.*mm; // top cone
0399     cone1RadiusMax   = 76./2.*mm;
0400     cone1HalfLengthZ = 13.7/2.*mm;
0401 
0402     cone2RadiusMin   = 8./2.*mm; // mid cone
0403     cone2RadiusMax   = 54./2.*mm;
0404     cone2HalfLengthZ = (44.3-13.7)/2.*mm;
0405 
0406     cone3RadiusMin   = 0.000001*mm; // bottom cone
0407     cone3RadiusMax   = 8./2.*mm;
0408     cone3HalfLengthZ = (46.8-44.3)/2.*mm;
0409     tubeFFFirstFaceZ   = 149.50*mm;
0410 
0411     ffZ = tubeFFFirstFaceZ + tube1HalfLengthZ; // the half point is located ad the centre of the tube1 because of the solids unions and translations
0412 
0413     G4Material* ffMaterial = GetMaterial("steel");
0414 
0415     // Physical volumes
0416     G4Tubs* tube1 = new G4Tubs("ffTube1", 0.*mm, tubeFFRadius, tube1HalfLengthZ, 0.*deg, 360.*deg);
0417     G4Cons* cone1 = new G4Cons("ffCone1", 0., cone1RadiusMax, 0., cone1RadiusMin, cone1HalfLengthZ, 0, 360.*deg);
0418     G4Cons* cone2 = new G4Cons("ffCone2", 0., cone2RadiusMax, 0., cone2RadiusMin, cone2HalfLengthZ, 0, 360.*deg);
0419     G4Cons* cone3 = new G4Cons("ffCone3", 0., cone3RadiusMax, 0., cone3RadiusMin, cone3HalfLengthZ, 0, 360.*deg);
0420 
0421     G4double pos = (cone1HalfLengthZ - tube1HalfLengthZ - 1.)*mm; // == (13.7/2. - 7.5/2. -1.)*mm
0422     G4UnionSolid *tube1AndCone1 = new G4UnionSolid("ffTube1AndCone1", tube1,
0423             cone1, 0, {0., 0., pos});
0424 
0425     pos += cone1HalfLengthZ + cone2HalfLengthZ;
0426     G4UnionSolid* tubeCone1AndCone2 = new G4UnionSolid("fftubeCone1AndCone2",
0427             tube1AndCone1, cone2, 0, {0., 0., pos});
0428 
0429     pos += cone2HalfLengthZ + cone3HalfLengthZ;
0430     G4UnionSolid* tubeCone12AndCone3 = new G4UnionSolid("fftubeCone12AndCone3",
0431             tubeCone1AndCone2, cone3, 0, {0., 0., pos});
0432 
0433     G4LogicalVolume* ffLV = new G4LogicalVolume(tubeCone12AndCone3, ffMaterial, "ffLV",0,0,0);
0434 
0435     new G4PVPlacement(0, {0,0,ffZ}, "ffPV", ffLV, accWorld_phys, false, 0);
0436 
0437     // Visualization
0438     G4VisAttributes* simpleAlSVisAtt= new G4VisAttributes(G4Colour::Magenta());
0439     simpleAlSVisAtt->SetVisibility(true);
0440     simpleAlSVisAtt->SetForceSolid(true);
0441     ffLV->SetVisAttributes(simpleAlSVisAtt);
0442 
0443     // Region for cuts
0444     G4Region* regVol;
0445     regVol= new G4Region("ffR");
0446     G4ProductionCuts* cuts = new G4ProductionCuts;
0447     cuts->SetProductionCut(0.1*cm);
0448     regVol->SetProductionCuts(cuts);
0449     ffLV->SetRegion(regVol);
0450     regVol->AddRootLogicalVolume(ffLV);
0451 }
0452 
0453 void DetectorConstruction::ConstructIonizationChamber()
0454 {
0455     G4bool checkOverlap = true;
0456 
0457     G4double tubeRadius, tubeA1Z, tubeA2Z, tubeA3Z, tubeA4Z, tubeA5Z, tubeA6Z;
0458     G4double kaptonThickness, AlThickness1, AlThickness8;
0459 
0460     kaptonThickness = 0.025/2. *mm;
0461     AlThickness1    = 1.6e-5/2. *cm;
0462     AlThickness8    = 8e-6/2. *cm;
0463 
0464     tubeRadius = 110./2.*mm; // upper principal tube
0465     tubeA1Z = (202.5 + AlThickness8)*mm;
0466     tubeA2Z = (204.5 + AlThickness1)*mm;
0467     tubeA3Z = (206.5 + AlThickness8)*mm;
0468     tubeA4Z = (207.5 + AlThickness8)*mm;
0469     tubeA5Z = (209.5 + AlThickness1)*mm;
0470     tubeA6Z = (211.5 + AlThickness8)*mm;
0471 
0472     G4double pos1 = (AlThickness1 + kaptonThickness)*mm;
0473     G4double pos8 = (AlThickness8 + kaptonThickness)*mm;
0474 
0475     G4Material* el_Al = G4NistManager::Instance()->FindOrBuildMaterial("G4_Al");
0476     G4Material* Kapton = GetMaterial("kapton");
0477 
0478     // Physical volumes
0479     G4Tubs* tubeAl8 = new G4Tubs("ICTube1A", 0.*mm, tubeRadius, AlThickness8, 0.*deg, 360.*deg);
0480     G4Tubs* tubeK = new G4Tubs("ICTube1B", 0.*mm, tubeRadius, kaptonThickness, 0.*deg, 360.*deg);
0481     G4Tubs* tubeAl1 = new G4Tubs("ICTube2A", 0.*mm, tubeRadius, AlThickness1, 0.*deg, 360.*deg);
0482 
0483     G4LogicalVolume* IC_Al8_LV = new G4LogicalVolume(tubeAl8, el_Al, "IC_Al8_LV", 0, 0, 0);
0484     G4LogicalVolume* IC_Al1_LV = new G4LogicalVolume(tubeAl1, el_Al, "IC_Al1_LV", 0, 0, 0);
0485     G4LogicalVolume* IC_K_LV   = new G4LogicalVolume(tubeK, Kapton, "IC_Al_LV",  0, 0, 0);
0486 
0487     new G4PVPlacement(0, {0,0,tubeA1Z}, "IC_AL8_PV1", IC_Al8_LV, accWorld_phys, false, checkOverlap);
0488     new G4PVPlacement(0, {0,0,tubeA1Z + pos8}, "IC_K_PV1", IC_K_LV, accWorld_phys, false, checkOverlap);
0489 
0490     new G4PVPlacement(0, {0,0,tubeA2Z}, "IC_AL1_PV2", IC_Al1_LV, accWorld_phys, false, checkOverlap);
0491     new G4PVPlacement(0, {0,0,tubeA2Z + pos1}, "IC_K_PV2", IC_K_LV, accWorld_phys, false, checkOverlap);
0492 
0493     new G4PVPlacement(0, {0,0,tubeA3Z}, "IC_AL8_PV3", IC_Al8_LV, accWorld_phys, false, checkOverlap);
0494     new G4PVPlacement(0, {0,0,tubeA3Z + pos8}, "IC_K_PV3", IC_K_LV, accWorld_phys, false, checkOverlap);
0495 
0496     new G4PVPlacement(0, {0,0,tubeA4Z}, "IC_AL8_PV4", IC_Al8_LV, accWorld_phys, false, checkOverlap);
0497     new G4PVPlacement(0, {0,0,tubeA4Z + pos8}, "IC_K_PV4", IC_K_LV, accWorld_phys, false, checkOverlap);
0498 
0499     new G4PVPlacement(0, {0,0,tubeA5Z}, "IC_AL1_PV5", IC_Al1_LV, accWorld_phys, false, checkOverlap);
0500     new G4PVPlacement(0, {0,0,tubeA5Z + pos1}, "IC_K_PV5", IC_K_LV, accWorld_phys, false, checkOverlap);
0501 
0502     new G4PVPlacement(0, {0,0,tubeA6Z}, "IC_AL8_PV6", IC_Al8_LV, accWorld_phys, false, checkOverlap);
0503     new G4PVPlacement(0, {0,0,tubeA6Z + pos8}, "IC_K_PV6", IC_K_LV, accWorld_phys, false, checkOverlap);
0504 
0505     // Visualization
0506     G4VisAttributes* simpleAlSVisAttAl1 = new G4VisAttributes(G4Colour::Red());
0507     simpleAlSVisAttAl1 -> SetVisibility(true);
0508     simpleAlSVisAttAl1 -> SetForceSolid(true);
0509     IC_Al1_LV -> SetVisAttributes(simpleAlSVisAttAl1);
0510 
0511     G4VisAttributes* simpleAlSVisAttAl8 = new G4VisAttributes(G4Colour::Magenta());
0512     simpleAlSVisAttAl8 -> SetVisibility(true);
0513     simpleAlSVisAttAl8 -> SetForceSolid(true);
0514     IC_Al8_LV -> SetVisAttributes(simpleAlSVisAttAl8);
0515 
0516     G4VisAttributes* simpleAlSVisAttK = new G4VisAttributes(G4Colour::Blue());
0517     simpleAlSVisAttK -> SetVisibility(true);
0518     simpleAlSVisAttK -> SetForceSolid(true);
0519     IC_K_LV -> SetVisAttributes(simpleAlSVisAttK);
0520 
0521     // Region for cuts
0522     G4Region * regVol = new G4Region("IonChamberR");
0523     G4ProductionCuts* cuts = new G4ProductionCuts;
0524     cuts -> SetProductionCut(0.1*mm);
0525     regVol -> SetProductionCuts(cuts);
0526 
0527     IC_Al1_LV -> SetRegion(regVol);
0528     regVol -> AddRootLogicalVolume(IC_Al1_LV);
0529     IC_Al8_LV -> SetRegion(regVol);
0530     regVol -> AddRootLogicalVolume(IC_Al8_LV);
0531     IC_K_LV -> SetRegion(regVol);
0532     regVol -> AddRootLogicalVolume(IC_K_LV);
0533 }
0534 
0535 void DetectorConstruction::ConstructVacuumWindow()
0536 {
0537     G4double tubeRadius, tubeHalfLengthZ, tubeZ;
0538 
0539     tubeRadius      = 116.53/2.*mm; // upper principal tube
0540     tubeHalfLengthZ = 2./2.*mm;
0541     tubeZ           = 215.75*mm + tubeHalfLengthZ;
0542 
0543     G4Material* el_Al= G4NistManager::Instance()->FindOrBuildMaterial("G4_Al");
0544 
0545     // Physical volumes
0546     G4Tubs* tube = new G4Tubs("VacWindowTube", 0.*mm, tubeRadius, tubeHalfLengthZ, 0.*deg, 360.*deg);
0547     G4LogicalVolume* tubeLV = new G4LogicalVolume(tube, el_Al, "VacWindowTubeLV",0,0,0);
0548 
0549     new G4PVPlacement(0, {0,0,tubeZ}, "VacWindowTubePV", tubeLV, accWorld_phys, false, 0);
0550 
0551     // Visualization
0552     G4VisAttributes* simpleAlSVisAtt = new G4VisAttributes(G4Colour::Green());
0553     simpleAlSVisAtt -> SetVisibility(true);
0554     simpleAlSVisAtt -> SetForceSolid(true);
0555     tubeLV -> SetVisAttributes(simpleAlSVisAtt);
0556 
0557     // Region for cuts
0558     G4Region* regVol= new G4Region("VacWindowR");
0559     G4ProductionCuts* cuts = new G4ProductionCuts;
0560     cuts->SetProductionCut(0.1*cm);
0561     regVol->SetProductionCuts(cuts);
0562     tubeLV -> SetRegion(regVol);
0563     regVol -> AddRootLogicalVolume(tubeLV);
0564 }
0565 void DetectorConstruction::ConstructJawsX()
0566 {
0567     G4bool checkOverlap = false;
0568 
0569     G4double boxSide         = 101.*mm;
0570     G4double box1HalfLengthZ =   3.*mm;
0571     G4double box2HalfLengthZ =  35.*mm;
0572     G4double box3HalfLengthZ =  35.*mm;
0573     G4double box4HalfLengthZ =  27.*mm;
0574     G4double box5HalfLengthZ =  10.*mm;
0575     G4double box6HalfLengthZ =   5.*mm;
0576 
0577     G4Material* el_Pb = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
0578     G4Material* XC10 = GetMaterial("xc10");
0579     G4Material* WNICU = GetMaterial("wnicu");
0580 
0581     G4ThreeVector boxJawXSide = {boxSide, boxSide, 212.*mm};
0582 
0583     G4Box* boxJaw1X = new G4Box("Jaw1Xbox", boxJawXSide.getX()/2., boxJawXSide.getY()/2., boxJawXSide.getZ()/2.);
0584     G4LogicalVolume* boxJaw1XLV = new G4LogicalVolume(boxJaw1X,
0585             G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw1XboxLV", 0, 0, 0);
0586     G4Box* boxJaw2X = new G4Box("Jaw2Xbox", boxJawXSide.getX()/2., boxJawXSide.getY()/2., boxJawXSide.getZ()/2.);
0587     G4LogicalVolume* boxJaw2XLV = new G4LogicalVolume(boxJaw2X,
0588                 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw2XboxLV", 0, 0, 0);
0589 
0590     jaw1XInitialPos.setX(boxJawXSide.getX()/2.);
0591     jaw1XInitialPos.setY(0.);
0592     jaw1XInitialPos.setZ(275.5*mm + boxJawXSide.getZ()/2.);
0593 
0594     G4RotationMatrix* cRotationJaw1X = new G4RotationMatrix();
0595     cRotationJaw1X -> rotateY(std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
0596 
0597     G4ThreeVector position = *cRotationJaw1X * jaw1XInitialPos;
0598     *cRotationJaw1X = cRotationJaw1X -> inverse();
0599 
0600     boxJaw1X_phys =
0601         new G4PVPlacement (cRotationJaw1X, position,
0602                 "Jaw1XPV", boxJaw1XLV, accWorld_phys, false, 0, checkOverlap);
0603 
0604     jaw2XInitialPos.setX(-jaw1XInitialPos.getX());
0605     jaw2XInitialPos.setY(jaw1XInitialPos.getY());
0606     jaw2XInitialPos.setZ(jaw1XInitialPos.getZ());
0607 
0608     G4RotationMatrix* cRotationJaw2X = new G4RotationMatrix();
0609     cRotationJaw2X -> rotateY(-std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
0610 
0611     position = *cRotationJaw2X * jaw2XInitialPos;
0612     *cRotationJaw2X = cRotationJaw2X -> inverse();
0613 
0614     boxJaw2X_phys =
0615                 new G4PVPlacement (cRotationJaw2X, position,
0616                         "Jaw2XPV", boxJaw2XLV, accWorld_phys, false, 0, 0);
0617 
0618     // Physical volumes
0619     G4Box* box1 = new G4Box("Jaws1XBox1", boxSide/2., boxSide/2., box1HalfLengthZ/2.);
0620     G4Box* box2 = new G4Box("Jaws1XBox2", boxSide/2., boxSide/2., box2HalfLengthZ/2.);
0621     G4Box* box3 = new G4Box("Jaws1XBox3", boxSide/2., boxSide/2., box3HalfLengthZ/2.);
0622     G4Box* box4 = new G4Box("Jaws1XBox4", boxSide/2., boxSide/2., box4HalfLengthZ/2.);
0623     G4Box* box5 = new G4Box("Jaws1XBox5", boxSide/2., boxSide/2., box5HalfLengthZ/2.);
0624     G4Box* box6 = new G4Box("Jaws1XBox6", boxSide/2., boxSide/2., box6HalfLengthZ/2.);
0625     G4LogicalVolume* box1LV1 = new G4LogicalVolume(box1, XC10,  "Jaws1XLV1", 0, 0, 0);
0626     G4LogicalVolume* box2LV1 = new G4LogicalVolume(box2, el_Pb, "Jaws1XLV2", 0, 0, 0);
0627     G4LogicalVolume* box3LV1 = new G4LogicalVolume(box3, WNICU, "Jaws1XLV3", 0, 0, 0);
0628     G4LogicalVolume* box4LV1 = new G4LogicalVolume(box4, el_Pb, "Jaws1XLV4", 0, 0, 0);
0629     G4LogicalVolume* box5LV1 = new G4LogicalVolume(box5, el_Pb, "Jaws1XLV5", 0, 0, 0);
0630     G4LogicalVolume* box6LV1 = new G4LogicalVolume(box6, WNICU, "Jaws1XLV6", 0, 0, 0);
0631 
0632     G4LogicalVolume* box1LV2 = new G4LogicalVolume(box1, XC10,  "Jaws2XLV1", 0, 0, 0);
0633     G4LogicalVolume* box2LV2 = new G4LogicalVolume(box2, el_Pb, "Jaws2XLV2", 0, 0, 0);
0634     G4LogicalVolume* box3LV2 = new G4LogicalVolume(box3, WNICU, "Jaws2XLV3", 0, 0, 0);
0635     G4LogicalVolume* box4LV2 = new G4LogicalVolume(box4, el_Pb, "Jaws2XLV4", 0, 0, 0);
0636     G4LogicalVolume* box5LV2 = new G4LogicalVolume(box5, el_Pb, "Jaws2XLV5", 0, 0, 0);
0637     G4LogicalVolume* box6LV2 = new G4LogicalVolume(box6, WNICU, "Jaws2XLV6", 0, 0, 0);
0638 
0639     G4ThreeVector centre = {0., 0., 0.};
0640     G4double zCentreCurrentBox = -boxJawXSide.getZ()/2.+box1HalfLengthZ/2.*mm;
0641 
0642     centre.setZ(zCentreCurrentBox);
0643     new G4PVPlacement(nullptr, centre, "Jaws1XPV1", box1LV1, boxJaw1X_phys, false, 0, checkOverlap);
0644     new G4PVPlacement(nullptr, centre, "Jaws2XPV1", box1LV2, boxJaw2X_phys, false, 0, checkOverlap);
0645 
0646     zCentreCurrentBox += (box1HalfLengthZ+box2HalfLengthZ)/2.*mm;
0647     centre.setZ(zCentreCurrentBox);
0648     new G4PVPlacement(nullptr, centre, "Jaws1XPV2", box2LV1, boxJaw1X_phys, false, 0, checkOverlap);
0649     new G4PVPlacement(nullptr, centre, "Jaws2XPV2", box2LV2, boxJaw2X_phys, false, 0, checkOverlap);
0650 
0651     zCentreCurrentBox += (box2HalfLengthZ+box3HalfLengthZ)/2.*mm;
0652     centre.setZ(zCentreCurrentBox);
0653     new G4PVPlacement(nullptr, centre, "Jaws1XPV3", box3LV1, boxJaw1X_phys, false, 0, checkOverlap);
0654     new G4PVPlacement(nullptr, centre, "Jaws2XPV3", box3LV2, boxJaw2X_phys, false, 0, checkOverlap);
0655 
0656     zCentreCurrentBox += (box3HalfLengthZ+box4HalfLengthZ)/2.*mm;
0657     centre.setZ(zCentreCurrentBox);
0658     new G4PVPlacement(nullptr, centre, "Jaws1XPV4", box4LV1, boxJaw1X_phys, false, 0, checkOverlap);
0659     new G4PVPlacement(nullptr, centre, "Jaws2XPV4", box4LV2, boxJaw2X_phys, false, 0, checkOverlap);
0660 
0661     zCentreCurrentBox += (box4HalfLengthZ+box5HalfLengthZ)/2.+97.*mm;
0662     centre.setZ(zCentreCurrentBox);
0663     new G4PVPlacement(nullptr, centre, "Jaws1XPV5", box5LV1, boxJaw1X_phys, false, 0, checkOverlap);
0664     new G4PVPlacement(nullptr, centre, "Jaws2XPV5", box5LV2, boxJaw2X_phys, false, 0, checkOverlap);
0665 
0666     zCentreCurrentBox += (box5HalfLengthZ+box6HalfLengthZ)/2.;
0667     centre.setZ(zCentreCurrentBox);
0668     new G4PVPlacement(nullptr, centre, "Jaws1XPV6", box6LV1, boxJaw1X_phys, false, 0, checkOverlap);
0669     new G4PVPlacement(nullptr, centre, "Jaws2XPV6", box6LV2, boxJaw2X_phys, false, 0, checkOverlap);
0670 
0671     // Visualization
0672     G4VisAttributes* visAtt = new G4VisAttributes(G4Colour::Yellow());
0673     visAtt -> SetVisibility(false);
0674     visAtt -> SetForceSolid(false);
0675 
0676     G4VisAttributes* simpleAlSVisAttPb = new G4VisAttributes(G4Colour::Blue());
0677     simpleAlSVisAttPb -> SetVisibility(true);
0678     simpleAlSVisAttPb -> SetForceSolid(true);
0679 
0680     G4VisAttributes* simpleAlSVisAttXC10 = new G4VisAttributes(G4Colour::Green());
0681     simpleAlSVisAttXC10 -> SetVisibility(true);
0682     simpleAlSVisAttXC10 ->SetForceSolid(true);
0683 
0684     G4VisAttributes* simpleAlSVisAttWNICU = new G4VisAttributes(G4Colour::Red());
0685     simpleAlSVisAttWNICU ->SetVisibility(true);
0686     simpleAlSVisAttWNICU ->SetForceSolid(true);
0687 
0688     box1LV1 -> SetVisAttributes(simpleAlSVisAttXC10);
0689     box2LV1 -> SetVisAttributes(simpleAlSVisAttPb);
0690     box3LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
0691     box4LV1 -> SetVisAttributes(simpleAlSVisAttPb);
0692     box5LV1 -> SetVisAttributes(simpleAlSVisAttPb);
0693     box6LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
0694 
0695     box1LV2 -> SetVisAttributes(simpleAlSVisAttXC10);
0696     box2LV2 -> SetVisAttributes(simpleAlSVisAttPb);
0697     box3LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
0698     box4LV2 -> SetVisAttributes(simpleAlSVisAttPb);
0699     box5LV2 -> SetVisAttributes(simpleAlSVisAttPb);
0700     box6LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
0701 
0702     boxJaw1XLV -> SetVisAttributes(visAtt);
0703     boxJaw2XLV -> SetVisAttributes(visAtt);
0704 
0705     // Region for cuts
0706     G4Region* regVol = new G4Region("JawsXR");
0707     G4ProductionCuts* cuts = new G4ProductionCuts;
0708     cuts -> SetProductionCut(2.*cm);
0709     regVol -> SetProductionCuts(cuts);
0710     box1LV1 -> SetRegion(regVol);
0711     regVol -> AddRootLogicalVolume(box1LV1);
0712     box2LV1 -> SetRegion(regVol);
0713     regVol -> AddRootLogicalVolume(box2LV1);
0714     box3LV1 -> SetRegion(regVol);
0715     regVol -> AddRootLogicalVolume(box3LV1);
0716     box4LV1 -> SetRegion(regVol);
0717     regVol -> AddRootLogicalVolume(box4LV1);
0718     box5LV1 -> SetRegion(regVol);
0719     regVol -> AddRootLogicalVolume(box5LV1);
0720     box6LV1 -> SetRegion(regVol);
0721     regVol -> AddRootLogicalVolume(box6LV1);
0722 
0723     box1LV2 -> SetRegion(regVol);
0724     regVol -> AddRootLogicalVolume(box1LV2);
0725     box2LV2 -> SetRegion(regVol);
0726     regVol -> AddRootLogicalVolume(box2LV2);
0727     box3LV2 -> SetRegion(regVol);
0728     regVol -> AddRootLogicalVolume(box3LV2);
0729     box4LV2 -> SetRegion(regVol);
0730     regVol -> AddRootLogicalVolume(box4LV2);
0731     box5LV2 -> SetRegion(regVol);
0732     regVol -> AddRootLogicalVolume(box5LV2);
0733     box6LV2 -> SetRegion(regVol);
0734     regVol -> AddRootLogicalVolume(box6LV2);
0735 }
0736 
0737 void DetectorConstruction::ConstructJawsY()
0738 {
0739     G4bool checkOverlap = false;
0740 
0741     G4double boxSide         = 101.*mm;
0742     G4double box1HalfLengthZ =  15.*mm;
0743     G4double box2HalfLengthZ =  31.*mm;
0744     G4double box3HalfLengthZ =  35.*mm;
0745     G4double box4HalfLengthZ =  21.*mm;
0746 
0747     G4Material* el_Pb = G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb");
0748     G4Material* XC10 = GetMaterial("xc10");
0749     G4Material* WNICU = GetMaterial("wnicu");
0750 
0751     G4ThreeVector boxJawYSide = {boxSide, boxSide, 219.5*mm}; // 219, perché 219.5??
0752 
0753     G4Box* boxJaw1Y = new G4Box("Jaw1Ybox", boxJawYSide.getX()/2., boxJawYSide.getY()/2., boxJawYSide.getZ()/2.);
0754     G4LogicalVolume* boxJaw1YLV = new G4LogicalVolume(boxJaw1Y,
0755             G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw1YboxLV", 0, 0, 0);
0756     G4Box* boxJaw2Y = new G4Box("Jaw2Ybox", boxJawYSide.getX()/2., boxJawYSide.getY()/2., boxJawYSide.getZ()/2.);
0757         G4LogicalVolume* boxJaw2YLV = new G4LogicalVolume(boxJaw2Y,
0758                 G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "Jaw2YboxLV", 0, 0, 0);
0759 
0760     jaw1YInitialPos.setX(0.);
0761     jaw1YInitialPos.setY(boxJawYSide.getX()/2.);
0762     jaw1YInitialPos.setZ(248.5*mm + boxJawYSide.getZ()/2.);
0763 
0764     G4RotationMatrix* cRotationJaw1Y = new G4RotationMatrix();
0765     cRotationJaw1Y -> rotateX(-std::fabs(std::atan(jawAperture/2/-(-(sourceToSkinDistance+10*cm)))));
0766 
0767     G4ThreeVector position = *cRotationJaw1Y * jaw1YInitialPos;
0768     *cRotationJaw1Y = cRotationJaw1Y -> inverse();
0769 
0770     boxJaw1Y_phys =
0771             new G4PVPlacement (cRotationJaw1Y, position,
0772                     "Jaw1YPV", boxJaw1YLV, accWorld_phys, false, 0, checkOverlap);
0773 
0774     jaw2YInitialPos.setX(jaw1YInitialPos.getX());
0775     jaw2YInitialPos.setY(-jaw1YInitialPos.getY());
0776     jaw2YInitialPos.setZ(jaw1YInitialPos.getZ());
0777 
0778     G4RotationMatrix* cRotationJaw2Y = new G4RotationMatrix();
0779     cRotationJaw2Y -> rotateX(std::fabs(std::atan(jawAperture/2/(-(sourceToSkinDistance+10*cm)))));
0780 
0781     position = *cRotationJaw2Y * jaw2YInitialPos;
0782     *cRotationJaw2Y = cRotationJaw2Y -> inverse();
0783 
0784     boxJaw2Y_phys =
0785             new G4PVPlacement (cRotationJaw2Y, position,
0786                     "Jaw2YPV", boxJaw2YLV, accWorld_phys, false, 0, checkOverlap);
0787 
0788     // Physical volumes
0789     G4Box *box1 = new G4Box("Jaws1YBox1", boxSide/2., boxSide/2., box1HalfLengthZ/2.);
0790     G4Box *box2 = new G4Box("Jaws1YBox2", boxSide/2., boxSide/2., box2HalfLengthZ/2.);
0791     G4Box *box3 = new G4Box("Jaws1YBox3", boxSide/2., boxSide/2., box3HalfLengthZ/2.);
0792     G4Box *box4 = new G4Box("Jaws1YBox4", boxSide/2., boxSide/2., box4HalfLengthZ/2.);
0793     G4LogicalVolume* box1LV1 = new G4LogicalVolume(box1, XC10,  "Jaws1YLV1", 0,0,0);
0794     G4LogicalVolume* box2LV1 = new G4LogicalVolume(box2, el_Pb, "Jaws1YLV2", 0,0,0);
0795     G4LogicalVolume* box3LV1 = new G4LogicalVolume(box3, WNICU, "Jaws1YLV3", 0,0,0);
0796     G4LogicalVolume* box4LV1 = new G4LogicalVolume(box4, el_Pb, "Jaws1YLV4", 0,0,0);
0797 
0798     G4LogicalVolume* box1LV2 = new G4LogicalVolume(box1, XC10,  "Jaws2YLV1", 0,0,0);
0799     G4LogicalVolume* box2LV2 = new G4LogicalVolume(box2, el_Pb, "Jaws2YLV2", 0,0,0);
0800     G4LogicalVolume* box3LV2 = new G4LogicalVolume(box3, WNICU, "Jaws2YLV3", 0,0,0);
0801     G4LogicalVolume* box4LV2 = new G4LogicalVolume(box4, el_Pb, "Jaws2YLV4", 0,0,0);
0802 
0803     G4ThreeVector centre = {0., 0., 0.};
0804     G4double zCentreCurrentBox = -boxJawYSide.getZ()/2.+box1HalfLengthZ/2.*mm;
0805 
0806     centre.setZ(zCentreCurrentBox);
0807     new G4PVPlacement(nullptr, centre, "Jaws1YPV1", box1LV1, boxJaw1Y_phys, false, 0, checkOverlap);
0808     new G4PVPlacement(nullptr, centre, "Jaws2YPV1", box1LV2, boxJaw2Y_phys, false, 0, checkOverlap);
0809 
0810     zCentreCurrentBox += (box1HalfLengthZ+box2HalfLengthZ)/2.*mm + 117.5*mm;
0811     centre.setZ(zCentreCurrentBox);
0812     new G4PVPlacement(nullptr, centre, "Jaws1YPV2", box2LV1, boxJaw1Y_phys, false, 0, checkOverlap);
0813     new G4PVPlacement(nullptr, centre, "Jaws2YPV2", box2LV2, boxJaw2Y_phys, false, 0, checkOverlap);
0814 
0815     zCentreCurrentBox += (box2HalfLengthZ+box3HalfLengthZ)/2.*mm;
0816     centre.setZ(zCentreCurrentBox);
0817     new G4PVPlacement(nullptr, centre, "Jaws1YPV3", box3LV1, boxJaw1Y_phys, false, 0, checkOverlap);
0818     new G4PVPlacement(nullptr, centre, "Jaws2YPV3", box3LV2, boxJaw2Y_phys, false, 0, checkOverlap);
0819 
0820     zCentreCurrentBox += (box3HalfLengthZ+box4HalfLengthZ)/2.*mm;
0821     centre.setZ(zCentreCurrentBox);
0822     new G4PVPlacement(nullptr, centre, "Jaws1YPV4", box4LV1, boxJaw1Y_phys, false, 0, checkOverlap);
0823     new G4PVPlacement(nullptr, centre, "Jaws2YPV4", box4LV2, boxJaw2Y_phys, false, 0, checkOverlap);
0824 
0825     // Visualization
0826     G4VisAttributes* visAtt = new G4VisAttributes(G4Colour::Yellow());
0827     visAtt -> SetVisibility(false);
0828     visAtt -> SetForceSolid(false);
0829 
0830     G4VisAttributes* simpleAlSVisAttPb = new G4VisAttributes(G4Colour::Blue());
0831     simpleAlSVisAttPb -> SetVisibility(true);
0832     simpleAlSVisAttPb -> SetForceSolid(true);
0833 
0834     G4VisAttributes* simpleAlSVisAttXC10 = new G4VisAttributes(G4Colour::Green());
0835     simpleAlSVisAttXC10 -> SetVisibility(true);
0836     simpleAlSVisAttXC10 -> SetForceSolid(true);
0837 
0838     G4VisAttributes* simpleAlSVisAttWNICU = new G4VisAttributes(G4Colour::Red());
0839     simpleAlSVisAttWNICU -> SetVisibility(true);
0840     simpleAlSVisAttWNICU -> SetForceSolid(true);
0841 
0842     box1LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
0843     box2LV1 -> SetVisAttributes(simpleAlSVisAttPb);
0844     box3LV1 -> SetVisAttributes(simpleAlSVisAttWNICU);
0845     box4LV1 -> SetVisAttributes(simpleAlSVisAttPb);
0846 
0847     box1LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
0848     box2LV2 -> SetVisAttributes(simpleAlSVisAttPb);
0849     box3LV2 -> SetVisAttributes(simpleAlSVisAttWNICU);
0850     box4LV2 -> SetVisAttributes(simpleAlSVisAttPb);
0851 
0852     boxJaw1YLV -> SetVisAttributes(visAtt);
0853     boxJaw2YLV -> SetVisAttributes(visAtt);
0854 
0855     // Region for cuts
0856     G4Region* regVol = new G4Region("JawsYR");
0857     G4ProductionCuts* cuts = new G4ProductionCuts;
0858     cuts -> SetProductionCut(2.*cm);
0859     regVol -> SetProductionCuts(cuts);
0860 
0861     box1LV1 -> SetRegion(regVol);
0862     regVol -> AddRootLogicalVolume(box1LV1);
0863     box2LV1 -> SetRegion(regVol);
0864     regVol -> AddRootLogicalVolume(box2LV1);
0865     box3LV1 -> SetRegion(regVol);
0866     regVol -> AddRootLogicalVolume(box3LV1);
0867     box4LV1 -> SetRegion(regVol);
0868     regVol -> AddRootLogicalVolume(box4LV1);
0869 
0870     box1LV2 -> SetRegion(regVol);
0871     regVol -> AddRootLogicalVolume(box1LV2);
0872     box2LV2 -> SetRegion(regVol);
0873     regVol -> AddRootLogicalVolume(box2LV2);
0874     box3LV2 -> SetRegion(regVol);
0875     regVol -> AddRootLogicalVolume(box3LV2);
0876     box4LV2 -> SetRegion(regVol);
0877     regVol -> AddRootLogicalVolume(box4LV2);
0878 }
0879 
0880 void DetectorConstruction::PhantomSegmentation(G4LogicalVolume* Phantom_log)
0881 {
0882 
0883     G4ThreeVector phantomHalfDim = {phantomSideDim/2., phantomSideDim/2., phantomSideDim/2.};
0884 
0885     nSideCells = CheckPhantomSegmentation(phantomSideDim/voxelSideDim);
0886     nDepthCells = CheckPhantomSegmentation(phantomSideDim/voxelDepthDim);
0887 
0888     G4Material* water = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
0889 
0890     G4String yRepName("RepY");
0891     G4VSolid* solYRep = new G4Box(yRepName, phantomHalfDim.getX(),
0892             phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ());
0893     G4LogicalVolume* logYRep = new G4LogicalVolume(solYRep, water, yRepName);
0894     new G4PVReplica(yRepName, logYRep, Phantom_log, kYAxis, nSideCells,
0895             2.*phantomHalfDim.getY()/nSideCells);
0896 
0897     G4String xRepName("RepX");
0898     G4VSolid* solXRep = new G4Box(xRepName, phantomHalfDim.getX()/nSideCells,
0899             phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ());
0900     G4LogicalVolume* logXRep = new G4LogicalVolume(solXRep, water, xRepName);
0901     new G4PVReplica(xRepName, logXRep, logYRep, kXAxis, nSideCells,
0902             2.*phantomHalfDim.getX()/nSideCells);
0903 
0904     G4String zVoxName("phantomSens");
0905     G4VSolid* solVoxel = new G4Box(zVoxName, phantomHalfDim.getX()/nSideCells,
0906             phantomHalfDim.getY()/nSideCells, phantomHalfDim.getZ()/nDepthCells);
0907 
0908     LVPhantomSens = new G4LogicalVolume(solVoxel, water, zVoxName); // This is the Sensitive Volume
0909     new G4PVReplica(zVoxName, LVPhantomSens, logXRep, kZAxis, nDepthCells,
0910             2.*phantomHalfDim.getZ()/nDepthCells);
0911 
0912 }
0913 
0914 G4int DetectorConstruction::CheckPhantomSegmentation(G4int nCells)
0915 {
0916     // I want an odd number of voxels
0917     // check whether nCells is integer
0918     if ( nCells != (int)nCells )
0919     {
0920         nCells = std::ceil(nCells);
0921     }
0922     // check whether nCells is even
0923     if ( std::fmod(nCells,2) == 0)
0924     {
0925         nCells += 1;
0926     }
0927 
0928     return nCells;
0929 }
0930 
0931 G4double DetectorConstruction::GetAccOriginPosition()
0932 {
0933     return sourceToSkinDistance;
0934 }
0935 
0936 G4double DetectorConstruction::GetFFilterRadius()
0937 {
0938     return tubeFFRadius;
0939 }
0940 
0941 G4double DetectorConstruction::GetFFilterZ()
0942 {
0943     return tubeFFFirstFaceZ;
0944 }
0945 
0946 G4double DetectorConstruction::GetVoxelDepthDim()
0947 {
0948     return voxelDepthDim;
0949 }
0950 
0951 G4int DetectorConstruction::GetNumberSideCells() const
0952 {
0953     return nSideCells;
0954 }
0955 
0956 G4int DetectorConstruction::GetNumberDepthCells() const
0957 {
0958     return nDepthCells;
0959 }
0960 
0961 G4int DetectorConstruction::GetPhantomDepth() const
0962 {
0963     return phantomSideDim;
0964 }
0965 
0966 void DetectorConstruction::SetJaws(G4double value)
0967 {
0968     jawAperture = value;
0969 }
0970 
0971 void DetectorConstruction::SetTargetPosition(G4double value)
0972 {
0973     sourceToSkinDistance = value;
0974 }
0975 
0976 void DetectorConstruction::SetPhantomSide(G4double value)
0977 {
0978     phantomSideDim = value;
0979 }
0980 
0981 void DetectorConstruction::SetVoxelSide(G4double value)
0982 {
0983     voxelSideDim = value;
0984 }
0985 
0986 void DetectorConstruction::SetVoxelDepth(G4double value)
0987 {
0988     voxelDepthDim = value;
0989 }
0990 
0991 void DetectorConstruction::UpdateGeometry(G4String string, G4double value)
0992 {
0993 
0994     if ( string == "fieldSide" )
0995     {
0996         //
0997         G4double halfFieldSide = value/2;
0998         G4ThreeVector position;
0999         G4RotationMatrix* rMatrix1X = new G4RotationMatrix();
1000 
1001         rMatrix1X -> rotateY(std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1002         position = *rMatrix1X * jaw1XInitialPos;
1003         boxJaw1X_phys -> SetTranslation(position);
1004         *rMatrix1X = rMatrix1X -> inverse();
1005         boxJaw1X_phys -> SetRotation(rMatrix1X);
1006 
1007         G4RotationMatrix* rMatrix2X = new G4RotationMatrix();
1008         rMatrix2X -> rotateY(-std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1009 
1010         position = *rMatrix2X * jaw2XInitialPos;
1011         boxJaw2X_phys -> SetTranslation(position);
1012         *rMatrix2X = rMatrix2X -> inverse();
1013         boxJaw2X_phys -> SetRotation(rMatrix2X);
1014 
1015         G4RotationMatrix* rMatrix1Y = new G4RotationMatrix();
1016         rMatrix1Y -> rotateX(-std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1017         position = *rMatrix1Y * jaw1YInitialPos;
1018         boxJaw1Y_phys -> SetTranslation(position);
1019         *rMatrix1Y = rMatrix1Y -> inverse();
1020         boxJaw1Y_phys -> SetRotation(rMatrix1Y);
1021 
1022         G4RotationMatrix* rMatrix2Y = new G4RotationMatrix();
1023         rMatrix2Y -> rotateX(std::fabs(std::atan(halfFieldSide/-sourceToSkinDistance)));
1024 
1025         position = *rMatrix2Y * jaw2YInitialPos;
1026         boxJaw2Y_phys -> SetTranslation(position);
1027         *rMatrix2Y = rMatrix2Y -> inverse();
1028         boxJaw2Y_phys -> SetRotation(rMatrix2Y);
1029     }
1030     else if ( string == "sourceToSkinDistance")
1031     {
1032         accWorld_phys -> SetTranslation({0., 0., value});
1033     }
1034     else if (string == "phantomSide")
1035     {
1036         G4Box* box = (G4Box *) phantom_phys -> GetLogicalVolume() -> GetSolid();
1037 
1038         // change phantom side
1039         box -> SetXHalfLength(phantomSideDim/2.);
1040         box -> SetYHalfLength(phantomSideDim/2.);
1041         box -> SetZHalfLength(phantomSideDim/2.);
1042         phantom_phys -> SetTranslation({0., 0., phantomSideDim/2.});
1043 
1044         // clear daughters and rebuild them
1045         phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1046         PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1047     }
1048     else if (string == "voxelSide")
1049     {
1050         phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1051         PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1052     }
1053     else if (string == "voxelDepth")
1054     {
1055         phantom_phys -> GetLogicalVolume() -> ClearDaughters();
1056         PhantomSegmentation(phantom_phys -> GetLogicalVolume());
1057     }
1058     else
1059         G4cerr << "*** DetectorMessenger::UpdateGeometry: Command not found" << G4endl;
1060 
1061     G4RunManager::GetRunManager()->GeometryHasBeenModified();
1062 }
1063 
1064
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!