Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-03-29 07:51:20

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 DetectorConstruction class
0028 
0029 #include "DetectorConstruction.hh"
0030 
0031 #include "DetectorMessenger.hh"
0032 
0033 #include "G4GeometryManager.hh"
0034 #include "G4LogicalVolume.hh"
0035 #include "G4LogicalVolumeStore.hh"
0036 #include "G4Material.hh"
0037 #include "G4NistManager.hh"
0038 #include "G4PVPlacement.hh"
0039 #include "G4PhysicalConstants.hh"
0040 #include "G4PhysicalVolumeStore.hh"
0041 #include "G4ReflectionFactory.hh"
0042 #include "G4RotationMatrix.hh"
0043 #include "G4SolidStore.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4Transform3D.hh"
0046 #include "G4Trd.hh"
0047 #include "G4Tubs.hh"
0048 
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050 
0051 DetectorConstruction::DetectorConstruction()
0052 {
0053   fMessenger = new DetectorMessenger(this);
0054 }
0055 
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057 
0058 DetectorConstruction::~DetectorConstruction()
0059 {
0060   delete fMessenger;
0061 }
0062 
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064 
0065 G4VPhysicalVolume* DetectorConstruction::Construct()
0066 {
0067   // Materials
0068   G4NistManager* nist = G4NistManager::Instance();
0069   G4Material* material = nist->FindOrBuildMaterial("G4_AIR");
0070 
0071   // Clean old geometry, if any
0072   //
0073   G4GeometryManager::GetInstance()->OpenGeometry();
0074   G4PhysicalVolumeStore::GetInstance()->Clean();
0075   G4LogicalVolumeStore::GetInstance()->Clean();
0076   G4SolidStore::GetInstance()->Clean();
0077   G4ReflectionFactory::Instance()->Clean();
0078 
0079   // World
0080   //
0081   G4double rmin = 0.;
0082   G4double rmax = 5 * cm;
0083   G4double hz = 5 * cm;
0084   G4double phiMin = 0.;
0085   G4double deltaPhi = 360 * degree;
0086 
0087   auto solidWorld = new G4Tubs("World",  // name
0088                                rmin, rmax, hz, phiMin, deltaPhi);  // size
0089 
0090   fWorldVolume = new G4LogicalVolume(solidWorld,  // solid
0091                                      material,  // material
0092                                      "World");  // name
0093 
0094   G4VPhysicalVolume* physiWorld = new G4PVPlacement(nullptr,  // no rotation
0095                                                     G4ThreeVector(),  // at (0,0,0)
0096                                                     fWorldVolume,  // logical volume
0097                                                     "World",  // name
0098                                                     nullptr,  // mother volume
0099                                                     false,  // no boolean operation
0100                                                     0);  // copy number
0101 
0102   // Trd volume
0103   //
0104   G4double dX1 = 1 * cm;
0105   G4double dX2 = 1 * cm;
0106   G4double dY1 = 1 * cm;
0107   G4double dY2 = 2 * cm;
0108   G4double dZ = 3 * cm;
0109 
0110   auto solidTrd = new G4Trd("trd",  // name
0111                             dX1 / 2, dX2 / 2, dY1 / 2, dY2 / 2, dZ / 2);  // size
0112 
0113   fTrdVolume = new G4LogicalVolume(solidTrd,  // solid
0114                                    material,  // material
0115                                    "trd");  // name
0116 
0117   // Place Volume1 and Volume2 according to selected methods
0118   //
0119   switch (fMethod) {
0120     case kWithDirectMatrix:
0121       PlaceWithDirectMatrix();
0122       break;
0123     case kWithInverseMatrix:
0124       PlaceWithInverseMatrix();
0125       break;
0126     case kWithAxialRotations:
0127       PlaceWithAxialRotations();
0128       break;
0129     case kWithEulerAngles:
0130       PlaceWithEulerAngles();
0131       break;
0132     case kWithReflections:
0133       PlaceWithReflections();
0134       break;
0135     default:;
0136       ;
0137   }
0138 
0139   // Return the root volume
0140   //
0141   return physiWorld;
0142 }
0143 
0144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0145 
0146 void DetectorConstruction::PlaceWithDirectMatrix()
0147 {
0148   G4double og = 3 * cm;
0149 
0150   // 1st position
0151   //
0152   G4double phi = 30 * deg;
0153   // u, v, w are the daughter axes, projected on the mother frame
0154   G4ThreeVector u = G4ThreeVector(0, 0, -1);
0155   G4ThreeVector v = G4ThreeVector(-std::sin(phi), std::cos(phi), 0.);
0156   G4ThreeVector w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0157   G4RotationMatrix rotm1 = G4RotationMatrix(u, v, w);
0158   G4cout << "\n --> phi = " << phi / deg << " deg;  direct rotation matrix : ";
0159   rotm1.print(G4cout);
0160   G4ThreeVector position1 = og * w;
0161   G4Transform3D transform1 = G4Transform3D(rotm1, position1);
0162 
0163   new G4PVPlacement(transform1,  // position, rotation
0164                     fTrdVolume,  // logical volume
0165                     "Trd",  // name
0166                     fWorldVolume,  // mother volume
0167                     false,  // no boolean operation
0168                     1);  // copy number
0169 
0170   // 2nd position
0171   //
0172   phi = phi + 90 * deg;
0173   v = G4ThreeVector(-std::sin(phi), std::cos(phi), 0.);
0174   w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0175   G4RotationMatrix rotm2 = G4RotationMatrix(u, v, w);
0176   G4ThreeVector position2 = og * w;
0177   G4Transform3D transform2 = G4Transform3D(rotm2, position2);
0178   new G4PVPlacement(transform2,  // position, rotation
0179                     fTrdVolume,  // logical volume
0180                     "Trd",  // name
0181                     fWorldVolume,  // mother volume
0182                     false,  // no boolean operation
0183                     2);  // copy number
0184 }
0185 
0186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0187 
0188 void DetectorConstruction::PlaceWithInverseMatrix()
0189 {
0190   G4double og = 3 * cm;
0191 
0192   // 1st position
0193   //
0194   G4double phi = 30 * deg;
0195   // u, v, w are the daughter axes, projected on the mother frame
0196   G4ThreeVector u = G4ThreeVector(0, 0, -1);
0197   G4ThreeVector v = G4ThreeVector(-std::sin(phi), std::cos(phi), 0.);
0198   G4ThreeVector w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0199   G4RotationMatrix rotm1 = G4RotationMatrix(u, v, w);
0200   auto rotm1Inv = new G4RotationMatrix(rotm1.inverse());
0201   G4cout << "\n --> phi = " << phi / deg << " deg;  inverse rotation matrix : ";
0202   rotm1Inv->print(G4cout);
0203   G4ThreeVector position1 = og * w;
0204 
0205   new G4PVPlacement(rotm1Inv, position1,
0206                     fTrdVolume,  // logical volume
0207                     "Trd",  // name
0208                     fWorldVolume,  // mother volume
0209                     false,  // no boolean operation
0210                     1);  // copy number
0211 
0212   // 2nd position
0213   //
0214   phi = phi + 90 * deg;
0215   v = G4ThreeVector(-std::sin(phi), std::cos(phi), 0.);
0216   w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0217   G4RotationMatrix rotm2 = G4RotationMatrix(u, v, w);
0218   auto rotm2Inv = new G4RotationMatrix(rotm2.inverse());
0219   G4ThreeVector position2 = og * w;
0220 
0221   new G4PVPlacement(rotm2Inv,  // rotation
0222                     position2,  // position
0223                     fTrdVolume,  // logical volume
0224                     "Trd",  // name
0225                     fWorldVolume,  // mother volume
0226                     false,  // no boolean operation
0227                     2);  // copy number
0228 }
0229 
0230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0231 
0232 void DetectorConstruction::PlaceWithAxialRotations()
0233 {
0234   G4double og = 3 * cm;
0235 
0236   // 1st position (with first G4PVPlacement constructor)
0237   //
0238   G4double phi = 30 * deg, theta = 90 * deg;
0239   G4ThreeVector rotAxis = G4ThreeVector(std::sin(theta - pi / 2), 0., std::cos(theta - pi / 2));
0240   G4RotationMatrix rotm1 = G4RotationMatrix();
0241   rotm1.rotateY(theta);
0242   rotm1.rotate(phi, rotAxis);
0243   G4cout << "\n --> direct rotation matrix : "
0244          << " theta = " << theta / deg << " deg;"
0245          << " phi  = " << phi / deg << " deg;";
0246   rotm1.print(G4cout);
0247   G4ThreeVector w =
0248     G4ThreeVector(std::sin(theta) * std::cos(phi), std::sin(phi), std::cos(theta) * std::cos(phi));
0249   G4ThreeVector position1 = og * w;
0250   G4Transform3D transform1(rotm1, position1);
0251 
0252   new G4PVPlacement(transform1,  // rotation,position
0253                     fTrdVolume,  // logical volume
0254                     "Trd",  // name
0255                     fWorldVolume,  // mother volume
0256                     false,  // no boolean operation
0257                     1);  // copy number
0258 
0259   // 2nd position (with second G4PVPlacement constructor)
0260   //
0261   phi = phi + 90 * deg;
0262   // rotm2Inv could be calculated with rotm2.inverse()
0263   // but also by the following :
0264   auto rotm2Inv = new G4RotationMatrix();
0265   rotm2Inv->rotate(-phi, rotAxis);
0266   rotm2Inv->rotateY(-theta);
0267   w =
0268     G4ThreeVector(std::sin(theta) * std::cos(phi), std::sin(phi), std::cos(theta) * std::cos(phi));
0269   G4ThreeVector position2 = og * w;
0270 
0271   new G4PVPlacement(rotm2Inv,  // rotation
0272                     position2,  // position
0273                     fTrdVolume,  // logical volume
0274                     "Trd",  // name
0275                     fWorldVolume,  // mother volume
0276                     false,  // no boolean operation
0277                     2);  // copy number
0278 }
0279 
0280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0281 
0282 void DetectorConstruction::PlaceWithEulerAngles()
0283 {
0284   // definitions : mother frame = {x,y,z} ; daughter frame = {u,v,w}
0285   //  n = node line = intercept of xy and uv planes
0286   //  phi_euler   = (x,n) : precession
0287   //  theta_euler = (z,w) : nutation
0288   //  psi_euler   = (n,u) : proper rotation
0289 
0290   G4double og = 3 * cm;
0291 
0292   // 1st position (with first G4PVPlacement constructor)
0293   //
0294   G4double phi = 30 * deg;
0295   G4double phi_euler = phi + pi / 2;
0296   G4double theta_euler = 90 * deg;
0297   G4double psi_euler = -90 * deg;
0298   // attention : clhep Euler constructor build inverse matrix !
0299   G4RotationMatrix rotm1Inv = G4RotationMatrix(phi_euler, theta_euler, psi_euler);
0300   G4RotationMatrix rotm1 = rotm1Inv.inverse();
0301   // remark : could be built as rotm1 = G4RotationMatrix(-psi, -theta, -phi)
0302   G4cout << "\n --> phi = " << phi / deg << " deg;  direct rotation matrix : ";
0303   rotm1.print(G4cout);
0304   G4ThreeVector w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0305   G4ThreeVector position1 = og * w;
0306   G4Transform3D transform1 = G4Transform3D(rotm1, position1);
0307 
0308   new G4PVPlacement(transform1,  // position, rotation
0309                     fTrdVolume,  // logical volume
0310                     "Trd",  // name
0311                     fWorldVolume,  // mother volume
0312                     false,  // no boolean operation
0313                     1);  // copy number
0314 
0315   // 2nd position (with second G4PVPlacement constructor)
0316   //
0317   phi = phi + 90 * deg;
0318 
0319   phi_euler = phi + pi / 2;
0320   auto rotm2Inv = new G4RotationMatrix(phi_euler, theta_euler, psi_euler);
0321   w = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0322   G4ThreeVector position2 = og * w;
0323 
0324   new G4PVPlacement(rotm2Inv,  // rotation
0325                     position2,  // position
0326                     fTrdVolume,  // logical volume
0327                     "Trd",  // name
0328                     fWorldVolume,  // mother volume
0329                     false,  // no boolean operation
0330                     2);  // copy number
0331 }
0332 
0333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0334 
0335 void DetectorConstruction::PlaceWithReflections()
0336 {
0337   /// Placement with reflections.
0338   /// In order to better show the reflection symmetry we do not apply
0339   /// the rotation along Y axis.
0340 
0341   G4double og = 3 * cm;
0342 
0343   // Place first two positionz in z = + 3cm
0344   //
0345 
0346   // 1st position
0347   G4double phi = 30 * deg;
0348   G4RotationMatrix rotm1;
0349   // rotm1.rotateY(90*deg);
0350   rotm1.rotateZ(phi);
0351   G4ThreeVector uz = G4ThreeVector(std::cos(phi), std::sin(phi), 0);
0352   G4ThreeVector position = og * uz;
0353   G4Transform3D transform1(rotm1, position);
0354   G4Transform3D translateZ = HepGeom::Translate3D(0, 0, 3. * cm);
0355 
0356   new G4PVPlacement(translateZ * transform1,  // rotation,position
0357                     fTrdVolume,  // logical volume
0358                     "Trd",  // name
0359                     fWorldVolume,  // mother volume
0360                     false,  // no boolean operation
0361                     1);  // copy number
0362 
0363   // 2nd position
0364   phi = phi + pi / 2;
0365   G4RotationMatrix rotm2;
0366   // rotm2.rotateY(90*deg);
0367   rotm2.rotateZ(phi);
0368   uz = G4ThreeVector(std::cos(phi), std::sin(phi), 0.);
0369   position = og * uz;
0370   G4Transform3D transform2 = G4Transform3D(rotm2, position);
0371 
0372   new G4PVPlacement(translateZ * transform2,  // rotation, position
0373                     fTrdVolume,  // logical volume
0374                     "Trd",  // name
0375                     fWorldVolume,  // mother volume
0376                     false,  // no boolean operation
0377                     2);  // copy number
0378 
0379   // Place next two positionz in z = - 3cm with reflection
0380   //
0381 
0382   // 3rd position
0383   translateZ = HepGeom::Translate3D(0, 0, -3. * cm);
0384   G4Transform3D reflect3D = HepGeom::ReflectZ3D();
0385 
0386   G4ReflectionFactory::Instance()->Place(translateZ * transform1 * reflect3D,  // rotation,position
0387                                          "Trd",  // name
0388                                          fTrdVolume,  // logical volume
0389                                          fWorldVolume,  // mother volume
0390                                          false,  // no boolean operation
0391                                          3);  // copy number
0392 
0393   // 4rd position
0394   G4ReflectionFactory::Instance()->Place(translateZ * transform2 * reflect3D,  // rotation,position
0395                                          "Trd",  // name
0396                                          fTrdVolume,  // logical volume
0397                                          fWorldVolume,  // mother volume
0398                                          false,  // no boolean operation
0399                                          4);  // copy number
0400 }
0401 
0402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0403 
0404 #include "G4RunManager.hh"
0405 
0406 void DetectorConstruction::SetMethod(EMethod method)
0407 {
0408   fMethod = method;
0409   G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
0410 }
0411 
0412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......