[807] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 27 | // |
---|
| 28 | #include "exrdmMaterial.hh" |
---|
| 29 | #include "exrdmMaterialData.hh" |
---|
| 30 | #include "exrdmMaterialMessenger.hh" |
---|
| 31 | |
---|
| 32 | #include "globals.hh" |
---|
| 33 | #include "G4UnitsTable.hh" |
---|
| 34 | #include "G4ios.hh" |
---|
| 35 | #include <vector> |
---|
| 36 | #include <iomanip> |
---|
| 37 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 38 | // |
---|
| 39 | exrdmMaterial::exrdmMaterial () |
---|
| 40 | { |
---|
| 41 | Material.clear(); |
---|
| 42 | Element.clear(); |
---|
| 43 | Isotope.clear(); |
---|
| 44 | // some default materials vacuum (0), air (1) and aluminium (2) defined here |
---|
| 45 | // examples of vacuum |
---|
| 46 | // |
---|
| 47 | // G4double a,z; |
---|
| 48 | |
---|
| 49 | static G4bool bmat = false ; |
---|
| 50 | |
---|
| 51 | if (!bmat) { |
---|
| 52 | |
---|
| 53 | // vacuum |
---|
| 54 | G4double density = universe_mean_density; //from PhysicalConstants.h |
---|
| 55 | G4double pressure = 3.e-18*pascal; |
---|
| 56 | G4double temperature = 2.73*kelvin; |
---|
| 57 | AddMaterial("Vacuum", "H", density,"gas",temperature,pressure); |
---|
| 58 | // air |
---|
| 59 | density = 1.290*mg/cm3; |
---|
| 60 | AddMaterial("Air", "N0.78-O0.22", density, "gas"); |
---|
| 61 | |
---|
| 62 | // aluminium |
---|
| 63 | density=2.700*g/cm3 ; |
---|
| 64 | AddMaterial ("Aluminium", "Al", density,""); |
---|
| 65 | |
---|
| 66 | //silicon |
---|
| 67 | density=2.3290*g/cm3 ; |
---|
| 68 | AddMaterial ("Silicon", "Si", density,""); |
---|
| 69 | |
---|
| 70 | bmat = true; |
---|
| 71 | } |
---|
| 72 | // create commands for interactive definition of the geometry |
---|
| 73 | materialMessenger = new exrdmMaterialMessenger(this); |
---|
| 74 | } |
---|
| 75 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 76 | // |
---|
| 77 | exrdmMaterial::~exrdmMaterial () |
---|
| 78 | { |
---|
| 79 | delete materialMessenger; |
---|
| 80 | } |
---|
| 81 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 82 | // |
---|
| 83 | void exrdmMaterial::AddMaterial (G4String name, G4String formula, G4double density, |
---|
| 84 | G4String state, G4double tem, G4double pres) |
---|
| 85 | { |
---|
| 86 | G4int isotope, Z; |
---|
| 87 | size_t i; |
---|
| 88 | for (i = 0; i<Material.size(); i++) { |
---|
| 89 | if (Material[i]->GetName() == name) { |
---|
| 90 | G4cerr <<" AddMaterial : material " <<name |
---|
| 91 | <<" already exists." <<G4endl; |
---|
| 92 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 93 | return; |
---|
| 94 | } |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | char *tokenPtr1 = NULL; |
---|
| 98 | char *sname = NULL; |
---|
| 99 | G4String s, s1("0123456789"); |
---|
| 100 | G4String element, isotopename; |
---|
| 101 | G4int ncomponents, natoms; |
---|
| 102 | G4double fatoms = 0.; |
---|
| 103 | size_t ls, id=0, ll, lr; |
---|
| 104 | ncomponents = 0; |
---|
| 105 | |
---|
| 106 | sname = new char[strlen(formula)+1]; |
---|
| 107 | strcpy(sname,formula); |
---|
| 108 | tokenPtr1 = strtok(sname,"-"); |
---|
| 109 | |
---|
| 110 | while (tokenPtr1 != NULL) { |
---|
| 111 | ncomponents++; |
---|
| 112 | tokenPtr1 = strtok( NULL, "-"); |
---|
| 113 | } |
---|
| 114 | delete[] sname; |
---|
| 115 | |
---|
| 116 | G4Material* aMaterial = 0; |
---|
| 117 | // G4cout << name <<" "<< formula << " " << density/(g/cm3) << " " << tem <<" " <<pres << G4endl; |
---|
| 118 | |
---|
| 119 | if (state == "") { |
---|
| 120 | aMaterial = new G4Material(name, density, ncomponents); |
---|
| 121 | } else if (state == "solid" && tem > 0.) { |
---|
| 122 | aMaterial = new G4Material(name, density, ncomponents, |
---|
| 123 | kStateSolid, tem ); |
---|
| 124 | } else if (state == "gas" && pres > 0.) { |
---|
| 125 | aMaterial = new G4Material(name, density, ncomponents, |
---|
| 126 | kStateGas, tem, pres ); |
---|
| 127 | } |
---|
| 128 | if (aMaterial == 0) { |
---|
| 129 | G4cerr <<" AddMaterial : Name " <<name <<"." <<G4endl; |
---|
| 130 | G4cerr <<"--> Command failed." <<G4endl; |
---|
| 131 | return; |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | sname=new char[strlen(formula)+1]; |
---|
| 135 | strcpy(sname,formula); |
---|
| 136 | tokenPtr1 = strtok(sname,"-"); |
---|
| 137 | |
---|
| 138 | while (tokenPtr1 != NULL) { |
---|
| 139 | isotope = 0; |
---|
| 140 | // G4cout << tokenPtr1 << G4endl; |
---|
| 141 | s = G4String(tokenPtr1); |
---|
| 142 | ls = s.length(); |
---|
| 143 | ll = s.find("("); |
---|
| 144 | lr = s.find(")"); |
---|
| 145 | if (ll == lr) { |
---|
| 146 | id = s.find_first_of(s1); |
---|
| 147 | element = s.substr(0,id); |
---|
| 148 | |
---|
| 149 | if (element.length() == 1) element.insert(0," "); |
---|
| 150 | for (i = 0; i<110; i++) { |
---|
| 151 | if (element == ELU[i]) break; |
---|
| 152 | } |
---|
| 153 | if (i == 110) { |
---|
| 154 | for (i = 0; i<110; i++) { |
---|
| 155 | if (element == ELL[i]) break; |
---|
| 156 | } |
---|
| 157 | if (i == 110) { |
---|
| 158 | for (i = 0; i<110; i++) { |
---|
| 159 | if (element == EUU[i]) break; |
---|
| 160 | } |
---|
| 161 | } |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | if (i == 110) { |
---|
| 165 | G4cerr <<"AddMaterial : Invalid element in material formula." |
---|
| 166 | <<element <<G4endl; |
---|
| 167 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 168 | // delete aMaterial; |
---|
| 169 | // Material[NbMat] = NULL; |
---|
| 170 | return; |
---|
| 171 | } |
---|
| 172 | |
---|
| 173 | Z = i+1; |
---|
| 174 | element = ELU[i]; |
---|
| 175 | if (id == std::string::npos) { |
---|
| 176 | natoms = 1; |
---|
| 177 | } else { |
---|
| 178 | natoms = atoi((s.substr(id,ls-id)).c_str()); |
---|
| 179 | } |
---|
| 180 | if (natoms < 1) fatoms = atof((s.substr(id,ls-id)).c_str()); |
---|
| 181 | // G4cout << " Elements = " << element << G4endl; |
---|
| 182 | //G4cout << " Nb of atoms = " << natoms << G4endl; |
---|
| 183 | } else { |
---|
| 184 | element = s.substr(0,ll); |
---|
| 185 | isotope = atoi((s.substr(ll+1,lr-ll)).c_str()); |
---|
| 186 | if (element.length() == 1) element.insert(0," "); |
---|
| 187 | for (i = 0; i<110; i++) { |
---|
| 188 | if (element == ELU[i]) break; |
---|
| 189 | } |
---|
| 190 | if (i == 110) { |
---|
| 191 | for (i = 0; i<110; i++) { |
---|
| 192 | if (element == ELL[i]) break; |
---|
| 193 | } |
---|
| 194 | if (i == 110) { |
---|
| 195 | for (i = 0; i<110; i++) { |
---|
| 196 | if (element == EUU[i]) break; |
---|
| 197 | } |
---|
| 198 | } |
---|
| 199 | } |
---|
| 200 | if (i == 110) { |
---|
| 201 | G4cerr <<"AddMaterial : Invalid element in material formula." |
---|
| 202 | <<element <<G4endl; |
---|
| 203 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 204 | // delete aMaterial; |
---|
| 205 | // Material[NbMat] = NULL; |
---|
| 206 | return; |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | Z = i+1; |
---|
| 210 | element = ELU[i]; |
---|
| 211 | isotopename = element+s.substr(ll+1,lr-ll-1); |
---|
| 212 | if (lr == std::string::npos ) { |
---|
| 213 | natoms = 1; |
---|
| 214 | } else { |
---|
| 215 | natoms = atoi((s.substr(lr+1,ls-lr)).c_str()); |
---|
| 216 | } |
---|
| 217 | if (natoms < 1) fatoms = atof((s.substr(id,ls-id)).c_str()); |
---|
| 218 | if (fatoms == 0.) natoms = 1; |
---|
| 219 | // |
---|
| 220 | // G4cout << " Elements = " << element << G4endl; |
---|
| 221 | // G4cout << " Isotope Nb = " << isotope << G4endl; |
---|
| 222 | // G4cout << " Nb of atoms = " << natoms << G4endl; |
---|
| 223 | } |
---|
| 224 | if (isotope != 0) { |
---|
| 225 | if (G4Isotope::GetIsotope(isotopename) == NULL) { |
---|
| 226 | // G4Isotope* aIsotope = new G4Isotope(isotopename, Z, isotope, A[Z-1]*g/mole); |
---|
| 227 | G4Isotope* aIsotope = new G4Isotope(isotopename, Z, isotope, isotope*g/mole); |
---|
| 228 | G4Element* aElement = new G4Element(isotopename, element, 1); |
---|
| 229 | aElement->AddIsotope(aIsotope, 100.*perCent); |
---|
| 230 | Isotope.push_back(aIsotope); |
---|
| 231 | if (natoms>0) { |
---|
| 232 | aMaterial->AddElement(aElement, natoms); |
---|
| 233 | } else { |
---|
| 234 | aMaterial->AddElement(aElement, fatoms); |
---|
| 235 | } |
---|
| 236 | Element.push_back(aElement); |
---|
| 237 | } else { |
---|
| 238 | if (natoms>0) { |
---|
| 239 | aMaterial->AddElement( G4Element::GetElement(isotopename) , natoms); |
---|
| 240 | } else { |
---|
| 241 | aMaterial->AddElement( G4Element::GetElement(isotopename) , fatoms); |
---|
| 242 | } |
---|
| 243 | } |
---|
| 244 | } else { |
---|
| 245 | if ( G4Element::GetElement(element) == NULL) { |
---|
| 246 | G4Element* aElement = new G4Element(element, element, Z, A[Z-1]*g/mole); |
---|
| 247 | if (natoms>0) { |
---|
| 248 | aMaterial->AddElement(aElement, natoms); |
---|
| 249 | } else { |
---|
| 250 | aMaterial->AddElement(aElement, fatoms); |
---|
| 251 | } |
---|
| 252 | Element.push_back(aElement); |
---|
| 253 | } else { |
---|
| 254 | if (natoms>0) { |
---|
| 255 | aMaterial->AddElement( G4Element::GetElement(element) , natoms); |
---|
| 256 | } else { |
---|
| 257 | aMaterial->AddElement( G4Element::GetElement(element) , fatoms); |
---|
| 258 | } |
---|
| 259 | } |
---|
| 260 | } |
---|
| 261 | tokenPtr1 = strtok( NULL, "-"); |
---|
| 262 | // s.empty(); |
---|
| 263 | //element.erase(); |
---|
| 264 | // |
---|
| 265 | } |
---|
| 266 | |
---|
| 267 | delete[] sname; |
---|
| 268 | Material.push_back(aMaterial); |
---|
| 269 | G4cout <<" Material:" <<name <<" with formula: " <<formula <<" added! " |
---|
| 270 | <<G4endl; |
---|
| 271 | G4cout <<" Nb of Material = " <<Material.size() <<G4endl; |
---|
| 272 | G4cout <<" Nb of Isotope = " <<Isotope.size() <<G4endl; |
---|
| 273 | G4cout <<" Nb of Element = " <<Element.size() <<G4endl; |
---|
| 274 | } |
---|
| 275 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 276 | // |
---|
| 277 | void exrdmMaterial::DeleteMaterial (G4int j) |
---|
| 278 | { |
---|
| 279 | size_t i(j-1); |
---|
| 280 | if (i > Material.size()) { |
---|
| 281 | G4cerr <<"DeleteMaterial : Invalid material index " <<j <<"." <<G4endl; |
---|
| 282 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 283 | } else { |
---|
| 284 | G4cerr <<"It seems there is no mechanism in G4 for deleting a material yet!" |
---|
| 285 | <<G4endl; |
---|
| 286 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 287 | } |
---|
| 288 | } |
---|
| 289 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 290 | // |
---|
| 291 | void exrdmMaterial::DeleteMaterial (G4String ) |
---|
| 292 | { |
---|
| 293 | G4cerr <<"It seems there is no mechanism in G4 for deleting a material yet!" |
---|
| 294 | <<G4endl; |
---|
| 295 | G4cerr <<"--> Command rejected." <<G4endl; |
---|
| 296 | } |
---|
| 297 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 298 | // |
---|
| 299 | G4int exrdmMaterial::GetMaterialIndex (G4String name) |
---|
| 300 | { |
---|
| 301 | size_t i ; |
---|
| 302 | for (i = 0; i < Material.size(); i++) { |
---|
| 303 | if (Material[i]->GetName() == name) break; |
---|
| 304 | } |
---|
| 305 | G4int k = G4int(i); |
---|
| 306 | if (i == Material.size()) k = -1; |
---|
| 307 | return k; |
---|
| 308 | } |
---|
| 309 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 310 | // |
---|
| 311 | void exrdmMaterial::ListMaterial () |
---|
| 312 | { |
---|
| 313 | G4cout <<" There are" <<std::setw(3) <<Material.size() |
---|
| 314 | <<" materials defined." <<G4endl; |
---|
| 315 | for (size_t i = 0; i< Material.size(); i++) |
---|
| 316 | G4cout <<" Material Index " <<std::setw(3) <<i+1 <<" " |
---|
| 317 | <<std::setw(14) <<Material[i]->GetName() |
---|
| 318 | <<" density: " <<std::setw(6) <<std::setprecision(3) |
---|
| 319 | <<G4BestUnit(Material[i]->GetDensity(),"Volumic Mass") <<G4endl; |
---|
| 320 | G4cout <<G4endl; |
---|
| 321 | |
---|
| 322 | } |
---|
| 323 | //////////////////////////////////////////////////////////////////////////////// |
---|