TFluka/TFlukaMCGeometry.cxx

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