moleculardna/include/DNAGeometry.hh

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: DNAGeometry.hh
0028 /// brief: This file handls MolecularDNA Geometry
0029
0030 #ifndef MOLECULAR_DNA_GEOMETRY_HH
0031 #define MOLECULAR_DNA_GEOMETRY_HH
0032
0033 #include "ChromosomeMapper.hh"
0034 #include "DNAWorld.hh"
0035 #include "DamageModel.hh"
0036 #include "MoleculeList.hh"
0037 #include "UtilityFunctions.hh"
0038
0039 #include "G4LogicalVolume.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4ThreeVector.hh"
0042 #include "G4VDNAMolecularGeometry.hh"
0043 #include "G4VPhysicalVolume.hh"
0044 #include "G4VTouchable.hh"
0045 #include "globals.hh"
0046
0047 #include <array>
0048 #include <map>
0049 #include <unordered_map>
0050 #include <vector>
0051
0052 class G4Material;
0053
0054 class G4MolecularConfiguration;
0055
0056 class G4VisAttributes;
0057
0058 class OctreeNode;
0059
0060 class DNAGeometryMessenger;
0061
0062 class PlacementVolumeInfo;
0063
0064 using placement = std::tuple<G4LogicalVolume*, int64_t, int64_t, G4ThreeVector,
0065                              G4ThreeVector>;  // dousatsu
0066 // G4int, G4int, G4ThreeVector, G4ThreeVector> placement;//ORG
0067
0068 // Holder to describe the global co-ordinates of each DNA chain
0069 // in a given placement volume
0070 // Format << array of global chain indices for local indices 0, 1, 2, 3>
0071 //         < array of start indices for the base pairs along each chain,
0072 //           given by GLOBAL chain index>
0073 //         < array of end indices for the base pairs along each chain,
0074 //           given by GLOBAL chain index>
0075 //         < Boolean to indicate if the strand ID should be flipped due
0076 //           to a 180 degree rotation of the DNA >
0077 //         < Boolean to indicate if volume has actually been placed >
0078 using placement_transform =
0079   std::tuple<std::array<int, 8>, std::array<int64_t, 8>, std::array<int64_t, 8>, G4bool, G4bool>;
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081
0082 struct molecule_t
0083 {
0084     G4String fname;
0085     G4String fshape;
0086     int64_t fchain_id;
0087     int64_t fstrand_id;
0088     int64_t fbase_idx;
0089     G4ThreeVector fposition;
0090     G4ThreeVector frotation;
0091     G4ThreeVector fsize;
0092 };
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095
0096 struct compareLVByName
0097 {
0098     bool operator()(const G4LogicalVolume* lhs, const G4LogicalVolume* rhs) const;
0099 };
0100
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0102
0103 class DNAGeometry : public G4VDNAMolecularGeometry
0104 {
0105   public:
0106     DNAGeometry();
0107
0108     ~DNAGeometry() override;
0109
0110     void BuildDNA(G4LogicalVolume*);
0111
0112     inline DNAWorld* GetDNAWorld();
0113
0114     void AddVoxelFile(const G4String& vname, const G4String& fname, const G4bool& twist)
0115     {
0116       fVoxelNames[vname] = fname;
0117       fVoxelTwist[vname] = twist;
0118     }
0119
0120     void SetFractalFilename(const G4String& fname)
0121     {
0122       fFractalCurveFile = fname;
0123       FillVoxelVectors();
0124     }
0125
0126     void SetFractalAnglesAsPi(const G4bool& value) { fAnglesAsPi = value; }
0127
0128     void EnableCustomMoleculeSizes(const G4bool& value) { fEnableCustomMoleculeSizes = value; }
0129
0130     void AddChangeMoleculeSize(const G4String& name, const G4ThreeVector& size);
0131
0132     inline void SetVoxelSideLength(const G4ThreeVector& length) { fVoxelSideLength = length; }
0133
0134     inline void SetFractalScaling(const G4ThreeVector& length) { fFractalScaling = length; }
0135
0136     OctreeNode* GetTopOctreeNode(G4LogicalVolume*) const;
0137
0138     const PlacementVolumeInfo* GetPVInfo(const G4LogicalVolume*) const;
0139
0140     inline auto GetChromosomeMapper() const { return fpChromosomeMapper; };
0141
0142     inline void SetOverlaps(G4bool overlap) { fCheckOverlaps = overlap; };
0143
0144     inline G4bool GetOverlaps() const { return fCheckOverlaps; };
0145
0146     inline void SetDrawCellVolumes(G4bool b) { fDrawCellVolumes = b; };
0147
0148     inline G4bool GetDrawCellVolumes() const { return fDrawCellVolumes; };
0149
0150     inline void SetVerbosity(G4int i) { fVerbosity = i; };
0151
0152     inline G4int GetVerbosity() const { return fVerbosity; };
0153
0154     inline void SetSmartless(G4int i) { fSmartless = i; };
0155
0156     inline G4int GetSmartless() const { return fSmartless; };
0157
0158     inline void SetMediumMaterial(G4Material* mat) { fpMediumMaterial = mat; };
0159
0160     inline G4LogicalVolume* GetDNAChemVolumePointer() const;
0161
0162     inline G4LogicalVolume* GetDNAPhysicsVolumePointer() const;
0163
0164     G4LogicalVolume* GetMatchingChemVolume(G4LogicalVolume*) const;
0165
0166     G4LogicalVolume* GetMatchingPhysVolume(G4LogicalVolume*) const;
0167
0168     void PrintOctreeStats();
0169
0170     inline void SetDirectInteractionRange(G4double r) { fDirectInteractionRange = r; };
0171
0172     void SetRadicalKillDistance(G4double r) { fRadicalKillDistance = r; };
0173
0174     void SetHistoneScav(G4bool hf) { fUseHistoneScav = hf; };
0175
0176     inline G4double GetDirectInteractionRange() const { return fDirectInteractionRange; };
0177
0178     G4double GetRadicalKillDistance() const { return fRadicalKillDistance; };
0179
0180     const DamageModel* GetDamageModel() { return fpDamageModel; };
0181
0182     // Tests
0183     void ChromosomeTest();
0184
0185     void BasePairIndexTest();
0186
0187     void UniqueIDTest();
0188
0189     int64_t GetGlobalUniqueIDTest(  // dousatsu
0190       int64_t, int64_t, int64_t, int64_t, int64_t) const;  // dousatsu
0191
0192   private:
0193     DNAGeometryMessenger* fpMessenger;
0194     G4String fFractalCurveFile;
0195     G4bool fAnglesAsPi = false;
0196     G4bool fEnableCustomMoleculeSizes = false;
0197     G4bool fCheckOverlaps = false;
0198     G4bool fDrawCellVolumes = false;
0199     G4int fVerbosity = 0;
0200     G4bool fUseHistoneScav = true;
0201     // vis attrs to draw cells
0202     G4VisAttributes* fpDrawCellAttrs;
0203     // Materials
0204     G4Material* fpMediumMaterial;
0205     G4Material *fpSugarMaterial, *fpPhosphateMaterial;
0206     G4Material *fpGuanineMaterial, *fpAdenineMaterial, *fpThymineMaterial;
0207     G4Material *fpCytosineMaterial, *fpHistoneMaterial;
0208
0209     // Placement of molecules
0210     void FillParameterisedSpace();
0211
0212     G4VPhysicalVolume* PlacePhosphate(G4LogicalVolume*, const molecule_t&, const molecule_t&);
0213
0214     G4VPhysicalVolume* PlaceSugar(G4LogicalVolume*, const molecule_t&, const molecule_t&);
0215
0216     G4VPhysicalVolume* PlaceBase(G4LogicalVolume*, const molecule_t&, const molecule_t&,
0217                                  const molecule_t&, const molecule_t&);
0218
0219     G4VPhysicalVolume* PlaceHistone(G4LogicalVolume*,  // dousatsu
0220                                     const molecule_t&);
0221
0222     // defining voxels
0223     std::pair<G4LogicalVolume*, PlacementVolumeInfo*> LoadVoxelVolume(const G4String&,
0224                                                                       const G4String&);
0225
0226     void FillVoxelVectors();
0227
0228     G4LogicalVolume *fpDNAVolumePhys{}, *fpDNAVolumeChem{};
0229     std::map<const G4LogicalVolume*, G4LogicalVolume*> fPhysToChem;
0230     std::map<const G4LogicalVolume*, G4LogicalVolume*> fChemToPhys;
0231     std::vector<G4ThreeVector> fVoxelPositions;
0232     std::vector<G4ThreeVector> fVoxelRotations;
0233     std::vector<G4String> fVoxelTypes;
0234     std::vector<G4int> fVoxelIndices;
0235     // std::vector<G4LogicalVolume *> fVoxelLogicalVolumes;
0236     G4ThreeVector fVoxelSideLength = G4ThreeVector(75 * nm, 75 * nm, 75 * nm);
0237     G4ThreeVector fFractalScaling = G4ThreeVector(1, 1, 1);
0238     std::map<G4String, G4String> fVoxelNames;  // For these, the first el.
0239     std::map<G4String, G4bool> fVoxelTwist;  // is also lv->GetName()
0240
0241     std::map<const G4LogicalVolume*, PlacementVolumeInfo*> fInfoMap;
0242     std::map<G4String, G4ThreeVector> fMoleculeSizes;
0243     G4double fSmartless = 2.0;
0244     G4double fRadicalKillDistance = 10 * nm, fDirectInteractionRange = 6 * angstrom;
0245     G4double fHistoneInteractionRadius = 25 * angstrom;
0246
0247     ChromosomeMapper* fpChromosomeMapper;
0248     DamageModel* fpDamageModel;
0249     DNAWorld* fpDNAPhysicsWorld;
0250
0251     // Members for placement transformations
0252   public:
0253     inline int64_t GetNumberOfPlacements() const
0254     {
0255       return (int64_t)fPlacementTransformations.size();
0256     }
0257
0258     G4int GetGlobalChain(G4int vol_idx, G4int local_chain) const
0259     {
0260       return std::get<0>(fPlacementTransformations[vol_idx])[local_chain];
0261     };
0262
0263     G4int GetLocalChain(G4int vol_idx, G4int global_chain) const;
0264
0265     inline int64_t GetStartIdx(int64_t vol_idx,
0266                                int64_t global_chn) const  // dousatsu
0267     {
0268       return std::get<1>(fPlacementTransformations[vol_idx])[global_chn];
0269     };
0270
0271     inline int64_t GetEndIdx(int64_t vol_idx,
0272                              int64_t global_chn) const  // dousatsu
0273     {
0274       return std::get<2>(fPlacementTransformations[vol_idx])[global_chn];
0275     };
0276
0277     int64_t GetMaxBPIdx() const;
0278
0279     inline G4int GetNumberOfChains() const
0280     {
0281       return std::get<0>(fPlacementTransformations[0]).size();
0282     };
0283
0284     inline G4bool GetStrandsFlipped(G4int vol_idx) const
0285     {
0286       return std::get<3>(fPlacementTransformations[vol_idx]);
0287     };
0288
0289   private:
0290     std::vector<placement_transform> fPlacementTransformations;
0291
0292     void AddNewPlacement(const G4LogicalVolume*, std::array<int, 8>, G4bool, G4bool);
0293
0294     void AddFourChainPlacement(std::vector<placement>::iterator, std::vector<placement>::iterator,
0295                                G4bool);
0296
0297     void AddSingleChainPlacement(std::vector<placement>::iterator, std::vector<placement>::iterator,
0298                                  G4bool);
0299
0300     // Lookup methods for chemistry handling
0301   public:
0302     // Unique ID generators
0303     inline int64_t GetGlobalPairID(G4int, G4int, int64_t) const;
0304
0305     int64_t GetGlobalUniqueID(G4VPhysicalVolume*,  // dousatsu
0306                               const G4VTouchable*) const;
0307
0308     // And reconstructors from unique ID
0309     int64_t GetBasePairFromUniqueID(int64_t) const;
0310
0311     static inline molecule GetMoleculeFromUniqueID(int64_t);
0312
0313     G4Material* GetMaterialFromUniqueID(int64_t) const;
0314
0315     G4int GetChainIDFromUniqueID(int64_t) const;
0316
0317     G4int GetStrandIDFromUniqueID(int64_t) const;
0318
0319     G4int GetPlacementIndexFromUniqueID(int64_t) const;
0320
0321     void FindNearbyMolecules(const G4LogicalVolume* motherLogical,
0322                              const G4ThreeVector& localPosition,
0323                              std::vector<G4VPhysicalVolume*>& out, double searchRange) override;
0324
0325     G4bool IsInsideHistone(const G4LogicalVolume* motherLogical,
0326                            const G4ThreeVector& localPosition) const;
0327 };
0328
0329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0330
0331 inline int64_t DNAGeometry::GetGlobalPairID(G4int place_idx, G4int chain_idx,
0332                                             int64_t base_idx) const
0333 {
0334   return base_idx + this->GetStartIdx(place_idx, chain_idx);
0335 }
0336
0337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0338
0339 inline molecule DNAGeometry::GetMoleculeFromUniqueID(int64_t idx)
0340 {
0341   return (molecule)(idx % ::molecule::Count);
0342 }
0343
0344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0345
0346 inline G4LogicalVolume* DNAGeometry::GetDNAPhysicsVolumePointer() const
0347 {
0348   return fpDNAVolumePhys;
0349 }
0350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0351
0352 inline DNAWorld* DNAGeometry::GetDNAWorld()
0353 {
0354   return fpDNAPhysicsWorld;
0355 }
0356
0357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0358
0359 inline G4LogicalVolume* DNAGeometry::GetDNAChemVolumePointer() const
0360 {
0361   return fpDNAVolumeChem;
0362 }
0363 #endif  // MOLECULAR_DNA_GEOMETRY_HH