Correction in CompareElement (to prevent from crash when vacuum (z<1) is passed.
[u/mrichter/AliRoot.git] / TGeant4 / TG4GeometryServices.cxx
1 // $Id$
2 // Category: geometry
3 //
4 // Author: I. Hrivnacova
5 //
6 // Class TG4GeometryServices
7 // -------------------------
8 // See the class description in the header file.
9
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"
17
18 #include <G4LogicalVolumeStore.hh>
19 #include <G4LogicalVolume.hh>
20 #include <G4Material.hh>
21 #include <G4Element.hh>
22 #include <G4UserLimits.hh>
23 #include <G3toG4.hh> 
24 #include <G3EleTable.hh> 
25 #include <g4std/vector>
26
27 #include <math.h>
28
29 TG4GeometryServices* TG4GeometryServices::fgInstance = 0;
30 const G4double       TG4GeometryServices::fgkAZTolerance = 0.001;      
31 const G4double       TG4GeometryServices::fgkDensityTolerance = 0.005; 
32
33 //_____________________________________________________________________________
34 TG4GeometryServices::TG4GeometryServices(TG4IntMap* mediumMap, 
35                                          TG4NameMap* nameMap) 
36   : fMediumMap(mediumMap),
37     fNameMap(nameMap)                            
38 {
39 //
40   if (fgInstance) {
41     TG4Globals::Exception(
42       "TG4GeometryServices: attempt to create two instances of singleton.");
43   }
44
45   fgInstance = this;
46 }
47
48 //_____________________________________________________________________________
49 TG4GeometryServices::TG4GeometryServices() {
50 // 
51   TG4Globals::Exception(
52     "TG4GeometryServices default constructor is protected.");
53 }
54
55
56 //_____________________________________________________________________________
57 TG4GeometryServices::TG4GeometryServices(const TG4GeometryServices& right) {
58 // 
59   TG4Globals::Exception(
60     "Attempt to copy TG4GeometryServices singleton.");
61 }
62
63
64 //_____________________________________________________________________________
65 TG4GeometryServices::~TG4GeometryServices() {
66 //
67 }
68
69 //=============================================================================
70 //
71 // operators
72 //
73 //=============================================================================
74
75 //_____________________________________________________________________________
76 TG4GeometryServices& 
77 TG4GeometryServices::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           
88 //=============================================================================
89 //
90 // private methods
91 //
92 //=============================================================================
93
94 //_____________________________________________________________________________
95 G4bool 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 //_____________________________________________________________________________
109 G4bool 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
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();
127
128   if (abs(ax - ae) < fgkAZTolerance && abs(z  - ze) < fgkAZTolerance) 
129     return true;
130   else  
131     return false;   
132 }                                            
133
134 //_____________________________________________________________________________
135 G4bool 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 //_____________________________________________________________________________
155 G4double* 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 }
188
189 //=============================================================================
190 //
191 // public methods
192 //
193 //=============================================================================
194
195 //_____________________________________________________________________________
196 G4double* 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
214 //_____________________________________________________________________________
215 G4String 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
227 //_____________________________________________________________________________
228 G4String 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 //_____________________________________________________________________________
240 G4String  TG4GeometryServices::G4ToG3VolumeName(const G4String& name) const
241 {
242 // Cuts _copyNo extension added to logical volume name in case 
243 // the logical volume was created by Gsposp method.
244 // ---
245
246   G4String cutName = name;
247   if (cutName.contains(gSeparator)) 
248   cutName = cutName(0,cutName.first(gSeparator));
249   
250   return cutName;
251 }
252
253 //_____________________________________________________________________________
254 G4String  TG4GeometryServices::GenerateLimitsName(G4int id, 
255                                        const G4String& medName,
256                                        const G4String& matName) const
257 {
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;
264   
265   return name;
266 }
267
268 //_____________________________________________________________________________
269 G4Material* TG4GeometryServices::MixMaterials(G4String name, G4double density, 
270                                     const TG4StringVector& matNames, 
271                                     const TG4doubleVector& matWeights)
272 {
273 // Creates a mixture of selected materials
274 // ---
275
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);
282   }    
283   // add verbose
284   // G4cout << "Nof of materials to be mixed: " << nofMaterials << G4endl;
285
286   // fill vector of materials
287   G4std::vector<G4Material*> matVector;  
288   G4int im;
289   for (im=0; im< nofMaterials; im++) {
290     // material
291     G4Material* material = G4Material::GetMaterial(matNames[im]);
292     matVector.push_back(material);
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];
299     G4double fraction = matWeights[im];
300     mixture->AddMaterial(material, fraction);
301   }
302
303   return mixture;
304 }  
305
306 //_____________________________________________________________________________
307 void TG4GeometryServices::PrintNameMap() const
308 {
309 // Prints the map of volumes names to second names.
310 // ---
311
312   fNameMap->PrintAll();
313 }
314  
315 //_____________________________________________________________________________
316 void 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 //_____________________________________________________________________________
327 void 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 //_____________________________________________________________________________
346 void 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 //_____________________________________________________________________________
371 Int_t TG4GeometryServices::NofG3Volumes() const
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
379   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
380
381   G4int counter = 0;  
382   for (G4int i=0; i<lvStore->size(); i++) {
383     G4LogicalVolume* lv = (*lvStore)[i];
384     if (IsG3Volume(lv->GetName())) counter++;
385   }
386   
387   return counter;  
388 }
389
390 //_____________________________________________________________________________
391 Int_t TG4GeometryServices::NofG4LogicalVolumes() const
392 {
393 // Returns the total number of logical volumes in the geometry.
394 // ---
395
396   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
397   return lvStore->size();
398 }
399
400 //_____________________________________________________________________________
401 Int_t TG4GeometryServices::NofG4PhysicalVolumes() const
402 {
403 // Returns the total number of physical volumes in the geometry.
404 // ---
405
406   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
407
408   G4int counter = 0;
409   for (G4int i=0; i<lvStore->size(); i++) {
410     counter += ((*lvStore)[i])->GetNoDaughters();
411   }
412   
413   return counter;  
414 }
415
416
417 //______________________________________________________________________________
418 G4bool TG4GeometryServices::IsSpecialControls()  const
419 {
420 // Returns true if a process control in some limits instance is set. 
421 // ---
422
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;    
431 }
432
433 //_____________________________________________________________________________
434 TG4Limits* 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 //_____________________________________________________________________________
452 const 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
461 //_____________________________________________________________________________
462 G4LogicalVolume* 
463 TG4GeometryServices::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
483
484 //_____________________________________________________________________________
485 TG4Limits* 
486 TG4GeometryServices::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 //_____________________________________________________________________________
508 G4int TG4GeometryServices::GetMediumId(G4LogicalVolume* lv) const
509 {
510 // Returns the second index for materials (having its origin in
511 // G4 tracking media concept)
512 // ---
513
514   return fMediumMap->GetSecond(lv->GetName());
515 }  
516
517 //_____________________________________________________________________________
518 G4double 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];      
535       a += aOfElement*massFraction /(TG4G3Units::AtomicWeight());
536     }
537   }
538   else { 
539     a = material->GetA();
540     a /= TG4G3Units::AtomicWeight();
541   }
542   return a;
543 }
544
545 //_____________________________________________________________________________
546 G4double 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 }
571
572 //_____________________________________________________________________________
573 G4Material* 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 //_____________________________________________________________________________
596 G4Material* 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