TFluka/TFlukaMCGeometry.cxx

   1 // @(#):$Name$:$Id$
   2 // Author: Andrei Gheata 10/07/2003
   3
   4 #include "TObjString.h"
   5 #include "TFluka.h"
   6 //#include "TVirtualMCApplication.h"
   7 #include "TFlukaMCGeometry.h"
   8 #include "TGeoManager.h"
   9 #include "TGeoVolume.h"
  10
  11 #include "TCallf77.h"
  12
  13 #ifndef WIN32
  14 # define idnrwr idnrwr_
  15 # define g1wr   g1wr_
  16 # define g1rtwr g1rtwr_
  17 # define conhwr conhwr_
  18 # define inihwr inihwr_
  19 # define jomiwr jomiwr_
  20 # define lkdbwr lkdbwr_
  21 # define lkfxwr lkfxwr_
  22 # define lkmgwr lkmgwr_
  23 # define lkwr lkwr_
  24 # define magfld magfld_
  25 # define nrmlwr nrmlwr_
  26 # define rgrpwr rgrpwr_
  27 # define isvhwr isvhwr_
  28
  29 #else
  30
  31 # define idnrwr IDNRWR
  32 # define g1wr   G1WR
  33 # define g1rtwr G1RTWR
  34 # define conhwr CONHWR
  35 # define inihwr INIHWR
  36 # define jomiwr JOMIWR
  37 # define lkdbwr LKDBWR
  38 # define lkfxwr LKFXWR
  39 # define lkmgwr LKMGWR
  40 # define lkwr   LKWR
  41 # define magfld MAGFLD
  42 # define nrmlwr NRMLWR
  43 # define rgrpwr RGRPWR
  44 # define isvhwr ISVHWR
  45
  46 #endif
  47
  48 //____________________________________________________________________________
  49 extern "C"
  50 {
  51    //
  52    // Prototypes for FLUKA navigation methods
  53    //
  54    Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
  55    void  type_of_call   g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  56                              Double_t * /*pV*/,  Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
  57                              Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
  58                                   Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
  59                              Double_t *s /*Lt*/, Int_t * /*jrLt*/);
  60
  61    void  type_of_call g1rtwr();
  62    void  type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/);
  63    void  type_of_call inihwr(Int_t & /*intHist*/);
  64    void  type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
  65                              Int_t & /*flukaReg*/);
  66    void  type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  67                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  68                              Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
  69    void  type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  70                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  71                              Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
  72    void  type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  73                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  74                                        Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
  75    void  type_of_call   lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  76                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  77                                   Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
  78 //   void  type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
  79 //                             Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
  80 //                             Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
  81    void  type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  82                              Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
  83                                   Double_t * /*norml*/, const Int_t & /*oldReg*/,
  84                                   const Int_t & /*newReg*/, Int_t & /*flagErr*/);
  85    void  type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
  86                              Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
  87    Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
  88 };
  89
  90 // TFluka global pointer
  91 TFluka *fluka = 0;
  92 TFlukaMCGeometry *mcgeom = 0;
  93
  94 ClassImp(TFlukaMCGeometry)
  95
  96 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
  97
  98 //_____________________________________________________________________________
  99 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
 100   : TVirtualMCGeometry(name, title)
 101 {
 102   //
 103   // Standard constructor
 104   //
 105   fDebug        = kFALSE;
 106   fLastMaterial = 0;
 107   fNextRegion   = 0;
 108   fNextLattice  = 0;
 109   fRegionList   = 0;
 110   fluka = (TFluka*)gMC;
 111   mcgeom = this;
 112 }
 113
 114 //_____________________________________________________________________________
 115 TFlukaMCGeometry::TFlukaMCGeometry()
 116   : TVirtualMCGeometry()
 117 {
 118   //
 119   // Default constructor
 120   //
 121   fDebug        = kFALSE;
 122   fLastMaterial = 0;
 123   fNextRegion   = 0;
 124   fNextLattice  = 0;
 125   fRegionList   = 0;
 126   fluka = (TFluka*)gMC;
 127   mcgeom = this;
 128 }
 129
 130 //_____________________________________________________________________________
 131 TFlukaMCGeometry::~TFlukaMCGeometry()
 132 {
 133   //
 134   // Destructor
 135   //
 136   fgInstance=0;
 137   if (fRegionList) delete [] fRegionList;
 138   if (gGeoManager) delete gGeoManager;
 139 }
 140
 141 //
 142 // private methods
 143 //
 144 //_____________________________________________________________________________
 145 TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
 146   : TVirtualMCGeometry()
 147 {
 148   //
 149   // Copy constructor
 150   //
 151 }
 152
 153 //_____________________________________________________________________________
 154 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
 155 {
 156 // Converts Float_t* array to Double_t*,
 157 // !! The new array has to be deleted by user.
 158 // ---
 159
 160   Double_t* doubleArray;
 161   if (size>0) {
 162     doubleArray = new Double_t[size];
 163     for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
 164   }
 165   else {
 166     //doubleArray = 0;
 167     doubleArray = new Double_t[1];
 168   }
 169   return doubleArray;
 170 }
 171 //
 172 // public methods
 173 //_____________________________________________________________________________
 174 void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 175                        Float_t &dens, Float_t &radl, Float_t &absl,
 176                        Float_t* /*ubuf*/, Int_t& /*nbuf*/)
 177 {
 178    if (fDebug) printf("Gfmate %i\n", imat);
 179    TGeoMaterial *mat;
 180    TIter next (gGeoManager->GetListOfMaterials());
 181    while ((mat = (TGeoMaterial*)next())) {
 182      if (mat->GetUniqueID() == (UInt_t)imat) break;
 183    }
 184    if (!mat) {
 185       Error("Gfmate", "no material with index %i found", imat);
 186       return;
 187    }
 188    sprintf(name, "%s", mat->GetName());
 189    a = mat->GetA();
 190    z = mat->GetZ();
 191    dens = mat->GetDensity();
 192    radl = mat->GetRadLen();
 193    absl = mat->GetIntLen();
 194    if (fDebug) printf("   ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
 195 }
 196
 197 //_____________________________________________________________________________
 198 void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 199                        Double_t &dens, Double_t &radl, Double_t &absl,
 200                        Double_t* /*ubuf*/, Int_t& /*nbuf*/)
 201 {
 202    if (fDebug) printf("Gfmate %i\n", imat);
 203    TGeoMaterial *mat;
 204    TIter next (gGeoManager->GetListOfMaterials());
 205    while ((mat = (TGeoMaterial*)next())) {
 206      if (mat->GetUniqueID() == (UInt_t)imat) break;
 207    }
 208    if (!mat) {
 209       Error("Gfmate", "no material with index %i found", imat);
 210       return;
 211    }
 212    sprintf(name, "%s", mat->GetName());
 213    a = mat->GetA();
 214    z = mat->GetZ();
 215    dens = mat->GetDensity();
 216    radl = mat->GetRadLen();
 217    absl = mat->GetIntLen();
 218    if (fDebug) printf("   ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
 219 }
 220
 221 //_____________________________________________________________________________
 222 void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z,
 223                        Double_t dens, Double_t radl, Double_t absl, Float_t* buf,
 224                        Int_t nwbuf)
 225 {
 226   //
 227   // Defines a Material
 228   //
 229   //  kmat               number assigned to the material
 230   //  name               material name
 231   //  a                  atomic mass in au
 232   //  z                  atomic number
 233   //  dens               density in g/cm3
 234   //  absl               absorbtion length in cm
 235   //                     if >=0 it is ignored and the program
 236   //                     calculates it, if <0. -absl is taken
 237   //  radl               radiation length in cm
 238   //                     if >=0 it is ignored and the program
 239   //                     calculates it, if <0. -radl is taken
 240   //  buf                pointer to an array of user words
 241   //  nbuf               number of user words
 242   //
 243
 244   Double_t* dbuf = CreateDoubleArray(buf, nwbuf);
 245   Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 246   delete [] dbuf;
 247 }
 248
 249 //_____________________________________________________________________________
 250 void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z,
 251                        Double_t dens, Double_t radl, Double_t absl, Double_t* /*buf*/,
 252                        Int_t /*nwbuf*/)
 253 {
 254   //
 255   // Defines a Material
 256   //
 257   //  kmat               number assigned to the material
 258   //  name               material name
 259   //  a                  atomic mass in au
 260   //  z                  atomic number
 261   //  dens               density in g/cm3
 262   //  absl               absorbtion length in cm
 263   //                     if >=0 it is ignored and the program
 264   //                     calculates it, if <0. -absl is taken
 265   //  radl               radiation length in cm
 266   //                     if >=0 it is ignored and the program
 267   //                     calculates it, if <0. -radl is taken
 268   //  buf                pointer to an array of user words
 269   //  nbuf               number of user words
 270   //
 271
 272   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 273   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 274   if (fDebug) printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens);
 275 }
 276
 277 //_____________________________________________________________________________
 278 void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Float_t* a, Float_t* z,
 279                       Double_t dens, Int_t nlmat, Float_t* wmat)
 280 {
 281   //
 282   // Defines mixture OR COMPOUND IMAT as composed by
 283   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 284   //
 285   // If NLMAT > 0 then wmat contains the proportion by
 286   // weights of each basic material in the mixture.
 287   //
 288   // If nlmat < 0 then WMAT contains the number of atoms
 289   // of a given kind into the molecule of the COMPOUND
 290   // In this case, WMAT in output is changed to relative
 291   // weigths.
 292   //
 293
 294   Double_t* da = CreateDoubleArray(a, TMath::Abs(nlmat));
 295   Double_t* dz = CreateDoubleArray(z, TMath::Abs(nlmat));
 296   Double_t* dwmat = CreateDoubleArray(wmat, TMath::Abs(nlmat));
 297
 298   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 299   for (Int_t i=0; i<nlmat; i++) {
 300     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 301   }
 302
 303   delete [] da;
 304   delete [] dz;
 305   delete [] dwmat;
 306 }
 307
 308 //_____________________________________________________________________________
 309 void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Double_t* a, Double_t* z,
 310                       Double_t dens, Int_t nlmat, Double_t* wmat)
 311 {
 312   //
 313   // Defines mixture OR COMPOUND IMAT as composed by
 314   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 315   //
 316   // If NLMAT > 0 then wmat contains the proportion by
 317   // weights of each basic material in the mixture.
 318   //
 319   // If nlmat < 0 then WMAT contains the number of atoms
 320   // of a given kind into the molecule of the COMPOUND
 321   // In this case, WMAT in output is changed to relative
 322   // weigths.
 323   //
 324
 325   if (nlmat < 0) {
 326      nlmat = - nlmat;
 327      Double_t amol = 0;
 328      Int_t i;
 329      for (i=0;i<nlmat;i++) {
 330         amol += a[i]*wmat[i];
 331      }
 332      for (i=0;i<nlmat;i++) {
 333         wmat[i] *= a[i]/amol;
 334      }
 335   }
 336   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 337   if (fDebug) {
 338      printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens);
 339      for (Int_t j=0; j<nlmat; j++) printf("  Elem %i: z=%g  a=%g  w=%g\n",j,z[j],a[j],wmat[j]);
 340   }
 341   gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 342 }
 343 //_____________________________________________________________________________
 344 Int_t TFlukaMCGeometry::GetMedium() const
 345 {
 346 // Get current medium number
 347    Int_t imed = 0;
 348    TGeoNode *node = gGeoManager->GetCurrentNode();
 349    if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
 350    else       imed = node->GetVolume()->GetMedium()->GetId();
 351    if (fDebug) printf("GetMedium=%i\n", imed);
 352    return imed;
 353 }
 354
 355 //_____________________________________________________________________________
 356 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
 357 {
 358 // Returns FLUKA material index for medium IMED
 359    TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
 360    if (!med) {
 361       Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
 362       return -1;
 363    }
 364    Int_t imatfl = med->GetMaterial()->GetIndex();
 365    return imatfl;
 366 }
 367
 368 //_____________________________________________________________________________
 369 Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
 370 {
 371 // Get an ordered list of regions matching a given medium number
 372    nreg = 0;
 373    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
 374    TIter next(gGeoManager->GetListOfUVolumes());
 375    TGeoVolume *vol;
 376    Int_t imedium, ireg;
 377    while ((vol = (TGeoVolume*)next())) {
 378       imedium = vol->GetMedium()->GetId();
 379       if (imedium == imed) {
 380          ireg = vol->GetNumber();
 381          fRegionList[nreg++] = ireg;
 382       }
 383    }
 384    return fRegionList;
 385 }
 386
 387 //_____________________________________________________________________________
 388 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
 389 {
 390 // Get an ordered list of regions matching a given medium number
 391    nreg = 0;
 392    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
 393    TIter next(gGeoManager->GetListOfUVolumes());
 394    TGeoVolume *vol;
 395    Int_t imaterial, ireg;
 396    while ((vol = (TGeoVolume*)next())) {
 397       imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
 398       if (imaterial == imat) {
 399          ireg = vol->GetNumber();
 400          fRegionList[nreg++] = ireg;
 401       }
 402    }
 403    return fRegionList;
 404 }
 405 //_____________________________________________________________________________
 406 void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol,
 407                      Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 408                      Double_t stemax, Double_t deemax, Double_t epsil,
 409                      Double_t stmin, Float_t* ubuf, Int_t nbuf)
 410 {
 411   //
 412   //  kmed      tracking medium number assigned
 413   //  name      tracking medium name
 414   //  nmat      material number
 415   //  isvol     sensitive volume flag
 416   //  ifield    magnetic field
 417   //  fieldm    max. field value (kilogauss)
 418   //  tmaxfd    max. angle due to field (deg/step)
 419   //  stemax    max. step allowed
 420   //  deemax    max. fraction of energy lost in a step
 421   //  epsil     tracking precision (cm)
 422   //  stmin     min. step due to continuous processes (cm)
 423   //
 424   //  ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim;
 425   //  ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking
 426   //  performed with g3helix; ifield = 3 if tracking performed with g3helx3.
 427   //
 428
 429   //printf("Creating mediuma: %s, numed=%d, nmat=%d\n",name,kmed,nmat);
 430   Double_t* dubuf = CreateDoubleArray(ubuf, nbuf);
 431   Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil,
 432          stmin, dubuf, nbuf);
 433   delete [] dubuf;
 434 }
 435
 436 //_____________________________________________________________________________
 437 void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol,
 438                      Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 439                      Double_t stemax, Double_t deemax, Double_t epsil,
 440                      Double_t stmin, Double_t* /*ubuf*/, Int_t /*nbuf*/)
 441 {
 442   //
 443   //  kmed      tracking medium number assigned
 444   //  name      tracking medium name
 445   //  nmat      material number
 446   //  isvol     sensitive volume flag
 447   //  ifield    magnetic field
 448   //  fieldm    max. field value (kilogauss)
 449   //  tmaxfd    max. angle due to field (deg/step)
 450   //  stemax    max. step allowed
 451   //  deemax    max. fraction of energy lost in a step
 452   //  epsil     tracking precision (cm)
 453   //  stmin     min. step due to continuos processes (cm)
 454   //
 455   //  ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim;
 456   //  ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking
 457   //  performed with g3helix; ifield = 3 if tracking performed with g3helx3.
 458   //
 459
 460   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 461   gGeoManager->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin);
 462   if (fDebug) printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol);
 463 }
 464
 465 //_____________________________________________________________________________
 466 void TFlukaMCGeometry::Matrix(Int_t& krot, Double_t thex, Double_t phix, Double_t they,
 467                      Double_t phiy, Double_t thez, Double_t phiz)
 468 {
 469   //
 470   //  krot     rotation matrix number assigned
 471   //  theta1   polar angle for axis i
 472   //  phi1     azimuthal angle for axis i
 473   //  theta2   polar angle for axis ii
 474   //  phi2     azimuthal angle for axis ii
 475   //  theta3   polar angle for axis iii
 476   //  phi3     azimuthal angle for axis iii
 477   //
 478   //  it defines the rotation matrix number irot.
 479   //
 480
 481   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 482   gGeoManager->Matrix(krot, thex, phix, they, phiy, thez, phiz);
 483   if (fDebug) printf("Rotation %i defined\n", krot);
 484 }
 485
 486 //_____________________________________________________________________________
 487 Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed,
 488                       Float_t *upar, Int_t npar)
 489 {
 490   //
 491   //  NAME   Volume name
 492   //  SHAPE  Volume type
 493   //  NUMED  Tracking medium number
 494   //  NPAR   Number of shape parameters
 495   //  UPAR   Vector containing shape parameters
 496   //
 497   //  It creates a new volume in the JVOLUM data structure.
 498   //
 499
 500   Double_t* dupar = CreateDoubleArray(upar, npar);
 501   Int_t id = Gsvolu(name, shape, nmed, dupar, npar);
 502   delete [] dupar;
 503   return id;
 504 }
 505
 506 //_____________________________________________________________________________
 507 Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed,
 508                       Double_t *upar, Int_t npar)
 509 {
 510   //
 511   //  NAME   Volume name
 512   //  SHAPE  Volume type
 513   //  NUMED  Tracking medium number
 514   //  NPAR   Number of shape parameters
 515   //  UPAR   Vector containing shape parameters
 516   //
 517   //  It creates a new volume in the JVOLUM data structure.
 518   //
 519   char vname[5];
 520   Vname(name,vname);
 521   char vshape[5];
 522   Vname(shape,vshape);
 523
 524   TGeoVolume* vol = gGeoManager->Volume(vname, shape, nmed, upar, npar);
 525   if (fDebug) printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed);
 526   return vol->GetNumber();
 527 }
 528
 529 //_____________________________________________________________________________
 530 void  TFlukaMCGeometry::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 531                      Int_t iaxis)
 532 {
 533   //
 534   // Create a new volume by dividing an existing one
 535   //
 536   //  NAME   Volume name
 537   //  MOTHER Mother volume name
 538   //  NDIV   Number of divisions
 539   //  IAXIS  Axis value
 540   //
 541   //  X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS.
 542   //  It divides a previously defined volume.
 543   //
 544   char vname[5];
 545   Vname(name,vname);
 546   char vmother[5];
 547   Vname(mother,vmother);
 548
 549   gGeoManager->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, "n");
 550   if (fDebug) printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv);
 551 }
 552
 553 //_____________________________________________________________________________
 554 void  TFlukaMCGeometry::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 555                       Int_t iaxis, Double_t c0i, Int_t numed)
 556 {
 557   //
 558   // Create a new volume by dividing an existing one
 559   //
 560   // Divides mother into ndiv divisions called name
 561   // along axis iaxis starting at coordinate value c0.
 562   // the new volume created will be medium number numed.
 563   //
 564   char vname[5];
 565   Vname(name,vname);
 566   char vmother[5];
 567   Vname(mother,vmother);
 568
 569   gGeoManager->Division(vname, vmother, iaxis, ndiv, c0i, 0, numed, "nx");
 570 }
 571 //_____________________________________________________________________________
 572 void  TFlukaMCGeometry::Gsdvt(const char *name, const char *mother, Double_t step,
 573                      Int_t iaxis, Int_t numed, Int_t /*ndvmx*/)
 574 {
 575   //
 576   // Create a new volume by dividing an existing one
 577   //
 578   //       Divides MOTHER into divisions called NAME along
 579   //       axis IAXIS in steps of STEP. If not exactly divisible
 580   //       will make as many as possible and will centre them
 581   //       with respect to the mother. Divisions will have medium
 582   //       number NUMED. If NUMED is 0, NUMED of MOTHER is taken.
 583   //       NDVMX is the expected maximum number of divisions
 584   //          (If 0, no protection tests are performed)
 585   //
 586   char vname[5];
 587   Vname(name,vname);
 588   char vmother[5];
 589   Vname(mother,vmother);
 590
 591   gGeoManager->Division(vname, vmother, iaxis, 0, 0, step, numed, "s");
 592 }
 593
 594 //_____________________________________________________________________________
 595 void  TFlukaMCGeometry::Gsdvt2(const char *name, const char *mother, Double_t step,
 596                       Int_t iaxis, Double_t c0, Int_t numed, Int_t /*ndvmx*/)
 597 {
 598   //
 599   // Create a new volume by dividing an existing one
 600   //
 601   //           Divides MOTHER into divisions called NAME along
 602   //            axis IAXIS starting at coordinate value C0 with step
 603   //            size STEP.
 604   //           The new volume created will have medium number NUMED.
 605   //           If NUMED is 0, NUMED of mother is taken.
 606   //           NDVMX is the expected maximum number of divisions
 607   //             (If 0, no protection tests are performed)
 608   //
 609   char vname[5];
 610   Vname(name,vname);
 611   char vmother[5];
 612   Vname(mother,vmother);
 613
 614   gGeoManager->Division(vname, vmother, iaxis, 0, c0, step, numed, "sx");
 615 }
 616
 617 //_____________________________________________________________________________
 618 void  TFlukaMCGeometry::Gsord(const char * /*name*/, Int_t /*iax*/)
 619 {
 620   //
 621   //    Flags volume CHNAME whose contents will have to be ordered
 622   //    along axis IAX, by setting the search flag to -IAX
 623   //           IAX = 1    X axis
 624   //           IAX = 2    Y axis
 625   //           IAX = 3    Z axis
 626   //           IAX = 4    Rxy (static ordering only  -> GTMEDI)
 627   //           IAX = 14   Rxy (also dynamic ordering -> GTNEXT)
 628   //           IAX = 5    Rxyz (static ordering only -> GTMEDI)
 629   //           IAX = 15   Rxyz (also dynamic ordering -> GTNEXT)
 630   //           IAX = 6    PHI   (PHI=0 => X axis)
 631   //           IAX = 7    THETA (THETA=0 => Z axis)
 632   //
 633
 634   // TBC - keep this function
 635   // nothing to be done for TGeo  //xx
 636 }
 637
 638 //_____________________________________________________________________________
 639 void  TFlukaMCGeometry::Gspos(const char *name, Int_t nr, const char *mother, Double_t x,
 640                      Double_t y, Double_t z, Int_t irot, const char *konly)
 641 {
 642   //
 643   // Position a volume into an existing one
 644   //
 645   //  NAME   Volume name
 646   //  NUMBER Copy number of the volume
 647   //  MOTHER Mother volume name
 648   //  X      X coord. of the volume in mother ref. sys.
 649   //  Y      Y coord. of the volume in mother ref. sys.
 650   //  Z      Z coord. of the volume in mother ref. sys.
 651   //  IROT   Rotation matrix number w.r.t. mother ref. sys.
 652   //  ONLY   ONLY/MANY flag
 653   //
 654   //  It positions a previously defined volume in the mother.
 655   //
 656
 657   TString only = konly;
 658   only.ToLower();
 659   Bool_t isOnly = kFALSE;
 660   if (only.Contains("only")) isOnly = kTRUE;
 661   char vname[5];
 662   Vname(name,vname);
 663   char vmother[5];
 664   Vname(mother,vmother);
 665
 666   Double_t *upar=0;
 667   gGeoManager->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar);
 668   if (fDebug) printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
 669 }
 670
 671 //_____________________________________________________________________________
 672 void  TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother,
 673                       Double_t x, Double_t y, Double_t z, Int_t irot,
 674                       const char *konly, Float_t *upar, Int_t np )
 675 {
 676   //
 677   //      Place a copy of generic volume NAME with user number
 678   //      NR inside MOTHER, with its parameters UPAR(1..NP)
 679   //
 680
 681   Double_t* dupar = CreateDoubleArray(upar, np);
 682   Gsposp(name, nr, mother, x, y, z, irot, konly, dupar, np);
 683   delete [] dupar;
 684 }
 685
 686 //_____________________________________________________________________________
 687 void  TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother,
 688                       Double_t x, Double_t y, Double_t z, Int_t irot,
 689                       const char *konly, Double_t *upar, Int_t np )
 690 {
 691   //
 692   //      Place a copy of generic volume NAME with user number
 693   //      NR inside MOTHER, with its parameters UPAR(1..NP)
 694   //
 695
 696   TString only = konly;
 697   only.ToLower();
 698   Bool_t isOnly = kFALSE;
 699   if (only.Contains("only")) isOnly = kTRUE;
 700   char vname[5];
 701   Vname(name,vname);
 702   char vmother[5];
 703   Vname(mother,vmother);
 704
 705   gGeoManager->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np);
 706   if (fDebug) printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
 707 }
 708
 709 //_____________________________________________________________________________
 710 Int_t TFlukaMCGeometry::VolId(const Text_t *name) const
 711 {
 712   //
 713   // Return the unique numeric identifier for volume name
 714   //
 715
 716   Int_t uid = gGeoManager->GetUID(name);
 717   if (uid<0) {
 718      printf("VolId: Volume %s not found\n",name);
 719      return 0;
 720   }
 721   if (fDebug) printf("VolId for %s: %i\n", name, uid);
 722   return uid;
 723 }
 724
 725 //_____________________________________________________________________________
 726 const char* TFlukaMCGeometry::VolName(Int_t id) const
 727 {
 728   //
 729   // Return the volume name given the volume identifier
 730   //
 731   TGeoVolume *volume = gGeoManager->GetVolume(id);
 732   if (!volume) {
 733      Error("VolName","volume with id=%d does not exist",id);
 734      return "NULL";
 735   }
 736   if (fDebug) printf("VolName for id=%i: %s\n", id, volume->GetName());
 737   return volume->GetName();
 738 }
 739
 740 //_____________________________________________________________________________
 741 Int_t TFlukaMCGeometry::NofVolumes() const
 742 {
 743   //
 744   // Return total number of volumes in the geometry
 745   //
 746
 747   return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
 748 }
 749
 750 //_____________________________________________________________________________
 751 Int_t TFlukaMCGeometry::VolId2Mate(Int_t id) const
 752 {
 753   //
 754   // Return material number for a given volume id
 755   //
 756   TGeoVolume *volume = gGeoManager->GetVolume(id);
 757   if (!volume) {
 758      Error("VolId2Mate","volume with id=%d does not exist",id);
 759      return 0;
 760   }
 761   TGeoMedium *med = volume->GetMedium();
 762   if (!med) return 0;
 763   if (fDebug) printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId());
 764   return med->GetId();
 765 }
 766
 767 //_____________________________________________________________________________
 768 Int_t TFlukaMCGeometry::CurrentVolID(Int_t& copyNo) const
 769 {
 770   // Returns the current volume ID and copy number
 771   if (gGeoManager->IsOutside()) return 0;
 772   TGeoNode *node = gGeoManager->GetCurrentNode();
 773   copyNo = node->GetNumber();
 774   Int_t id = node->GetVolume()->GetNumber();
 775   if (fDebug) printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id);
 776   return id;
 777 }
 778
 779 //_____________________________________________________________________________
 780 Int_t TFlukaMCGeometry::CurrentVolOffID(Int_t off, Int_t& copyNo) const
 781 {
 782   // Return the current volume "off" upward in the geometrical tree
 783   // ID and copy number
 784   if (off<0 || off>gGeoManager->GetLevel()) return 0;
 785   if (off==0) return CurrentVolID(copyNo);
 786   TGeoNode *node = gGeoManager->GetMother(off);
 787   if (!node) return 0;
 788   copyNo = node->GetNumber();
 789   if (fDebug) printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() );
 790   return node->GetVolume()->GetNumber();
 791 }
 792 // FLUKA specific
 793
 794 //_____________________________________________________________________________
 795 const char* TFlukaMCGeometry::CurrentVolName() const
 796 {
 797   //
 798   // Returns the current volume name
 799   //
 800   if (gGeoManager->IsOutside()) return 0;
 801   if (fDebug) printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName());
 802   return gGeoManager->GetCurrentVolume()->GetName();
 803 }
 804 //_____________________________________________________________________________
 805 const char* TFlukaMCGeometry::CurrentVolOffName(Int_t off) const
 806 {
 807   //
 808   // Return the current volume "off" upward in the geometrical tree
 809   // ID, name and copy number
 810   // if name=0 no name is returned
 811   //
 812   if (off<0 || off>gGeoManager->GetLevel()) return 0;
 813   if (off==0) return CurrentVolName();
 814   TGeoNode *node = gGeoManager->GetMother(off);
 815   if (!node) return 0;
 816   if (fDebug) printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName());
 817   return node->GetVolume()->GetName();
 818 }
 819
 820 //_____________________________________________________________________________
 821 void TFlukaMCGeometry::Gsatt(const char *name, const char *att, Int_t val)
 822 {
 823   //
 824   //  NAME   Volume name
 825   //  IOPT   Name of the attribute to be set
 826   //  IVAL   Value to which the attribute is to be set
 827   // see: TFluka::Gsatt
 828   char vname[5];
 829   Vname(name,vname);
 830   char vatt[5];
 831   Vname(att,vatt);
 832   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 833 }
 834
 835 //_____________________________________________________________________________
 836 void TFlukaMCGeometry::Gdtom(Float_t *xd, Float_t *xm, Int_t iflag)
 837 {
 838   //
 839   //  Computes coordinates XM (Master Reference System
 840   //  knowing the coordinates XD (Detector Ref System)
 841   //  The local reference system can be initialized by
 842   //    - the tracking routines and GDTOM used in GUSTEP
 843   //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
 844   //        (inverse routine is GMTOD)
 845   //
 846   //   If IFLAG=1  convert coordinates
 847   //      IFLAG=2  convert direction cosinus
 848   //
 849    Double_t XM[3], XD[3];
 850    Int_t i;
 851    for (i=0;i<3;i++) XD[i] = xd[i];
 852    if (iflag == 1) gGeoManager->LocalToMaster(XD,XM);
 853    else            gGeoManager->LocalToMasterVect(XD,XM);
 854    for (i=0;i<3;i++) xm[i]=XM[i];
 855 }
 856
 857 //_____________________________________________________________________________
 858 void TFlukaMCGeometry::Gdtom(Double_t *xd, Double_t *xm, Int_t iflag)
 859 {
 860    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
 861    else            gGeoManager->LocalToMasterVect(xd,xm);
 862 }
 863
 864 //_____________________________________________________________________________
 865 void TFlukaMCGeometry::Gmtod(Float_t *xm, Float_t *xd, Int_t iflag)
 866 {
 867   //
 868   //       Computes coordinates XD (in DRS)
 869   //       from known coordinates XM in MRS
 870   //       The local reference system can be initialized by
 871   //         - the tracking routines and GMTOD used in GUSTEP
 872   //         - a call to GMEDIA(XM,NUMED,CHECK)
 873   //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
 874   //             (inverse routine is GDTOM)
 875   //
 876   //        If IFLAG=1  convert coordinates
 877   //           IFLAG=2  convert direction cosinus
 878   //
 879    Double_t XM[3], XD[3];
 880    Int_t i;
 881    for (i=0;i<3;i++) XM[i]=xm[i];
 882    if (iflag == 1) gGeoManager->MasterToLocal(XM,XD);
 883    else            gGeoManager->MasterToLocalVect(XM,XD);
 884    for (i=0;i<3;i++) xd[i] = XD[i];
 885 }
 886
 887 //_____________________________________________________________________________
 888 void TFlukaMCGeometry::Gmtod(Double_t *xm, Double_t *xd, Int_t iflag)
 889 {
 890    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
 891    else            gGeoManager->MasterToLocalVect(xm,xd);
 892 }
 893
 894 //_____________________________________________________________________________
 895 void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
 896 {
 897   // ==== from FLUGG ====
 898   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
 899   // according to the fluka standard. In addition,. they must be equal to the
 900   // names of the fluka materials - see fluka manual - in order that the
 901   // program load the right cross sections, and equal to the names included in
 902   // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
 903   // own .pemf, in order to get the right cross sections loaded in memory.
 904
 905    Int_t zelem[128];
 906    static char elNames[220] = {
 907    //  1 ============================= 5 ==================================== 10 ===================================== 15 ===
 908       'H','_','H','E','L','I','B','E','B','_','C','_','N','_','O','_','F','_','N','E','N','A','M','G','A','L','S','I','P','_',
 909       'S','_','C','L','A','R','K','_','C','A','S','C','T','I','V','_','C','R','M','N','F','E','C','O','N','I','C','U','Z','N',
 910       'G','A','G','E','A','S','S','E','B','R','K','R','R','B','S','R','Y','_','Z','R','N','B','M','O','T','C','R','U','R','H',
 911       'P','D','A','G','C','D','I','N','S','N','S','B','T','E','I','_','X','E','C','S','B','A','L','A','C','E','P','R','N','D',
 912       'P','M','S','M','E','U','G','D','T','B','D','Y','H','O','E','R','T','M','Y','B','L','U','H','F','T','A','W','_','R','E',
 913       'O','S','I','R','P','T','A','U','H','G','T','L','P','B','B','I','P','O','A','T','R','N','F','R','R','A','A','C','T','H',
 914       'P','A','U','_','N','P','P','U','A','M','C','M','B','K','C','F','E','S','F','M','M','D','N','O','L','R','R','F','D','B',
 915       'S','G','B','H','H','S','M','T','D','S'};
 916    memset(zelem, 0, 128*sizeof(Int_t));
 917    TString sname;
 918    gGeoManager->Export("flgeom.root");
 919    if (fname) sname = fname;
 920    else       sname = "flukaMat.inp";
 921    ofstream out;
 922    out.open(sname.Data(), ios::out);
 923    if (!out.good()) {
 924       Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
 925       return;
 926    }
 927    PrintHeader(out, "MATERIALS AND COMPOUNDS");
 928    PrintHeader(out, "MATERIALS");
 929    TList *matlist = gGeoManager->GetListOfMaterials();
 930    TIter next(matlist);
 931    Int_t nmater = matlist->GetSize();
 932    Int_t nfmater = 0;
 933    TObjArray *listfluka = new TObjArray(nmater+50);
 934    TObjArray *listflukanames = new TObjArray(nmater+50);
 935    TGeoMaterial *mat, *matorig;
 936    TGeoMixture *mix = 0;
 937    TString matname;
 938    TObjString *objstr, *objstrother;
 939    Int_t i,j,k,idmat;
 940    Bool_t done;
 941    Int_t nelem, nidmat;
 942    Double_t amat,zmat,rhomat;
 943    Double_t  zel, ael, wel, rho;
 944    char elname[8] = {' ',' ','_', 'E','L','E','M','\0'};
 945    char digit[3];
 946    Bool_t found = kFALSE;
 947
 948    if (fDebug) printf("Creating materials and compounds\n");
 949    for (i=0; i<nmater; i++) {
 950       mat = (TGeoMaterial*)matlist->At(i);
 951       if (mat->GetZ()<1E-1) {
 952          mat->SetIndex(2); // vacuum, built-in inside FLUKA
 953          continue;
 954       }
 955 //      printf("material: %s index=%i: Z=%f A=%f rho=%f\n", mat->GetName(), mat->GetIndex(),mat->GetZ(),mat->GetA(),mat->GetDensity());
 956       matorig = gGeoManager->FindDuplicateMaterial(mat);
 957       if (matorig) {
 958          idmat = matorig->GetIndex();
 959          mat->SetIndex(idmat);
 960 //         printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat);
 961          matorig = 0;
 962       } else  {
 963 //         printf(" Adding to temp list with index %i\n", nfmater+3);
 964          listfluka->Add(mat);
 965          mat->SetIndex(nfmater+3);
 966          matorig = mat;
 967          objstr = new TObjString(mat->GetName());
 968          listflukanames->Add(objstr);
 969          nfmater++;
 970          // look if name is existing
 971          nidmat = 0;
 972          matname = objstr->GetString();
 973          ToFlukaString(matname);
 974          objstr->SetString(matname.Data());
 975          done = kFALSE;
 976          while (!done) {
 977             if (nfmater == 1) break;
 978             for (j=0; j<nfmater-1; j++) {
 979                objstrother = (TObjString*)listflukanames->At(j);
 980                if (objstr->IsEqual(objstrother)) {
 981                   // we have to change the name
 982                   if (nidmat>98) {
 983                      Error("CreateFlukaMatFile", "too many materials having same name");
 984                      return;
 985                   }
 986                   nidmat++;
 987                   k = matname.Index(" ");
 988                   if (k<0 || k>6) k=6;
 989                   if (nidmat>9) {
 990                      sprintf(digit, "%d", nidmat);
 991                   } else {
 992                      digit[0] = '0';
 993                      sprintf(&digit[1], "%d", nidmat);
 994                   }
 995                   matname.Insert(k,digit);
 996                   matname.Remove(8);
 997                   objstr->SetString(matname.Data());
 998                   break;
 999                }
1000                if (j == nfmater-2) {
1001                   done = kTRUE;
1002                   break;
1003                }
1004             }
1005          }
1006 //         printf(" newmat name: %s\n", matname.Data());
1007       }
1008       // now we have unique materials with unique names in the lists
1009
1010       if (matorig && matorig->IsMixture()) {
1011       // create dummy materials for elements
1012          rho = 0.999;
1013          mix = (TGeoMixture*)matorig;
1014          nelem = mix->GetNelements();
1015 //         printf(" material is a MIXTURE with %i elements:\n", nelem);
1016          for (j=0; j<nelem; j++) {
1017             found = kFALSE;
1018             zel = (mix->GetZmixt())[j];
1019             ael = (mix->GetAmixt())[j];
1020 //            printf("   Zelem[%i] = %g\n",j,zel);
1021             if ((zel-Int_t(zel))>0.01) {
1022                TGeoMaterial *mat1;
1023                for (Int_t imat=0; imat<nfmater; imat++) {
1024                   mat1 = (TGeoMaterial*)listfluka->At(imat);
1025                   if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
1026                   if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
1027                   found = kTRUE;
1028                   break;
1029                }
1030                if (!found) Warning("CreateFlukaMatFile", "element with Z=%f\n", zel);
1031             }
1032             if (!zelem[Int_t(zel)] && !found) {
1033                // write fluka element
1034                memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
1035                zelem[Int_t(zel)] = 1;
1036                mat = new TGeoMaterial(elname, ael, zel, rho);
1037                mat->SetIndex(nfmater+3);
1038 //               printf("  element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n",
1039 //                       elname,nfmater+3,zel,ael,rho);
1040                listfluka->Add(mat);
1041                objstr = new TObjString(elname);
1042                listflukanames->Add(objstr);
1043                nfmater++;
1044             }
1045          }
1046       }
1047    }
1048    // now dump materials in the file
1049 //   printf("DUMPING %i materials\n", nfmater);
1050    for (i=0; i<nfmater; i++) {
1051       mat = (TGeoMaterial*)listfluka->At(i);
1052       out << setw(10) << "MATERIAL  ";
1053       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1054 //      matname = mat->GetName();
1055       objstr = (TObjString*)listflukanames->At(i);
1056       matname = objstr->GetString();
1057       ToFlukaString(matname);
1058       zmat = mat->GetZ();
1059       if (zmat-Int_t(zmat)>0.01) {
1060          if (zmat-Int_t(zmat)>0.5) zmat = Int_t(zmat)+1.;
1061          else zmat = Int_t(zmat);
1062       }
1063       amat = mat->GetA();
1064       rhomat = mat->GetDensity();
1065       // write material card
1066       if (mat->IsMixture()) {
1067          out << setw(10) << " ";
1068          out << setw(10) << " ";
1069          mix = (TGeoMixture*)mat;
1070       } else {
1071          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << zmat;
1072          out << setw(10) << setprecision(3) << amat;
1073       }
1074       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1075       out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rhomat;
1076       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1077       out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3);
1078       out << setw(10) << " ";
1079       out << setw(10) << " ";
1080       out << setw(8) << matname.Data() << endl;
1081    }
1082    // write mixture header
1083    PrintHeader(out, "COMPOUNDS");
1084    Int_t counttothree;
1085    TGeoMaterial *element;
1086    for (i=0; i<nfmater; i++) {
1087       mat = (TGeoMaterial*)listfluka->At(i);
1088       if (!mat->IsMixture()) continue;
1089       mix = (TGeoMixture*)mat;
1090       counttothree = 0;
1091       out << setw(10) << "COMPOUND  ";
1092       nelem = mix->GetNelements();
1093       objstr = (TObjString*)listflukanames->At(i);
1094       matname = objstr->GetString();
1095 //      printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem);
1096       for (j=0; j<nelem; j++) {
1097          // dump mixture cards
1098 //         printf(" #elem %i: Z=%g, A=%g, W=%g\n", j, (mix->GetZmixt())[j],
1099 //                (mix->GetAmixt())[j],(mix->GetWmixt())[j]);
1100          wel = (mix->GetWmixt())[j];
1101          zel = (mix->GetZmixt())[j];
1102          ael = (mix->GetAmixt())[j];
1103          if (zel-Int_t(zel)>0.01) {
1104             // loop the temporary list
1105             element = 0;
1106             TGeoMaterial *mat1;
1107             for (Int_t imat=0; imat<i; imat++) {
1108                mat1 = (TGeoMaterial*)listfluka->At(imat);
1109                if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
1110                if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
1111                element = mat1;
1112                break;
1113             }
1114          } else {
1115             memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
1116             element = (TGeoMaterial*)listfluka->FindObject(elname);
1117          }
1118          if (!element) {
1119             Error("CreateFlukaMatFile", "Element Z=%g %s not found", zel, elname);
1120             return;
1121          }
1122          idmat = element->GetIndex();
1123 //         printf("element %s , index=%i\n", element->GetName(), idmat);
1124          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1125          out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -wel;
1126          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1127          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
1128          counttothree++;
1129          if (counttothree == 3) {
1130             out << matname.Data();
1131             out << endl;
1132             if ( (j+1) != nelem) out << setw(10) << "COMPOUND  ";
1133             counttothree = 0;
1134          }
1135       }
1136       //Unless we have 3, 6, 9... submaterials we need to put some empty
1137       //space and the compound name
1138       if (nelem%3) {
1139          for (j=0; j<(3-(nelem%3)); j++)
1140             out << setw(10) << " " << setw(10) << " ";
1141          out << matname.Data();
1142          out << endl;
1143       }
1144    }
1145
1146    // Now print the list of regions (volumes in TGeo)
1147    Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
1148    TGeoVolume *vol;
1149 /*
1150    PrintHeader(out, "TGEO VOLUMES");
1151    for (i=1; i<=nvols; i++) {
1152       vol = gGeoManager->GetVolume(i);
1153       out.setf(std::ios::left, std::ios::adjustfield);
1154       out << setw(10) << i;
1155       out << setw(20) << vol->GetName() << endl;
1156    }
1157 */
1158    // Now print the material assignments
1159    Double_t flagfield;
1160    PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
1161    for (i=1; i<=nvols; i++) {
1162       vol = gGeoManager->GetVolume(i);
1163       mat = vol->GetMedium()->GetMaterial();
1164       idmat = mat->GetIndex();
1165 //      flagfield = (vol->GetField())?1.:0.;
1166       flagfield = 1.;
1167       out << setw(10) << "ASSIGNMAT ";
1168       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1169       out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
1170       out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
1171       out << setw(10) << "0.0";
1172       out << setw(10) << "0.0";
1173       out << setw(10) << setiosflags(ios::fixed) << flagfield;
1174       out << setw(10) << "0.0";
1175       out << endl;
1176    }
1177    delete listfluka;
1178    listflukanames->Delete();
1179    delete listflukanames;
1180    out.close();
1181    fLastMaterial = nfmater+2;
1182 }
1183
1184 //_____________________________________________________________________________
1185 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
1186 {
1187 // Print a FLUKA header.
1188   out << "*\n" << "*\n" << "*\n";
1189   out << "*********************  " << text << " *********************\n"
1190      << "*\n";
1191   out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
1192      << endl;
1193   out << "*" << endl;
1194 }
1195
1196 //_____________________________________________________________________________
1197 Int_t TFlukaMCGeometry::RegionId() const
1198 {
1199 // Returns current region id <-> TGeo node id
1200    if (gGeoManager->IsOutside()) return 0;
1201    return gGeoManager->GetCurrentNode()->GetUniqueID();
1202 }
1203 //_____________________________________________________________________________
1204 void TFlukaMCGeometry::SetMreg(Int_t mreg)
1205 {
1206 // Update if needed next history;
1207    Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1208    if (mreg==curreg) return;
1209    if (mreg==fNextRegion) {
1210       if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1);
1211       return;
1212    }
1213    if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg);
1214 }
1215
1216 //_____________________________________________________________________________
1217 void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
1218 {
1219 // Set index/history for next entered region
1220    fNextRegion = mreg;
1221    fNextLattice = latt;
1222 }
1223
1224 //_____________________________________________________________________________
1225 void TFlukaMCGeometry::ToFlukaString(TString &str) const
1226 {
1227 // ToFlukaString converts an string to something usefull in FLUKA:
1228 // * Capital letters
1229 // * Only 8 letters
1230 // * Replace ' ' by '_'
1231    if (str.Length()<8) {
1232       str += "        ";
1233    }
1234    str.Remove(8);
1235    Int_t ilast;
1236    for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1237    str.ToUpper();
1238    for (Int_t pos=0; pos<ilast; pos++)
1239       if (str(pos)==' ') str.Replace(pos,1,"_",1);
1240    return;
1241 }
1242 //______________________________________________________________________________
1243 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
1244 {
1245   //
1246   //  convert name to upper case. Make vname at least 4 chars
1247   //
1248   Int_t l = strlen(name);
1249   Int_t i;
1250   l = l < 4 ? l : 4;
1251   for (i=0;i<l;i++) vname[i] = toupper(name[i]);
1252   for (i=l;i<4;i++) vname[i] = ' ';
1253   vname[4] = 0;
1254 }
1255
1256
1257 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1258
1259 //_____________________________________________________________________________
1260 Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1261 {
1262 //   from FLUGG:
1263 // Wrapper for setting DNEAR option on fluka side. Must return 0
1264 // if user doesn't want Fluka to use DNEAR to compute the
1265 // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1266 // card in fluka input), returns 1 if user wants Fluka always to
1267 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1268 // coming from all directions!!!)
1269    if (mcgeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
1270    return 0;
1271 }
1272
1273 //_____________________________________________________________________________
1274 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1275           Double_t *pV,  Int_t &oldReg , const Int_t &oldLttc, Double_t & propStep,
1276           Int_t & /*nascFlag*/, Double_t &retStep, Int_t &newReg,
1277                Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
1278           Double_t *sLt, Int_t *jrLt)
1279 {
1280 //   from FLUGG:
1281 // Wrapper for geometry tracking: returns approved step of
1282 // particle and all variables that fluka G1 computes.
1283
1284    // Initialize current point/direction
1285    if (mcgeom->IsDebugging()) {
1286       printf("========== Inside G1WR\n");
1287       printf("   point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
1288    }
1289    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1290    gGeoManager->SetCurrentDirection(pV);
1291    if (mcgeom->IsDebugging()) printf("   oldReg=%i  oldLttc=%i  pstep=%f\n",oldReg, oldLttc, propStep);
1292    if (oldLttc==999999999) printf("WOOPS - wrong old lattice\n");
1293    if (gGeoManager->IsOutside()) {
1294       gGeoManager->SetOutside(kFALSE);
1295       gGeoManager->CdTop();
1296    }
1297    Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
1298    Int_t curreg = gGeoManager->GetCurrentVolume()->GetNumber();
1299    if (mcgeom->IsDebugging()) printf("   curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1300    Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
1301    if (!regsame && mcgeom->IsDebugging()) printf("   REGIONS DOES NOT MATCH\n");
1302    if (oldLttc != curLttc) {
1303       if (mcgeom->IsDebugging()) printf("   HISTORIES DOES NOT MATCH\n");
1304       gGeoManager->CdNode(oldLttc-1);
1305       curLttc = gGeoManager->GetCurrentNodeId()+1;
1306       curreg  = gGeoManager->GetCurrentVolume()->GetNumber();
1307       if (mcgeom->IsDebugging()) printf("   re-initialized point: curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1308    }
1309    lttcFlag = 0;
1310    sLt[lttcFlag] = 0.;
1311    jrLt[lttcFlag] = curLttc;
1312    // now 'oldregion' contains the real region, matching or not the old history
1313
1314    // Compute geometry step/safety within physical step limit
1315 //   newReg = oldregion;
1316    Double_t *point = gGeoManager->GetCurrentPoint();
1317    Double_t *dir = gGeoManager->GetCurrentDirection();
1318    Double_t steptot = 0.;
1319    Double_t snext = 0.;
1320    Int_t istep = 0;
1321    Bool_t done = kFALSE;
1322    Double_t pst;
1323    Int_t i;
1324    while (!done) {
1325       gGeoManager->FindNextBoundary(-propStep);
1326       snext = gGeoManager->GetStep();
1327       if (mcgeom->IsDebugging()) printf("   FindNextBoundary(%g) snext=%g\n", propStep, snext);
1328       if (steptot == 0) {
1329          saf = gGeoManager->GetSafeDistance();
1330          if (mcgeom->IsDebugging()) printf("   Safety: %g\n", saf);
1331       }
1332       sLt[lttcFlag] = propStep;
1333       jrLt[lttcFlag] = gGeoManager->GetCurrentNodeId()+1;
1334       lttcFlag++; //1
1335       sLt[lttcFlag] = 0.;
1336       jrLt[lttcFlag] = -1;
1337       newReg = curreg;
1338       newLttc = oldLttc;
1339       if (snext<propStep) {
1340          // There is a boundary on the way.
1341          // Make a step=snext+1E-6 to force boundary crossing
1342          lttcFlag--; // 0
1343          steptot += snext;
1344          sLt[lttcFlag] = snext;
1345          retStep = snext;
1346 //         lttcFlag++;
1347          // make the step to get into the next region
1348          for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
1349          gGeoManager->FindNode();
1350          istep = 0;
1351          if (mcgeom->IsDebugging()) printf("   boundary: step made %g\n", snext);
1352          while (gGeoManager->IsSameLocation() && steptot<propStep) {
1353             if (istep>1E3) {
1354                printf("Geometry error: could not cross boundary after extra 10 microns\n");
1355                return;
1356             }
1357             for (i=0;i<3;i++) point[i]+=1E-6*dir[i];
1358             gGeoManager->FindNode();
1359             sLt[lttcFlag] += 1E-6;
1360             retStep = sLt[lttcFlag];
1361             steptot += 1E-6;
1362             istep++;
1363          }
1364          if (steptot>propStep) {printf("Error\n");return;}
1365          // we managed to cross the boundary -> in which region
1366          newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1367          lttcFlag++; //1
1368          newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
1369          sLt[lttcFlag] = snext; // at 1
1370          jrLt[lttcFlag] = newLttc;
1371          sLt[lttcFlag+1] = 0.;
1372          jrLt[lttcFlag+1] = -1;
1373          // !!!!!!!!!!
1374
1375          while (newReg==oldReg && steptot<propStep) {
1376             if (mcgeom->IsDebugging()) printf("   Entered SAME region... continue\n");
1377             pst = propStep-steptot;
1378             gGeoManager->FindNextBoundary(-pst);
1379             snext = gGeoManager->GetStep();
1380             steptot += snext;
1381             if (snext<pst) {
1382                printf("Found new boundary\n");
1383                sLt[lttcFlag] = snext;
1384                retStep = steptot; // ???
1385                for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
1386                steptot += 1E-6;
1387                gGeoManager->FindNode();
1388                if (gGeoManager->IsSameLocation()) {
1389                   printf("Cannot cross boundary\n");
1390                   break;
1391                }
1392                newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1393                newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
1394                if (mcgeom->IsDebugging()) printf("Found newreg=%i, newLttc=%i, lttFlag is: %i\n", newReg, newLttc, lttcFlag);
1395                sLt[lttcFlag-1] += snext; // correct step in old region
1396                sLt[lttcFlag] = propStep-snext;
1397                jrLt[lttcFlag] = newLttc;
1398                sLt[lttcFlag+1] = 0.;
1399                jrLt[lttcFlag+1] = -1;
1400                if (newReg != oldReg) break; // lttcFlag=1
1401                lttcFlag++;
1402             } else {
1403                 if (mcgeom->IsDebugging()) printf("Not crossing next\n");
1404                 lttcFlag--; //0
1405                 retStep=steptot;
1406                 sLt[lttcFlag] = retStep;
1407                 sLt[lttcFlag+1] = 0.;
1408                 jrLt[lttcFlag+1] = -1;
1409                 done = kTRUE;
1410             }
1411          }
1412
1413          lttcFlag++; //2
1414          if (mcgeom->IsDebugging()) {
1415             if (!gGeoManager->IsOutside()) {
1416                printf("   ENTERED region %i, newLttc=%i in: %s\n", newReg,newLttc,gGeoManager->GetPath());
1417             } else printf("   EXIT GEOMETRY: BLKHOLE reg=%i\n", newReg);
1418          }
1419       }
1420       // no boundary within proposed step
1421       lttcFlag--;
1422       done = kTRUE;
1423    }
1424    if (mcgeom->IsDebugging()) printf("=> newReg=%i newLttc=%i lttcFlag=%i\n", newReg, newLttc, lttcFlag);
1425    mcgeom->SetNextRegion(newReg, newLttc);
1426    if (mcgeom->IsDebugging()) {
1427       printf("=> snext=%g safe=%g\n", snext, saf);
1428       for (Int_t i=0; i<lttcFlag+1; i++) printf("   jrLt[%i]=%i  sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1429    }
1430    if (newLttc!=oldLttc) {
1431       if (gGeoManager->IsOutside()) {
1432          gGeoManager->SetOutside(kFALSE);
1433          gGeoManager->CdTop();
1434       }
1435       gGeoManager->CdNode(oldLttc-1);
1436    }
1437    if (mcgeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
1438 }
1439
1440 //_____________________________________________________________________________
1441 void g1rtwr()
1442 {
1443    if (mcgeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
1444 }
1445
1446 //_____________________________________________________________________________
1447 void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
1448 {
1449    if (mcgeom->IsDebugging()) printf("========== Dummy CONHWR\n");
1450 }
1451
1452 //_____________________________________________________________________________
1453 void inihwr(Int_t &intHist)
1454 {
1455    if (mcgeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
1456    if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1457    if (intHist<=0) {
1458 //      printf("=== wrong history number\n");
1459       return;
1460    }
1461    if (intHist==0) gGeoManager->CdTop();
1462    else gGeoManager->CdNode(intHist-1);
1463    if (mcgeom->IsDebugging()) {
1464       printf(" --- current path: %s\n", gGeoManager->GetPath());
1465       printf("<= INIHWR\n");
1466    }
1467 }
1468
1469 //_____________________________________________________________________________
1470 void  jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1471              Int_t &flukaReg)
1472 {
1473 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1474 // number of regions (volumes in TGeo)
1475    // build application geometry
1476    if (mcgeom->IsDebugging()) printf("========== Inside JOMIWR\n");
1477    flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
1478    if (mcgeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
1479 }
1480
1481 //_____________________________________________________________________________
1482 void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1483             Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1484             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1485 {
1486    if (mcgeom->IsDebugging()) {
1487       printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
1488 //      printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1489       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1490    }
1491    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1492    if (gGeoManager->IsOutside()) {
1493       newReg = mcgeom->NofVolumes()+1;
1494 //      newLttc = gGeoManager->GetCurrentNodeId();
1495       newLttc = 999999999;
1496       if (mcgeom->IsDebugging()) {
1497          printf("OUTSIDE\n");
1498          printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1499          printf("<= LKMGWR\n");
1500       }
1501       flagErr = newReg;
1502       return;
1503    }
1504    newReg = node->GetVolume()->GetNumber();
1505    newLttc = gGeoManager->GetCurrentNodeId()+1;
1506    flagErr = newReg;
1507    if (mcgeom->IsDebugging()) {
1508       printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1509       printf("<= LKDBWR\n");
1510    }
1511 }
1512
1513 //_____________________________________________________________________________
1514 void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1515             Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1516             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1517 {
1518    if (mcgeom->IsDebugging()) {
1519       printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
1520 //      printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1521       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1522    }
1523    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1524    if (gGeoManager->IsOutside()) {
1525       newReg = mcgeom->NofVolumes()+1;
1526 //      newLttc = gGeoManager->GetCurrentNodeId();
1527       newLttc = 999999999;
1528       if (mcgeom->IsDebugging()) {
1529          printf("OUTSIDE\n");
1530          printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1531          printf("<= LKMGWR\n");
1532       }
1533       flagErr = newReg;
1534       return;
1535    }
1536    newReg = node->GetVolume()->GetNumber();
1537    newLttc = gGeoManager->GetCurrentNodeId()+1;
1538    flagErr = newReg;
1539    if (mcgeom->IsDebugging()) {
1540       printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1541       printf("<= LKFXWR\n");
1542    }
1543 }
1544
1545 //_____________________________________________________________________________
1546 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1547             Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1548                       Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1549 {
1550    if (mcgeom->IsDebugging()) {
1551       printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
1552 //      printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1553       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1554    }
1555    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1556    if (gGeoManager->IsOutside()) {
1557       newReg = mcgeom->NofVolumes()+1;
1558 //      newLttc = gGeoManager->GetCurrentNodeId();
1559       newLttc = 999999999;
1560       if (mcgeom->IsDebugging()) {
1561          printf("OUTSIDE\n");
1562          printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1563          printf("<= LKMGWR\n");
1564       }
1565       flagErr = newReg;
1566       return;
1567    }
1568    newReg = node->GetVolume()->GetNumber();
1569    newLttc = gGeoManager->GetCurrentNodeId()+1;
1570    flagErr = newReg;
1571    if (mcgeom->IsDebugging()) {
1572       printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1573       printf("<= LKMGWR\n");
1574    }
1575 }
1576
1577 //_____________________________________________________________________________
1578 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1579           Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1580                Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1581 {
1582    if (mcgeom->IsDebugging()) {
1583       printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
1584 //      printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1585       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1586    }
1587    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1588    if (gGeoManager->IsOutside()) {
1589       newReg = mcgeom->NofVolumes()+1;
1590 //      newLttc = gGeoManager->GetCurrentNodeId();
1591       newLttc = 999999999;
1592       if (mcgeom->IsDebugging()) {
1593          printf("OUTSIDE\n");
1594          printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1595          printf("<= LKMGWR\n");
1596       }
1597       flagErr = newReg;
1598       return;
1599    }
1600    newReg = node->GetVolume()->GetNumber();
1601    newLttc = gGeoManager->GetCurrentNodeId()+1;
1602    flagErr = newReg;
1603    if (mcgeom->IsDebugging()) {
1604       printf("  out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
1605       printf("<= LKWR\n");
1606    }
1607 }
1608
1609 //_____________________________________________________________________________
1610 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1611             Double_t &pVx, Double_t &pVy, Double_t &pVz,
1612                  Double_t *norml, const Int_t &oldReg,
1613                  const Int_t &newReg, Int_t &flagErr)
1614 {
1615    if (mcgeom->IsDebugging()) {
1616       printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1617       printf("   oldReg=%i, newReg=%i\n", oldReg,newReg);
1618    }
1619    Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1620    Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
1621    if (mcgeom->IsDebugging()) printf("   curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath());
1622    Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
1623    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1624    gGeoManager->SetCurrentDirection(pVx,pVy,pVz);
1625    if (!regsame) {
1626       if (mcgeom->IsDebugging()) printf("   REGIONS DOEN NOT MATCH\n");
1627       gGeoManager->FindNode();
1628       curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1629       curLttc = gGeoManager->GetCurrentNodeId()+1;
1630       if (mcgeom->IsDebugging()) printf("   re-initialized point: curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1631    }
1632    Double_t *dnorm = gGeoManager->FindNormalFast();
1633    flagErr = 0;
1634    if (!dnorm) {
1635       printf("   ERROR: Cannot compute fast normal\n");
1636       flagErr = 1;
1637       norml[0] = -pVx;
1638       norml[1] = -pVy;
1639       norml[2] = -pVz;
1640    }
1641    norml[0] = -dnorm[0];
1642    norml[1] = -dnorm[1];
1643    norml[2] = -dnorm[2];
1644    if (mcgeom->IsDebugging()) printf("   normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
1645    curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1646    curLttc = gGeoManager->GetCurrentNodeId()+1;
1647    if (mcgeom->IsDebugging()) {
1648       printf("   final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
1649       printf("<= NRMLWR\n");
1650    }
1651 }
1652
1653 //_____________________________________________________________________________
1654 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1655             Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1656 {
1657    if (mcgeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
1658 }
1659
1660 //_____________________________________________________________________________
1661 Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1662 {
1663 //   from FLUGG:
1664 // Wrapper for saving current navigation history (fCheck=default)
1665 // and returning its pointer. If fCheck=-1 copy of history pointed
1666 // by intHist is made in NavHistWithCount object, and its pointer
1667 // is returned. fCheck=1 and fCheck=2 cases are only in debugging
1668 // version: an array is created by means of FGeometryInit functions
1669 // (but could be a static int * ptrArray = new int[10000] with
1670 // file scope as well) that stores a flag for deleted/undeleted
1671 // histories and at the end of event is checked to verify that
1672 // all saved history objects have been deleted.
1673
1674 // For TGeo, just return the current node ID. No copy need to be made.
1675
1676    if (mcgeom->IsDebugging()) printf("=> Inside ISVHWR\n");
1677    if (check<0) return intHist;
1678    Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1679    if (mcgeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
1680    return histInt;
1681 }
1682
1683
1684
1685