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