File indexing completed on 2025-02-23 09:21:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
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
0051
0052 DetectorConstruction::DetectorConstruction()
0053 {
0054 fMessenger = new DetectorMessenger(this);
0055 }
0056
0057
0058
0059 DetectorConstruction::~DetectorConstruction()
0060 {
0061 delete fMessenger;
0062 }
0063
0064
0065
0066 G4VPhysicalVolume* DetectorConstruction::Construct()
0067 {
0068
0069 G4NistManager* nist = G4NistManager::Instance();
0070 G4Material* material = nist->FindOrBuildMaterial("G4_AIR");
0071
0072
0073
0074 G4GeometryManager::GetInstance()->OpenGeometry();
0075 G4PhysicalVolumeStore::GetInstance()->Clean();
0076 G4LogicalVolumeStore::GetInstance()->Clean();
0077 G4SolidStore::GetInstance()->Clean();
0078 G4ReflectionFactory::Instance()->Clean();
0079
0080
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",
0089 rmin, rmax, hz, phiMin, deltaPhi);
0090
0091 fWorldVolume = new G4LogicalVolume(solidWorld,
0092 material,
0093 "World");
0094
0095 G4VPhysicalVolume* physiWorld = new G4PVPlacement(nullptr,
0096 G4ThreeVector(),
0097 fWorldVolume,
0098 "World",
0099 nullptr,
0100 false,
0101 0);
0102
0103
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",
0112 dX1 / 2, dX2 / 2, dY1 / 2, dY2 / 2, dZ / 2);
0113
0114 fTrdVolume = new G4LogicalVolume(solidTrd,
0115 material,
0116 "trd");
0117
0118
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
0141
0142 return physiWorld;
0143 }
0144
0145
0146
0147 void DetectorConstruction::PlaceWithDirectMatrix()
0148 {
0149 G4double og = 3 * cm;
0150
0151
0152
0153 G4double phi = 30 * deg;
0154
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,
0165 fTrdVolume,
0166 "Trd",
0167 fWorldVolume,
0168 false,
0169 1);
0170
0171
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,
0180 fTrdVolume,
0181 "Trd",
0182 fWorldVolume,
0183 false,
0184 2);
0185 }
0186
0187
0188
0189 void DetectorConstruction::PlaceWithInverseMatrix()
0190 {
0191 G4double og = 3 * cm;
0192
0193
0194
0195 G4double phi = 30 * deg;
0196
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,
0208 "Trd",
0209 fWorldVolume,
0210 false,
0211 1);
0212
0213
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,
0223 position2,
0224 fTrdVolume,
0225 "Trd",
0226 fWorldVolume,
0227 false,
0228 2);
0229 }
0230
0231
0232
0233 void DetectorConstruction::PlaceWithAxialRotations()
0234 {
0235 G4double og = 3 * cm;
0236
0237
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,
0254 fTrdVolume,
0255 "Trd",
0256 fWorldVolume,
0257 false,
0258 1);
0259
0260
0261
0262 phi = phi + 90 * deg;
0263
0264
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,
0273 position2,
0274 fTrdVolume,
0275 "Trd",
0276 fWorldVolume,
0277 false,
0278 2);
0279 }
0280
0281
0282
0283 void DetectorConstruction::PlaceWithEulerAngles()
0284 {
0285
0286
0287
0288
0289
0290
0291 G4double og = 3 * cm;
0292
0293
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
0300 G4RotationMatrix rotm1Inv = G4RotationMatrix(phi_euler, theta_euler, psi_euler);
0301 G4RotationMatrix rotm1 = rotm1Inv.inverse();
0302
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,
0310 fTrdVolume,
0311 "Trd",
0312 fWorldVolume,
0313 false,
0314 1);
0315
0316
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,
0326 position2,
0327 fTrdVolume,
0328 "Trd",
0329 fWorldVolume,
0330 false,
0331 2);
0332 }
0333
0334
0335
0336 void DetectorConstruction::PlaceWithReflections()
0337 {
0338
0339
0340
0341
0342 G4double og = 3 * cm;
0343
0344
0345
0346
0347
0348 G4double phi = 30 * deg;
0349 G4RotationMatrix rotm1;
0350
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,
0358 fTrdVolume,
0359 "Trd",
0360 fWorldVolume,
0361 false,
0362 1);
0363
0364
0365 phi = phi + pi / 2;
0366 G4RotationMatrix rotm2;
0367
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,
0374 fTrdVolume,
0375 "Trd",
0376 fWorldVolume,
0377 false,
0378 2);
0379
0380
0381
0382
0383
0384 translateZ = HepGeom::Translate3D(0, 0, -3. * cm);
0385 G4Transform3D reflect3D = HepGeom::ReflectZ3D();
0386
0387 G4ReflectionFactory::Instance()->Place(translateZ * transform1 * reflect3D,
0388 "Trd",
0389 fTrdVolume,
0390 fWorldVolume,
0391 false,
0392 3);
0393
0394
0395 G4ReflectionFactory::Instance()->Place(translateZ * transform2 * reflect3D,
0396 "Trd",
0397 fTrdVolume,
0398 fWorldVolume,
0399 false,
0400 4);
0401 }
0402
0403
0404
0405 #include "G4RunManager.hh"
0406
0407 void DetectorConstruction::SetMethod(EMethod method)
0408 {
0409 fMethod = method;
0410 G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
0411 }
0412
0413