File indexing completed on 2026-03-30 07:50:54
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
0031
0032
0033
0034
0035
0036
0037
0038 #include "DetectorConstruction.hh"
0039
0040
0041 #include "CLHEP/Units/SystemOfUnits.h"
0042 #include "ChromosomeParameterisation.hh"
0043
0044 #include "G4Box.hh"
0045 #include "G4Ellipsoid.hh"
0046 #include "G4LogicalVolume.hh"
0047 #include "G4NistManager.hh"
0048 #include "G4Orb.hh"
0049 #include "G4PVParameterised.hh"
0050 #include "G4PVPlacement.hh"
0051 #include "G4RotationMatrix.hh"
0052 #include "G4Tubs.hh"
0053 #include "G4UnionSolid.hh"
0054 #include "G4VisAttributes.hh"
0055 #include "globals.hh"
0056
0057 #define countof(x) (sizeof(x) / sizeof(x[0]))
0058
0059 using namespace std;
0060 using CLHEP::degree;
0061 using CLHEP::micrometer;
0062 using CLHEP::mm;
0063 using CLHEP::nanometer;
0064
0065 static G4VisAttributes visInvBlue(false, G4Colour(0.0, 0.0, 1.0));
0066 static G4VisAttributes visInvWhite(false, G4Colour(1.0, 1.0, 1.0));
0067 static G4VisAttributes visInvPink(false, G4Colour(1.0, 0.0, 1.0));
0068 static G4VisAttributes visInvCyan(false, G4Colour(0.0, 1.0, 1.0));
0069 static G4VisAttributes visInvRed(false, G4Colour(1.0, 0.0, 0.0));
0070 static G4VisAttributes visInvGreen(false, G4Colour(0.0, 1.0, 0.0));
0071 static G4VisAttributes visBlue(true, G4Colour(0.0, 0.0, 1.0));
0072 static G4VisAttributes visWhite(true, G4Colour(1.0, 1.0, 1.0));
0073 static G4VisAttributes visPink(true, G4Colour(1.0, 0.0, 1.0));
0074 static G4VisAttributes visCyan(true, G4Colour(0.0, 1.0, 1.0));
0075 static G4VisAttributes visRed(true, G4Colour(1.0, 0.0, 0.0));
0076 static G4VisAttributes visGreen(true, G4Colour(0.0, 1.0, 0.0));
0077
0078
0079
0080 DetectorConstruction::DetectorConstruction()
0081 : G4VUserDetectorConstruction(), fBuildChromatineFiber(true), fBuildBases(false)
0082 {}
0083
0084
0085
0086 DetectorConstruction::~DetectorConstruction() {}
0087
0088
0089
0090 G4VPhysicalVolume* DetectorConstruction::Construct()
0091 {
0092 DefineMaterials();
0093 return ConstructDetector();
0094 }
0095
0096
0097
0098 void DetectorConstruction::DefineMaterials()
0099 {
0100
0101 G4NistManager* man = G4NistManager::Instance();
0102 man->FindOrBuildMaterial("G4_WATER");
0103 }
0104
0105
0106
0107 void DetectorConstruction::LoadChromosome(const char* filename, G4VPhysicalVolume* chromBox,
0108 G4LogicalVolume* logicBoxros)
0109 {
0110 ChromosomeParameterisation* cp = new ChromosomeParameterisation(filename);
0111 new G4PVParameterised("box ros", logicBoxros, chromBox, kUndefined, cp->GetNumRosettes(), cp);
0112
0113 G4cout << filename << " done" << G4endl;
0114 }
0115
0116
0117
0118 G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
0119 {
0120 if (fBuildBases == false && fBuildChromatineFiber == false) {
0121 G4cout << "======================================================" << G4endl;
0122 G4cout << "WARNING from DetectorConstruction::ConstructDetector:" << G4endl;
0123 G4cout << "As long as the flags fBuildBases and fBuildChromatineFiber are "
0124 "false, the output root file will be empty"
0125 << G4endl;
0126 G4cout << "This is intended for fast computation, display, testing ..." << G4endl;
0127 G4cout << "======================================================" << G4endl;
0128 }
0129
0130 G4String name;
0131
0132
0133
0134
0135
0136 DefineMaterials();
0137 G4Material* waterMaterial = G4Material::GetMaterial("G4_WATER");
0138
0139 G4Box* solidWorld = new G4Box("world", 10.0 * mm, 10.0 * mm, 10.0 * mm);
0140 G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, waterMaterial, "world");
0141 G4PVPlacement* physiWorld =
0142 new G4PVPlacement(0, G4ThreeVector(), "world", logicWorld, 0, false, 0);
0143 logicWorld->SetVisAttributes(&visInvWhite);
0144
0145
0146
0147
0148
0149 G4Box* solidTin = new G4Box("tin", 13 * micrometer, 10 * micrometer, 5 * micrometer);
0150 G4LogicalVolume* logicTin = new G4LogicalVolume(solidTin, waterMaterial, "tin");
0151 G4VPhysicalVolume* physiTin =
0152 new G4PVPlacement(0, G4ThreeVector(), "tin", logicTin, physiWorld, false, 0);
0153 logicTin->SetVisAttributes(&visInvWhite);
0154
0155
0156
0157
0158
0159 G4Ellipsoid* solidNucleus =
0160 new G4Ellipsoid("nucleus", 11.82 * micrometer, 8.52 * micrometer, 3 * micrometer, 0, 0);
0161 G4LogicalVolume* logicNucleus = new G4LogicalVolume(solidNucleus, waterMaterial, "logic nucleus");
0162 G4VPhysicalVolume* physiNucleus =
0163 new G4PVPlacement(0, G4ThreeVector(), "physi nucleus", logicNucleus, physiTin, false, 0);
0164 logicNucleus->SetVisAttributes(&visPink);
0165
0166
0167
0168
0169
0170
0171 G4double chromosomePositionSizeRotation[][7] = {
0172 {4.467, 2.835, 0, 1.557, 1.557, 1.557, 90},
0173 {-4.467, 2.835, 0, 1.557, 1.557, 1.557, 0},
0174 {4.423, -2.831, 0, 1.553, 1.553, 1.553, 90},
0175 {-4.423, -2.831, 0, 1.553, 1.553, 1.553, 0},
0176 {1.455, 5.63, 0, 1.455, 1.455, 1.455, 0},
0177 {-1.455, 5.63, 0, 1.455, 1.455, 1.455, 90},
0178 {1.435, 0, 1.392, 1.435, 1.435, 1.435, 0},
0179 {-1.435, 0, 1.392, 1.435, 1.435, 1.435, 90},
0180 {1.407, 0, -1.450, 1.407, 1.407, 1.407, 90},
0181 {-1.407, 0, -1.450, 1.407, 1.407, 1.407, 0},
0182 {1.380, -5.437, 0, 1.380, 1.380, 1.380, 0},
0183 {-1.380, -5.437, 0, 1.380, 1.380, 1.380, 90},
0184 {1.347, 2.782, -1.150, 1.347, 1.347, 1.347, 90},
0185 {-1.347, 2.782, -1.150, 1.347, 1.347, 1.347, 0},
0186 {1.311, -2.746, -1.220, 1.311, 1.311, 1.311, 90},
0187 {-1.311, -2.746, -1.220, 1.311, 1.311, 1.311, 0},
0188 {7.251, -2.541, 0, 1.275, 1.275, 1.275, 0},
0189 {-6.701, 0, -0.85, 1.275, 1.275, 1.275, 90},
0190 {4.148, 0, 1.278, 1.278, 1.278, 1.278, 90},
0191 {-4.148, 0, 1.278, 1.278, 1.278, 1.278, 0},
0192 {4.147, 0, -1.277, 1.277, 1.277, 1.277, 0},
0193 {-4.147, 0, -1.277, 1.277, 1.277, 1.277, 90},
0194 {8.930, 0.006, 0, 1.272, 1.272, 1.272, 90},
0195 {-7.296, 2.547, 0, 1.272, 1.272, 1.272, 90},
0196 {1.207, -2.642, 1.298, 1.207, 1.207, 1.207, 0},
0197 {-1.207, -2.642, 1.298, 1.207, 1.207, 1.207, 90},
0198 {1.176, 2.611, 1.368, 1.176, 1.176, 1.176, 0},
0199 {-1.176, 2.611, 1.368, 1.176, 1.176, 1.176, 90},
0200 {4.065, 5.547, 0, 1.155, 1.155, 1.155, 90},
0201 {-4.065, 5.547, 0, 1.155, 1.155, 1.155, 0},
0202 {6.542, 0.159, 1.116, 1.116, 1.116, 1.116, 0},
0203 {-9.092, 0, 0, 1.116, 1.116, 1.116, 0},
0204 {6.507, 0.159, -1.081, 1.081, 1.081, 1.081, 90},
0205 {-7.057, -2.356, 0, 1.081, 1.081, 1.081, 90},
0206 {3.824, -5.448, 0, 1.064, 1.064, 1.064, 90},
0207 {-3.824, -5.448, 0, 1.064, 1.064, 1.064, 0},
0208 {5.883, -5.379, 0, 0.995, 0.995, 0.995, 0},
0209 {-9.133, -2.111, 0, 0.995, 0.995, 0.995, 0},
0210 {6.215, 5.387, 0, 0.995, 0.995, 0.995, 0},
0211 {-6.971, -4.432, 0, 0.995, 0.995, 0.995, 90},
0212 {9.583, 2.177, 0, 0.899, 0.899, 0.899, 90},
0213 {-9.467, 2.03, 0, 0.899, 0.899, 0.899, 0},
0214 {9.440, -2.180, 0, 0.914, 0.914, 0.914, 90},
0215 {-6.34, 0, 1.339, 0.914, 0.914, 0.914, 0},
0216 {-6.947, 4.742, 0, 0.923, 0.923, 0.923, 90},
0217 {7.354, 2.605, 0, 1.330, 1.330, 1.330, 0}
0218 };
0219
0220 G4RotationMatrix* rotch = new G4RotationMatrix;
0221 rotch->rotateY(90 * degree);
0222
0223 vector<G4VPhysicalVolume*> physiBox(48);
0224
0225 for (unsigned int i = 0; i < countof(chromosomePositionSizeRotation); i++) {
0226 G4double* p = &chromosomePositionSizeRotation[i][0];
0227 G4double* size = &chromosomePositionSizeRotation[i][3];
0228 G4double rotation = chromosomePositionSizeRotation[i][6];
0229 G4ThreeVector pos(p[0] * micrometer, p[1] * micrometer, p[2] * micrometer);
0230 G4RotationMatrix* rot = rotation == 0 ? 0 : rotch;
0231
0232 ostringstream ss;
0233 ss << "box" << (i / 2) + 1 << (i % 2 ? 'l' : 'r');
0234 name = ss.str();
0235 ss.str("");
0236 ss.clear();
0237
0238
0239
0240
0241
0242 G4Box* solidBox =
0243 new G4Box(name, size[0] * micrometer, size[1] * micrometer, size[2] * micrometer);
0244 G4LogicalVolume* logicBox = new G4LogicalVolume(solidBox, waterMaterial, name);
0245 physiBox[i] = new G4PVPlacement(rot, pos, "chromo", logicBox, physiNucleus, false, 0);
0246 logicBox->SetVisAttributes(&visBlue);
0247 }
0248
0249
0250
0251
0252
0253 G4Tubs* solidBoxros = new G4Tubs("solid box ros", 0 * nanometer, 399 * nanometer, 20 * nanometer,
0254 0 * degree, 360 * degree);
0255 G4LogicalVolume* logicBoxros = new G4LogicalVolume(solidBoxros, waterMaterial, "box ros");
0256 logicBoxros->SetVisAttributes(&visInvBlue);
0257
0258
0259
0260 for (int k = 0; k < 22; k++) {
0261 ostringstream oss;
0262 oss << "chromo" << k + 1 << ".dat";
0263 name = oss.str();
0264 oss.str("");
0265 oss.clear();
0266
0267 LoadChromosome(name.c_str(), physiBox[k * 2], logicBoxros);
0268 LoadChromosome(name.c_str(), physiBox[k * 2 + 1], logicBoxros);
0269 }
0270
0271 LoadChromosome("chromoY.dat", physiBox[44], logicBoxros);
0272 LoadChromosome("chromoX.dat", physiBox[45], logicBoxros);
0273
0274
0275 if (fBuildChromatineFiber) {
0276
0277 G4Tubs* solidEnv = new G4Tubs("chromatin fiber", 0, 15.4 * nanometer, 80.5 * nanometer,
0278 0 * degree, 360 * degree);
0279 G4LogicalVolume* logicEnv = new G4LogicalVolume(solidEnv, waterMaterial, "LV chromatin fiber");
0280 logicEnv->SetVisAttributes(&visInvPink);
0281
0282
0283 for (G4int i = 0; i < 7; i++) {
0284 G4RotationMatrix* rotFiber = new G4RotationMatrix;
0285 rotFiber->rotateX(90 * degree);
0286 rotFiber->rotateY(i * 25.72 * degree);
0287 G4ThreeVector posFiber = G4ThreeVector(0, 152 * nanometer, 0);
0288 posFiber.rotateZ(i * 25.72 * degree);
0289 new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
0290
0291 rotFiber = new G4RotationMatrix;
0292 rotFiber->rotateX(90 * degree);
0293 rotFiber->rotateY((7 + i) * 25.72 * degree);
0294 posFiber = G4ThreeVector(0, 152 * nanometer, 0);
0295 posFiber.rotateZ((7 + i) * 25.72 * degree);
0296 new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
0297
0298 rotFiber = new G4RotationMatrix;
0299 rotFiber->rotateX(90 * degree);
0300 rotFiber->rotateY((25.72 + (i - 14) * 51.43) * degree);
0301 posFiber = G4ThreeVector(-36.5 * nanometer, 312 * nanometer, 0);
0302 posFiber.rotateZ((i - 14) * 51.43 * degree);
0303 new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
0304
0305 rotFiber = new G4RotationMatrix;
0306 rotFiber->rotateX(90 * degree);
0307 rotFiber->rotateY(180 * degree);
0308 rotFiber->rotateY((i - 21) * 51.43 * degree);
0309 posFiber = G4ThreeVector(-103 * nanometer, 297 * nanometer, 0);
0310 posFiber.rotateZ((i - 21) * 51.43 * degree);
0311 new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
0312 }
0313
0314 if (fBuildBases) {
0315
0316 G4Tubs* solidHistone = new G4Tubs("solid histone", 0, 3.25 * nanometer, 2.85 * nanometer,
0317 0 * degree, 360 * degree);
0318 G4LogicalVolume* logicHistone =
0319 new G4LogicalVolume(solidHistone, waterMaterial, "logic histone");
0320
0321
0322 G4Orb* solidBp1 = new G4Orb("blue sphere", 0.17 * nanometer);
0323 G4LogicalVolume* logicBp1 = new G4LogicalVolume(solidBp1, waterMaterial, "logic blue sphere");
0324 G4Orb* solidBp2 = new G4Orb("pink sphere", 0.17 * nanometer);
0325 G4LogicalVolume* logicBp2 = new G4LogicalVolume(solidBp2, waterMaterial, "logic pink sphere");
0326
0327
0328
0329 G4Orb* solidSugar_48em1_nm = new G4Orb("sugar", 0.48 * nanometer);
0330
0331 G4ThreeVector posi(0.180248 * nanometer, 0.32422 * nanometer, 0.00784 * nanometer);
0332 G4UnionSolid* uniDNA =
0333 new G4UnionSolid("move", solidSugar_48em1_nm, solidSugar_48em1_nm, 0, posi);
0334
0335 G4ThreeVector posi2(-0.128248 * nanometer, 0.41227 * nanometer, 0.03584 * nanometer);
0336 G4UnionSolid* uniDNA2 =
0337 new G4UnionSolid("move2", solidSugar_48em1_nm, solidSugar_48em1_nm, 0, posi2);
0338
0339
0340
0341
0342
0343 for (G4int n = 2; n < 200; n++) {
0344 G4double SP1[2][3] = {
0345 {(-0.6 * nanometer) * cos(n * 0.26), 0, (0.6 * nanometer) * sin(n * 0.26)},
0346 {(0.6 * nanometer) * cos(n * 0.26), 0, (-0.6 * nanometer) * sin(0.26 * n)}};
0347 G4double matriceSP1[3][3] = {
0348 {cos(n * 0.076), -sin(n * 0.076), 0}, {sin(n * 0.076), cos(n * 0.076), 0}, {0, 0, 1}};
0349 G4double matriceSP2[2][3];
0350
0351 for (G4int i = 0; i < 3; i++) {
0352 G4double sumSP1 = 0;
0353 G4double sumSP2 = 0;
0354 for (G4int j = 0; j < 3; j++) {
0355 sumSP1 += matriceSP1[i][j] * SP1[0][j];
0356 sumSP2 += matriceSP1[i][j] * SP1[1][j];
0357 }
0358 matriceSP2[0][i] = sumSP1;
0359 matriceSP2[1][i] = sumSP2;
0360 }
0361
0362 G4double heliceSP[3] = {(4.85 * nanometer) * cos(n * 0.076),
0363 (4.85 * nanometer) * sin(n * 0.076), (n * 0.026 * nanometer)};
0364
0365 for (G4int i = 0; i < 3; i++) {
0366 matriceSP2[0][i] += heliceSP[i];
0367 matriceSP2[1][i] += heliceSP[i];
0368 }
0369 G4ThreeVector posSugar1(matriceSP2[0][2], matriceSP2[0][1],
0370 (matriceSP2[0][0]) - (4.25 * nanometer));
0371 G4ThreeVector posSugar2(matriceSP2[1][2], matriceSP2[1][1],
0372 (matriceSP2[1][0]) - (5.45 * nanometer));
0373
0374 ostringstream ss;
0375 ss << "sugar_" << n;
0376 name = ss.str().c_str();
0377 ss.str("");
0378 ss.clear();
0379
0380
0381 uniDNA = new G4UnionSolid(name, uniDNA, solidSugar_48em1_nm, 0, posSugar1);
0382
0383 ss << "sugar_" << n;
0384 name = ss.str().c_str();
0385 ss.str("");
0386 ss.clear();
0387
0388
0389 uniDNA2 = new G4UnionSolid(name, uniDNA2, solidSugar_48em1_nm, 0, posSugar2);
0390 }
0391 G4LogicalVolume* logicSphere3 = new G4LogicalVolume(uniDNA, waterMaterial, "logic sugar 2");
0392 G4LogicalVolume* logicSphere4 = new G4LogicalVolume(uniDNA2, waterMaterial, "logic sugar 4");
0393
0394
0395
0396
0397 for (G4int n = 0; n < 200; n++) {
0398 G4double bp1[2][3] = {
0399 {(-0.34 * nanometer) * cos(n * 0.26), 0, (0.34 * nanometer) * sin(n * 0.26)},
0400 {(0.34 * nanometer) * cos(n * 0.26), 0, (-0.34 * nanometer) * sin(0.26 * n)}};
0401 G4double matriceBP1[3][3] = {
0402 {cos(n * 0.076), -sin(n * 0.076), 0}, {sin(n * 0.076), cos(n * 0.076), 0}, {0, 0, 1}};
0403 G4double matriceBP2[2][3];
0404
0405 for (G4int i = 0; i < 3; i++) {
0406 G4double sumBP1 = 0;
0407 G4double sumBP2 = 0;
0408 for (G4int j = 0; j < 3; j++) {
0409 sumBP1 += matriceBP1[i][j] * bp1[0][j];
0410 sumBP2 += matriceBP1[i][j] * bp1[1][j];
0411 }
0412 matriceBP2[0][i] = sumBP1;
0413 matriceBP2[1][i] = sumBP2;
0414 }
0415 G4double heliceBP[3] = {(4.8 * nanometer) * cos(n * 0.076),
0416 (4.8 * nanometer) * sin(n * 0.076), n * 0.026 * nanometer};
0417
0418 for (G4int i = 0; i < 3; i++) {
0419 matriceBP2[0][i] += heliceBP[i];
0420 matriceBP2[1][i] += heliceBP[i];
0421 }
0422 G4ThreeVector position1(matriceBP2[0][2], matriceBP2[0][1],
0423 matriceBP2[0][0] - (4.25 * nanometer));
0424 G4ThreeVector position2(matriceBP2[1][2], matriceBP2[1][1],
0425 matriceBP2[1][0] - (5.45 * nanometer));
0426
0427 new G4PVPlacement(0, position1, logicBp1, "physi blue sphere", logicSphere3, false, 0);
0428 new G4PVPlacement(0, position2, logicBp2, "physi pink sphere", logicSphere4, false, 0);
0429 }
0430
0431
0432
0433
0434
0435 for (int j = 0; j < 90; j++) {
0436
0437 G4RotationMatrix* rotStrand1 = new G4RotationMatrix;
0438 rotStrand1->rotateZ(j * -51.43 * degree);
0439 G4ThreeVector posStrand1(-2.7 * nanometer, 9.35 * nanometer,
0440 (-69.9 * nanometer) + (j * 1.67 * nanometer));
0441 posStrand1.rotateZ(j * 51.43 * degree);
0442 new G4PVPlacement(rotStrand1, posStrand1, logicSphere3, "physi sugar 2", logicEnv, false,
0443 0);
0444
0445 G4RotationMatrix* rotStrand2 = new G4RotationMatrix;
0446 rotStrand2->rotateZ(j * -51.43 * degree);
0447 G4ThreeVector posStrand2(-2.7 * nanometer, 9.35 * nanometer,
0448 (-68.7 * nanometer) + (j * 1.67 * nanometer));
0449 posStrand2.rotateZ(j * 51.43 * degree);
0450 new G4PVPlacement(rotStrand2, posStrand2, logicSphere4, "physi sugar 4", logicEnv, false,
0451 0);
0452
0453
0454 G4RotationMatrix* rotHistone = new G4RotationMatrix;
0455 rotHistone->rotateY(90 * degree);
0456 rotHistone->rotateX(j * (-51.43 * degree));
0457 G4ThreeVector posHistone(0.0, 9.35 * nanometer, (-74.15 + j * 1.67) * nanometer);
0458 posHistone.rotateZ(j * 51.43 * degree);
0459 new G4PVPlacement(rotHistone, posHistone, logicHistone, "PV histone", logicEnv, false, 0);
0460 }
0461
0462
0463
0464
0465 logicBp1->SetVisAttributes(&visInvCyan);
0466 logicBp2->SetVisAttributes(&visInvPink);
0467
0468 logicSphere3->SetVisAttributes(&visInvWhite);
0469 logicSphere4->SetVisAttributes(&visInvRed);
0470
0471 logicHistone->SetVisAttributes(&visInvBlue);
0472 }
0473 }
0474
0475 G4cout << "Geometry has been loaded" << G4endl;
0476 return physiWorld;
0477 }