4 // Author: I. Hrivnacova
6 // Class TG4GeometryServices
7 // -------------------------
8 // See the class description in the header file.
10 #include "TG4GeometryServices.h"
11 #include "TG4Limits.h"
12 #include "TG4IntMap.h"
13 #include "TG4NameMap.h"
14 #include "TG4G3Units.h"
15 #include "TG4G3ControlVector.h"
16 #include "TG4Globals.h"
18 #include <G4LogicalVolumeStore.hh>
19 #include <G4LogicalVolume.hh>
20 #include <G4Material.hh>
21 #include <G4Element.hh>
22 #include <G4UserLimits.hh>
24 #include <G3EleTable.hh>
25 #include <g4std/vector>
29 TG4GeometryServices* TG4GeometryServices::fgInstance = 0;
30 const G4double TG4GeometryServices::fgkAZTolerance = 0.001;
31 const G4double TG4GeometryServices::fgkDensityTolerance = 0.005;
33 //_____________________________________________________________________________
34 TG4GeometryServices::TG4GeometryServices(TG4IntMap* mediumMap,
36 : fMediumMap(mediumMap),
41 TG4Globals::Exception(
42 "TG4GeometryServices: attempt to create two instances of singleton.");
48 //_____________________________________________________________________________
49 TG4GeometryServices::TG4GeometryServices() {
51 TG4Globals::Exception(
52 "TG4GeometryServices default constructor is protected.");
56 //_____________________________________________________________________________
57 TG4GeometryServices::TG4GeometryServices(const TG4GeometryServices& right) {
59 TG4Globals::Exception(
60 "Attempt to copy TG4GeometryServices singleton.");
64 //_____________________________________________________________________________
65 TG4GeometryServices::~TG4GeometryServices() {
69 //=============================================================================
73 //=============================================================================
75 //_____________________________________________________________________________
77 TG4GeometryServices::operator=(const TG4GeometryServices& right)
79 // check assignement to self
80 if (this == &right) return *this;
82 TG4Globals::Exception(
83 "Attempt to assign TG4GeometryServices singleton.");
88 //=============================================================================
92 //=============================================================================
94 //_____________________________________________________________________________
95 G4bool TG4GeometryServices::IsG3Volume(const G4String& lvName) const
97 // Returns true if the logical volume of given volumeName
98 // was not created by Gsposp method with a generic name
102 if (lvName.contains(gSeparator))
108 //_____________________________________________________________________________
109 G4bool TG4GeometryServices::CompareElement(G4double a, G4double z,
110 const G4Element* element) const
112 // Compares given parameters with those of a given element,
113 // returns true if they are equal, false otherwise.
116 G4double ae = element->GetA()/TG4G3Units::AtomicWeight();
117 G4double ze = element->GetZ();
119 // g3tog4 can redefine A
126 ax = G3Ele.GetEle(z)->GetA()/TG4G3Units::AtomicWeight();
128 if (abs(ax - ae) < fgkAZTolerance && abs(z - ze) < fgkAZTolerance)
134 //_____________________________________________________________________________
135 G4bool TG4GeometryServices::CompareMaterial(G4int nofElements, G4double density,
136 const G4Material* material) const
138 // Compares given density with those of a given material,
139 // returns true if they are equal, false otherwise.
142 G4double dm = material->GetDensity()/TG4G3Units::MassDensity();
143 G4int ne = material->GetNumberOfElements();
145 // density percentual difference
146 G4double diff = abs(density - dm)/(density + dm)*2.;
148 if (nofElements == ne && diff < fgkDensityTolerance)
154 //_____________________________________________________________________________
155 G4double* TG4GeometryServices::ConvertAtomWeight(G4int nmat,
156 G4double* a, G4double* wmat) const
158 // In case of proportions given in atom counts (nmat<0),
159 // the wmat[i] are converted to weight fractions.
160 // (From g3tog4 G4gsmixt.)
161 // The new array has to be delete by client.
164 G4double* weight = new G4double[abs(nmat)];
169 for (i=0; i<abs(nmat); i++) {
170 // total molecular weight
171 aMol += wmat[i]*a[i];
174 G4String text = "TG4GeometryServices::ConvertAtomWeight:\n";
175 text = text + " Total molecular weight = 0";
176 TG4Globals::Warning(text);
178 for (G4int i=0; i<abs(nmat); i++) {
180 weight[i] = wmat[i]*a[i]/aMol;
184 for (G4int j=0; j<nmat; j++) weight[j] = wmat[j];
189 //=============================================================================
193 //=============================================================================
195 //_____________________________________________________________________________
196 G4double* TG4GeometryServices::CreateG4doubleArray(Float_t* array,
199 // Converts Float_t* array to G4double*,
200 // !! The new array has to be deleted by user.
203 G4double* doubleArray;
205 doubleArray = new G4double[size];
206 for (G4int i=0; i<size; i++) doubleArray[i] = array[i];
214 //_____________________________________________________________________________
215 G4String TG4GeometryServices::CutName(const char* name) const
217 // Removes spaces after the name if present.
220 G4String cutName = name;
221 G4int i = cutName.length();
222 while (cutName(--i) == ' ') cutName = cutName(0,i);
227 //_____________________________________________________________________________
228 G4String TG4GeometryServices::CutMaterialName(const char* name) const
230 // Removes the $ with precedent spaces at the name if present.
233 G4String cutName = name;
234 cutName = cutName.substr(0,cutName.find('$'));
236 return CutName(cutName);
239 //_____________________________________________________________________________
240 G4String TG4GeometryServices::G4ToG3VolumeName(const G4String& name) const
242 // Cuts _copyNo extension added to logical volume name in case
243 // the logical volume was created by Gsposp method.
246 G4String cutName = name;
247 if (cutName.contains(gSeparator))
248 cutName = cutName(0,cutName.first(gSeparator));
253 //_____________________________________________________________________________
254 G4String TG4GeometryServices::GenerateLimitsName(G4int id,
255 const G4String& medName,
256 const G4String& matName) const
258 // Generate unique name for user limits composed from the tracking medium id,
259 // name and its material name.
262 TG4Globals::AppendNumberToString(name, id);
263 name = name + "__med_" + medName + "__mat_" + matName;
268 //_____________________________________________________________________________
269 G4Material* TG4GeometryServices::MixMaterials(G4String name, G4double density,
270 const TG4StringVector& matNames,
271 const TG4doubleVector& matWeights)
273 // Creates a mixture of selected materials
276 // number of materials to be mixed
277 G4int nofMaterials = matNames.size();
278 if (nofMaterials != matWeights.size()) {
279 G4String text = "TG4GeometryServices::MixMaterials: ";
280 text = text + "different number of material names and weigths.";
281 TG4Globals::Exception(text);
284 // G4cout << "Nof of materials to be mixed: " << nofMaterials << G4endl;
286 // fill vector of materials
287 G4std::vector<G4Material*> matVector;
289 for (im=0; im< nofMaterials; im++) {
291 G4Material* material = G4Material::GetMaterial(matNames[im]);
292 matVector.push_back(material);
295 // create the mixed material
296 G4Material* mixture = new G4Material(name, density, nofMaterials);
297 for (im=0; im< nofMaterials; im++) {
298 G4Material* material = matVector[im];
299 G4double fraction = matWeights[im];
300 mixture->AddMaterial(material, fraction);
306 //_____________________________________________________________________________
307 void TG4GeometryServices::PrintNameMap() const
309 // Prints the map of volumes names to second names.
312 fNameMap->PrintAll();
315 //_____________________________________________________________________________
316 void TG4GeometryServices::PrintLimits(const G4String& name) const
318 // Finds the limits with the specified name and prints them.
321 TG4Limits* limits = FindLimits(name, true);
323 if (limits) limits->Print();
326 //_____________________________________________________________________________
327 void TG4GeometryServices::PrintVolumeLimits(const G4String& volumeName) const
329 // Finds a logical volume with the specified name and prints
333 G4LogicalVolume* lv = FindLogicalVolume(volumeName, false);
336 TG4Limits* limits = GetLimits(lv->GetUserLimits());
337 G4cout << volumeName << " ";
341 G4cout << "has not the limits set." << G4endl;
345 //_____________________________________________________________________________
346 void TG4GeometryServices::PrintStatistics(G4bool open, G4bool close) const
348 // Print G4 geometry statistics
352 if (open) TG4Globals::PrintStars(true);
354 G4cout << " GEANT4 Geometry statistics: " << G4endl
355 << " " << NofG4LogicalVolumes()
356 << " logical volumes" << G4endl
358 << NofG4PhysicalVolumes()
359 << " physical volumes" << G4endl
361 << G4Material::GetNumberOfMaterials()
362 << " materials" << G4endl
364 << TG4Limits::GetNofLimits()
365 << " user limits" << G4endl;
367 if (close) TG4Globals::PrintStars(false);
370 //_____________________________________________________________________________
371 Int_t TG4GeometryServices::NofG3Volumes() const
373 // Returns the total number of logical volumes corresponding
375 // The logical volume that were created by Gsposp method
376 // with a generic name (name_copyNo) are NOT included.
379 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
382 for (G4int i=0; i<lvStore->size(); i++) {
383 G4LogicalVolume* lv = (*lvStore)[i];
384 if (IsG3Volume(lv->GetName())) counter++;
390 //_____________________________________________________________________________
391 Int_t TG4GeometryServices::NofG4LogicalVolumes() const
393 // Returns the total number of logical volumes in the geometry.
396 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
397 return lvStore->size();
400 //_____________________________________________________________________________
401 Int_t TG4GeometryServices::NofG4PhysicalVolumes() const
403 // Returns the total number of physical volumes in the geometry.
406 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
409 for (G4int i=0; i<lvStore->size(); i++) {
410 counter += ((*lvStore)[i])->GetNoDaughters();
417 //______________________________________________________________________________
418 G4bool TG4GeometryServices::IsSpecialControls() const
420 // Returns true if a process control in some limits instance is set.
423 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
425 for (G4int i=0; i<lvStore->size(); i++) {
426 TG4Limits* limits = GetLimits((*lvStore)[i]->GetUserLimits());
427 if (limits && limits->IsControl()) return true;
433 //_____________________________________________________________________________
434 TG4Limits* TG4GeometryServices::GetLimits(G4UserLimits* limits) const
436 // Checks and converts type of the given limits.
439 if (!limits) return 0;
441 TG4Limits* tg4Limits = dynamic_cast<TG4Limits*> (limits);
444 G4Exception("TG4GeometryServices::GetLimits: Wrong limits type.");
451 //_____________________________________________________________________________
452 const G4String& TG4GeometryServices::GetMapSecond(const G4String& name)
454 // Returns the second string associated with the name in
458 return fNameMap->GetSecond(name);
461 //_____________________________________________________________________________
463 TG4GeometryServices::FindLogicalVolume(const G4String& name, G4bool silent) const
465 // Finds a logical volume with the specified name in G4LogicalVolumeStore.
468 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
470 for (G4int i=0; i<lvStore->size(); i++) {
471 G4LogicalVolume* lv = (*lvStore)[i];
472 if (lv->GetName() == name) return lv;
476 G4String text = "TG4GeometryServices:: FindLogicalVolume:\n";
477 text = text + " Logical volume " + name + " does not exist.";
478 TG4Globals::Warning(text);
484 //_____________________________________________________________________________
486 TG4GeometryServices::FindLimits(const G4String& name, G4bool silent) const
488 // Finds limits with the specified name.
491 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
493 for (G4int i=0; i<lvStore->size(); i++) {
494 G4LogicalVolume* lv = (*lvStore)[i];
495 TG4Limits* limits = GetLimits(lv->GetUserLimits());
496 if (limits && limits->GetName() == name) return limits;
500 G4String text = "TG4GeometryServices:: FindLimits:\n";
501 text = text + " Logical volume " + name + " does not exist.";
502 TG4Globals::Warning(text);
507 //_____________________________________________________________________________
508 G4int TG4GeometryServices::GetMediumId(G4LogicalVolume* lv) const
510 // Returns the second index for materials (having its origin in
511 // G4 tracking media concept)
514 return fMediumMap->GetSecond(lv->GetName());
517 //_____________________________________________________________________________
518 G4double TG4GeometryServices::GetEffA(G4Material* material) const
520 // Returns A or effective A=sum(pi*Ai) (if compound/mixture)
521 // of given material.
525 G4int nofElements = material->GetNumberOfElements();
526 if (nofElements > 1) {
527 G4String text = "Effective A for material mixture (";
528 text = text + material->GetName();
529 text = text + ") is used.";
530 //TG4Globals::Warning(text);
532 for (G4int i=0; i<nofElements; i++) {
533 G4double aOfElement = material->GetElement(i)->GetA();
534 G4double massFraction = material->GetFractionVector()[i];
535 a += aOfElement*massFraction /(TG4G3Units::AtomicWeight());
539 a = material->GetA();
540 a /= TG4G3Units::AtomicWeight();
545 //_____________________________________________________________________________
546 G4double TG4GeometryServices::GetEffZ(G4Material* material) const
548 // Returns Z or effective Z=sum(pi*Zi) (if compound/mixture)
549 // of given material.
553 G4int nofElements = material->GetNumberOfElements();
554 if (nofElements > 1) {
555 G4String text = "Effective Z for material mixture (";
556 text = text + material->GetName();
557 text = text + ") is used.";
558 //TG4Globals::Warning(text);
560 for (G4int i=0; i<nofElements; i++) {
561 G4double zOfElement = material->GetElement(i)->GetZ();
562 G4double massFraction = material->GetFractionVector()[i];
563 z += zOfElement*massFraction;
567 z = material->GetZ();
572 //_____________________________________________________________________________
573 G4Material* TG4GeometryServices::FindMaterial(G4double a, G4double z,
574 G4double density) const
576 // Finds material in G4MaterialTable with specified parameters,
577 // returns 0 if not found.
580 const G4MaterialTable* kpMatTable = G4Material::GetMaterialTable();
582 for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
584 G4Material* material = (*kpMatTable)[i];
586 if (CompareElement(a, z, material->GetElement(0)) &&
587 CompareMaterial(1, density, material))
595 //_____________________________________________________________________________
596 G4Material* TG4GeometryServices::FindMaterial(G4double* a, G4double* z,
598 G4int nmat, G4double* wmat) const
600 // Finds material in G4MaterialTable with specified parameters,
601 // returns 0 if not found.
604 G4double* weight = ConvertAtomWeight(nmat, a, wmat);
606 // loop over materials
607 G4Material* found = 0;
608 for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
610 G4Material* material = (*G4Material::GetMaterialTable())[i];
611 G4int nm = material->GetNumberOfElements();
613 if (CompareMaterial(nm, density, material)) {
615 // loop over elements
617 for (G4int ie=0; ie<nm; ie++) {
619 G4double we = (material->GetFractionVector())[ie];
621 if (!CompareElement(a[ie], z[ie], material->GetElement(ie)) ||
622 abs(weight[ie] - we) > fgkAZTolerance) {