From 1b9adf0fbb1dc6e56fd4c8d0ff9cdf486e90b455 Mon Sep 17 00:00:00 2001
From: John William Spencer <john.william.spencer@cern.ch>
Date: Tue, 2 Feb 2021 20:11:43 +0000
Subject: [PATCH] Replace DetectorConstruction.cc

---
 src/DetectorConstruction.cc | 365 ++++++++++++++++--------------------
 1 file changed, 160 insertions(+), 205 deletions(-)

diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index fb8ed7f..8d19fab 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -29,11 +29,11 @@
 
 #include "DetectorConstruction.hh"
 #include "CalorimeterSD.hh"
-
 #include "G4Material.hh"
 #include "G4NistManager.hh"
 
 #include "G4Box.hh"
+#include "G4Trap.hh"
 #include "G4LogicalVolume.hh"
 #include "G4PVPlacement.hh"
 #include "G4PVReplica.hh"
@@ -58,8 +58,7 @@ G4GlobalMagFieldMessenger* DetectorConstruction::fMagFieldMessenger = 0;
 DetectorConstruction::DetectorConstruction()
  : G4VUserDetectorConstruction(),
    fCheckOverlaps(true),
-   fNofLayers(-1),
-   fNofPlanes(-1)
+   fNofLayers(-1)
 {
 }
 
@@ -115,7 +114,7 @@ void DetectorConstruction::DefineMaterials()
   //G4Element* elPb = new G4Element("Lead"    , "Pb", z= 82.,a=207.19*g/mole);
   //G4Element* elFe = new G4Element("Iron"    , "Fe", z= 26.,a=55.85*g/mole);
   //G4Element* elTi = new G4Element("Titanium", "Ti", z= 22.,a=44.87*g/mole);
-  //G4Element* elCa = new G4Element("Calcium",  "Ca", z= 20.,a=40.078*g/mole);
+  G4Element* elCa = new G4Element("Calcium",  "Ca", z= 20.,a=40.078*g/mole);
 
   // //http://inspirehep.net/record/641997/files/641997.pdf
   // new G4Material("EmulsionMaterial", z=8.9, a=18.2*g/mole, density=2.40*g/cm3);
@@ -147,6 +146,26 @@ void DetectorConstruction::DefineMaterials()
   H2O->AddElement(elH, natoms=2);
   H2O->AddElement(elO, natoms=1);
 
+  G4Material* SiO2
+    = new G4Material("SiO2", 
+                     density     = 2.2*g/cm3, 
+                     ncomponents = 2);
+  SiO2->AddElement(elSi, natoms=1);
+  SiO2->AddElement(elO, natoms=2);
+
+  G4Material* CaCO3
+    = new G4Material("CaCO3",
+                     density     = 2.93*g/cm3, 
+                     ncomponents = 3);
+  CaCO3->AddElement(elCa, natoms=1);
+  CaCO3->AddElement(elC, natoms=1);
+  CaCO3->AddElement(elO, natoms=3);
+  
+  G4Material* RockMaterial
+    = new G4Material("RockMaterial", density=2.5*g/cm3, ncomponents=2);
+  RockMaterial->AddMaterial(SiO2, fractionmass=0.589); // Sandstone
+  RockMaterial->AddMaterial(CaCO3, fractionmass=0.411); // Marl
+  
   G4Material* TAC
     = new G4Material("TAC", 
                      density     = 1.35*g/cm3, 
@@ -186,47 +205,43 @@ void DetectorConstruction::DefineMaterials()
 G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
 {
   // Geometry parameters
-  fNofLayers = 1000; // temporary
-  G4double absoThickness = 1.*mm;
-  G4double emulThickness = 0.05*mm;
-  G4double baseThickness = 0.20*mm;
-  G4double calorSizeX = 25*cm;
-  G4double calorSizeY = 25*cm;
+  G4double rockSizeX = 15.0*cm;
+  G4double rockSizeY = 12.5*cm;
+  G4double rockThickness = 500.*cm;
+  G4double rockCaloGap = 315.4*cm;
   
-  fNofPlanes = 3;
-  G4double sctThickness = 0.285*mm;
-  G4double baseBoardThickness = 0.380*mm;
-  G4double glueThickness = 0.025*mm;
-  G4double planeSizeXY = 24*cm; // interface detector
-  G4double gapPlane = 35.*mm;
-  G4double gapCaloSCT = 20.*mm;
+  fNofLayers = 101;
+  G4double absoThickness = 1.00*mm;
+  // fNofLayers = 120;
+  // G4double absoThickness = 0.50*mm;
   
+  G4double emulThickness = 0.05*mm;
+  G4double baseThickness = 0.20*mm;
+  G4double calorSizeX = 12.5*cm;
+  G4double calorSizeY = 10.0*cm;
+
   auto layerThickness = absoThickness + emulThickness*2 + baseThickness;
   auto calorThickness = fNofLayers * layerThickness;
-  auto planeThickness = baseBoardThickness + glueThickness*2 + sctThickness*2;
-  auto stationThickness = gapCaloSCT+gapPlane*2+planeThickness*3;
-  
-  auto worldSizeX  = 1.2 * calorSizeX;
-  auto worldSizeY  = 1.2 * calorSizeY;
-  auto worldSizeZ  = 1.2 * (calorThickness+stationThickness);
+
+  auto worldSizeX  = 1.2 * rockSizeX;
+  auto worldSizeY  = 1.2 * rockSizeY;
+  auto worldSizeZ  = 1.2 * (rockThickness+rockCaloGap+calorThickness);
+  printf("Xin: rockThickness= %f mm, rockCaloGap = %f mm, calorThickness = %f mm\n",rockThickness/mm,rockCaloGap/mm,calorThickness/mm);
   
   // Get materials
   auto defaultMaterial = G4Material::GetMaterial("Galactic");
-  auto absorberMaterial = G4Material::GetMaterial("G4_W"); //Tungsten
-  //auto absorberMaterial = G4Material::GetMaterial("G4_Pb"); //Lead
+  auto absorberMaterial = G4Material::GetMaterial("G4_Pb"); //Lead
+  //auto absorberMaterial = G4Material::GetMaterial("G4_W"); //Tungsten
   auto emulsionMaterial = G4Material::GetMaterial("EmulsionMaterial");
   auto baseMaterial = G4Material::GetMaterial("BaseMaterial");
-  auto sctMaterial = G4Material::GetMaterial("G4_Si"); //Silicon
-  auto baseBoardMaterial = G4Material::GetMaterial("TpgMaterial");
-  auto glueMaterial = G4Material::GetMaterial("Epoxy");
-  
-  if ( ! defaultMaterial || ! absorberMaterial || ! emulsionMaterial || ! baseMaterial || ! sctMaterial || ! baseBoardMaterial || ! glueMaterial ) {
+  auto rockMaterial = G4Material::GetMaterial("RockMaterial");
+
+  if ( ! defaultMaterial || ! absorberMaterial || ! emulsionMaterial || ! baseMaterial || ! rockMaterial ) {
     G4ExceptionDescription msg;
     msg << "Cannot retrieve materials already defined."; 
-    G4Exception("DetectorConstruction::DefineVolumes()",
-      "MyCode0001", FatalException, msg);
+    G4Exception("DetectorConstruction::DefineVolumes()", "MyCode0001", FatalException, msg);
   }
-   
+  
   //     
   // World
   //
@@ -250,6 +265,95 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
                  false,            // no boolean operation
                  0,                // copy number
                  fCheckOverlaps);  // checking overlaps 
+
+  //                                 
+  // Station
+  //
+  auto stationS 
+    = new G4Box("Station",           // its name
+		rockSizeX/2, rockSizeY/2, (rockThickness+rockCaloGap)/2); //its size
+                         
+  auto stationLV
+    = new G4LogicalVolume(
+			  stationS,         // its solid
+			  defaultMaterial,  // its material
+			  "Station");       // its name
+
+  new G4PVPlacement(
+		    0,                // no rotation
+		    G4ThreeVector(0., 0., -calorThickness/2),
+		    stationLV,        // its logical volume                         
+		    "Station",        // its name
+		    worldLV,          // its mother  volume
+		    false,            // no boolean operation
+		    0,                // copy number
+		    fCheckOverlaps);  // checking overlaps
+
+  //                                 
+  // rock
+  //
+  // auto rockS 
+  //   = new G4Box("Rock",           // its name
+  //                rockSizeX/2, rockSizeY/2, rockThickness/2); //its size
+
+  //http://geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomSolids.html#constructed-solid-geometry-csg-solids
+  G4double scale = 0.895;
+  auto rockS
+    = new G4Trap("Rock",
+		 rockSizeY/2, //pDz
+		 atan(rockCaloGap*scale/2/rockSizeY), //pTheta
+		 acos(-1)/2, //pPhi
+		 (rockThickness+rockCaloGap*scale)/2, //pDy1
+		 rockSizeX/2, //pDx1
+		 rockSizeX/2, //pDx2
+		 0., //pAlp1
+		 rockThickness/2, //pDy2
+		 rockSizeX/2, //pDx3
+		 rockSizeX/2, //pDx4
+		 0.); //pAlp2
+
+  G4RotationMatrix* xRot = new G4RotationMatrix();
+  xRot->rotateX(acos(-1)/2);
+
+  auto rockLV
+    = new G4LogicalVolume(
+			  rockS,           // its solid
+			  rockMaterial,    // its material
+			  "Rock");         // its name
+
+  new G4PVPlacement(
+		    //0,                // no rotation
+		    xRot,             // rotation
+		    G4ThreeVector(0., 0., -rockCaloGap*(1-scale)/2-rockCaloGap*scale/4), // its position
+		    rockLV,           // its logical volume                         
+		    "Rock",           // its name
+		    stationLV,        // its mother  volume
+		    false,            // no boolean operation
+		    0,                // copy number
+		    fCheckOverlaps);  // checking overlaps 
+  
+  //                                 
+  // gap
+  //
+  auto gapS 
+    = new G4Box("Gap",           // its name
+		rockSizeX/2, rockSizeY/2, rockCaloGap*(1-scale)/2); //its size
+
+  auto gapLV
+    = new G4LogicalVolume(
+			  gapS,             // its solid
+			  defaultMaterial,  // its material
+			  "Gap");           // its name
+
+  new G4PVPlacement(
+		    0,                // no rotation
+		    G4ThreeVector(0., 0., (rockThickness+rockCaloGap*scale)/2), // its position
+		    gapLV,            // its logical volume                         
+		    "Gap",            // its name
+		    stationLV,        // its mother  volume
+		    false,            // no boolean operation
+		    0,                // copy number
+		    fCheckOverlaps);  // checking overlaps 
   
   // 
   // Calorimeter
@@ -266,7 +370,7 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
                                    
   new G4PVPlacement(
                  0,                // no rotation
-                 G4ThreeVector(0., 0., -stationThickness/2),
+                 G4ThreeVector(0., 0., (rockThickness+rockCaloGap)/2),
                  calorLV,          // its logical volume                         
                  "Calorimeter",    // its name
                  worldLV,          // its mother  volume
@@ -292,7 +396,7 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
                  layerLV,          // its logical volume
                  calorLV,          // its mother
                  kZAxis,           // axis of replication
-                 fNofLayers,        // number of replica
+                 fNofLayers,       // number of replica
                  layerThickness);  // witdth of replica
   
   //                               
@@ -374,161 +478,6 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
                  1,                // copy number
                  fCheckOverlaps);  // checking overlaps
 
-  //                                 
-  // Station
-  //
-  auto stationS 
-    = new G4Box("Station",           // its name
-                 planeSizeXY/2, planeSizeXY/2, stationThickness/2); //its size
-                         
-  auto stationLV
-    = new G4LogicalVolume(
-                 stationS,         // its solid
-                 defaultMaterial,  // its material
-                 "Station");       // its name
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., calorThickness/2),
-                 stationLV,        // its logical volume                         
-                 "Station",        // its name
-                 worldLV,          // its mother  volume
-                 false,            // no boolean operation
-                 0,                // copy number
-                 fCheckOverlaps);  // checking overlaps
-
-  //                                 
-  // Plane
-  //
-  auto planeS 
-    = new G4Box("Plane",           // its name
-                 planeSizeXY/2, planeSizeXY/2, planeThickness/2); //its size
-
-  auto planeLV
-    = new G4LogicalVolume(
-                 planeS,           // its solid
-                 defaultMaterial,  // its material
-                 "Plane");         // its name
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., -(planeThickness+gapPlane)+0.5*gapCaloSCT), // its position
-                 planeLV,          // its logical volume                         
-                 "Plane1",         // its name
-                 stationLV,        // its mother  volume
-                 false,            // no boolean operation
-                 0,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., 0.5*gapCaloSCT), // its position
-                 planeLV,          // its logical volume                         
-                 "Plane2",         // its name
-                 stationLV,        // its mother  volume
-                 false,            // no boolean operation
-                 1,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., (planeThickness+gapPlane)+0.5*gapCaloSCT), // its position
-                 planeLV,          // its logical volume                         
-                 "Plane3",         // its name
-                 stationLV,        // its mother  volume
-                 false,            // no boolean operation
-                 2,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  //                               
-  // BaseBoard
-  //
-  auto baseBoardS 
-    = new G4Box("BaseBoard",        // its name
-                 planeSizeXY/2, planeSizeXY/2, baseBoardThickness/2); // its size
-
-  auto baseBoardLV
-    = new G4LogicalVolume(
-                 baseBoardS,        // its solid
-                 baseBoardMaterial, // its material
-                 "BaseBoardLV");    // its name
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., 0.), // its position
-                 baseBoardLV,      // its logical volume                         
-                 "BaseBoard",      // its name
-                 planeLV,          // its mother  volume
-                 false,            // no boolean operation
-                 0,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  //                               
-  // Glue
-  //
-  auto glueS 
-    = new G4Box("Glue",        // its name
-                 planeSizeXY/2, planeSizeXY/2, glueThickness/2); // its size
-                         
-  auto glueLV
-    = new G4LogicalVolume(
-                 glueS,        // its solid
-                 glueMaterial, // its material
-                 "GlueLV");    // its name
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., -(baseBoardThickness+glueThickness)/2), // its position
-                 glueLV,           // its logical volume                         
-                 "Glue1",          // its name
-                 planeLV,          // its mother  volume
-                 false,            // no boolean operation
-                 0,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., (baseBoardThickness+glueThickness)/2), // its position
-                 glueLV,           // its logical volume                         
-                 "Glue2",          // its name
-                 planeLV,          // its mother  volume
-                 false,            // no boolean operation
-                 1,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  //                               
-  // SCT
-  //
-  auto sctS 
-    = new G4Box("SCT",        // its name
-                 planeSizeXY/2, planeSizeXY/2, sctThickness/2); // its size
-                         
-  auto sctLV
-    = new G4LogicalVolume(
-                 sctS,        // its solid
-                 sctMaterial, // its material
-                 "SCTLV");    // its name
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., -(baseBoardThickness+sctThickness)/2-glueThickness), // its position
-                 sctLV,            // its logical volume                         
-                 "SCT1",           // its name
-                 planeLV,          // its mother  volume
-                 false,            // no boolean operation
-                 0,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
-  new G4PVPlacement(
-                 0,                // no rotation
-                 G4ThreeVector(0., 0., (baseBoardThickness+sctThickness)/2+glueThickness), // its position
-                 sctLV,            // its logical volume                         
-                 "SCT2",           // its name
-                 planeLV,          // its mother  volume
-                 false,            // no boolean operation
-                 1,                // copy number
-                 fCheckOverlaps);  // checking overlaps 
-
   //
   // print parameters
   //
@@ -551,30 +500,41 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
   worldLV->SetVisAttributes(G4VisAttributes::GetInvisible());
 
   //auto visAttributes = new G4VisAttributes(G4Colour(1.0,1.0,1.0)); // white
-  auto visAttributes = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)); // grey
+  auto visAttributes = new G4VisAttributes(G4Colour::White());
   //visAttributes->SetVisibility(false);
   calorLV->SetVisAttributes(visAttributes);
 
   layerLV->SetVisAttributes(G4VisAttributes::GetInvisible());
   
-  visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); // blue
-  visAttributes->SetVisibility(false);
+  //http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4Colour.html
+  visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.0));
+  visAttributes->SetColour(0.5,0.5,0.5,0.1); // Grey
+  //visAttributes->SetVisibility(false);
+  visAttributes->SetForceSolid(true);
   absorberLV->SetVisAttributes(visAttributes);
 
-  visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.0)); // yellow
-  visAttributes->SetVisibility(false);
+  visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.0));
+  visAttributes->SetColour(0.9,0.9,0.0,0.1); // Yellow
+  //visAttributes->SetVisibility(false);
+  visAttributes->SetForceSolid(true);
   baseLV->SetVisAttributes(visAttributes);
 
   visAttributes = new G4VisAttributes(G4Colour(0.8888,0.0,0.0)); // red
   visAttributes->SetVisibility(false);
   emulsionLV->SetVisAttributes(visAttributes);
 
-  stationLV->SetVisAttributes(G4VisAttributes::GetInvisible());
+  visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.0));
+  visAttributes->SetColour(0.5,0.5,0.5,0.1); // Grey
+  //visAttributes->SetVisibility(false);
+  visAttributes->SetForceSolid(true);
+  rockLV->SetVisAttributes(visAttributes);
+
+  // visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.0));
+  // visAttributes->SetColour(0.9,0.9,0.0,0.1); // Yellow
+  // //visAttributes->SetVisibility(false);
+  // visAttributes->SetForceSolid(true);
+  // gapLV->SetVisAttributes(visAttributes);
 
-  visAttributes = new G4VisAttributes(G4Colour(0.0,1.0,0.0)); // green
-  //visAttributes = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9));
-  sctLV->SetVisAttributes(visAttributes);
-  
   //
   // Always return the physical World
   //
@@ -595,11 +555,6 @@ void DetectorConstruction::ConstructSDandField()
   G4SDManager::GetSDMpointer()->AddNewDetector(emulsionSD);
   SetSensitiveDetector("EmulsionLV",emulsionSD);
 
-
-
-
-
-
   // 
   // Magnetic field
   //
-- 
GitLab