B5/src/DetectorConstruction.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 //
0027 /// \file B5/src/DetectorConstruction.cc
0028 /// \brief Implementation of the B5::DetectorConstruction class
0029
0030 #include "DetectorConstruction.hh"
0031
0032 #include "CellParameterisation.hh"
0033 #include "Constants.hh"
0034 #include "DriftChamberSD.hh"
0035 #include "EmCalorimeterSD.hh"
0036 #include "HadCalorimeterSD.hh"
0037 #include "HodoscopeSD.hh"
0038 #include "MagneticField.hh"
0039
0040 #include "G4Box.hh"
0041 #include "G4Colour.hh"
0042 #include "G4FieldManager.hh"
0043 #include "G4GenericMessenger.hh"
0044 #include "G4LogicalVolume.hh"
0045 #include "G4Material.hh"
0046 #include "G4NistManager.hh"
0047 #include "G4PVParameterised.hh"
0048 #include "G4PVPlacement.hh"
0049 #include "G4PVReplica.hh"
0050 #include "G4RotationMatrix.hh"
0051 #include "G4RunManager.hh"
0052 #include "G4SDManager.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "G4ThreeVector.hh"
0055 #include "G4Tubs.hh"
0056 #include "G4UserLimits.hh"
0057 #include "G4VisAttributes.hh"
0058
0059 namespace B5
0060 {
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 G4ThreadLocal MagneticField* DetectorConstruction::fMagneticField = nullptr;
0065 G4ThreadLocal G4FieldManager* DetectorConstruction::fFieldMgr = nullptr;
0066
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068
0069 DetectorConstruction::DetectorConstruction()
0070 {
0071   fArmRotation = new G4RotationMatrix();
0072   fArmRotation->rotateY(fArmAngle);
0073
0074   // define commands for this class
0075   DefineCommands();
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079
0080 DetectorConstruction::~DetectorConstruction()
0081 {
0082   delete fArmRotation;
0083   delete fMessenger;
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 G4VPhysicalVolume* DetectorConstruction::Construct()
0089 {
0090   // Construct materials
0091   ConstructMaterials();
0092   auto air = G4Material::GetMaterial("G4_AIR");
0093   // auto argonGas = G4Material::GetMaterial("_Ar");
0094   auto argonGas = G4Material::GetMaterial("G4_Ar");
0095   auto scintillator = G4Material::GetMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
0096   auto csI = G4Material::GetMaterial("G4_CESIUM_IODIDE");
0097   auto lead = G4Material::GetMaterial("G4_Pb");
0098
0099   // Option to switch on/off checking of volumes overlaps
0100   //
0101   G4bool checkOverlaps = true;
0102
0103   // geometries --------------------------------------------------------------
0104   // experimental hall (world volume)
0105   auto worldSolid = new G4Box("worldBox", 10. * m, 3. * m, 10. * m);
0106   auto worldLogical = new G4LogicalVolume(worldSolid, air, "worldLogical");
0107   auto worldPhysical = new G4PVPlacement(nullptr, G4ThreeVector(), worldLogical, "worldPhysical",
0108                                          nullptr, false, 0, checkOverlaps);
0109
0110   // Tube with Local Magnetic field
0111
0112   auto magneticSolid = new G4Tubs("magneticTubs", 0., 1. * m, 1. * m, 0., 360. * deg);
0113
0114   fMagneticLogical = new G4LogicalVolume(magneticSolid, air, "magneticLogical");
0115
0116   // placement of Tube
0117
0118   auto fieldRot = new G4RotationMatrix();
0119   fieldRot->rotateX(90. * deg);
0120   new G4PVPlacement(fieldRot, G4ThreeVector(), fMagneticLogical, "magneticPhysical", worldLogical,
0121                     false, 0, checkOverlaps);
0122
0123   // set step limit in tube with magnetic field
0124   auto userLimits = new G4UserLimits(1 * m);
0125   fMagneticLogical->SetUserLimits(userLimits);
0126
0127   // first arm
0128   auto firstArmSolid = new G4Box("firstArmBox", 1.5 * m, 1. * m, 3. * m);
0129   auto firstArmLogical = new G4LogicalVolume(firstArmSolid, air, "firstArmLogical");
0130   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -5. * m), firstArmLogical, "firstArmPhysical",
0131                     worldLogical, false, 0, checkOverlaps);
0132
0133   // second arm
0134   auto secondArmSolid = new G4Box("secondArmBox", 2. * m, 2. * m, 3.5 * m);
0135   auto secondArmLogical = new G4LogicalVolume(secondArmSolid, air, "secondArmLogical");
0136   auto x = -5. * m * std::sin(fArmAngle);
0137   auto z = 5. * m * std::cos(fArmAngle);
0138   fSecondArmPhys = new G4PVPlacement(fArmRotation, G4ThreeVector(x, 0., z), secondArmLogical,
0139                                      "fSecondArmPhys", worldLogical, false, 0, checkOverlaps);
0140
0141   // hodoscopes in first arm
0142   auto hodoscope1Solid = new G4Box("hodoscope1Box", 5. * cm, 20. * cm, 0.5 * cm);
0143   fHodoscope1Logical = new G4LogicalVolume(hodoscope1Solid, scintillator, "hodoscope1Logical");
0144
0145   for (auto i = 0; i < kNofHodoscopes1; i++) {
0146     G4double x1 = (i - kNofHodoscopes1 / 2) * 10. * cm;
0147     new G4PVPlacement(nullptr, G4ThreeVector(x1, 0., -1.5 * m), fHodoscope1Logical,
0148                       "hodoscope1Physical", firstArmLogical, false, i, checkOverlaps);
0149   }
0150
0151   // drift chambers in first arm
0152   auto chamber1Solid = new G4Box("chamber1Box", 1. * m, 30. * cm, 1. * cm);
0153   auto chamber1Logical = new G4LogicalVolume(chamber1Solid, argonGas, "chamber1Logical");
0154
0155   for (auto i = 0; i < kNofChambers; i++) {
0156     G4double z1 = (i - kNofChambers / 2) * 0.5 * m;
0157     new G4PVPlacement(nullptr, G4ThreeVector(0., 0., z1), chamber1Logical, "chamber1Physical",
0158                       firstArmLogical, false, i, checkOverlaps);
0159   }
0160
0161   // "virtual" wire plane
0162   auto wirePlane1Solid = new G4Box("wirePlane1Box", 1. * m, 30. * cm, 0.1 * mm);
0163   fWirePlane1Logical = new G4LogicalVolume(wirePlane1Solid, argonGas, "wirePlane1Logical");
0164   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fWirePlane1Logical, "wirePlane1Physical",
0165                     chamber1Logical, false, 0, checkOverlaps);
0166
0167   // hodoscopes in second arm
0168   auto hodoscope2Solid = new G4Box("hodoscope2Box", 5. * cm, 20. * cm, 0.5 * cm);
0169   fHodoscope2Logical = new G4LogicalVolume(hodoscope2Solid, scintillator, "hodoscope2Logical");
0170
0171   for (auto i = 0; i < kNofHodoscopes2; i++) {
0172     G4double x2 = (i - kNofHodoscopes2 / 2) * 10. * cm;
0173     new G4PVPlacement(nullptr, G4ThreeVector(x2, 0., 0.), fHodoscope2Logical, "hodoscope2Physical",
0174                       secondArmLogical, false, i, checkOverlaps);
0175   }
0176
0177   // drift chambers in second arm
0178   auto chamber2Solid = new G4Box("chamber2Box", 1.5 * m, 30. * cm, 1. * cm);
0179   auto chamber2Logical = new G4LogicalVolume(chamber2Solid, argonGas, "chamber2Logical");
0180
0181   for (auto i = 0; i < kNofChambers; i++) {
0182     G4double z2 = (i - kNofChambers / 2) * 0.5 * m - 1.5 * m;
0183     new G4PVPlacement(nullptr, G4ThreeVector(0., 0., z2), chamber2Logical, "chamber2Physical",
0184                       secondArmLogical, false, i, checkOverlaps);
0185   }
0186
0187   // "virtual" wire plane
0188   auto wirePlane2Solid = new G4Box("wirePlane2Box", 1.5 * m, 30. * cm, 0.1 * mm);
0189   fWirePlane2Logical = new G4LogicalVolume(wirePlane2Solid, argonGas, "wirePlane2Logical");
0190   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fWirePlane2Logical, "wirePlane2Physical",
0191                     chamber2Logical, false, 0, checkOverlaps);
0192
0193   // CsI calorimeter
0194   auto emCalorimeterSolid = new G4Box("EMcalorimeterBox", 1.5 * m, 30. * cm, 15. * cm);
0195   auto emCalorimeterLogical = new G4LogicalVolume(emCalorimeterSolid, csI, "EMcalorimeterLogical");
0196   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 2. * m), emCalorimeterLogical,
0197                     "EMcalorimeterPhysical", secondArmLogical, false, 0, checkOverlaps);
0198
0199   // EMcalorimeter cells
0200   auto cellSolid = new G4Box("cellBox", 7.5 * cm, 7.5 * cm, 15. * cm);
0201   fCellLogical = new G4LogicalVolume(cellSolid, csI, "cellLogical");
0202   G4VPVParameterisation* cellParam = new CellParameterisation();
0203   new G4PVParameterised("cellPhysical", fCellLogical, emCalorimeterLogical, kXAxis, kNofEmCells,
0204                         cellParam);
0205
0206   // hadron calorimeter
0207   auto hadCalorimeterSolid = new G4Box("HadCalorimeterBox", 1.5 * m, 30. * cm, 50. * cm);
0208   auto hadCalorimeterLogical =
0209     new G4LogicalVolume(hadCalorimeterSolid, lead, "HadCalorimeterLogical");
0210   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 3. * m), hadCalorimeterLogical,
0211                     "HadCalorimeterPhysical", secondArmLogical, false, 0, checkOverlaps);
0212
0213   // hadron calorimeter column
0214   auto HadCalColumnSolid = new G4Box("HadCalColumnBox", 15. * cm, 30. * cm, 50. * cm);
0215   auto HadCalColumnLogical = new G4LogicalVolume(HadCalColumnSolid, lead, "HadCalColumnLogical");
0216   new G4PVReplica("HadCalColumnPhysical", HadCalColumnLogical, hadCalorimeterLogical, kXAxis,
0217                   kNofHadColumns, 30. * cm);
0218
0219   // hadron calorimeter cell
0220   auto HadCalCellSolid = new G4Box("HadCalCellBox", 15. * cm, 15. * cm, 50. * cm);
0221   auto HadCalCellLogical = new G4LogicalVolume(HadCalCellSolid, lead, "HadCalCellLogical");
0222   new G4PVReplica("HadCalCellPhysical", HadCalCellLogical, HadCalColumnLogical, kYAxis, kNofHadRows,
0223                   30. * cm);
0224
0225   // hadron calorimeter layers
0226   auto HadCalLayerSolid = new G4Box("HadCalLayerBox", 15. * cm, 15. * cm, 2.5 * cm);
0227   auto HadCalLayerLogical = new G4LogicalVolume(HadCalLayerSolid, lead, "HadCalLayerLogical");
0228   new G4PVReplica("HadCalLayerPhysical", HadCalLayerLogical, HadCalCellLogical, kZAxis,
0229                   kNofHadCells, 5. * cm);
0230
0231   // scintillator plates
0232   auto HadCalScintiSolid = new G4Box("HadCalScintiBox", 15. * cm, 15. * cm, 0.5 * cm);
0233   fHadCalScintiLogical =
0234     new G4LogicalVolume(HadCalScintiSolid, scintillator, "HadCalScintiLogical");
0235   new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 2. * cm), fHadCalScintiLogical,
0236                     "HadCalScintiPhysical", HadCalLayerLogical, false, 0, checkOverlaps);
0237
0238   // visualization attributes ------------------------------------------------
0239
0240   G4VisAttributes invisible(G4VisAttributes::GetInvisible());
0241   G4VisAttributes invisibleBlue(false, G4Colour::Blue());
0242   G4VisAttributes invisibleGreen(false, G4Colour::Green());
0243   G4VisAttributes invisibleYellow(false, G4Colour::Yellow());
0244   G4VisAttributes blue(G4Colour::Blue());
0245   G4VisAttributes cgray(G4Colour::Gray());
0246   G4VisAttributes green(G4Colour::Green());
0247   G4VisAttributes red(G4Colour::Red());
0248   G4VisAttributes yellow(G4Colour::Yellow());
0249
0250   worldLogical->SetVisAttributes(invisible);
0251   firstArmLogical->SetVisAttributes(invisible);
0252   secondArmLogical->SetVisAttributes(invisible);
0253
0254   fMagneticLogical->SetVisAttributes(cgray);
0255   fHodoscope1Logical->SetVisAttributes(red);
0256   fHodoscope2Logical->SetVisAttributes(red);
0257
0258   chamber1Logical->SetVisAttributes(green);
0259   chamber2Logical->SetVisAttributes(green);
0260   fWirePlane1Logical->SetVisAttributes(invisibleGreen);
0261   fWirePlane2Logical->SetVisAttributes(invisibleGreen);
0262
0263   emCalorimeterLogical->SetVisAttributes(invisibleYellow);
0264   fCellLogical->SetVisAttributes(yellow);
0265
0266   hadCalorimeterLogical->SetVisAttributes(blue);
0267   HadCalColumnLogical->SetVisAttributes(invisibleBlue);
0268   HadCalCellLogical->SetVisAttributes(invisibleBlue);
0269   HadCalLayerLogical->SetVisAttributes(invisibleBlue);
0270   fHadCalScintiLogical->SetVisAttributes(invisibleBlue);
0271
0272   // return the world physical volume ----------------------------------------
0273
0274   return worldPhysical;
0275 }
0276
0277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0278
0279 void DetectorConstruction::ConstructSDandField()
0280 {
0281   // sensitive detectors -----------------------------------------------------
0282   auto sdManager = G4SDManager::GetSDMpointer();
0283   G4String SDname;
0284
0285   auto hodoscope1 = new HodoscopeSD(SDname = "/hodoscope1");
0286   sdManager->AddNewDetector(hodoscope1);
0287   fHodoscope1Logical->SetSensitiveDetector(hodoscope1);
0288
0289   auto hodoscope2 = new HodoscopeSD(SDname = "/hodoscope2");
0290   sdManager->AddNewDetector(hodoscope2);
0291   fHodoscope2Logical->SetSensitiveDetector(hodoscope2);
0292
0293   auto chamber1 = new DriftChamberSD(SDname = "/chamber1");
0294   sdManager->AddNewDetector(chamber1);
0295   fWirePlane1Logical->SetSensitiveDetector(chamber1);
0296
0297   auto chamber2 = new DriftChamberSD(SDname = "/chamber2");
0298   sdManager->AddNewDetector(chamber2);
0299   fWirePlane2Logical->SetSensitiveDetector(chamber2);
0300
0301   auto emCalorimeter = new EmCalorimeterSD(SDname = "/EMcalorimeter");
0302   sdManager->AddNewDetector(emCalorimeter);
0303   fCellLogical->SetSensitiveDetector(emCalorimeter);
0304
0305   auto hadCalorimeter = new HadCalorimeterSD(SDname = "/HadCalorimeter");
0306   sdManager->AddNewDetector(hadCalorimeter);
0307   fHadCalScintiLogical->SetSensitiveDetector(hadCalorimeter);
0308
0309   // magnetic field ----------------------------------------------------------
0310   fMagneticField = new MagneticField();
0311   fFieldMgr = new G4FieldManager();
0312   fFieldMgr->SetDetectorField(fMagneticField);
0313   fFieldMgr->CreateChordFinder(fMagneticField);
0314   G4bool forceToAllDaughters = true;
0315   fMagneticLogical->SetFieldManager(fFieldMgr, forceToAllDaughters);
0316 }
0317
0318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0319
0320 void DetectorConstruction::ConstructMaterials()
0321 {
0322   auto nistManager = G4NistManager::Instance();
0323
0324   // Air
0325   nistManager->FindOrBuildMaterial("G4_AIR");
0326
0327   // Argon gas
0328   nistManager->FindOrBuildMaterial("G4_Ar");
0329   // With a density different from the one defined in NIST
0330   // G4double density = 1.782e-03*g/cm3;
0331   // nistManager->BuildMaterialWithNewDensity("_Ar","G4_Ar",density);
0332   // !! cases segmentation fault
0333
0334   // Scintillator
0335   // (PolyVinylToluene, C_9H_10)
0336   nistManager->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
0337
0338   // CsI
0339   nistManager->FindOrBuildMaterial("G4_CESIUM_IODIDE");
0340
0341   // Lead
0342   nistManager->FindOrBuildMaterial("G4_Pb");
0343
0344   // Vacuum "Galactic"
0345   // nistManager->FindOrBuildMaterial("G4_Galactic");
0346
0347   // Vacuum "Air with low density"
0348   // auto air = G4Material::GetMaterial("G4_AIR");
0349   // G4double density = 1.0e-5*air->GetDensity();
0350   // nistManager
0351   //   ->BuildMaterialWithNewDensity("Air_lowDensity", "G4_AIR", density);
0352
0353   G4cout << G4endl << "The materials defined are : " << G4endl << G4endl;
0354   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0355 }
0356
0357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0358
0359 void DetectorConstruction::SetArmAngle(G4double val)
0360 {
0361   if (!fSecondArmPhys) {
0362     G4cerr << "Detector has not yet been constructed." << G4endl;
0363     return;
0364   }
0365
0366   fArmAngle = val;
0367   *fArmRotation = G4RotationMatrix();  // make it unit vector
0368   fArmRotation->rotateY(fArmAngle);
0369   auto x = -5. * m * std::sin(fArmAngle);
0370   auto z = 5. * m * std::cos(fArmAngle);
0371   fSecondArmPhys->SetTranslation(G4ThreeVector(x, 0., z));
0372
0373   // tell G4RunManager that we change the geometry
0374   G4RunManager::GetRunManager()->GeometryHasBeenModified();
0375 }
0376
0377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0378
0379 void DetectorConstruction::DefineCommands()
0380 {
0381   // Define /B5/detector command directory using generic messenger class
0382   fMessenger = new G4GenericMessenger(this, "/B5/detector/", "Detector control");
0383
0384   // armAngle command
0385   auto& armAngleCmd = fMessenger->DeclareMethodWithUnit(
0386     "armAngle", "deg", &DetectorConstruction::SetArmAngle, "Set rotation angle of the second arm.");
0387   armAngleCmd.SetParameterName("angle", true);
0388   armAngleCmd.SetRange("angle>=0. && angle<180.");
0389   armAngleCmd.SetDefaultValue("30.");
0390 }
0391
0392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0393
0394 }  // namespace B5