TFluka/TFlukaMCGeometry.cxx

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