Correction in CompareElement (to prevent from crash when vacuum (z<1) is passed.
[u/mrichter/AliRoot.git] / TGeant4 / TG4GeometryServices.cxx
CommitLineData
b68f7176 1// $Id$
2// Category: geometry
3//
e5967ab3 4// Author: I. Hrivnacova
5//
6// Class TG4GeometryServices
7// -------------------------
b68f7176 8// See the class description in the header file.
9
10#include "TG4GeometryServices.h"
e5967ab3 11#include "TG4Limits.h"
12#include "TG4IntMap.h"
13#include "TG4NameMap.h"
c4294cc7 14#include "TG4G3Units.h"
e5967ab3 15#include "TG4G3ControlVector.h"
16#include "TG4Globals.h"
b68f7176 17
b68f7176 18#include <G4LogicalVolumeStore.hh>
19#include <G4LogicalVolume.hh>
20#include <G4Material.hh>
e5967ab3 21#include <G4Element.hh>
e4a64a3a 22#include <G4UserLimits.hh>
b68f7176 23#include <G3toG4.hh>
e5967ab3 24#include <G3EleTable.hh>
25#include <g4std/vector>
26
27#include <math.h>
b68f7176 28
29TG4GeometryServices* TG4GeometryServices::fgInstance = 0;
e5967ab3 30const G4double TG4GeometryServices::fgkAZTolerance = 0.001;
31const G4double TG4GeometryServices::fgkDensityTolerance = 0.005;
b68f7176 32
fed0483a 33//_____________________________________________________________________________
e5967ab3 34TG4GeometryServices::TG4GeometryServices(TG4IntMap* mediumMap,
b68f7176 35 TG4NameMap* nameMap)
e5967ab3 36 : fMediumMap(mediumMap),
37 fNameMap(nameMap)
b68f7176 38{
39//
40 if (fgInstance) {
41 TG4Globals::Exception(
42 "TG4GeometryServices: attempt to create two instances of singleton.");
43 }
44
b68f7176 45 fgInstance = this;
46}
47
fed0483a 48//_____________________________________________________________________________
b68f7176 49TG4GeometryServices::TG4GeometryServices() {
50//
51 TG4Globals::Exception(
e5967ab3 52 "TG4GeometryServices default constructor is protected.");
b68f7176 53}
54
55
fed0483a 56//_____________________________________________________________________________
b68f7176 57TG4GeometryServices::TG4GeometryServices(const TG4GeometryServices& right) {
58//
59 TG4Globals::Exception(
60 "Attempt to copy TG4GeometryServices singleton.");
61}
62
63
fed0483a 64//_____________________________________________________________________________
b68f7176 65TG4GeometryServices::~TG4GeometryServices() {
66//
67}
68
69//=============================================================================
70//
71// operators
72//
73//=============================================================================
74
fed0483a 75//_____________________________________________________________________________
b68f7176 76TG4GeometryServices&
77TG4GeometryServices::operator=(const TG4GeometryServices& right)
78{
79 // check assignement to self
80 if (this == &right) return *this;
81
82 TG4Globals::Exception(
83 "Attempt to assign TG4GeometryServices singleton.");
84
85 return *this;
86}
87
e5967ab3 88//=============================================================================
89//
90// private methods
91//
92//=============================================================================
93
94//_____________________________________________________________________________
95G4bool TG4GeometryServices::IsG3Volume(const G4String& lvName) const
96{
97// Returns true if the logical volume of given volumeName
98// was not created by Gsposp method with a generic name
99// (name_copyNo).
100// ---
101
102 if (lvName.contains(gSeparator))
103 return false;
104 else
105 return true;
106}
107
108//_____________________________________________________________________________
109G4bool TG4GeometryServices::CompareElement(G4double a, G4double z,
110 const G4Element* element) const
111{
112// Compares given parameters with those of a given element,
113// returns true if they are equal, false otherwise.
114// ---
115
116 G4double ae = element->GetA()/TG4G3Units::AtomicWeight();
117 G4double ze = element->GetZ();
118
119 // g3tog4 can redefine A
7c00fc32 120 G4double ax;
121 if (z<1) {
122 // vacuum
123 ax = 1.01*g/mole;
124 }
125 else
126 ax = G3Ele.GetEle(z)->GetA()/TG4G3Units::AtomicWeight();
e5967ab3 127
128 if (abs(ax - ae) < fgkAZTolerance && abs(z - ze) < fgkAZTolerance)
129 return true;
130 else
131 return false;
132}
133
134//_____________________________________________________________________________
135G4bool TG4GeometryServices::CompareMaterial(G4int nofElements, G4double density,
136 const G4Material* material) const
137{
138// Compares given density with those of a given material,
139// returns true if they are equal, false otherwise.
140// ---
141
142 G4double dm = material->GetDensity()/TG4G3Units::MassDensity();
143 G4int ne = material->GetNumberOfElements();
144
145 // density percentual difference
146 G4double diff = abs(density - dm)/(density + dm)*2.;
147
148 if (nofElements == ne && diff < fgkDensityTolerance)
149 return true;
150 else
151 return false;
152}
153
154//_____________________________________________________________________________
155G4double* TG4GeometryServices::ConvertAtomWeight(G4int nmat,
156 G4double* a, G4double* wmat) const
157{
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.
162// ---
163
164 G4double* weight = new G4double[abs(nmat)];
165
166 if (nmat<0) {
167 G4double aMol = 0.;
168 G4int i;
169 for (i=0; i<abs(nmat); i++) {
170 // total molecular weight
171 aMol += wmat[i]*a[i];
172 }
173 if (aMol == 0.) {
174 G4String text = "TG4GeometryServices::ConvertAtomWeight:\n";
175 text = text + " Total molecular weight = 0";
176 TG4Globals::Warning(text);
177 }
178 for (G4int i=0; i<abs(nmat); i++) {
179 // weight fractions
180 weight[i] = wmat[i]*a[i]/aMol;
181 }
182 }
183 else
184 for (G4int j=0; j<nmat; j++) weight[j] = wmat[j];
185
186 return weight;
187}
b68f7176 188
189//=============================================================================
190//
191// public methods
192//
193//=============================================================================
194
fed0483a 195//_____________________________________________________________________________
b68f7176 196G4double* TG4GeometryServices::CreateG4doubleArray(Float_t* array,
197 G4int size) const
198{
199// Converts Float_t* array to G4double*,
200// !! The new array has to be deleted by user.
201// ---
202
203 G4double* doubleArray;
204 if (size>0) {
205 doubleArray = new G4double[size];
206 for (G4int i=0; i<size; i++) doubleArray[i] = array[i];
207 }
208 else {
209 doubleArray = 0;
210 }
211 return doubleArray;
212}
213
fed0483a 214//_____________________________________________________________________________
b68f7176 215G4String TG4GeometryServices::CutName(const char* name) const
216{
217// Removes spaces after the name if present.
218// ---
219
220 G4String cutName = name;
221 G4int i = cutName.length();
222 while (cutName(--i) == ' ') cutName = cutName(0,i);
223
224 return cutName;
225}
226
fed0483a 227//_____________________________________________________________________________
e5967ab3 228G4String TG4GeometryServices::CutMaterialName(const char* name) const
229{
230// Removes the $ with precedent spaces at the name if present.
231// ---
232
233 G4String cutName = name;
234 cutName = cutName.substr(0,cutName.find('$'));
235
236 return CutName(cutName);
237}
238
239//_____________________________________________________________________________
240G4String TG4GeometryServices::G4ToG3VolumeName(const G4String& name) const
b68f7176 241{
242// Cuts _copyNo extension added to logical volume name in case
243// the logical volume was created by Gsposp method.
244// ---
245
e5967ab3 246 G4String cutName = name;
247 if (cutName.contains(gSeparator))
248 cutName = cutName(0,cutName.first(gSeparator));
249
250 return cutName;
b68f7176 251}
252
fed0483a 253//_____________________________________________________________________________
e5967ab3 254G4String TG4GeometryServices::GenerateLimitsName(G4int id,
255 const G4String& medName,
256 const G4String& matName) const
e4a64a3a 257{
e5967ab3 258// Generate unique name for user limits composed from the tracking medium id,
259// name and its material name.
260//
261 G4String name = "";
262 TG4Globals::AppendNumberToString(name, id);
263 name = name + "__med_" + medName + "__mat_" + matName;
e4a64a3a 264
e5967ab3 265 return name;
e4a64a3a 266}
b68f7176 267
fed0483a 268//_____________________________________________________________________________
b68f7176 269G4Material* TG4GeometryServices::MixMaterials(G4String name, G4double density,
e5967ab3 270 const TG4StringVector& matNames,
271 const TG4doubleVector& matWeights)
b68f7176 272{
273// Creates a mixture of selected materials
274// ---
275
276 // number of materials to be mixed
e5967ab3 277 G4int nofMaterials = matNames.size();
278 if (nofMaterials != matWeights.size()) {
b68f7176 279 G4String text = "TG4GeometryServices::MixMaterials: ";
280 text = text + "different number of material names and weigths.";
281 TG4Globals::Exception(text);
282 }
283 // add verbose
284 // G4cout << "Nof of materials to be mixed: " << nofMaterials << G4endl;
285
286 // fill vector of materials
e5967ab3 287 G4std::vector<G4Material*> matVector;
b68f7176 288 G4int im;
289 for (im=0; im< nofMaterials; im++) {
290 // material
e5967ab3 291 G4Material* material = G4Material::GetMaterial(matNames[im]);
292 matVector.push_back(material);
b68f7176 293 }
294
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];
e5967ab3 299 G4double fraction = matWeights[im];
b68f7176 300 mixture->AddMaterial(material, fraction);
301 }
302
303 return mixture;
304}
305
fed0483a 306//_____________________________________________________________________________
e5967ab3 307void TG4GeometryServices::PrintNameMap() const
308{
309// Prints the map of volumes names to second names.
310// ---
311
312 fNameMap->PrintAll();
313}
314
315//_____________________________________________________________________________
316void TG4GeometryServices::PrintLimits(const G4String& name) const
317{
318// Finds the limits with the specified name and prints them.
319// ---
320
321 TG4Limits* limits = FindLimits(name, true);
322
323 if (limits) limits->Print();
324}
325
326//_____________________________________________________________________________
327void TG4GeometryServices::PrintVolumeLimits(const G4String& volumeName) const
328{
329// Finds a logical volume with the specified name and prints
330// its limits.
331// ---
332
333 G4LogicalVolume* lv = FindLogicalVolume(volumeName, false);
334
335 if (lv) {
336 TG4Limits* limits = GetLimits(lv->GetUserLimits());
337 G4cout << volumeName << " ";
338 if (limits)
339 limits->Print();
340 else
341 G4cout << "has not the limits set." << G4endl;
342 }
343}
344
345//_____________________________________________________________________________
346void TG4GeometryServices::PrintStatistics(G4bool open, G4bool close) const
347{
348// Print G4 geometry statistics
349// ---
350
351
352 if (open) TG4Globals::PrintStars(true);
353
354 G4cout << " GEANT4 Geometry statistics: " << G4endl
355 << " " << NofG4LogicalVolumes()
356 << " logical volumes" << G4endl
357 << " "
358 << NofG4PhysicalVolumes()
359 << " physical volumes" << G4endl
360 << " "
361 << G4Material::GetNumberOfMaterials()
362 << " materials" << G4endl
363 << " "
364 << TG4Limits::GetNofLimits()
365 << " user limits" << G4endl;
366
367 if (close) TG4Globals::PrintStars(false);
368}
369
370//_____________________________________________________________________________
fed0483a 371Int_t TG4GeometryServices::NofG3Volumes() const
b68f7176 372{
373// Returns the total number of logical volumes corresponding
374// to G3 volumes. (
375// The logical volume that were created by Gsposp method
376// with a generic name (name_copyNo) are NOT included.
377// ---
378
e5967ab3 379 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
b68f7176 380
381 G4int counter = 0;
e5967ab3 382 for (G4int i=0; i<lvStore->size(); i++) {
383 G4LogicalVolume* lv = (*lvStore)[i];
b68f7176 384 if (IsG3Volume(lv->GetName())) counter++;
385 }
386
387 return counter;
388}
389
fed0483a 390//_____________________________________________________________________________
b68f7176 391Int_t TG4GeometryServices::NofG4LogicalVolumes() const
392{
393// Returns the total number of logical volumes in the geometry.
394// ---
395
e5967ab3 396 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
397 return lvStore->size();
b68f7176 398}
399
fed0483a 400//_____________________________________________________________________________
b68f7176 401Int_t TG4GeometryServices::NofG4PhysicalVolumes() const
402{
403// Returns the total number of physical volumes in the geometry.
404// ---
405
e5967ab3 406 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
b68f7176 407
408 G4int counter = 0;
e5967ab3 409 for (G4int i=0; i<lvStore->size(); i++) {
410 counter += ((*lvStore)[i])->GetNoDaughters();
b68f7176 411 }
412
413 return counter;
414}
415
e5967ab3 416
417//______________________________________________________________________________
418G4bool TG4GeometryServices::IsSpecialControls() const
b68f7176 419{
e5967ab3 420// Returns true if a process control in some limits instance is set.
b68f7176 421// ---
422
e5967ab3 423 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
424
425 for (G4int i=0; i<lvStore->size(); i++) {
426 TG4Limits* limits = GetLimits((*lvStore)[i]->GetUserLimits());
427 if (limits && limits->IsControl()) return true;
428 }
429
430 return false;
b68f7176 431}
432
fed0483a 433//_____________________________________________________________________________
e5967ab3 434TG4Limits* TG4GeometryServices::GetLimits(G4UserLimits* limits) const
435{
436// Checks and converts type of the given limits.
437// ---
438
439 if (!limits) return 0;
440
441 TG4Limits* tg4Limits = dynamic_cast<TG4Limits*> (limits);
442
443 if (!tg4Limits) {
444 G4Exception("TG4GeometryServices::GetLimits: Wrong limits type.");
445 return 0;
446 }
447 else
448 return tg4Limits;
449}
450
451//_____________________________________________________________________________
b68f7176 452const G4String& TG4GeometryServices::GetMapSecond(const G4String& name)
453{
454// Returns the second string associated with the name in
455// the name map.
456// ---
457
458 return fNameMap->GetSecond(name);
459}
460
e5967ab3 461//_____________________________________________________________________________
462G4LogicalVolume*
463TG4GeometryServices::FindLogicalVolume(const G4String& name, G4bool silent) const
464{
465// Finds a logical volume with the specified name in G4LogicalVolumeStore.
466// ---
467
468 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
469
470 for (G4int i=0; i<lvStore->size(); i++) {
471 G4LogicalVolume* lv = (*lvStore)[i];
472 if (lv->GetName() == name) return lv;
473 }
474
475 if (!silent) {
476 G4String text = "TG4GeometryServices:: FindLogicalVolume:\n";
477 text = text + " Logical volume " + name + " does not exist.";
478 TG4Globals::Warning(text);
479 }
480 return 0;
481}
482
b68f7176 483
fed0483a 484//_____________________________________________________________________________
e5967ab3 485TG4Limits*
486TG4GeometryServices::FindLimits(const G4String& name, G4bool silent) const
487{
488// Finds limits with the specified name.
489// ---
490
491 G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
492
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;
497 }
498
499 if (!silent) {
500 G4String text = "TG4GeometryServices:: FindLimits:\n";
501 text = text + " Logical volume " + name + " does not exist.";
502 TG4Globals::Warning(text);
503 }
504 return 0;
505}
506
507//_____________________________________________________________________________
508G4int TG4GeometryServices::GetMediumId(G4LogicalVolume* lv) const
b68f7176 509{
510// Returns the second index for materials (having its origin in
511// G4 tracking media concept)
512// ---
513
e5967ab3 514 return fMediumMap->GetSecond(lv->GetName());
b68f7176 515}
516
fed0483a 517//_____________________________________________________________________________
b68f7176 518G4double TG4GeometryServices::GetEffA(G4Material* material) const
519{
520// Returns A or effective A=sum(pi*Ai) (if compound/mixture)
521// of given material.
522// ---
523
524 G4double a = 0.;
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);
531
532 for (G4int i=0; i<nofElements; i++) {
533 G4double aOfElement = material->GetElement(i)->GetA();
534 G4double massFraction = material->GetFractionVector()[i];
c4294cc7 535 a += aOfElement*massFraction /(TG4G3Units::AtomicWeight());
b68f7176 536 }
537 }
538 else {
539 a = material->GetA();
c4294cc7 540 a /= TG4G3Units::AtomicWeight();
b68f7176 541 }
542 return a;
543}
544
fed0483a 545//_____________________________________________________________________________
b68f7176 546G4double TG4GeometryServices::GetEffZ(G4Material* material) const
547{
548// Returns Z or effective Z=sum(pi*Zi) (if compound/mixture)
549// of given material.
550// ---
551
552 G4double z = 0.;
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);
559
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;
564 }
565 }
566 else {
567 z = material->GetZ();
568 }
569 return z;
570}
e5967ab3 571
572//_____________________________________________________________________________
573G4Material* TG4GeometryServices::FindMaterial(G4double a, G4double z,
574 G4double density) const
575{
576// Finds material in G4MaterialTable with specified parameters,
577// returns 0 if not found.
578// ---
579
580 const G4MaterialTable* kpMatTable = G4Material::GetMaterialTable();
581
582 for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
583
584 G4Material* material = (*kpMatTable)[i];
585
586 if (CompareElement(a, z, material->GetElement(0)) &&
587 CompareMaterial(1, density, material))
588
589 return material;
590 }
591
592 return 0;
593}
594
595//_____________________________________________________________________________
596G4Material* TG4GeometryServices::FindMaterial(G4double* a, G4double* z,
597 G4double density,
598 G4int nmat, G4double* wmat) const
599{
600// Finds material in G4MaterialTable with specified parameters,
601// returns 0 if not found.
602// ---
603
604 G4double* weight = ConvertAtomWeight(nmat, a, wmat);
605
606 // loop over materials
607 G4Material* found = 0;
608 for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
609
610 G4Material* material = (*G4Material::GetMaterialTable())[i];
611 G4int nm = material->GetNumberOfElements();
612
613 if (CompareMaterial(nm, density, material)) {
614
615 // loop over elements
616 G4bool equal = true;
617 for (G4int ie=0; ie<nm; ie++) {
618
619 G4double we = (material->GetFractionVector())[ie];
620
621 if (!CompareElement(a[ie], z[ie], material->GetElement(ie)) ||
622 abs(weight[ie] - we) > fgkAZTolerance) {
623
624 equal = false;
625 break;
626 }
627 }
628 if (equal) {
629 found = material;
630 break;
631 }
632 }
633 }
634
635 delete [] weight;
636 return found;
637}
638