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