TGeant4/TG4GeometryManager.cxx

   1 // $Id$
   2 // Category: geometry
   3 //
   4 // C++ interface to Geant3 basic routines
   5 // for building Geant4 geometry
   6 //
   7 // by V. Berejnoi, 25.2.1999
   8 // materials, tracking media support
   9 // added by I.Hrivnacova, 27.5.1999
  10
  11 #include "TG4GeometryManager.h"
  12 #include "TG4GeometryOutputManager.h"
  13 #include "TG4PhysicsManager.h"
  14 #include "TG4VSensitiveDetector.h"
  15 #include "TG4Limits.h"
  16 #include "TG4Globals.h"
  17 #include "TG3Units.h"
  18
  19 #include <G3toG4.hh>
  20 #include <G3toG4BuildTree.hh>
  21 #include <G3VolTable.hh>
  22 #include <G3RotTable.hh>
  23 #include <G3EleTable.hh>
  24 #include <G3MatTable.hh>
  25 #include <G3MedTable.hh>
  26 #include <G3SensVolVector.hh>
  27
  28 #include <G4SDManager.hh>
  29 #include <G4VSensitiveDetector.hh>
  30 #include <G4LogicalVolumeStore.hh>
  31 #include <G4PVPlacement.hh>
  32 #include <G4Material.hh>
  33 #include <G4MaterialPropertiesTable.hh>
  34 #include <G4Element.hh>
  35
  36 #include <math.h>
  37 #include <ctype.h>
  38
  39 // extern global method from g3tog4
  40 void G3CLRead(G4String &, char *);
  41
  42 TG4GeometryManager* TG4GeometryManager::fgInstance = 0;
  43
  44 TG4GeometryManager::TG4GeometryManager()
  45   : fMediumCounter(0),
  46     fMaterialCounter(0),
  47     fMatrixCounter(0),
  48     fUseG3TMLimits(false),
  49     fWriteGeometry(true)
  50 {
  51 //
  52   if (fgInstance) {
  53     TG4Globals::Exception(
  54       "TG4GeometryManager: attempt to create two instances of singleton.");
  55   }
  56
  57   fOutputManager = new TG4GeometryOutputManager();
  58
  59   fgInstance = this;
  60
  61   // instantiate the default element table
  62   //TG4ElementTable::Instance();
  63 }
  64
  65 TG4GeometryManager::TG4GeometryManager(const TG4GeometryManager& right) {
  66 //
  67   TG4Globals::Exception(
  68     "Attempt to copy TG4GeometryManager singleton.");
  69 }
  70
  71
  72 TG4GeometryManager::~TG4GeometryManager() {
  73 //
  74   delete fOutputManager;
  75 }
  76
  77 //==================================================================== =========
  78 //
  79 // operators
  80 //
  81 //==================================================================== =========
  82
  83 TG4GeometryManager&
  84 TG4GeometryManager::operator=(const TG4GeometryManager& right)
  85 {
  86   // check assignement to self
  87   if (this == &right) return *this;
  88
  89   TG4Globals::Exception(
  90     "Attempt to assign TG4GeometryManager singleton.");
  91
  92   return *this;
  93 }
  94
  95
  96 //=============================================================================
  97 //
  98 // private methods
  99 //
 100 //=============================================================================
 101
 102
 103 G4double* TG4GeometryManager::CreateG4doubleArray(Float_t* array,
 104                G4int size) const
 105 {
 106 // Converts Float_t* array to G4double*,
 107 // !! The new array has to be deleted by user.
 108 // ---
 109
 110   G4double* doubleArray;
 111   if (size>0) {
 112     doubleArray = new G4double[size];
 113     for (G4int i=0; i<size; i++) doubleArray[i] = array[i];
 114   }
 115   else {
 116     doubleArray = 0;
 117   }
 118   return doubleArray;
 119 }
 120
 121
 122 G4String TG4GeometryManager::CutName(const char* name) const
 123 {
 124 // Removes spaces after the name if present.
 125 // ---
 126
 127   G4String cutName = name;
 128   G4int i = cutName.length();
 129   while (cutName(--i) == ' ') cutName = cutName(0,i);
 130
 131   return cutName;
 132 }
 133
 134
 135 void TG4GeometryManager::GstparCut(G4int itmed, TG3Cut par, G4double parval)
 136 {
 137 // Sets special tracking medium parameter.
 138 // It is applied to all logical volumes that use the specified
 139 // tracking medium.
 140 // ---
 141
 142   // get medium from table
 143   G3MedTableEntry* medium = G3Med.get(itmed);
 144   if (!medium) {
 145     G4String text = "TG4GeometryManager::GstparCut: \n";
 146     text = text + "    Medium not found.";
 147     G4Exception(text);
 148   }
 149
 150   // get/create user limits
 151   G4UserLimits* limits = medium->GetLimits();
 152   TG4Limits* tg4Limits;
 153   if (limits) {
 154     tg4Limits = dynamic_cast<TG4Limits*> (limits);
 155     if (!tg4Limits)
 156       G4Exception("TG4GeometryManager::GstparCut: Wrong limits type.");
 157   }
 158   else {
 159     tg4Limits = new TG4Limits();
 160     medium->SetLimits(tg4Limits);
 161     // add verbose
 162     G4cout << "TG4GeometryManager::GstparCut: new TG4Limits() for medium "
 163            << itmed << " has been created." << endl;
 164   }
 165   // set parameter
 166   tg4Limits->SetG3Cut(par, parval*GeV);
 167 }
 168
 169
 170 void TG4GeometryManager::GstparFlag(G4int itmed, TG3Flag par, G4double parval)
 171 {
 172 // Sets special tracking medium parameter.
 173 // It is applied to all logical volumes that use the specified
 174 // tracking medium.
 175 // ---
 176
 177   // get medium from table
 178   G3MedTableEntry* medium = G3Med.get(itmed);
 179   if (!medium) {
 180     G4String text = "TG4GeometryManager::GstparFlag: \n";
 181     text = text + "    Medium not found.";
 182     G4Exception(text);
 183   }
 184
 185   // get/create user limits
 186   G4UserLimits* limits = medium->GetLimits();
 187   TG4Limits* tg4Limits;
 188   if (limits) {
 189     tg4Limits = dynamic_cast<TG4Limits*> (limits);
 190     if (!tg4Limits)
 191       G4Exception("TG4GeometryManager::GstparFlag: Wrong limits type.");
 192   }
 193   else {
 194     tg4Limits = new TG4Limits();
 195     medium->SetLimits(tg4Limits);
 196     // add verbose
 197     G4cout << "TG4GeometryManager::GstparFlag: new TG4Limits() for medium "
 198            << itmed << " has been created." << endl;
 199   }
 200   // set parameter
 201   tg4Limits->SetG3Flag(par, parval);
 202 }
 203
 204
 205 void TG4GeometryManager::FillMediumIdVector()
 206 {
 207 // The second index for materials (having its origin in
 208 // G4 tracking media concept) is stored in a vector of G4int
 209 // parallel to G4MaterialTable.
 210 // ---
 211
 212   // initialize vector
 213   G4int nofMaterials = G4Material::GetNumberOfMaterials();
 214   G4int i;
 215   for (i=0; i<nofMaterials; i++)
 216     fMediumIdVector.push_back(0);
 217
 218   // fill vector
 219   for (i=0; i<fMediumCounter; i++) {
 220     // material index (index in G4Material table)
 221     G3MedTableEntry* mte = G3Med.get(i+1);
 222     G4int materialIndex = mte->GetMaterial()->GetIndex();
 223
 224     // medium index (ID of G3MedTableEntry)
 225     G4int mediumIndex = mte->GetID();
 226
 227     // store medium index in the vector
 228     fMediumIdVector[materialIndex] = mediumIndex;
 229   }
 230
 231   // add verbose
 232   G4cout << "Total nof materials: " << nofMaterials << endl;
 233   G4cout << "Total nof tracking medias: " << fMediumCounter << endl;
 234 }
 235
 236
 237 //=============================================================================
 238 //
 239 // public methods - AliMC implementation
 240 //
 241 //=============================================================================
 242
 243
 244 void TG4GeometryManager::Material(Int_t& kmat, const char* name, Float_t a,
 245           Float_t z, Float_t dens, Float_t radl, Float_t absl, Float_t* buf,
 246           Int_t nwbuf)
 247 {
 248 // Creates G4Material.
 249 // !! Parameters radl, absl, buf, nwbuf are ignored in G4gsmate
 250 // Comment:
 251 // absl - this parameter is ignored by GEANT3, too
 252 // ---
 253
 254     kmat = ++fMaterialCounter;
 255     G4double* bufin = CreateG4doubleArray(buf, nwbuf);
 256
 257     // write token to the output file
 258     if (fWriteGeometry)
 259       fOutputManager->WriteGsmate(kmat, name, a, z, dens, radl, nwbuf, bufin);
 260
 261     G4gsmate(kmat, name, a, z, dens, radl, nwbuf, bufin);
 262
 263     delete [] bufin;
 264
 265     if (nwbuf > 0) {
 266       G4String matName = name;
 267       G4String text
 268         = "TG4GeometryManager: user defined parameters for material ";
 269       text = text + matName;
 270       text = text + " are ignored by Geant4.";
 271       TG4Globals::Warning(text);
 272     }
 273 }
 274
 275
 276 void TG4GeometryManager::Mixture(Int_t& kmat, const char *name, Float_t *a,
 277           Float_t *z, Float_t dens, Int_t nlmat, Float_t *wmat)
 278 {
 279 // Creates G4Material composed of more elements.
 280 // !! Parameters radl, absl, buf, nwbuf are ignored in G4gsmate
 281 // Comment:
 282 // absl - this parameter is ignored by GEANT3, too
 283 // ---
 284
 285    Int_t npar = abs(nlmat);
 286    G4double *ain = CreateG4doubleArray(a, npar);
 287    G4double *zin = CreateG4doubleArray(z, npar);
 288    G4double *wmatin = CreateG4doubleArray(wmat, npar);
 289
 290    kmat = ++fMaterialCounter;
 291
 292    // write token to the output file
 293    if (fWriteGeometry)
 294      fOutputManager->WriteGsmixt(kmat, name, ain, zin, dens, nlmat, wmatin);
 295
 296    G4gsmixt(kmat, name, ain, zin, dens, nlmat, wmatin);
 297
 298    // !!! in Geant3:
 299    // After a call with ratios by number (negative number of elements),
 300    // the ratio array is changed to the ratio by weight, so all successive
 301    // calls with the same array must specify the number of elements as
 302    // positive
 303
 304    // wmatin may be modified
 305    for (G4int i=0; i<npar; i++) wmat[i] = wmatin[i];
 306
 307    delete [] ain;
 308    delete [] zin;
 309    delete [] wmatin;
 310 }
 311
 312 void TG4GeometryManager::Medium(Int_t& kmed, const char *name, Int_t nmat,
 313           Int_t isvol, Int_t ifield, Float_t fieldm, Float_t tmaxfd,
 314           Float_t stemax, Float_t deemax, Float_t epsil,
 315           Float_t stmin, Float_t* ubuf, Int_t nbuf)
 316 {
 317 // Creates a temporary "medium" that is used for
 318 // assigning corresponding parameters to G4 objects:
 319 // NTMED is stored as a second material index;
 320 // ISVOL is used for builing TG3SensVolVector;
 321 // STEMAX is passed in TG4Limits (if fUseG3TMLimits is set true);
 322 // !! The other parameters (IFIELD, FIELDM, TMAXFD, DEEMAX, EPSIL, STMIN)
 323 // are ignored by Geant4.
 324 // ---
 325
 326 //  Geant3 desription:
 327 //  ==================
 328 //  NTMED  Tracking medium number
 329 //  NAME   Tracking medium name
 330 //  NMAT   Material number
 331 //  ISVOL  Sensitive volume flag
 332 //  IFIELD Magnetic field
 333 //  FIELDM Max. field value (Kilogauss)
 334 //  TMAXFD Max. angle due to field (deg/step)
 335 //  STEMAX Max. step allowed
 336 //  DEEMAX Max. fraction of energy lost in a step
 337 //  EPSIL  Tracking precision (cm)
 338 //  STMIN  Min. step due to continuos processes (cm)
 339 //
 340 //  IFIELD = 0 if no magnetic field; IFIELD = -1 if user decision in GUSWIM;
 341 //  IFIELD = 1 if tracking performed with GRKUTA; IFIELD = 2 if tracking
 342 //  performed with GHELIX; IFIELD = 3 if tracking performed with GHELX3.
 343 // ---
 344
 345   kmed = ++fMediumCounter;
 346
 347   // write token to the output file
 348   if (fWriteGeometry)
 349     fOutputManager->WriteGstmed(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd,
 350         stemax, deemax, epsil, stmin, 0, 0);
 351
 352   G4gstmed(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 353        epsil, stmin, 0, fUseG3TMLimits);
 354      // !! instead of the nbuf argument the bool fIsG3Default is passed
 355
 356   if (nbuf > 0) {
 357     G4String medName = name;
 358     G4String text
 359       = "TG4GeometryManager: user defined parameters for medium ";
 360     text = text + medName;
 361     text = text + " are ignored by Geant4.";
 362     TG4Globals::Warning(text);
 363   }
 364 }
 365
 366
 367 void TG4GeometryManager::Matrix(Int_t& krot, Float_t thetaX, Float_t phiX,
 368            Float_t thetaY, Float_t phiY, Float_t thetaZ, Float_t phiZ)
 369 {
 370 // Creates G4RotationMatrix.
 371 // ---
 372
 373   krot = ++fMatrixCounter;
 374
 375   // write token to the output file
 376   if (fWriteGeometry)
 377     fOutputManager->WriteGsrotm(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 378
 379   G4gsrotm(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 380 }
 381
 382
 383 G4Material* TG4GeometryManager::MixMaterials(G4String name, G4double density,
 384         TG4StringVector* matNames, TG4doubleVector* matWeights)
 385 {
 386 // Creates a mixture of selected materials
 387 // ---
 388
 389   // number of materials to be mixed
 390   G4int nofMaterials = matNames->entries();
 391   if (nofMaterials != matWeights->entries()) {
 392     G4String text = "TG4GeometryManager::MixMaterials: ";
 393     text = text +  "different number of material names and weigths.";
 394     TG4Globals::Exception(text);
 395   }
 396   // add verbose
 397   // G4cout << "Nof of materials to be mixed: " << nofMaterials << endl;
 398
 399   // fill vector of materials
 400   TG4MaterialVector matVector;
 401   G4int im;
 402   for (im=0; im< nofMaterials; im++) {
 403     // material
 404     G4Material* material = G4Material::GetMaterial((*matNames)[im]);
 405     matVector.insert(material);
 406   }
 407
 408   // create the mixed material
 409   G4Material* mixture = new G4Material(name, density, nofMaterials);
 410   for (im=0; im< nofMaterials; im++) {
 411     G4Material* material = matVector[im];
 412     G4double fraction = (*matWeights)[im];
 413     mixture->AddMaterial(material, fraction);
 414   }
 415
 416   return mixture;
 417 }
 418
 419
 420 void TG4GeometryManager::Ggclos()
 421 {
 422 // Sets the top VTE in temporary G3 volume table.
 423 // Close geometry output file (if fWriteGeometry is set true).
 424 //
 425 //  Geant3 desription:
 426 //  ==================
 427 // close out the geometry
 428 // ---
 429
 430   if (fWriteGeometry) fOutputManager->WriteGgclos();
 431
 432   G4ggclos();
 433 }
 434
 435
 436 void TG4GeometryManager::Gfmate(Int_t imat, char *name, Float_t &a,
 437           Float_t &z, Float_t &dens, Float_t &radl, Float_t &absl,
 438           Float_t* ubuf, Int_t& nbuf)
 439 {
 440 //  Geant3 desription:
 441 //  ==================
 442 // Return parameters for material IMAT
 443 // ---
 444
 445   G4Material* material = G3Mat.get(imat);
 446
 447   if (material) {
 448     // to do: change this correctly
 449     // !! unsafe conversion
 450     const char* chName = material->GetName();
 451     name = (char*)chName;
 452     a = GetEffA(material);
 453     z = GetEffZ(material);
 454
 455     dens = material->GetDensity();
 456     dens /= TG3Units::MassDensity();
 457
 458     radl = material->GetRadlen();
 459     radl /= TG3Units::Length();
 460
 461     // the following parameters are not defined in Geant4
 462     absl = 0.;
 463     ubuf = 0;
 464     nbuf = 0;
 465   }
 466   else {
 467     TG4Globals::Exception(
 468      "TG4GeometryManager::Gfmate: material has not been found.");
 469   }
 470 }
 471
 472
 473 void  TG4GeometryManager::Gstpar(Int_t itmed, const char *param,
 474            Float_t parval)
 475 {
 476 // Passes the tracking medium parameter to TG4Limits.
 477 // The tracking medium parameter is set only in case
 478 // its value is different from the "global" physics setup.
 479 // (If: CheckCut/FlagWithG3Defaults is used checking
 480 //  is performed with respect to G3 default values.)
 481 // When any cut/flag parameter is set in limits
 482 // the physics manager is locked and the physics setup
 483 // cannot be changed.
 484 // Applying this TG4Limits to particles and physical
 485 // processes is still in development.
 486 //
 487 //  Geant3 desription:
 488 //  ==================
 489 //  To change the value of cut  or mechanism "CHPAR"
 490 //  to a new value PARVAL  for tracking medium ITMED
 491 //  The  data   structure  JTMED   contains  the   standard  tracking
 492 //  parameters (CUTS and flags to control the physics processes)  which
 493 //  are used  by default  for all  tracking media.   It is  possible to
 494 //  redefine individually  with GSTPAR  any of  these parameters  for a
 495 //  given tracking medium.
 496 //  ITMED     tracking medium number
 497 //  CHPAR     is a character string (variable name)
 498 //  PARVAL    must be given as a floating point.
 499 // ---
 500
 501   // write token to the output file
 502   if (fWriteGeometry)
 503     fOutputManager->WriteGstpar(itmed, param, parval);
 504
 505   // get physics setup
 506   TG4PhysicsManager* physicsManager = TG4PhysicsManager::Instance();
 507   //TG4CutVector* cutVector = physicsManager->GetCutVector();
 508   //TG4FlagVector* flagVector = physicsManager->GetFlagVector();
 509
 510   G4String name = CutName(param);
 511   TG3Cut g3Cut;
 512   if (physicsManager->CheckCutWithCutVector(name, parval, g3Cut)) {
 513       GstparCut(itmed, g3Cut, parval);
 514       physicsManager->Lock();
 515   }
 516   else {
 517     TG3Flag g3Flag;
 518     if (physicsManager->CheckFlagWithFlagVector(name, parval, g3Flag)) {
 519       GstparFlag(itmed, g3Flag, parval);
 520       physicsManager->Lock();
 521     }
 522     else if (g3Cut==kNoG3Cuts && g3Flag==kNoG3Flags) {
 523       G4String text = "TG4GeometryManager::Gstpar:";
 524       text = text + name;
 525       text = text + " parameter is not yet implemented.";
 526       TG4Globals::Warning(text);
 527     }
 528   }
 529 }
 530
 531
 532 void  TG4GeometryManager::Gsckov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 533                          Float_t* absco, Float_t* effic, Float_t* rindex)
 534 {
 535 //
 536 //  Geant3 desription:
 537 //  ==================
 538 //
 539 //    Stores the tables for UV photon tracking in medium ITMED
 540 //    Please note that it is the user's responsability to
 541 //    provide all the coefficients:
 542 //
 543 //
 544 //       ITMED       Tracking medium number
 545 //       NPCKOV      Number of bins of each table
 546 //       PPCKOV      Value of photon momentum (in GeV)
 547 //       ABSCO       Absorbtion coefficients
 548 //                   dielectric: absorbtion length in cm
 549 //                   metals    : absorbtion fraction (0<=x<=1)
 550 //       EFFIC       Detection efficiency for UV photons
 551 //       RINDEX      Refraction index (if=0 metal)
 552 // ---
 553
 554   G4double* ppckovDbl = CreateG4doubleArray(ppckov, npckov);
 555   G4double* abscoDbl  = CreateG4doubleArray(absco, npckov);
 556   G4double* efficDbl  = CreateG4doubleArray(effic, npckov);
 557   G4double* rindexDbl = CreateG4doubleArray(rindex, npckov);
 558
 559   // add units
 560   G4int i;
 561   for (i=0; i<npckov; i++) {
 562     ppckovDbl[i] = ppckovDbl[i]*TG3Units::Energy();
 563     abscoDbl[i]  = abscoDbl[i]*TG3Units::Length();
 564   }
 565
 566   // create material properties table
 567   G4MaterialPropertiesTable* table = new G4MaterialPropertiesTable();
 568   table->AddProperty("ABSLENGTH", ppckovDbl, abscoDbl, npckov);
 569                     // used in G4OpAbsorption process
 570   table->AddProperty("EFFICIENCY", ppckovDbl, efficDbl, npckov);
 571                     // used in G4OpBoundary process
 572   table->AddProperty("RINDEX", ppckovDbl, rindexDbl, npckov);
 573                     // used in G4Cerenkov, G4OpRayleigh, G4OpBoundary
 574
 575   // get material of medium from table
 576   G3MedTableEntry* medium = G3Med.get(itmed);
 577   if (!medium) {
 578     G4String text = "TG4GeometryManager::Gsckov: \n";
 579     text = text + "    Medium not found.";
 580     G4Exception(text);
 581   }
 582   G4Material* material = medium->GetMaterial();
 583
 584   // set material properties table
 585   material->SetMaterialPropertiesTable(table);
 586
 587   G4cout << "The tables for UV photon tracking set for "
 588          << material->GetName() << endl;
 589   for (i=0; i<npckov; i++)
 590     G4cout << ppckovDbl[i] << " " << rindexDbl[i] << endl;
 591
 592   delete ppckovDbl;
 593   delete abscoDbl;
 594   delete efficDbl;
 595   delete rindexDbl;
 596 }
 597
 598
 599 void  TG4GeometryManager::Gsdvn(const char *name, const char *mother,
 600            Int_t ndiv, Int_t iaxis)
 601 {
 602 //  Geant3 desription:
 603 //  ==================
 604 //  NAME   Volume name
 605 //  MOTHER Mother volume name
 606 //  NDIV   Number of divisions
 607 //  IAXIS  Axis value
 608 //
 609 //  X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS.
 610 //  It divides a previously defined volume.
 611 //  ---
 612
 613     // write token to the output file
 614     if (fWriteGeometry)
 615       fOutputManager->WriteGsdvn(name, mother, ndiv, iaxis);
 616
 617     G4gsdvn(CutName(name), CutName(mother), ndiv, iaxis);
 618
 619     // register name in name map
 620     fNameMap.AddName(CutName(name));
 621 }
 622
 623
 624 void  TG4GeometryManager::Gsdvn2(const char *name, const char *mother,
 625            Int_t ndiv, Int_t iaxis, Float_t c0i, Int_t numed)
 626 {
 627 //  Geant3 desription:
 628 //  ==================
 629 //  DIVIDES MOTHER INTO NDIV DIVISIONS CALLED NAME
 630 //  ALONG AXIS IAXIS STARTING AT COORDINATE VALUE C0.
 631 //  THE NEW VOLUME CREATED WILL BE MEDIUM NUMBER NUMED.
 632 // ---
 633
 634     // write token to the output file
 635     if (fWriteGeometry)
 636       fOutputManager->WriteGsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 637
 638     G4gsdvn2(CutName(name),CutName(mother), ndiv, iaxis, c0i, numed);
 639
 640     // register name in name map
 641     fNameMap.AddName(CutName(name));
 642 }
 643
 644
 645 void  TG4GeometryManager::Gsdvt(const char *name, const char *mother,
 646            Float_t step, Int_t iaxis, Int_t numed, Int_t ndvmx)
 647 {
 648 //  Geant3 desription:
 649 //  ==================
 650 //       Divides MOTHER into divisions called NAME along
 651 //       axis IAXIS in steps of STEP. If not exactly divisible
 652 //       will make as many as possible and will centre them
 653 //       with respect to the mother. Divisions will have medium
 654 //       number NUMED. If NUMED is 0, NUMED of MOTHER is taken.
 655 //       NDVMX is the expected maximum number of divisions
 656 //          (If 0, no protection tests are performed)
 657 // ---
 658
 659     // write token to the output file
 660     if (fWriteGeometry)
 661       fOutputManager->WriteGsdvt(name, mother, step, iaxis, numed, ndvmx);
 662
 663     G4gsdvt(CutName(name), CutName(mother), step, iaxis, numed, ndvmx);
 664
 665     // register name in name map
 666     fNameMap.AddName(CutName(name));
 667 }
 668
 669
 670 void  TG4GeometryManager::Gsdvt2(const char *name, const char *mother,
 671            Float_t step, Int_t iaxis, Float_t c0, Int_t numed, Int_t ndvmx)
 672 {
 673 //  Geant3 desription:
 674 //  ==================
 675 // Create a new volume by dividing an existing one
 676 //
 677 //           Divides MOTHER into divisions called NAME along
 678 //            axis IAXIS starting at coordinate value C0 with step
 679 //            size STEP.
 680 //           The new volume created will have medium number NUMED.
 681 //           If NUMED is 0, NUMED of mother is taken.
 682 //           NDVMX is the expected maximum number of divisions
 683 //             (If 0, no protection tests are performed)
 684 // ---
 685
 686     // write token to the output file
 687     if (fWriteGeometry)
 688       fOutputManager->WriteGsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 689
 690     G4gsdvt2(CutName(name), CutName(mother), step, iaxis, c0, numed, ndvmx);
 691
 692     // register name in name map
 693     fNameMap.AddName(CutName(name));
 694 }
 695
 696
 697 void  TG4GeometryManager::Gsord(const char *name, Int_t iax)
 698 {
 699 // No corresponding action in G4.
 700 //
 701 //  Geant3 desription:
 702 //  ==================
 703 //    Flags volume CHNAME whose contents will have to be ordered
 704 //    along axis IAX, by setting the search flag to -IAX
 705 //           IAX = 1    X axis
 706 //           IAX = 2    Y axis
 707 //           IAX = 3    Z axis
 708 //           IAX = 4    Rxy (static ordering only  -> GTMEDI)
 709 //           IAX = 14   Rxy (also dynamic ordering -> GTNEXT)
 710 //           IAX = 5    Rxyz (static ordering only -> GTMEDI)
 711 //           IAX = 15   Rxyz (also dynamic ordering -> GTNEXT)
 712 //           IAX = 6    PHI   (PHI=0 => X axis)
 713 //           IAX = 7    THETA (THETA=0 => Z axis)
 714 // ---
 715
 716   TG4Globals::Warning("TG4GeometryManager::Gsord: dummy method.");
 717 }
 718
 719
 720 void  TG4GeometryManager::Gspos(const char *vname, Int_t num,
 721           const char *vmoth, Float_t x, Float_t y, Float_t z, Int_t irot,
 722           const char *vonly)
 723 {
 724 //  Geant3 desription:
 725 //  ==================
 726 // Position a volume into an existing one
 727 //
 728 //  NAME   Volume name
 729 //  NUMBER Copy number of the volume
 730 //  MOTHER Mother volume name
 731 //  X      X coord. of the volume in mother ref. sys.
 732 //  Y      Y coord. of the volume in mother ref. sys.
 733 //  Z      Z coord. of the volume in mother ref. sys.
 734 //  IROT   Rotation matrix number w.r.t. mother ref. sys.
 735 //  ONLY   ONLY/MANY flag
 736 //
 737 //  It positions a previously defined volume in the mother.
 738 // ---
 739
 740    // write token to the output file
 741    if (fWriteGeometry)
 742      fOutputManager->WriteGspos(vname, num, vmoth, x, y, z, irot, vonly);
 743
 744    G4gspos(CutName(vname), num, CutName(vmoth), x, y, z, irot, vonly);
 745
 746    // register name in name map
 747    fNameMap.AddName(CutName(vname));
 748 }
 749
 750
 751 void  TG4GeometryManager::Gsposp(const char *name, Int_t nr,
 752            const char *mother, Float_t x, Float_t y, Float_t z, Int_t irot,
 753            const char *konly, Float_t *upar, Int_t np )
 754 {
 755 //  Geant3 desription:
 756 //  ==================
 757 //      Place a copy of generic volume NAME with user number
 758 //      NR inside MOTHER, with its parameters UPAR(1..NP)
 759 // ---
 760
 761    G4double* parin = CreateG4doubleArray(upar, np);
 762
 763    // write token to the output file
 764    if (fWriteGeometry)
 765      fOutputManager->WriteGsposp(name, nr, mother, x, y, z, irot, konly, parin, np);
 766
 767    G4gsposp(CutName(name), nr, CutName(mother), x, y, z, irot, konly,
 768              parin, np);
 769
 770    delete [] parin;
 771
 772    // register name in name map
 773    fNameMap.AddName(CutName(name));
 774 }
 775
 776
 777 Int_t TG4GeometryManager::Gsvolu(const char *name, const char *shape,
 778           Int_t nmed, Float_t *upar, Int_t npar)
 779 {
 780 //  Geant3 desription:
 781 //  ==================
 782 //  NAME   Volume name
 783 //  SHAPE  Volume type
 784 //  NUMED  Tracking medium number
 785 //  NPAR   Number of shape parameters
 786 //  UPAR   Vector containing shape parameters
 787 //
 788 //  It creates a new volume in the JVOLUM data structure.
 789 // ---
 790
 791   G4double* parin = CreateG4doubleArray(upar, npar);
 792
 793   // write token to the output file
 794   if (fWriteGeometry)
 795     fOutputManager->WriteGsvolu(name, shape, nmed, parin, npar);
 796
 797   G4gsvolu(CutName(name), CutName(shape), nmed, parin, npar);
 798
 799   delete [] parin;
 800
 801   // register name in name map
 802   fNameMap.AddName(CutName(name));
 803
 804   return 0;
 805 }
 806
 807
 808 void TG4GeometryManager::WriteEuclid(const char* fileName,
 809           const char* topVolName, Int_t number, Int_t nlevel)
 810 {
 811 //  Geant3 desription:
 812 //  ==================
 813 //
 814 //     ******************************************************************
 815 //     *                                                                *
 816 //     *  Write out the geometry of the detector in EUCLID file format  *
 817 //     *                                                                *
 818 //     *       filnam : will be with the extension .euc                 *
 819 //     *       topvol : volume name of the starting node                *
 820 //     *       number : copy number of topvol (relevant for gsposp)     *
 821 //     *       nlevel : number of  levels in the tree structure         *
 822 //     *                to be written out, starting from topvol         *
 823 //     *                                                                *
 824 //     *       Author : M. Maire                                        *
 825 //     *                                                                *
 826 //     ******************************************************************
 827 //
 828 //     File filnam.tme is written out with the definitions of tracking
 829 //     medias and materials.
 830 //     As to restore original numbers for materials and medias, program
 831 //     searches in the file euc_medi.dat and comparing main parameters of
 832 //     the mat. defined inside geant and the one in file recognizes them
 833 //     and is able to take number from file. If for any material or medium,
 834 //     this procedure fails, ordering starts from 1.
 835 //     Arrays IOTMED and IOMATE are used for this procedure
 836 // ---
 837
 838   TG4Globals::Warning(
 839     "TG4GeometryManager::WriteEuclid(..) is not yet implemented.");
 840 }
 841
 842
 843 Int_t TG4GeometryManager::VolId(const Text_t* volName) const
 844 {
 845 // Returns the sensitive detector identifier.
 846 // !! Gives exception in case logical volume is not associated with
 847 // a sensitive detector.
 848 // ---
 849
 850   G4String g4VolName = CutName(volName);
 851   G4LogicalVolumeStore* pLVStore = G4LogicalVolumeStore::GetInstance();
 852
 853   for (G4int i=0; i<pLVStore->entries(); i++) {
 854     G4LogicalVolume* lv = pLVStore->at(i);
 855     G4VSensitiveDetector* sd = lv->GetSensitiveDetector();
 856
 857     if ((sd) && (sd->GetName()==g4VolName)) {
 858       TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
 859       if (tsd)
 860         return tsd->GetID();
 861       else {
 862         TG4Globals::Exception(
 863           "TG4GeometryManager::VolId: Unknown sensitive detector type");
 864         return 0;
 865       }
 866     }
 867   }
 868
 869   G4String text = "TG4GeometryManager::VolId: Sensitive detector ";
 870   text = text + g4VolName;
 871   text = text + " is not defined.\n";
 872   text = text + "    Set /alDet/setAllSensitive true in PreInit.";
 873   TG4Globals::Exception(text);
 874   return 0;
 875 }
 876
 877
 878 const char* TG4GeometryManager::VolName(Int_t id) const
 879 {
 880 // Returns the name of the sensitive detector with the given identifier.
 881 // !! Gives exception in case logical volume is not associated with
 882 // a sensitive detector.
 883 // ---
 884
 885   G4LogicalVolumeStore* pLVStore = G4LogicalVolumeStore::GetInstance();
 886
 887   for (G4int i=0; i<pLVStore->entries(); i++) {
 888     G4LogicalVolume* lv = pLVStore->at(i);
 889     G4VSensitiveDetector* sd = lv->GetSensitiveDetector();
 890
 891     if (sd) {
 892       G4int sdID;
 893       TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
 894       if (tsd)
 895         sdID = tsd->GetID();
 896       else {
 897         TG4Globals::Exception(
 898           "TG4GeometryManager::VolId: Unknown sensitive detector type");
 899         return 0;
 900       }
 901       if (sdID == id) return sd->GetName();
 902     }
 903   }
 904
 905   G4String text = "TG4GeometryManager::VolName:\n";
 906   text = text + "    Sensitive detector with given id is not defined. \n";
 907   text = text + "    Set /alDet/setAllSensitive true in PreInit.";
 908   TG4Globals::Exception(text);
 909   return "";
 910 }
 911
 912
 913 Int_t TG4GeometryManager::NofVolumes() const
 914 {
 915 // Returns the total number of sensitive detectors.
 916 // ---
 917
 918   return NofSensitiveDetectors();
 919 }
 920
 921
 922 void TG4GeometryManager::ReadG3Geometry(G4String filePath)
 923 {
 924 // Processes g3calls.dat file and fills G3 tables.
 925 // ---
 926
 927   // add verbose
 928   G4cout << "Reading the call list file " << filePath << "..." << endl;
 929   G3CLRead(filePath, NULL);
 930   G4cout << "Call list file read completed. Build geometry" << endl;
 931 }
 932
 933
 934 //=============================================================================
 935 //
 936 // public methods - Geant4 only
 937 //
 938 //=============================================================================
 939
 940
 941 G4VPhysicalVolume* TG4GeometryManager::CreateG4Geometry()
 942 {
 943 // Creates G4 geometry objects according to the G3VolTable
 944 // and returns the top physical volume in case it was created
 945 // (return zero otherwise).
 946 // ---
 947
 948   // set the first entry in the G3Vol table
 949   Ggclos();
 950   G3VolTableEntry* first = G3Vol.GetFirstVTE();
 951
 952   // close g3calls.dat
 953   if (fWriteGeometry) fOutputManager->CloseFile();
 954
 955   // create G4 geometry
 956   G3toG4BuildTree(first,0);
 957
 958   // print G3 volume table statistics
 959   G3Vol.VTEStat();
 960
 961   // print G4 geometry statistics
 962   G4cout << "G4 Stat: instantiated "
 963          << NofG4LogicalVolumes()  << " logical volumes \n"
 964          << "                      "
 965          << NofG4PhysicalVolumes() << " physical volumes" << endl;
 966
 967   // position the first entry
 968   // (in Geant3 the top volume cannot be positioned)
 969   //
 970   G4VPhysicalVolume* top = 0;
 971   if (first->GetLV()->GetNoDaughters() == 0) {
 972     top = new G4PVPlacement(0, G4ThreeVector(), first->GetName(),
 973                             first->GetLV(), 0, false, 0);
 974   }
 975   return top;
 976 }
 977
 978
 979 void TG4GeometryManager::UseG3TrackingMediaLimits()
 980 {
 981 // Sets fUseG3TMLimits option.
 982 // !! This method has to be called only before starting
 983 // creating geometry.
 984 // ---
 985
 986   if (fMediumCounter == 0) {
 987     fUseG3TMLimits = true;
 988   }
 989   else {
 990     G4String text = "TG4GeometryManager::UseG3TMLimits: \n";
 991     text = text + "    It is too late to set G3 defaults. \n";
 992     text = text + "    Some media has been already processed.";
 993     TG4Globals::Exception(text);
 994   }
 995 }
 996
 997
 998 void TG4GeometryManager::ClearG3Tables()
 999 {
1000 // Clears G3 volumes, materials, rotations(?) tables
1001 // and sensitive volumes vector.
1002 // The top volume is kept in the vol table.
1003 // ---
1004
1005   // clear volume table
1006   // but keep the top volume in the table
1007   G3VolTableEntry* top = G3Vol.GetFirstVTE();
1008   G4String name = top->GetName();
1009   G4String shape = top->GetShape();
1010   G3VolTableEntry* keep
1011     = new G3VolTableEntry(name, shape, top->GetRpar(), top->GetNpar(),
1012                            top->GetNmed(), top->GetSolid(), false);
1013   keep->SetLV(top->GetLV());
1014   G3Vol.Clear();
1015   G3Vol.PutVTE(keep);
1016
1017   // clear other tables
1018   G3Mat.Clear();
1019   //G3Rot.Clear();
1020   G3SensVol.clear();
1021 }
1022
1023
1024 void TG4GeometryManager::ClearG3TablesFinal()
1025 {
1026 // Clears G3 medias and volumes tables
1027 // (the top volume is removed from the vol table)
1028 // ---
1029
1030   // fill medium id vector
1031   FillMediumIdVector();
1032
1033   G3Med.Clear();
1034   G3Vol.Clear();
1035
1036   // reset medium counter
1037   //fMaterialCounter = 0;
1038   fMediumCounter = 0;
1039 }
1040
1041
1042 void TG4GeometryManager::OpenOutFile(G4String filePath)
1043 {
1044 // Opens output file.
1045 // ---
1046
1047   fOutputManager->OpenFile(filePath);
1048 }
1049
1050
1051 void TG4GeometryManager::PrintNameMap()
1052 {
1053 // Prints the map of volumes names to second names.
1054 // ---
1055
1056   fNameMap.PrintAll();
1057 }
1058
1059
1060 void TG4GeometryManager::SetWriteGeometry(G4bool writeGeometry)
1061 {
1062 // Controls geometry output.
1063 // ---
1064
1065   fWriteGeometry = writeGeometry;
1066 }
1067
1068
1069 void TG4GeometryManager::SetMapSecond(const G4String& name)
1070 {
1071 // Sets the second name for the map of volumes names.
1072 // ---
1073
1074   fNameMap.SetSecond(name);
1075 }
1076
1077
1078 Int_t TG4GeometryManager::NofG3Volumes() const
1079 {
1080 // Returns the total number of logical volumes corresponding
1081 // to G3 volumes. (
1082 // The logical volume that were created by Gsposp method
1083 // with a generic name (name_copyNo) are NOT included.
1084 // ---
1085
1086   G4LogicalVolumeStore* pLVStore = G4LogicalVolumeStore::GetInstance();
1087
1088   G4int counter = 0;
1089   for (G4int i=0; i<pLVStore->entries(); i++) {
1090     G4LogicalVolume* lv = (*pLVStore)[i];
1091     if (IsG3Volume(lv->GetName())) counter++;
1092   }
1093
1094   return counter;
1095 }
1096
1097
1098 Int_t TG4GeometryManager::NofG4LogicalVolumes() const
1099 {
1100 // Returns the total number of logical volumes in the geometry.
1101 // ---
1102
1103   G4LogicalVolumeStore* pLVStore = G4LogicalVolumeStore::GetInstance();
1104   return pLVStore->entries();
1105 }
1106
1107
1108 Int_t TG4GeometryManager::NofG4PhysicalVolumes() const
1109 {
1110 // Returns the total number of physical volumes in the geometry.
1111 // ---
1112
1113   G4LogicalVolumeStore* pLVStore = G4LogicalVolumeStore::GetInstance();
1114
1115   G4int counter = 0;
1116   for (G4int i=0; i<pLVStore->entries(); i++) {
1117     counter += ((*pLVStore)[i])->GetNoDaughters();
1118   }
1119
1120   return counter;
1121 }
1122
1123
1124 Int_t TG4GeometryManager::NofSensitiveDetectors() const
1125 {
1126 // Returns the total number of sensitive detectors.
1127 // ---
1128
1129   return TG4VSensitiveDetector::GetTotalNofSensitiveDetectors();
1130 }
1131
1132
1133 G4bool TG4GeometryManager::IsG3Volume(G4String lvName) const
1134 {
1135 // Returns true if the logical volume of given volumeName
1136 // was not created by Gsposp method with a generic name
1137 // (name_copyNo).
1138 // ---
1139
1140   if (lvName.contains(gSeparator))
1141     return false;
1142   else
1143     return true;
1144 }
1145
1146
1147 void TG4GeometryManager::G4ToG3VolumeName(G4String& name) const
1148 {
1149 // Cuts _copyNo extension added to logical volume name in case
1150 // the logical volume was created by Gsposp method.
1151 // ---
1152
1153   if (name.contains(gSeparator))
1154   name = name(0,name.first(gSeparator));
1155 }
1156
1157
1158 const G4String& TG4GeometryManager::GetMapSecond(const G4String& name)
1159 {
1160 // Returns the second string associated with the name in
1161 // the name map.
1162 // ---
1163
1164   return fNameMap.GetSecond(name);
1165 }
1166
1167
1168 G4int TG4GeometryManager::GetMediumId(G4Material* material) const
1169 {
1170 // Returns the second index for materials (having its origin in
1171 // G4 tracking media concept)
1172 // ---
1173
1174   return fMediumIdVector[material->GetIndex()];
1175 }
1176
1177
1178 G4double TG4GeometryManager::GetEffA(G4Material* material) const
1179 {
1180 // Returns A or effective A=sum(pi*Ai) (if compound/mixture)
1181 // of given material.
1182 // ---
1183
1184   G4double a = 0.;
1185   G4int nofElements = material->GetNumberOfElements();
1186   if (nofElements > 1) {
1187     G4String text = "Effective A for material mixture (";
1188     text = text + material->GetName();
1189     text = text + ") is used.";
1190     //TG4Globals::Warning(text);
1191
1192     for (G4int i=0; i<nofElements; i++) {
1193       G4double aOfElement = material->GetElement(i)->GetA();
1194       G4double massFraction = material->GetFractionVector()[i];
1195       a += aOfElement*massFraction /(TG3Units::AtomicWeight());
1196     }
1197   }
1198   else {
1199     a = material->GetA();
1200     a /= TG3Units::AtomicWeight();
1201   }
1202   return a;
1203 }
1204
1205
1206 G4double TG4GeometryManager::GetEffZ(G4Material* material) const
1207 {
1208 // Returns Z or effective Z=sum(pi*Zi) (if compound/mixture)
1209 // of given material.
1210 // ---
1211
1212   G4double z = 0.;
1213   G4int nofElements = material->GetNumberOfElements();
1214   if (nofElements > 1) {
1215     G4String text = "Effective Z for material mixture (";
1216     text = text + material->GetName();
1217     text = text + ") is used.";
1218     //TG4Globals::Warning(text);
1219
1220     for (G4int i=0; i<nofElements; i++) {
1221       G4double zOfElement = material->GetElement(i)->GetZ();
1222       G4double massFraction = material->GetFractionVector()[i];
1223       z += zOfElement*massFraction;
1224     }
1225   }
1226   else {
1227     z = material->GetZ();
1228   }
1229   return z;
1230 }