TFluka/TFlukaMCGeometry.cxx

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