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 = G3Ele.GetEle(z)->GetA()/TG4G3Units::AtomicWeight();
 121
 122   if (abs(ax - ae) < fgkAZTolerance && abs(z  - ze) < fgkAZTolerance)
 123     return true;
 124   else
 125     return false;
 126 }
 127
 128 //_____________________________________________________________________________
 129 G4bool TG4GeometryServices::CompareMaterial(G4int nofElements, G4double density,
 130                                             const G4Material* material) const
 131 {
 132 // Compares given density with those of a given material,
 133 // returns true if they are equal, false otherwise.
 134 // ---
 135
 136   G4double dm = material->GetDensity()/TG4G3Units::MassDensity();
 137   G4int ne = material->GetNumberOfElements();
 138
 139   // density percentual difference
 140   G4double diff = abs(density - dm)/(density + dm)*2.;
 141
 142   if (nofElements == ne && diff < fgkDensityTolerance)
 143     return true;
 144   else
 145     return false;
 146 }
 147
 148 //_____________________________________________________________________________
 149 G4double* TG4GeometryServices::ConvertAtomWeight(G4int nmat,
 150                                                  G4double* a, G4double* wmat) const
 151 {
 152 // In case of proportions given in atom counts (nmat<0),
 153 // the wmat[i] are converted to weight fractions.
 154 // (From g3tog4 G4gsmixt.)
 155 // The new array has to be delete by client.
 156 // ---
 157
 158   G4double* weight = new G4double[abs(nmat)];
 159
 160   if (nmat<0) {
 161     G4double aMol = 0.;
 162     G4int i;
 163     for (i=0; i<abs(nmat); i++) {
 164       // total molecular weight
 165       aMol += wmat[i]*a[i];
 166     }
 167     if (aMol == 0.) {
 168       G4String text = "TG4GeometryServices::ConvertAtomWeight:\n";
 169       text = text + " Total molecular weight = 0";
 170       TG4Globals::Warning(text);
 171     }
 172     for (G4int i=0; i<abs(nmat); i++) {
 173       // weight fractions
 174       weight[i] = wmat[i]*a[i]/aMol;
 175     }
 176   }
 177   else
 178     for (G4int j=0; j<nmat; j++) weight[j] = wmat[j];
 179
 180   return weight;
 181 }
 182
 183 //=============================================================================
 184 //
 185 // public methods
 186 //
 187 //=============================================================================
 188
 189 //_____________________________________________________________________________
 190 G4double* TG4GeometryServices::CreateG4doubleArray(Float_t* array,
 191                G4int size) const
 192 {
 193 // Converts Float_t* array to G4double*,
 194 // !! The new array has to be deleted by user.
 195 // ---
 196
 197   G4double* doubleArray;
 198   if (size>0) {
 199     doubleArray = new G4double[size];
 200     for (G4int i=0; i<size; i++) doubleArray[i] = array[i];
 201   }
 202   else {
 203     doubleArray = 0;
 204   }
 205   return doubleArray;
 206 }
 207
 208 //_____________________________________________________________________________
 209 G4String TG4GeometryServices::CutName(const char* name) const
 210 {
 211 // Removes spaces after the name if present.
 212 // ---
 213
 214   G4String cutName = name;
 215   G4int i = cutName.length();
 216   while (cutName(--i) == ' ') cutName = cutName(0,i);
 217
 218   return cutName;
 219 }
 220
 221 //_____________________________________________________________________________
 222 G4String TG4GeometryServices::CutMaterialName(const char* name) const
 223 {
 224 // Removes the $ with precedent spaces at the name if present.
 225 // ---
 226
 227   G4String cutName = name;
 228   cutName = cutName.substr(0,cutName.find('$'));
 229
 230   return CutName(cutName);
 231 }
 232
 233 //_____________________________________________________________________________
 234 G4String  TG4GeometryServices::G4ToG3VolumeName(const G4String& name) const
 235 {
 236 // Cuts _copyNo extension added to logical volume name in case
 237 // the logical volume was created by Gsposp method.
 238 // ---
 239
 240   G4String cutName = name;
 241   if (cutName.contains(gSeparator))
 242   cutName = cutName(0,cutName.first(gSeparator));
 243
 244   return cutName;
 245 }
 246
 247 //_____________________________________________________________________________
 248 G4String  TG4GeometryServices::GenerateLimitsName(G4int id,
 249                                        const G4String& medName,
 250                                        const G4String& matName) const
 251 {
 252 // Generate unique name for user limits composed from the tracking medium id,
 253 // name and its material name.
 254 //
 255   G4String name = "";
 256   TG4Globals::AppendNumberToString(name, id);
 257   name = name + "__med_" + medName + "__mat_" + matName;
 258
 259   return name;
 260 }
 261
 262 //_____________________________________________________________________________
 263 G4Material* TG4GeometryServices::MixMaterials(G4String name, G4double density,
 264                                     const TG4StringVector& matNames,
 265                                     const TG4doubleVector& matWeights)
 266 {
 267 // Creates a mixture of selected materials
 268 // ---
 269
 270   // number of materials to be mixed
 271   G4int nofMaterials = matNames.size();
 272   if (nofMaterials != matWeights.size()) {
 273     G4String text = "TG4GeometryServices::MixMaterials: ";
 274     text = text +  "different number of material names and weigths.";
 275     TG4Globals::Exception(text);
 276   }
 277   // add verbose
 278   // G4cout << "Nof of materials to be mixed: " << nofMaterials << G4endl;
 279
 280   // fill vector of materials
 281   G4std::vector<G4Material*> matVector;
 282   G4int im;
 283   for (im=0; im< nofMaterials; im++) {
 284     // material
 285     G4Material* material = G4Material::GetMaterial(matNames[im]);
 286     matVector.push_back(material);
 287   }
 288
 289   // create the mixed material
 290   G4Material* mixture = new G4Material(name, density, nofMaterials);
 291   for (im=0; im< nofMaterials; im++) {
 292     G4Material* material = matVector[im];
 293     G4double fraction = matWeights[im];
 294     mixture->AddMaterial(material, fraction);
 295   }
 296
 297   return mixture;
 298 }
 299
 300 //_____________________________________________________________________________
 301 void TG4GeometryServices::PrintNameMap() const
 302 {
 303 // Prints the map of volumes names to second names.
 304 // ---
 305
 306   fNameMap->PrintAll();
 307 }
 308
 309 //_____________________________________________________________________________
 310 void TG4GeometryServices::PrintLimits(const G4String& name) const
 311 {
 312 // Finds the limits with the specified name and prints them.
 313 // ---
 314
 315   TG4Limits* limits = FindLimits(name, true);
 316
 317   if (limits)  limits->Print();
 318 }
 319
 320 //_____________________________________________________________________________
 321 void TG4GeometryServices::PrintVolumeLimits(const G4String& volumeName) const
 322 {
 323 // Finds a logical volume with the specified name and prints
 324 // its limits.
 325 // ---
 326
 327   G4LogicalVolume* lv = FindLogicalVolume(volumeName, false);
 328
 329   if (lv) {
 330     TG4Limits* limits = GetLimits(lv->GetUserLimits());
 331     G4cout << volumeName << "  ";
 332     if (limits)
 333       limits->Print();
 334     else
 335       G4cout << "has not the limits set." << G4endl;
 336   }
 337 }
 338
 339 //_____________________________________________________________________________
 340 void TG4GeometryServices::PrintStatistics(G4bool open, G4bool close) const
 341 {
 342 // Print G4 geometry statistics
 343 // ---
 344
 345
 346   if (open)  TG4Globals::PrintStars(true);
 347
 348   G4cout << "    GEANT4 Geometry statistics: " << G4endl
 349          << "          " << NofG4LogicalVolumes()
 350                          << " logical volumes" << G4endl
 351          << "          "
 352                          << NofG4PhysicalVolumes()
 353                          << " physical volumes" << G4endl
 354          << "          "
 355                          << G4Material::GetNumberOfMaterials()
 356                          << " materials"        << G4endl
 357          << "          "
 358                          << TG4Limits::GetNofLimits()
 359                          << " user limits"      << G4endl;
 360
 361   if (close) TG4Globals::PrintStars(false);
 362 }
 363
 364 //_____________________________________________________________________________
 365 Int_t TG4GeometryServices::NofG3Volumes() const
 366 {
 367 // Returns the total number of logical volumes corresponding
 368 // to G3 volumes. (
 369 // The logical volume that were created by Gsposp method
 370 // with a generic name (name_copyNo) are NOT included.
 371 // ---
 372
 373   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 374
 375   G4int counter = 0;
 376   for (G4int i=0; i<lvStore->size(); i++) {
 377     G4LogicalVolume* lv = (*lvStore)[i];
 378     if (IsG3Volume(lv->GetName())) counter++;
 379   }
 380
 381   return counter;
 382 }
 383
 384 //_____________________________________________________________________________
 385 Int_t TG4GeometryServices::NofG4LogicalVolumes() const
 386 {
 387 // Returns the total number of logical volumes in the geometry.
 388 // ---
 389
 390   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 391   return lvStore->size();
 392 }
 393
 394 //_____________________________________________________________________________
 395 Int_t TG4GeometryServices::NofG4PhysicalVolumes() const
 396 {
 397 // Returns the total number of physical volumes in the geometry.
 398 // ---
 399
 400   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 401
 402   G4int counter = 0;
 403   for (G4int i=0; i<lvStore->size(); i++) {
 404     counter += ((*lvStore)[i])->GetNoDaughters();
 405   }
 406
 407   return counter;
 408 }
 409
 410
 411 //______________________________________________________________________________
 412 G4bool TG4GeometryServices::IsSpecialControls()  const
 413 {
 414 // Returns true if a process control in some limits instance is set.
 415 // ---
 416
 417   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 418
 419   for (G4int i=0; i<lvStore->size(); i++) {
 420     TG4Limits* limits = GetLimits((*lvStore)[i]->GetUserLimits());
 421     if (limits && limits->IsControl()) return true;
 422   }
 423
 424   return false;
 425 }
 426
 427 //_____________________________________________________________________________
 428 TG4Limits* TG4GeometryServices::GetLimits(G4UserLimits* limits) const
 429 {
 430 // Checks and converts type of the given limits.
 431 // ---
 432
 433   if (!limits) return 0;
 434
 435   TG4Limits* tg4Limits = dynamic_cast<TG4Limits*> (limits);
 436
 437   if (!tg4Limits) {
 438     G4Exception("TG4GeometryServices::GetLimits: Wrong limits type.");
 439     return 0;
 440   }
 441   else
 442     return tg4Limits;
 443 }
 444
 445 //_____________________________________________________________________________
 446 const G4String& TG4GeometryServices::GetMapSecond(const G4String& name)
 447 {
 448 // Returns the second string associated with the name in
 449 // the name map.
 450 // ---
 451
 452   return fNameMap->GetSecond(name);
 453 }
 454
 455 //_____________________________________________________________________________
 456 G4LogicalVolume*
 457 TG4GeometryServices::FindLogicalVolume(const G4String& name, G4bool silent) const
 458 {
 459 // Finds a logical volume with the specified name in G4LogicalVolumeStore.
 460 // ---
 461
 462   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 463
 464   for (G4int i=0; i<lvStore->size(); i++) {
 465     G4LogicalVolume* lv = (*lvStore)[i];
 466     if (lv->GetName() == name) return lv;
 467   }
 468
 469   if (!silent) {
 470     G4String text = "TG4GeometryServices:: FindLogicalVolume:\n";
 471     text = text + "    Logical volume " + name + " does not exist.";
 472     TG4Globals::Warning(text);
 473   }
 474   return 0;
 475 }
 476
 477
 478 //_____________________________________________________________________________
 479 TG4Limits*
 480 TG4GeometryServices::FindLimits(const G4String& name, G4bool silent) const
 481 {
 482 // Finds limits with the specified name.
 483 // ---
 484
 485   G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
 486
 487   for (G4int i=0; i<lvStore->size(); i++) {
 488     G4LogicalVolume* lv = (*lvStore)[i];
 489     TG4Limits* limits = GetLimits(lv->GetUserLimits());
 490     if (limits && limits->GetName() == name) return limits;
 491   }
 492
 493   if (!silent) {
 494     G4String text = "TG4GeometryServices:: FindLimits:\n";
 495     text = text + "    Logical volume " + name + " does not exist.";
 496     TG4Globals::Warning(text);
 497   }
 498   return 0;
 499 }
 500
 501 //_____________________________________________________________________________
 502 G4int TG4GeometryServices::GetMediumId(G4LogicalVolume* lv) const
 503 {
 504 // Returns the second index for materials (having its origin in
 505 // G4 tracking media concept)
 506 // ---
 507
 508   return fMediumMap->GetSecond(lv->GetName());
 509 }
 510
 511 //_____________________________________________________________________________
 512 G4double TG4GeometryServices::GetEffA(G4Material* material) const
 513 {
 514 // Returns A or effective A=sum(pi*Ai) (if compound/mixture)
 515 // of given material.
 516 // ---
 517
 518   G4double a = 0.;
 519   G4int nofElements = material->GetNumberOfElements();
 520   if (nofElements > 1) {
 521     G4String text = "Effective A for material mixture (";
 522     text = text + material->GetName();
 523     text = text + ") is used.";
 524     //TG4Globals::Warning(text);
 525
 526     for (G4int i=0; i<nofElements; i++) {
 527       G4double aOfElement = material->GetElement(i)->GetA();
 528       G4double massFraction = material->GetFractionVector()[i];
 529       a += aOfElement*massFraction /(TG4G3Units::AtomicWeight());
 530     }
 531   }
 532   else {
 533     a = material->GetA();
 534     a /= TG4G3Units::AtomicWeight();
 535   }
 536   return a;
 537 }
 538
 539 //_____________________________________________________________________________
 540 G4double TG4GeometryServices::GetEffZ(G4Material* material) const
 541 {
 542 // Returns Z or effective Z=sum(pi*Zi) (if compound/mixture)
 543 // of given material.
 544 // ---
 545
 546   G4double z = 0.;
 547   G4int nofElements = material->GetNumberOfElements();
 548   if (nofElements > 1) {
 549     G4String text = "Effective Z for material mixture (";
 550     text = text + material->GetName();
 551     text = text + ") is used.";
 552     //TG4Globals::Warning(text);
 553
 554     for (G4int i=0; i<nofElements; i++) {
 555       G4double zOfElement = material->GetElement(i)->GetZ();
 556       G4double massFraction = material->GetFractionVector()[i];
 557       z += zOfElement*massFraction;
 558     }
 559   }
 560   else {
 561     z = material->GetZ();
 562   }
 563   return z;
 564 }
 565
 566 //_____________________________________________________________________________
 567 G4Material* TG4GeometryServices::FindMaterial(G4double a, G4double z,
 568                                              G4double density) const
 569 {
 570 // Finds material in G4MaterialTable with specified parameters,
 571 // returns 0 if not found.
 572 // ---
 573
 574   const G4MaterialTable* kpMatTable = G4Material::GetMaterialTable();
 575
 576   for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
 577
 578     G4Material* material = (*kpMatTable)[i];
 579
 580     if (CompareElement(a, z, material->GetElement(0)) &&
 581         CompareMaterial(1, density, material))
 582
 583         return material;
 584   }
 585
 586   return 0;
 587 }
 588
 589 //_____________________________________________________________________________
 590 G4Material* TG4GeometryServices::FindMaterial(G4double* a, G4double* z,
 591                                               G4double density,
 592                                               G4int nmat, G4double* wmat) const
 593 {
 594 // Finds material in G4MaterialTable with specified parameters,
 595 // returns 0 if not found.
 596 // ---
 597
 598   G4double* weight = ConvertAtomWeight(nmat, a, wmat);
 599
 600   // loop over materials
 601   G4Material* found = 0;
 602   for (G4int i=0; i<G4Material::GetNumberOfMaterials(); i++) {
 603
 604     G4Material* material = (*G4Material::GetMaterialTable())[i];
 605     G4int nm = material->GetNumberOfElements();
 606
 607     if (CompareMaterial(nm, density, material)) {
 608
 609       // loop over elements
 610       G4bool equal = true;
 611       for (G4int ie=0; ie<nm; ie++) {
 612
 613         G4double we = (material->GetFractionVector())[ie];
 614
 615         if (!CompareElement(a[ie], z[ie], material->GetElement(ie)) ||
 616             abs(weight[ie] - we) > fgkAZTolerance) {
 617
 618           equal = false;
 619           break;
 620         }
 621       }
 622       if (equal) {
 623         found = material;
 624         break;
 625       }
 626     }
 627   }
 628
 629   delete [] weight;
 630   return found;
 631 }
 632