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