Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-04 08:07:06 

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 /// \file PhysGeoImport.cc
0027 /// \brief Implementation of the PhysGeoImport class
0028 ///
0029 /// Implementation of the plasmid load methods for the geometry
0030 
0031 // Authors: J. Naoki D. Kondo (UCSF, US) : 10/10/2021
0032 //          J. Ramos-Mendez and B. Faddegon (UCSF, US)
0033 //
0034 
0035 #include "PhysGeoImport.hh"
0036 
0037 #include "G4DNAChemistryManager.hh"
0038 #include "G4Ellipsoid.hh"
0039 #include "G4VPhysicalVolume.hh"
0040 #include "Randomize.hh"
0041 
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0043 
0044 PhysGeoImport::PhysGeoImport() {}
0045 
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0047 
0048 G4LogicalVolume* PhysGeoImport::CreateLogicVolumeXYZ(G4String fileName)
0049 {
0050   G4NistManager* man = G4NistManager::Instance();
0051   fEnvelopeWater = man->FindOrBuildMaterial("G4_WATER");
0052   fpWater =
0053     man->BuildMaterialWithNewDensity("G4_WATER_MODIFIED", "G4_WATER", 1.0 * g / cm / cm / cm);
0054 
0055   ReadFile(fileName);
0056 
0057   G4double des1 = 1.1344640137963142;
0058   G4double des2 = des1 + (CLHEP::pi * .5);
0059   G4double ang = 0.6283185307179586;
0060   G4double bet1 = 0.6283185307179586 * 2;
0061   G4double posi = 1.0471975511965976;
0062   G4double sep = .1 * angstrom;
0063   // Geometries Sizes
0064   G4double PxRs = 2.9389169420478556 * angstrom;
0065   G4double PyRs = 2.9389169420478556 * angstrom;
0066   G4double PzRs = 2.9389169420478556 * angstrom;
0067   G4double PxYs = 2.7 * angstrom;
0068   G4double PyYs = 2.7 * angstrom;
0069   G4double PzYs = 2.7 * angstrom;
0070   G4double PxBp = 2.45 * angstrom;
0071   G4double PyBp = 2.45 * angstrom;
0072   G4double PzBp = 2.45 * angstrom;
0073 
0074   G4double xin = -170 * angstrom;
0075   G4double yin = -170 * angstrom;
0076   G4double zin = -170 * angstrom;
0077   G4double xfn = 170 * angstrom;
0078   G4double yfn = 170 * angstrom;
0079   G4double zfn = 170 * angstrom;
0080 
0081   G4int nVertex = fVertexes.size();
0082 
0083   // Envelope
0084 
0085   std::string boxNameSolid = fGeoName + "_solid";
0086   G4Box* box_solid =
0087     new G4Box(boxNameSolid, 0.5 * (fXMax - fXMin) + 0.5 * 3.4 * nm,
0088               0.5 * (fYMax - fYMin) + 0.5 * 3.4 * nm, 0.5 * (fZMax - fZMin) + 0.5 * 3.4 * nm);
0089 
0090   G4String boxNameLogic = fGeoName + "_logic";
0091   G4LogicalVolume* box_logic =
0092     new G4LogicalVolume(box_solid, fEnvelopeWater, boxNameLogic, 0, 0, 0);
0093 
0094   // Desoxyribose
0095 
0096   G4Ellipsoid* RSolidSugar = new G4Ellipsoid("sdeoxyribose", PxRs, PyRs, PzRs, -PzRs, .445 * PzRs);
0097   G4LogicalVolume* RSugar = new G4LogicalVolume(RSolidSugar, fpWater, "ldeoxyribose", 0, 0, 0);
0098   G4VisAttributes* MyVisAtt_Rs = new G4VisAttributes(G4Colour(G4Colour::Red()));
0099   MyVisAtt_Rs->SetForceSolid(true);
0100   RSugar->SetVisAttributes(MyVisAtt_Rs);
0101 
0102   // Phosphoric Acid
0103 
0104   G4Ellipsoid* YSolidSugar = new G4Ellipsoid("sphosphate", PxYs, PyYs, PzYs, -PzYs, .9 * angstrom);
0105 
0106   G4LogicalVolume* YSugar = new G4LogicalVolume(YSolidSugar, fpWater, "lphosphate", 0, 0, 0);
0107   G4VisAttributes* MyVisAtt_Ys = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
0108   MyVisAtt_Ys->SetForceSolid(true);
0109   YSugar->SetVisAttributes(MyVisAtt_Ys);
0110 
0111   // Base Pairs
0112 
0113   G4Ellipsoid* Base1a = new G4Ellipsoid("BasePair1a", PxBp, PyBp, PzBp, -PzBp, 1.15 * angstrom);
0114   G4LogicalVolume* BaseP1a = new G4LogicalVolume(Base1a, fpWater, "BasePair1a", 0, 0, 0);
0115   G4VisAttributes* MyVisAtt_Bp1a = new G4VisAttributes(G4Colour(G4Colour::Green()));
0116   MyVisAtt_Bp1a->SetForceSolid(true);
0117   BaseP1a->SetVisAttributes(MyVisAtt_Bp1a);
0118 
0119   G4Ellipsoid* Base1b = new G4Ellipsoid("BasePair1b", PxBp, PyBp, PzBp, 1.15 * angstrom, PzBp);
0120   G4LogicalVolume* BaseP1b = new G4LogicalVolume(Base1b, fpWater, "BasePair1b", 0, 0, 0);
0121 
0122   G4VisAttributes* MyVisAtt_Bp1b = new G4VisAttributes(G4Colour(G4Colour::Green()));
0123   MyVisAtt_Bp1b->SetForceSolid(true);
0124   BaseP1b->SetVisAttributes(MyVisAtt_Bp1b);
0125 
0126   G4int index = 0;
0127   G4double cAngle = 0;
0128   G4double pi = CLHEP::pi;
0129   for (int vertex = 0; vertex < nVertex - 1; vertex++) {
0130     xin = fVertexes[vertex][0] - fOffsetX;
0131     yin = fVertexes[vertex][1] - fOffsetY;
0132     zin = fVertexes[vertex][2] - fOffsetZ;
0133     xfn = fVertexes[vertex + 1][0] - fOffsetX;
0134     yfn = fVertexes[vertex + 1][1] - fOffsetY;
0135     zfn = fVertexes[vertex + 1][2] - fOffsetZ;
0136 
0137     G4double phi0 =
0138       std::atan2(zfn - zin, std::sqrt(((xfn - xin) * (xfn - xin)) + ((yfn - yin) * (yfn - yin))));
0139     G4double theta0 = std::atan2(xfn - xin, yfn - yin);
0140     G4double lenght = std::sqrt(((xfn - xin) * (xfn - xin)) + ((yfn - yin) * (yfn - yin))
0141                                 + ((zfn - zin) * (zfn - zin)));
0142     G4double dl = 1.0 / (lenght / (3.4 * angstrom));
0143     G4int nChain = (fVertexes[vertex] - fVertexes[vertex + 1]).mag() / (0.34 * nm);
0144 
0145     for (G4int nseg = 0; nseg < nChain; nseg++) {
0146       cAngle += ang;
0147 
0148       G4double theta = cAngle;
0149       G4double x1 = 0;
0150       G4double y1 = 0;
0151       G4double z1 = ((2 * PzBp) + PzRs + sep);
0152 
0153       G4double x2 = 0;
0154       G4double y2 = 0;
0155       G4double z2 = ((2 * PzBp) + PzRs + sep);
0156 
0157       G4ThreeVector plus2 = G4ThreeVector(0, 0, (.5 * PzRs) + PzYs);
0158       plus2.rotateX(-des1);
0159       plus2.rotateZ(-posi);
0160 
0161       G4ThreeVector plus2alt = G4ThreeVector(0, 0, (.5 * PzRs) + PzYs);
0162       plus2alt.rotateX(-des1);
0163       plus2alt.rotateZ(posi);
0164 
0165       G4double x3 = 0;
0166       G4double y3 = 0;
0167       G4double z3 = PzBp + sep;
0168 
0169       G4ThreeVector position1i = G4ThreeVector(x1, y1, z1);
0170       G4ThreeVector position2i = G4ThreeVector(x2, y2, z2) + plus2;
0171       G4ThreeVector position2ialt = G4ThreeVector(x2, y2, -z2) - plus2alt;
0172       G4ThreeVector position3i = G4ThreeVector(x3, y3, z3);
0173 
0174       position1i.rotateY(theta);
0175       position2i.rotateY(theta);
0176       position2ialt.rotateY(theta);
0177       position3i.rotateY(theta);
0178 
0179       G4double x = dl * nseg * (xfn - xin) + xin;
0180       G4double y = dl * nseg * (yfn - yin) + yin;
0181       G4double z = dl * nseg * (zfn - zin) + zin;
0182 
0183       position1i.rotateX(phi0);
0184       position2i.rotateX(phi0);
0185       position2ialt.rotateX(phi0);
0186       position3i.rotateX(phi0);
0187 
0188       position1i.rotateZ(-theta0);
0189       position2i.rotateZ(-theta0);
0190       position2ialt.rotateZ(-theta0);
0191       position3i.rotateZ(-theta0);
0192 
0193       G4double yrot1 = theta;
0194       G4double xrot1 = -des1;
0195       G4RotationMatrix rotm1 = G4RotationMatrix();
0196       rotm1.rotateX(xrot1);
0197       rotm1.rotateZ(-posi);
0198       rotm1.rotateY(yrot1);
0199       rotm1.rotateX(phi0);
0200       rotm1.rotateZ(-theta0);
0201       G4ThreeVector position1 = position1i + G4ThreeVector(x, y, z);
0202       G4Transform3D transform1(rotm1, position1);
0203 
0204       G4double yrot1alt = theta + pi;
0205       G4double xrot1alt = des1;
0206       G4RotationMatrix rotm1alt = G4RotationMatrix();
0207       rotm1alt.rotateX(xrot1alt);
0208       rotm1alt.rotateZ(-posi);
0209       rotm1alt.rotateY(yrot1alt);
0210       rotm1alt.rotateX(phi0);
0211       rotm1alt.rotateZ(-theta0);
0212       G4ThreeVector position1alt = -position1i + G4ThreeVector(x, y, z);
0213       G4Transform3D transform1alt(rotm1alt, position1alt);
0214 
0215       G4double yrot2 = theta;
0216       G4double xrot2 = -des2;
0217       G4RotationMatrix rotm2 = G4RotationMatrix();
0218       rotm2.rotateX(xrot2);
0219       rotm2.rotateY(yrot2 - bet1 + 0.8726646259971648);
0220       rotm2.rotateX(phi0);
0221       rotm2.rotateZ(-theta0);
0222       G4ThreeVector position2 = position2i + G4ThreeVector(x, y, z);
0223       G4Transform3D transform2(rotm2, position2);
0224 
0225       G4double yrot2alt = theta + pi;
0226       G4double xrot2alt = des2;
0227       G4RotationMatrix rotm2alt = G4RotationMatrix();
0228       rotm2alt.rotateX(xrot2alt);
0229       rotm2alt.rotateY(yrot2alt + bet1 - 0.8726646259971648);
0230       rotm2alt.rotateX(phi0);
0231       rotm2alt.rotateZ(-theta0);
0232       G4ThreeVector position2alt = position2ialt + G4ThreeVector(x, y, z);
0233       G4Transform3D transform2alt(rotm2alt, position2alt);
0234 
0235       G4double yrot3 = theta;
0236       G4RotationMatrix rotm3 = G4RotationMatrix();
0237       rotm3.rotateX(-pi / 2);
0238       rotm3.rotateZ(-ang);
0239       rotm3.rotateY(yrot3);
0240       rotm3.rotateX(phi0);
0241       rotm3.rotateZ(-theta0);
0242       G4ThreeVector position3 = position3i + G4ThreeVector(x, y, z);
0243       G4Transform3D transform3(rotm3, position3);
0244 
0245       G4double yrot3alt = theta + pi;
0246       G4RotationMatrix rotm3alt = G4RotationMatrix();
0247       rotm3alt.rotateX(pi / 2);
0248       rotm3alt.rotateZ(-ang);
0249       rotm3alt.rotateY(yrot3alt);
0250       rotm3alt.rotateX(phi0);
0251       rotm3alt.rotateZ(-theta0);
0252       G4ThreeVector position3alt = -position3i + G4ThreeVector(x, y, z);
0253       G4Transform3D transform3alt(rotm3alt, position3alt);
0254 
0255       new G4PVPlacement(transform1, RSugar, "deoxyribose1", box_logic, false, index, false);
0256 
0257       new G4PVPlacement(transform1alt, RSugar, "deoxyribose2", box_logic, false, index, false);
0258 
0259       new G4PVPlacement(transform3, BaseP1a, "BasePair1", box_logic, false, index, false);
0260 
0261       new G4PVPlacement(transform3alt, BaseP1a, "BasePair2", box_logic, false, index, false);
0262 
0263       new G4PVPlacement(transform2, YSugar, "phosphate1", box_logic, false, index, false);
0264 
0265       new G4PVPlacement(transform2alt, YSugar, "phosphate2", box_logic, false, index, false);
0266 
0267       G4ThreeVector Deoxy1 = position1;
0268       G4ThreeVector Deoxy2 = position1alt;
0269 
0270       fSampleDNAPositions.push_back(Deoxy1);
0271       fSampleDNAPositions.push_back(Deoxy2);
0272       fSampleDNANames.push_back("Deoxyribose");
0273       fSampleDNANames.push_back("Deoxyribose");
0274       fSampleDNADetails.push_back({-1, index, 1});
0275       fSampleDNADetails.push_back({-1, index, 2});
0276 
0277       index++;
0278     }
0279   }
0280 
0281   return box_logic;
0282 }
0283 
0284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0285 
0286 void PhysGeoImport::ReadFile(G4String fileName)
0287 {
0288   G4double x, y, z;
0289   fXMin = 1 * mm, fYMin = 1 * mm, fZMin = 1 * mm;
0290   fXMax = -1 * mm, fYMax = -1 * mm, fZMax = -1 * mm;
0291   std::ifstream plasmidFile(fileName);
0292   while (true) {
0293     plasmidFile >> x >> y >> z;
0294     if (!plasmidFile.good()) break;
0295     x *= nm;
0296     y *= nm;
0297     z *= nm;
0298     fVertexes.push_back(G4ThreeVector(x, y, z));
0299     if (fXMin > x) fXMin = x;
0300     if (fXMax < x) fXMax = x;
0301     if (fYMin > y) fYMin = y;
0302     if (fYMax < y) fYMax = y;
0303     if (fZMin > z) fZMin = z;
0304     if (fZMax < z) fZMax = z;
0305   }
0306   plasmidFile.close();
0307   fOffsetX = (fXMin + fXMax) * 0.5;
0308   fOffsetY = (fYMin + fYMax) * 0.5;
0309   fOffsetZ = (fZMin + fZMax) * 0.5;
0310 
0311   std::vector<G4ThreeVector> VertRed;
0312   for (size_t i = 0; i < fVertexes.size(); i++) {
0313     if (i % 15 == 0) VertRed.push_back(fVertexes[i]);
0314   }
0315 
0316   VertRed.push_back(fVertexes[0]);
0317 
0318   fVertexes = VertRed;
0319 }
0320 
0321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!