Back to home page

EIC code displayed by LXR

 
 

    


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