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