Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:29:08 

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