source: trunk/examples/extended/radioactivedecay/exrdm/src/exrdmMaterial.cc@ 1036

Last change on this file since 1036 was 807, checked in by garnier, 17 years ago

update
File size: 10.5 KB
Line 
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//
39exrdmMaterial::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//
77exrdmMaterial::~exrdmMaterial ()
78{
79 delete materialMessenger;
80}
81////////////////////////////////////////////////////////////////////////////////
82//
83void 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//
277void 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//
291void 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//
299G4int 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//
311void 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////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.