Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:22:12

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 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // and the DNA geometry given in the Geom_DNA example
0029 // shall cite the following Geant4-DNA collaboration publications:
0030 // [1] NIM B 298 (2013) 47-54
0031 // [2] Med. Phys. 37 (2010) 4692-4708
0032 // The Geant4-DNA web site is available at http://geant4-dna.org
0033 //
0034 /// \file DetectorConstruction.cc
0035 /// \brief Implementation of the DetectorConstruction class
0036 
0037 #include "DetectorConstruction.hh"
0038 
0039 // Geant4
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0078 
0079 DetectorConstruction::DetectorConstruction()
0080   : G4VUserDetectorConstruction(), fBuildChromatineFiber(true), fBuildBases(false)
0081 {}
0082 
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084 
0085 DetectorConstruction::~DetectorConstruction() {}
0086 
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0088 
0089 G4VPhysicalVolume* DetectorConstruction::Construct()
0090 {
0091   DefineMaterials();
0092   return ConstructDetector();
0093 }
0094 
0095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0096 
0097 void DetectorConstruction::DefineMaterials()
0098 {
0099   // Water is defined from NIST material database
0100   G4NistManager* man = G4NistManager::Instance();
0101   man->FindOrBuildMaterial("G4_WATER");
0102 }
0103 
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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   //                               World
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   //                             Box nucleus
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   //                            Cell nucleus
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   //                        Chromosomes territories
0167   /****************************************************************************/
0168   // NOTE: The only supported values for the rotation are
0169   // 0 and 90 degrees on the Y axis.
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},  // 5 right
0180     {-1.407, 0, -1.450, 1.407, 1.407, 1.407, 0},  // 5 left
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},  // 10 right
0190     {-4.148, 0, 1.278, 1.278, 1.278, 1.278, 0},  // 10 left
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},  // 15 right
0200     {-4.065, 5.547, 0, 1.155, 1.155, 1.155, 0},  // 15 left
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},  // 20 right
0210     {-6.971, -4.432, 0, 0.995, 0.995, 0.995, 90},  // 20 left
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},  // Y
0216     {7.354, 2.605, 0, 1.330, 1.330, 1.330, 0}  // X
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      snprintf(name, countof(name), "box%d%c",
0239      (i / 2) + 1, i % 2 ? 'l' : 'r');
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   //                 Box containing the chromatin flowers
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   // Loading flower box position for each chromosome territory
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     // snprintf(name, countof(name), "chromo%d.dat", k + 1);
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     // chromatin fiber envelope
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     // Chromatin fiber position
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       // Histones
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       // Base pair
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       // Phosphodiester group
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        Phosphodiester group Position
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         //  snprintf(name, countof(name), "sugar %d", n);
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         //  snprintf(name, countof(name), "sugar %d", n);
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        Base pair Position
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       //                 Initial position of different elements
0432       /****************************************************************************/
0433       // DNA and histone positions
0434       for (int j = 0; j < 90; j++) {
0435         // DNA (bp-SP)
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         // histones
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       //                        Visualisation colors
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 }