moleculardna/src/ChromosomeFactory.cc

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 #include "ChromosomeFactory.hh"
0028
0029 #include "CylindricalChromosome.hh"
0030 #include "EllipticalChromosome.hh"
0031 #include "RodChromosome.hh"
0032 #include "SphericalChromosome.hh"
0033 #include "VirtualChromosome.hh"
0034
0035 #include "G4PhysicalConstants.hh"
0036 #include "G4ThreeVector.hh"
0037 #include "G4UnitsTable.hh"
0038 #include "Randomize.hh"
0039
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0041
0042 VirtualChromosome* ChromosomeFactory::MakeChromosome(const G4String& name,
0043                                                      const std::vector<G4String>& commands)
0044 {
0045   VirtualChromosome* chromosome = nullptr;
0046   G4String chromosome_type = commands[0];
0047   if (chromosome_type == CylindricalChromosome::fShape) {
0048     // interpret the command for cylinder
0049     // expect cyl rad height x y z unit rx ry rz
0050     // rotations are in degrees
0051     if (commands.size() == 7) {
0052       G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
0053       G4double radius = std::stod(commands[1]) * unit;
0054       G4double hgt = std::stod(commands[2]) * unit;
0055       G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
0056                            std::stod(commands[5]) * unit);
0057       chromosome = new CylindricalChromosome(name, center, radius, hgt);
0058     }
0059     else if (commands.size() == 10)  // euler angles given
0060     {
0061       G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
0062       G4double radius = std::stod(commands[1]) * unit;
0063       G4double hgt = std::stod(commands[2]) * unit;
0064       G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
0065                            std::stod(commands[5]) * unit);
0066
0067       // Rotations are first X, then Y, Then Z
0068       G4RotationMatrix rot;
0069       rot.rotateX(std::stod(commands[7]) * pi / 180);
0070       rot.rotateY(std::stod(commands[8]) * pi / 180);
0071       rot.rotateZ(std::stod(commands[9]) * pi / 180);
0072
0073       chromosome = new CylindricalChromosome(name, center, radius, hgt, rot);
0074     }
0075     else {
0076       G4cout << "The arguments for a cylinder are:" << G4endl;
0077       G4cout << "1)    name cyl rad height x y z unit" << G4endl;
0078       G4cout << "2)    name cyl rad height x y z unit rx ry rz" << G4endl;
0079       G4cout << "Note that rotations are in degrees" << G4endl;
0080       InvalidReading(chromosome_type);
0081     }
0082   }
0083   else if (chromosome_type == RodChromosome::fShape) {
0084     // interpret the command for rod
0085     // expect cyl rad height x y z unit rx ry rz
0086     // rotations are in degrees
0087     if (commands.size() == 7) {
0088       G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
0089       G4double radius = std::stod(commands[1]) * unit;
0090       G4double hgt = std::stod(commands[2]) * unit;
0091       G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
0092                            std::stod(commands[5]) * unit);
0093       chromosome = new RodChromosome(name, center, radius, hgt);
0094     }
0095     else if (commands.size() == 10)  // euler angles given
0096     {
0097       G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
0098       G4double radius = std::stod(commands[1]) * unit;
0099       G4double hgt = std::stod(commands[2]) * unit;
0100       G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
0101                            std::stod(commands[5]) * unit);
0102
0103       // Rotations are first X, then Y, Then Z
0104       G4RotationMatrix rot;
0105       rot.rotateX(std::stod(commands[7]) * pi / 180);
0106       rot.rotateY(std::stod(commands[8]) * pi / 180);
0107       rot.rotateZ(std::stod(commands[9]) * pi / 180);
0108
0109       chromosome = new RodChromosome(name, center, radius, hgt, rot);
0110     }
0111     else {
0112       G4cout << "The arguments for a cylinder are:" << G4endl;
0113       G4cout << "1)    name rod rad height x y z unit" << G4endl;
0114       G4cout << "2)    name rod rad height x y z unit rx ry rz" << G4endl;
0115       G4cout << "Note that rotations are in degrees" << G4endl;
0116       InvalidReading(chromosome_type);
0117     }
0118   }
0119   else if (chromosome_type == SphericalChromosome::fShape) {
0120     // interpret the command for sphere
0121     // expect sphere rad x y z unit rx ry rz
0122     // rotations are in degrees
0123     if (commands.size() == 6) {
0124       G4double unit = G4UnitDefinition::GetValueOf(commands[5]);
0125       G4double radius = std::stod(commands[1]) * unit;
0126       G4ThreeVector center(std::stod(commands[2]) * unit, std::stod(commands[3]) * unit,
0127                            std::stod(commands[4]) * unit);
0128       chromosome = new SphericalChromosome(name, center, radius);
0129     }
0130     else if (commands.size() == 9)  // euler angles given
0131     {
0132       G4double unit = G4UnitDefinition::GetValueOf(commands[5]);
0133       G4double radius = std::stod(commands[1]) * unit;
0134       G4ThreeVector center(std::stod(commands[2]) * unit, std::stod(commands[3]) * unit,
0135                            std::stod(commands[4]) * unit);
0136
0137       // Rotations are first X, then Y, Then Z
0138       G4RotationMatrix rot;
0139       rot.rotateX(std::stod(commands[6]) * pi / 180);
0140       rot.rotateY(std::stod(commands[7]) * pi / 180);
0141       rot.rotateZ(std::stod(commands[8]) * pi / 180);
0142
0143       chromosome = new SphericalChromosome(name, center, radius, rot);
0144     }
0145     else {
0146       G4cout << "The arguments for a sphere are:" << G4endl;
0147       G4cout << "1)    name sphere rad x y z unit" << G4endl;
0148       G4cout << "2)    name sphere rad x y z unit rx ry rz" << G4endl;
0149       G4cout << "Note that rotations are in degrees" << G4endl;
0150       InvalidReading(chromosome_type);
0151     }
0152   }
0153   else if (chromosome_type == EllipticalChromosome::fShape) {
0154     // interpret the command for Ellipse
0155     // expect ellipse sx sy sz x y z unit rx ry rz
0156     // rotations are in degrees
0157     // sx sy sz are semi-major axes
0158     if (commands.size() == 8) {
0159       G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
0160       G4double sx = std::stod(commands[1]) * unit;
0161       G4double sy = std::stod(commands[2]) * unit;
0162       G4double sz = std::stod(commands[3]) * unit;
0163       G4ThreeVector center(std::stod(commands[4]) * unit, std::stod(commands[5]) * unit,
0164                            std::stod(commands[6]) * unit);
0165       chromosome = new EllipticalChromosome(name, center, sx, sy, sz);
0166     }
0167     else if (commands.size() == 11)  // euler angles given
0168     {
0169       G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
0170       G4double sx = std::stod(commands[1]) * unit;
0171       G4double sy = std::stod(commands[2]) * unit;
0172       G4double sz = std::stod(commands[3]) * unit;
0173       G4ThreeVector center(std::stod(commands[4]) * unit, std::stod(commands[5]) * unit,
0174                            std::stod(commands[6]) * unit);
0175
0176       // Rotations are first X, then Y, Then Z
0177       G4RotationMatrix rot;
0178       rot.rotateX(std::stod(commands[8]) * pi / 180);
0179       rot.rotateY(std::stod(commands[9]) * pi / 180);
0180       rot.rotateZ(std::stod(commands[10]) * pi / 180);
0181
0182       chromosome = new EllipticalChromosome(name, center, sx, sy, sz, rot);
0183     }
0184     else {
0185       G4cout << "The arguments for a ellipse are:" << G4endl;
0186       G4cout << "1)    name ellipse sx sy sz x y z unit" << G4endl;
0187       G4cout << "2)    name ellipse sx sy sz x y z unit rx ry rz" << G4endl;
0188       G4cout << "Note that rotations are in degrees" << G4endl;
0189       G4cout << "Note that dimensions (sx, sy, sz) are semi-major axes" << G4endl;
0190       InvalidReading(chromosome_type);
0191     }
0192   }
0193   else {
0194     chromosome = nullptr;
0195     InvalidReading(chromosome_type);
0196   }
0197   return chromosome;
0198 }
0199
0200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0201 void ChromosomeFactory::InvalidReading(const G4String& chromosome_type)
0202 {
0203   G4ExceptionDescription errmsg;
0204   errmsg << "Chromosome type: " << chromosome_type << " is not valid" << G4endl;
0205   G4Exception("ChromosomeFactory::MakeChromosome", "", FatalException, errmsg);
0206 }
0207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0208
0209 void ChromosomeFactory::Test()
0210 {
0211   G4cout << "------------------------------------------------------------"
0212          << "--------------------" << G4endl;
0213   G4cout << "Chromosome Test" << G4endl;
0214   // populate vector for tests
0215   std::vector<std::vector<G4String>*> tests;
0216   G4double r;
0217   for (G4int ii = 0; ii != 50; ii++) {
0218     // No rotation constructor test
0219     auto vec = new std::vector<G4String>();
0220     r = G4UniformRand();
0221     if (r < 0.25) {
0222       vec->push_back("rod");
0223       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0224       vec->push_back(std::to_string(5 * G4UniformRand() + 10));
0225     }
0226     else if (r < 0.5) {
0227       vec->push_back("sphere");
0228       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0229     }
0230     else if (r < 0.75) {
0231       vec->push_back("cyl");
0232       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0233       vec->push_back(std::to_string(5 * G4UniformRand() + 5));
0234     }
0235     else {
0236       vec->push_back("ellipse");
0237       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0238       vec->push_back(std::to_string(5 * G4UniformRand() + 5));
0239       vec->push_back(std::to_string(20 * G4UniformRand() - 6));
0240     }
0241
0242     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0243     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0244     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0245     vec->push_back("um");
0246     tests.push_back(vec);
0247   }
0248   for (G4int ii = 0; ii != 50; ++ii) {
0249     // Test with rotation
0250     auto vec = new std::vector<G4String>();
0251     r = G4UniformRand();
0252     if (r < 0.25) {
0253       vec->push_back("rod");
0254       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0255       vec->push_back(std::to_string(5 * G4UniformRand() + 10));
0256     }
0257     else if (r < 0.5) {
0258       vec->push_back("sphere");
0259       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0260     }
0261     else if (r < 0.75) {
0262       vec->push_back("cyl");
0263       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0264       vec->push_back(std::to_string(5 * G4UniformRand() + 5));
0265     }
0266     else {
0267       vec->push_back("ellipse");
0268       vec->push_back(std::to_string(10 * G4UniformRand() + 2));
0269       vec->push_back(std::to_string(5 * G4UniformRand() + 5));
0270       vec->push_back(std::to_string(20 * G4UniformRand() - 6));
0271     }
0272
0273     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0274     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0275     vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
0276
0277     vec->push_back("um");
0278
0279     vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
0280     vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
0281     vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
0282
0283     tests.push_back(vec);
0284   }
0285
0286   VirtualChromosome* chromo;
0287   for (auto& test : tests) {
0288     chromo = this->MakeChromosome("test", (*test));
0289     G4int passes = 0;
0290     G4int n = 1000;
0291     for (G4int jj = 0; jj != n; jj++) {
0292       G4ThreeVector point = chromo->RandomPointInChromosome();
0293       if (chromo->PointInChromosome(point)) {
0294         passes++;
0295       }
0296       else {
0297         G4cout << point << G4endl;
0298       }
0299     }
0300     if (passes == n) {
0301       G4cout << "Chromosome Test Passed for " << chromo->GetShape() << " based on " << n
0302              << " test points" << G4endl;
0303     }
0304     else {
0305       G4cout << "Chromosome Test Failed for " << chromo->GetShape() << " based on " << n
0306              << " test points (only " << passes << " pass)" << G4endl;
0307       chromo->Print();
0308     }
0309     delete chromo;
0310   }
0311
0312   for (auto& test : tests) {
0313     delete test;
0314   }
0315   G4cout << "------------------------------------------------------------"
0316          << "--------------------" << G4endl;
0317 }
0318
0319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......