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