TFluka/TFlukaMCGeometry.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // $Id$
  17
  18 // Class TFlukaMCGeometry
  19 // --------------------
  20 // Implementation of the TVirtualMCGeometry interface
  21 // for defining and using TGeo geometry with FLUKA.
  22 // This allows the FLUKA MonteCarlo to run with the TGeo
  23 // geometrical modeller
  24 // Author: Andrei Gheata 10/07/2003
  25
  26 #include "Riostream.h"
  27 #include "TList.h"
  28 #include "TCallf77.h"
  29 #include "TFluka.h"
  30 #include "TSystem.h"
  31 #include "TFlukaMCGeometry.h"
  32 #include "TFlukaConfigOption.h"
  33 #include "TGeoManager.h"
  34 #include "TGeoVolume.h"
  35 #include "TObjString.h"
  36 #include "Fsourcm.h"
  37 #include "Ftrackr.h"
  38 #include "Fstepsz.h"       //(STEPSZ) fluka common
  39
  40 #ifndef WIN32
  41 # define idnrwr idnrwr_
  42 # define g1wr   g1wr_
  43 # define g1rtwr g1rtwr_
  44 # define conhwr conhwr_
  45 # define inihwr inihwr_
  46 # define jomiwr jomiwr_
  47 # define lkdbwr lkdbwr_
  48 # define lkfxwr lkfxwr_
  49 # define lkmgwr lkmgwr_
  50 # define lkwr lkwr_
  51 # define magfld magfld_
  52 # define nrmlwr nrmlwr_
  53 # define rgrpwr rgrpwr_
  54 # define isvhwr isvhwr_
  55
  56 #else
  57
  58 # define idnrwr IDNRWR
  59 # define g1wr   G1WR
  60 # define g1rtwr G1RTWR
  61 # define conhwr CONHWR
  62 # define inihwr INIHWR
  63 # define jomiwr JOMIWR
  64 # define lkdbwr LKDBWR
  65 # define lkfxwr LKFXWR
  66 # define lkmgwr LKMGWR
  67 # define lkwr   LKWR
  68 # define magfld MAGFLD
  69 # define nrmlwr NRMLWR
  70 # define rgrpwr RGRPWR
  71 # define isvhwr ISVHWR
  72
  73 #endif
  74
  75 //____________________________________________________________________________
  76 extern "C"
  77 {
  78    //
  79    // Prototypes for FLUKA navigation methods
  80    //
  81    Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
  82    void  type_of_call   g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  83                              Double_t * /*pV*/,  Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
  84                              Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
  85                                   Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
  86                              Double_t *s /*Lt*/, Int_t * /*jrLt*/);
  87
  88    void  type_of_call g1rtwr();
  89    void  type_of_call conhwr(Int_t & /*intHist*/, Int_t & /*incrCount*/);
  90    void  type_of_call inihwr(Int_t & /*intHist*/);
  91    void  type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
  92                              Int_t & /*flukaReg*/);
  93    void  type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  94                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  95                              Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
  96    void  type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
  97                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
  98                              Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
  99    void  type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
 100                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
 101                                        Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
 102    void  type_of_call   lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
 103                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
 104                                   Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
 105 //   void  type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
 106 //                             Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
 107 //                             Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
 108    void  type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
 109                              Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
 110                                   Double_t * /*norml*/, const Int_t & /*oldReg*/,
 111                                   const Int_t & /*newReg*/, Int_t & /*flagErr*/);
 112    void  type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
 113                              Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
 114    Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
 115 }
 116
 117 // TFluka global pointer
 118 TFluka *gFluka = 0;
 119 TFlukaMCGeometry *gMCGeom = 0;
 120 Int_t gNstep = 0;
 121
 122 ClassImp(TFlukaMCGeometry)
 123
 124 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
 125
 126 //_____________________________________________________________________________
 127 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
 128    :TNamed(name, title),
 129    fDebug(kFALSE),
 130    fLastMaterial(0),
 131    fDummyRegion(0),
 132    fCurrentRegion(0),
 133    fCurrentLattice(0),
 134    fNextRegion(0),
 135    fNextLattice(0),
 136    fRegionList(0),
 137    fIndmat(0),
 138    fMatList(new TObjArray(256)),
 139    fMatNames(new TObjArray(256))
 140 {
 141   //
 142   // Standard constructor
 143   //
 144   gFluka = (TFluka*)gMC;
 145   gMCGeom = this;
 146   gNstep = 0;
 147 }
 148
 149 //_____________________________________________________________________________
 150 TFlukaMCGeometry::TFlukaMCGeometry()
 151   :TNamed(),
 152    fDebug(kFALSE),
 153    fLastMaterial(0),
 154    fDummyRegion(0),
 155    fCurrentRegion(0),
 156    fCurrentLattice(0),
 157    fNextRegion(0),
 158    fNextLattice(0),
 159    fRegionList(0),
 160    fIndmat(0),
 161    fMatList(0),
 162    fMatNames(0)
 163
 164 {
 165   //
 166   // Default constructor
 167   //
 168   gFluka = (TFluka*)gMC;
 169   gMCGeom = this;
 170   gNstep = 0;
 171 }
 172
 173 //_____________________________________________________________________________
 174 TFlukaMCGeometry::~TFlukaMCGeometry()
 175 {
 176   //
 177   // Destructor
 178   //
 179   fgInstance=0;
 180   if (fRegionList) delete [] fRegionList;
 181   if (fMatList) delete fMatList;
 182   if (fMatNames) {fMatNames->Delete(); delete fMatNames;}
 183   if (gGeoManager) delete gGeoManager;
 184 }
 185
 186 //
 187 // private methods
 188 //
 189
 190 //_____________________________________________________________________________
 191 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
 192 {
 193 // Converts Float_t* array to Double_t*,
 194 // !! The new array has to be deleted by user.
 195 // ---
 196
 197   Double_t* doubleArray;
 198   if (size>0) {
 199     doubleArray = new Double_t[size];
 200     for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
 201   }
 202   else {
 203     //doubleArray = 0;
 204     doubleArray = new Double_t[1];
 205   }
 206   return doubleArray;
 207 }
 208 //
 209 // public methods
 210
 211
 212 //_____________________________________________________________________________
 213 Int_t TFlukaMCGeometry::GetMedium() const
 214 {
 215 // Get current medium number
 216    Int_t imed = 0;
 217    TGeoNode *node = gGeoManager->GetCurrentNode();
 218    if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
 219    else       imed = node->GetVolume()->GetMedium()->GetId();
 220    if (fDebug) printf("GetMedium=%i\n", imed);
 221    return imed;
 222 }
 223
 224 //_____________________________________________________________________________
 225 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
 226 {
 227 // Returns FLUKA material index for medium IMED
 228      TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
 229    if (!med) {
 230       Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
 231       return -1;
 232    }
 233    TGeoMaterial* mat = med->GetMaterial();
 234    if (!mat->IsUsed()) return -1;
 235    Int_t imatfl = med->GetMaterial()->GetIndex();
 236    return imatfl;
 237 }
 238
 239 //_____________________________________________________________________________
 240 Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
 241 {
 242 // Get an ordered list of regions matching a given medium number
 243    nreg = 0;
 244    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
 245    TIter next(gGeoManager->GetListOfUVolumes());
 246    TGeoVolume *vol;
 247    Int_t imedium, ireg;
 248    while ((vol = (TGeoVolume*)next())) {
 249        TGeoMedium* med;
 250        if ((med = vol->GetMedium()) == 0) continue;
 251        imedium = med->GetId();
 252        if (imedium == imed) {
 253            ireg = vol->GetNumber();
 254            fRegionList[nreg++] = ireg;
 255        }
 256    }
 257    return fRegionList;
 258 }
 259
 260 //_____________________________________________________________________________
 261 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
 262 {
 263 // Get an ordered list of regions matching a given medium number
 264    nreg = 0;
 265    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
 266    TIter next(gGeoManager->GetListOfUVolumes());
 267    TGeoVolume *vol;
 268    Int_t imaterial, ireg;
 269    while ((vol = (TGeoVolume*)next())) {
 270        TGeoMedium* med;
 271        if ((med = vol->GetMedium()) == 0) continue;
 272        imaterial = med->GetMaterial()->GetIndex();
 273        if (imaterial == imat) {
 274            ireg = vol->GetNumber();
 275            fRegionList[nreg++] = ireg;
 276        }
 277    }
 278    return fRegionList;
 279 }
 280
 281 //_____________________________________________________________________________
 282 Int_t TFlukaMCGeometry::NofVolumes() const
 283 {
 284   //
 285   // Return total number of volumes in the geometry
 286   //
 287
 288   return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
 289 }
 290
 291 //_____________________________________________________________________________
 292 TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
 293 {
 294 // Try to replace a wrongly-defined material
 295    static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177,
 296       25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193,
 297       15.3383, 4.55,   9.6502, 6.4561, 21.7963, 29.8246, 15.4021};
 298
 299    Int_t ind;
 300    Double_t dz;
 301    for (ind=0; ind<23; ind++) {
 302       dz = TMath::Abs(z-kz[ind]);
 303       if (dz<1E-4) break;
 304    }
 305    if (ind>22) {
 306       printf("Cannot patch material with Z=%g\n", z);
 307       return 0;
 308    }
 309    TGeoMixture *mix = 0;
 310    TGeoElement *element;
 311    TGeoElementTable *table = gGeoManager->GetElementTable();
 312    switch (ind) {
 313       case 0: // AIR
 314          mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
 315          element = table->GetElement(6); // C
 316          mix->DefineElement(0, element, 0.000124);
 317          element = table->GetElement(7); // N
 318          mix->DefineElement(1, element, 0.755267);
 319          element = table->GetElement(8); // O
 320          mix->DefineElement(2, element, 0.231781);
 321          element = table->GetElement(18); // AR
 322          mix->DefineElement(3, element, 0.012827);
 323          break;
 324       case 1: //SDD SI CHIP
 325          mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
 326          element = table->GetElement(1);
 327          mix->DefineElement(0, element, 0.004367771);
 328          element = table->GetElement(6);
 329          mix->DefineElement(1, element, 0.039730642);
 330          element = table->GetElement(7);
 331          mix->DefineElement(2, element, 0.001396798);
 332          element = table->GetElement(8);
 333          mix->DefineElement(3, element, 0.01169634);
 334          element = table->GetElement(14);
 335          mix->DefineElement(4, element, 0.844665);
 336          element = table->GetElement(47);
 337          mix->DefineElement(5, element, 0.09814344903);
 338          break;
 339       case 2:  // WATER
 340          mix = new TGeoMixture("ITS_WATER", 2, 1.0);
 341          element = table->GetElement(1);
 342          mix->DefineElement(0, element, 0.111898344);
 343          element = table->GetElement(8);
 344          mix->DefineElement(1, element, 0.888101656);
 345          break;
 346       case 3: // CERAMICS
 347          mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
 348          element = table->GetElement(8);
 349          mix->DefineElement(0, element, 0.59956);
 350          element = table->GetElement(13);
 351          mix->DefineElement(1, element, 0.3776);
 352          element = table->GetElement(14);
 353          mix->DefineElement(2, element, 0.00933);
 354          element = table->GetElement(24);
 355          mix->DefineElement(3, element, 0.002);
 356          element = table->GetElement(25);
 357          mix->DefineElement(4, element, 0.0115);
 358          break;
 359       case 4: // EPOXY
 360          mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
 361          element = table->GetElement(1);
 362          mix->DefineElement(0, element, 0.19);
 363          element = table->GetElement(6);
 364          mix->DefineElement(1, element, 0.18);
 365          element = table->GetElement(8);
 366          mix->DefineElement(2, element, 0.35);
 367          element = table->GetElement(14);
 368          mix->DefineElement(3, element, 0.28);
 369          break;
 370       case 5: // EPOXY
 371          mix = new TGeoMixture("G10FR4", 4, 1.8);
 372          element = table->GetElement(1);
 373          mix->DefineElement(0, element, 0.19);
 374          element = table->GetElement(6);
 375          mix->DefineElement(1, element, 0.18);
 376          element = table->GetElement(8);
 377          mix->DefineElement(2, element, 0.35);
 378          element = table->GetElement(14);
 379          mix->DefineElement(3, element, 0.28);
 380          break;
 381       case 6: // KAPTON
 382          mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
 383          element = table->GetElement(1);
 384          mix->DefineElement(0, element, 0.026363415);
 385          element = table->GetElement(6);
 386          mix->DefineElement(1, element, 0.6911272);
 387          element = table->GetElement(7);
 388          mix->DefineElement(2, element, 0.073271325);
 389          element = table->GetElement(8);
 390          mix->DefineElement(3, element, 0.209238060);
 391          break;
 392       case 7: // INOX
 393          mix = new TGeoMixture("ITS_INOX", 9, 7.9);
 394          element = table->GetElement(6);
 395          mix->DefineElement(0, element, 0.0003);
 396          element = table->GetElement(14);
 397          mix->DefineElement(1, element, 0.01);
 398          element = table->GetElement(15);
 399          mix->DefineElement(2, element, 0.00045);
 400          element = table->GetElement(16);
 401          mix->DefineElement(3, element, 0.0003);
 402          element = table->GetElement(24);
 403          mix->DefineElement(4, element, 0.17);
 404          element = table->GetElement(25);
 405          mix->DefineElement(5, element, 0.02);
 406          element = table->GetElement(26);
 407          mix->DefineElement(6, element, 0.654);
 408          element = table->GetElement(28);
 409          mix->DefineElement(7, element, 0.12);
 410          element = table->GetElement(42);
 411          mix->DefineElement(8, element, 0.025);
 412          break;
 413       case 8: // ROHACELL
 414          mix = new TGeoMixture("ROHACELL", 4, 0.05);
 415          element = table->GetElement(1);
 416          mix->DefineElement(0, element, 0.07836617);
 417          element = table->GetElement(6);
 418          mix->DefineElement(1, element, 0.64648941);
 419          element = table->GetElement(7);
 420          mix->DefineElement(2, element, 0.08376983);
 421          element = table->GetElement(8);
 422          mix->DefineElement(3, element, 0.19137459);
 423          break;
 424       case 9: // SDD-C-AL
 425          mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
 426          element = table->GetElement(1);
 427          mix->DefineElement(0, element, 0.022632);
 428          element = table->GetElement(6);
 429          mix->DefineElement(1, element, 0.8176579);
 430          element = table->GetElement(7);
 431          mix->DefineElement(2, element, 0.0093488);
 432          element = table->GetElement(8);
 433          mix->DefineElement(3, element, 0.0503618);
 434          element = table->GetElement(13);
 435          mix->DefineElement(4, element, 0.1);
 436          break;
 437       case 10: // X7R-CAP
 438          mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
 439          element = table->GetElement(8);
 440          mix->DefineElement(0, element, 0.085975822);
 441          element = table->GetElement(22);
 442          mix->DefineElement(1, element, 0.084755042);
 443          element = table->GetElement(28);
 444          mix->DefineElement(2, element, 0.038244751);
 445          element = table->GetElement(29);
 446          mix->DefineElement(3, element, 0.009471271);
 447          element = table->GetElement(50);
 448          mix->DefineElement(4, element, 0.321736471);
 449          element = table->GetElement(56);
 450          mix->DefineElement(5, element, 0.251639432);
 451          element = table->GetElement(82);
 452          mix->DefineElement(6, element, 0.2081768);
 453          break;
 454       case 11: // SDD ruby sph. Al2O3
 455          mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
 456          element = table->GetElement(8);
 457          mix->DefineElement(0, element, 0.5293);
 458          element = table->GetElement(13);
 459          mix->DefineElement(1, element, 0.4707);
 460          break;
 461       case 12: // SDD HV microcable
 462          mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
 463          element = table->GetElement(1);
 464          mix->DefineElement(0, element, 0.01983871336);
 465          element = table->GetElement(6);
 466          mix->DefineElement(1, element, 0.520088819984);
 467          element = table->GetElement(7);
 468          mix->DefineElement(2, element, 0.0551367996);
 469          element = table->GetElement(8);
 470          mix->DefineElement(3, element, 0.157399667056);
 471          element = table->GetElement(13);
 472          mix->DefineElement(4, element, 0.247536);
 473          break;
 474       case 13: //SDD LV+signal cable
 475          mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
 476          element = table->GetElement(1);
 477          mix->DefineElement(0, element, 0.0082859922);
 478          element = table->GetElement(6);
 479          mix->DefineElement(1, element, 0.21722436468);
 480          element = table->GetElement(7);
 481          mix->DefineElement(2, element, 0.023028867);
 482          element = table->GetElement(8);
 483          mix->DefineElement(3, element, 0.06574077612);
 484          element = table->GetElement(13);
 485          mix->DefineElement(4, element, 0.68572);
 486          break;
 487       case 14: //SDD hybrid microcab
 488          mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
 489          element = table->GetElement(1);
 490          mix->DefineElement(0, element, 0.00926228815);
 491          element = table->GetElement(6);
 492          mix->DefineElement(1, element, 0.24281879711);
 493          element = table->GetElement(7);
 494          mix->DefineElement(2, element, 0.02574224025);
 495          element = table->GetElement(8);
 496          mix->DefineElement(3, element, 0.07348667449);
 497          element = table->GetElement(13);
 498          mix->DefineElement(4, element, 0.64869);
 499          break;
 500       case 15: //SDD anode microcab
 501          mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
 502          element = table->GetElement(1);
 503          mix->DefineElement(0, element, 0.0128595919215);
 504          element = table->GetElement(6);
 505          mix->DefineElement(1, element, 0.392653705471);
 506          element = table->GetElement(7);
 507          mix->DefineElement(2, element, 0.041626868025);
 508          element = table->GetElement(8);
 509          mix->DefineElement(3, element, 0.118832707289);
 510          element = table->GetElement(13);
 511          mix->DefineElement(4, element, 0.431909);
 512          break;
 513       case 16: // inox/alum
 514          mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
 515          element = table->GetElement(13);
 516          mix->DefineElement(0, element, 0.816164);
 517          element = table->GetElement(14);
 518          mix->DefineElement(1, element, 0.000919182);
 519          element = table->GetElement(24);
 520          mix->DefineElement(2, element, 0.0330906);
 521          element = table->GetElement(26);
 522          mix->DefineElement(3, element, 0.131443);
 523          element = table->GetElement(28);
 524          mix->DefineElement(4, element, 0.0183836);
 525       case 17: // MYLAR
 526          mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
 527          element = table->GetElement(1);
 528          mix->DefineElement(0, element, 0.0416667);
 529          element = table->GetElement(6);
 530          mix->DefineElement(1, element, 0.625);
 531          element = table->GetElement(8);
 532          mix->DefineElement(2, element, 0.333333);
 533          break;
 534       case 18: // SPDBUS(AL+KPT+EPOX)   - unknown composition
 535          mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
 536          element = table->GetElement(9);
 537          mix->DefineElement(0, element, 1.);
 538          z = element->Z();
 539          break;
 540       case 19: // SDD/SSD rings   - unknown composition
 541          mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
 542          element = table->GetElement(6);
 543          mix->DefineElement(0, element, 1.);
 544          z = element->Z();
 545          break;
 546       case 20: // SPD end ladder   - unknown composition
 547          mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
 548          element = table->GetElement(22);
 549          mix->DefineElement(0, element, 1.);
 550          z = element->Z();
 551          break;
 552       case 21: // SDD end ladder   - unknown composition
 553          mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
 554          element = table->GetElement(30);
 555          mix->DefineElement(0, element, 1.);
 556          z = element->Z();
 557          break;
 558       case 22: // SSD end ladder   - unknown composition
 559          mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
 560          element = table->GetElement(16);
 561          mix->DefineElement(0, element, 1.);
 562          z = element->Z();
 563          break;
 564    }
 565    mix->SetZ(z);
 566    printf("Patched with mixture %s\n", mix->GetName());
 567    return mix;
 568 }
 569
 570 //_____________________________________________________________________________
 571 void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
 572 {
 573   // ==== from FLUGG ====
 574   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
 575   // according to the fluka standard. In addition,. they must be equal to the
 576   // names of the fluka materials - see fluka manual - in order that the
 577   // program load the right cross sections, and equal to the names included in
 578   // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
 579   // own .pemf, in order to get the right cross sections loaded in memory.
 580   // Materials defined by FLUKA
 581    TString sname;
 582    gGeoManager->Export("geometry.root");
 583    if (fname) sname = fname;
 584    else       sname = "flukaMat.inp";
 585    ofstream out;
 586    out.open(sname.Data(), ios::out);
 587    if (!out.good()) {
 588       Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
 589       return;
 590    }
 591    PrintHeader(out, "MATERIALS AND COMPOUNDS");
 592    PrintHeader(out, "MATERIALS");
 593    Int_t i,j,idmat;
 594    Int_t counttothree, nelem;
 595    Double_t a,z,rho, w;
 596    TGeoElementTable *table = gGeoManager->GetElementTable();
 597    TGeoElement *element;
 598    element = table->GetElement(13);
 599    element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
 600    element = table->GetElement(15);
 601    element->SetTitle("PHOSPHO");  // same story ...
 602    Int_t nelements = table->GetNelements();
 603    TList *matlist = gGeoManager->GetListOfMaterials();
 604    TIter next(matlist);
 605    Int_t nmater = matlist->GetSize();
 606    Int_t nfmater =  0;
 607    Int_t newind  = 26;  // here non predefined materials start
 608
 609    TGeoMaterial *mat;
 610    TGeoMixture *mix = 0;
 611    TString matname;
 612    TObjString *objstr;
 613    // Create all needed elements
 614    for (Int_t i=1; i<nelements; i++) {
 615       element = table->GetElement(i);
 616       // skip elements which are not defined
 617       if (!element->IsUsed() && !element->IsDefined()) continue;
 618       matname = element->GetTitle();
 619       ToFlukaString(matname);
 620       rho = 0.999;
 621
 622       mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
 623       // Check if the element has been predefined by FLUKA
 624       Int_t pmid = GetPredefinedMaterialId(Int_t(element->Z()));
 625       if (pmid > 0) {
 626           mat->SetIndex(pmid);
 627       } else {
 628           mat->SetIndex(newind++);
 629       }
 630
 631       mat->SetUsed(kTRUE);
 632       fMatList->Add(mat);
 633       objstr = new TObjString(matname.Data());
 634       fMatNames->Add(objstr);
 635       nfmater++;
 636    }
 637
 638    fIndmat = nfmater;
 639    // Adjust material names and add them to FLUKA list
 640    for (i=0; i<nmater; i++) {
 641       mat = (TGeoMaterial*)matlist->At(i);
 642       if (!mat->IsUsed()) continue;
 643       z = mat->GetZ();
 644       a = mat->GetA();
 645       rho = mat->GetDensity();
 646       if (mat->GetZ()<0.001) {
 647          mat->SetIndex(2); // vacuum, built-in inside FLUKA
 648          continue;
 649       }
 650       matname = mat->GetName();
 651       FlukaMatName(matname);
 652
 653       mat->SetIndex(newind++);
 654       objstr = new TObjString(matname.Data());
 655       fMatList->Add(mat);
 656       fMatNames->Add(objstr);
 657       nfmater++;
 658    }
 659
 660    // Dump all elements with MATERIAL cards
 661    for (i=0; i<nfmater; i++) {
 662       mat = (TGeoMaterial*)fMatList->At(i);
 663 //      mat->SetUsed(kFALSE);
 664       mix = 0;
 665       out << setw(10) << "MATERIAL  ";
 666       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 667       objstr = (TObjString*)fMatNames->At(i);
 668       matname = objstr->GetString();
 669       z = mat->GetZ();
 670       a = mat->GetA();
 671       rho = mat->GetDensity();
 672       if (mat->IsMixture()) {
 673          out << setw(10) << " ";
 674          out << setw(10) << " ";
 675          mix = (TGeoMixture*)mat;
 676       } else {
 677          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
 678          out << setw(10) << setprecision(3) << a;
 679       }
 680       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 681       out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho;
 682       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 683       out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
 684       out << setw(10) << " ";
 685       out << setw(10) << " ";
 686       out << setw(8) << matname.Data() << endl;
 687       if (!mix) {
 688          // add LOW-MAT card for NEON to associate with ARGON neutron xsec
 689          if (z==10) {
 690             out << setw(10) << "LOW-MAT   ";
 691             out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 692             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
 693             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.;
 694             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.;
 695             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.;
 696             out << setw(10) << " ";
 697             out << setw(10) << " ";
 698 //            out << setw(8) << matname.Data() << endl;
 699             out << setw(8) << " " << endl;
 700          }
 701          else {
 702             element = table->GetElement((int)z);
 703             TString elename = element->GetTitle();
 704             ToFlukaString(elename);
 705             if( matname.CompareTo( elename ) != 0 ) {
 706                out << setw(10) << "LOW-MAT   ";
 707                out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 708                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
 709                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
 710                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
 711                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
 712                out << setw(10) << " ";
 713                out << setw(10) << " ";
 714                // missing material at Low Energy Cross Section Table
 715                if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 ||
 716                    (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 ||
 717                    (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 ||
 718                    (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 ||
 719                    (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 ||
 720                    (int)z==89 || (int)z==91 )
 721                   out << setw(8) << "UNKNOWN " << endl;
 722                else
 723                   out << setw(8) << elename.Data() << endl;
 724    //               out << setw(8) << " " << endl;
 725             }
 726          }
 727          continue;
 728       }
 729       counttothree = 0;
 730       out << setw(10) << "COMPOUND  ";
 731       nelem = mix->GetNelements();
 732       objstr = (TObjString*)fMatNames->At(i);
 733       matname = objstr->GetString();
 734       for (j=0; j<nelem; j++) {
 735          w = (mix->GetWmixt())[j];
 736          if (w<0.00001) w=0.00001;
 737          z = (mix->GetZmixt())[j];
 738          a = (mix->GetAmixt())[j];
 739          idmat = GetElementIndex(Int_t(z));
 740          if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z);
 741          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 742          out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w;
 743          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 744          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
 745          counttothree++;
 746          if (counttothree == 3) {
 747             out << matname.Data();
 748             out << endl;
 749             if ( (j+1) != nelem) out << setw(10) << "COMPOUND  ";
 750             counttothree = 0;
 751          }
 752       }
 753       if (nelem%3) {
 754          for (j=0; j<(3-(nelem%3)); j++)
 755             out << setw(10) << " " << setw(10) << " ";
 756          out << matname.Data();
 757          out << endl;
 758       }
 759    }
 760    Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
 761    TGeoVolume *vol;
 762    // Now print the material assignments
 763    Double_t flagfield = 0.;
 764    printf("#############################################################\n");
 765    if (gFluka->IsFieldEnabled()) {
 766       flagfield = 1.;
 767       printf("Magnetic field enabled\n");
 768    } else printf("Magnetic field disabled\n");
 769    printf("#############################################################\n");
 770
 771    PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
 772    for (i=1; i<=nvols; i++) {
 773       TGeoMedium* med;
 774       vol = gGeoManager->GetVolume(i);
 775       if ((med = vol->GetMedium()) == 0) continue;
 776       mat = med->GetMaterial();
 777       idmat = mat->GetIndex();
 778       for (Int_t j=0; j<nfmater; j++) {
 779          mat = (TGeoMaterial*)fMatList->At(j);
 780          if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE);
 781       }
 782
 783       Float_t hasfield  = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
 784       out << "* Assigning material:   " << vol->GetMedium()->GetMaterial()->GetName() << "   to Volume: " << vol->GetName();
 785       out << endl;
 786
 787       out << setw(10) << "ASSIGNMAT ";
 788       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 789       out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
 790       out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
 791       out << setw(10) << "0.0";
 792       out << setw(10) << "0.0";
 793       out << setw(10) << setiosflags(ios::fixed) <<  hasfield;
 794       out << setw(10) << "0.0";
 795       out << endl;
 796    }
 797    // dummy region
 798    idmat = 2; // vacuum
 799    fDummyRegion = nvols+1;
 800    out << "* Dummy region:   " << endl;
 801    out << setw(10) << "ASSIGNMAT ";
 802    out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 803    out << setw(10) << setiosflags(ios::fixed) << idmat;
 804    out << setw(10) << setiosflags(ios::fixed) << fDummyRegion;
 805    out << setw(10) << "0.0";
 806    out << setw(10) << "0.0";
 807    out << setw(10) << "0.0";
 808    out << setw(10) << "0.0" << endl;
 809    out.close();
 810 //   fLastMaterial = nfmater+2;
 811    fLastMaterial = newind;
 812 }
 813
 814 //_____________________________________________________________________________
 815 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
 816 {
 817 // Print a FLUKA header.
 818   out << "*\n" << "*\n" << "*\n";
 819   out << "*********************  " << text << " *********************\n"
 820      << "*\n";
 821   out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
 822      << endl;
 823   out << "*" << endl;
 824 }
 825
 826 //_____________________________________________________________________________
 827 Int_t TFlukaMCGeometry::RegionId() const
 828 {
 829 // Returns current region id <-> TGeo node id
 830    if (gGeoManager->IsOutside()) return 0;
 831    return gGeoManager->GetCurrentNode()->GetUniqueID();
 832 }
 833
 834 //_____________________________________________________________________________
 835 Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const
 836 {
 837 // Get index of a material having a given Z element.
 838    TIter next(fMatList);
 839    TGeoMaterial *mat;
 840    Int_t index = 0;
 841    while ((mat=(TGeoMaterial*)next())) {
 842       if (mat->IsMixture()) continue;
 843       if (mat->GetElement()->Z() == z) return mat->GetIndex();
 844    }
 845    return index;
 846 }
 847
 848 //_____________________________________________________________________________
 849 void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc)
 850 {
 851 // Update if needed next history;
 852 //   if (gFluka->GetDummyBoundary()==2) {
 853 //      gGeoManager->CdNode(fNextLattice-1);
 854 //      return;
 855 //   }
 856    if (lttc == TFlukaMCGeometry::kLttcOutside) {
 857       fCurrentRegion = NofVolumes()+2;
 858       fCurrentLattice = lttc;
 859       gGeoManager->CdTop();
 860       gGeoManager->SetOutside(kTRUE);
 861    }
 862    if (lttc == TFlukaMCGeometry::kLttcVirtual) return;
 863    if (lttc <=0) {
 864       Error("TFlukaMCGeometry::SetMreg","Invalide lattice %i",lttc);
 865       return;
 866    }
 867    fCurrentRegion = mreg;
 868    fCurrentLattice = lttc;
 869
 870    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
 871    if (crtlttc == lttc) return;
 872    gGeoManager->CdNode(lttc-1);
 873    while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp();
 874 }
 875
 876 //_____________________________________________________________________________
 877 void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
 878 {
 879 // Set index/history for next entered region
 880    fCurrentRegion = mreg;
 881    fCurrentLattice = latt;
 882 }
 883
 884 //_____________________________________________________________________________
 885 void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
 886 {
 887 // Set index/history for next entered region
 888    fNextRegion = mreg;
 889    fNextLattice = latt;
 890 }
 891
 892 //_____________________________________________________________________________
 893 void TFlukaMCGeometry::ToFlukaString(TString &str) const
 894 {
 895 // ToFlukaString converts an string to something usefull in FLUKA:
 896 // * Capital letters
 897 // * Only 8 letters
 898 // * Replace ' ' by '_'
 899    if (str.Length()<8) {
 900       str += "        ";
 901    }
 902    str.Remove(8);
 903    Int_t ilast;
 904    for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
 905    str.ToUpper();
 906    for (Int_t pos=0; pos<ilast; pos++)
 907       if (str(pos)==' ') str.Replace(pos,1,"_",1);
 908    return;
 909 }
 910
 911 //_____________________________________________________________________________
 912 void TFlukaMCGeometry::FlukaMatName(TString &str) const
 913 {
 914 // Strip the detector name
 915     TObjArray * tokens = str.Tokenize("_");
 916     Int_t ntok = tokens->GetEntries();
 917     if (ntok > 1) {
 918         TString head = ((TObjString*) tokens->At(0))->GetString();
 919         Int_t nhead = head.Length();
 920         str = str.Remove(0, nhead + 1);
 921     }
 922     tokens->Clear();
 923     delete tokens;
 924 // Convert a name to upper case 8 chars.
 925    ToFlukaString(str);
 926    Int_t ilast;
 927    for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
 928    if (ilast>5) ilast = 5;
 929    char number[3];
 930    TIter next(fMatNames);
 931    TObjString *objstr;
 932    TString matname;
 933    Int_t index = 0;
 934    while ((objstr=(TObjString*)next())) {
 935       matname = objstr->GetString();
 936       if (matname == str) {
 937          index++;
 938          if (index<10) {
 939             number[0] = '0';
 940             sprintf(&number[1], "%d", index);
 941          } else if (index<100) {
 942             sprintf(number, "%d", index);
 943          } else {
 944             Error("FlukaMatName", "Too many materials %s", str.Data());
 945             return;
 946          }
 947          str.Replace(ilast+1, 2, number);
 948          str.Remove(8);
 949       }
 950    }
 951 }
 952
 953 //______________________________________________________________________________
 954 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
 955 {
 956   //
 957   //  convert name to upper case. Make vname at least 4 chars
 958   //
 959   Int_t l = strlen(name);
 960   Int_t i;
 961   l = l < 4 ? l : 4;
 962   for (i=0;i<l;i++) vname[i] = toupper(name[i]);
 963   for (i=l;i<4;i++) vname[i] = ' ';
 964   vname[4] = 0;
 965 }
 966
 967
 968 //______________________________________________________________________________
 969 Int_t TFlukaMCGeometry::GetNstep()
 970 {
 971    // return gNstep for debug propose
 972    return gNstep;
 973 }
 974
 975
 976 Int_t TFlukaMCGeometry::GetPredefinedMaterialId(Int_t z) const
 977 {
 978 // Get predifined material id from Z if present in list
 979     const Int_t kMax = 25;
 980
 981     static Int_t idFluka[kMax] =
 982         {-1, -1,  1,  2,  4,  6,  7,  8,  12,  13,
 983          26, 29, 47, 14, 79, 80, 82, 73,  11,  18,
 984          20, 50, 74, 22, 28};
 985
 986     Int_t id = -1;
 987
 988     for (Int_t i = 0; i < kMax; i++)
 989     {
 990         if (z == idFluka[i]) {
 991             id = i + 1;
 992             break;
 993         }
 994
 995   }
 996
 997     return id;
 998 }
 999
1000 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1001 //
1002 //_____________________________________________________________________________
1003 Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1004 {
1005 //   from FLUGG:
1006 // Wrapper for setting DNEAR option on fluka side. Must return 0
1007 // if user doesn't want Fluka to use DNEAR to compute the
1008 // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1009 // card in fluka input), returns 1 if user wants Fluka always to
1010 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1011 // coming from all directions!!!)
1012    if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
1013    return 0;
1014 }
1015
1016 //_____________________________________________________________________________
1017 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1018           Double_t *pV,  Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
1019           Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
1020           Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
1021           Double_t *sLt, Int_t *jrLt)
1022
1023 {
1024    // Initialize FLUKa point and direction;
1025    static Int_t ierr = 0;
1026    gNstep++;
1027 //      gMCGeom->SetDebugMode(kTRUE);
1028
1029    NORLAT.xn[0] = pSx;
1030    NORLAT.xn[1] = pSy;
1031    NORLAT.xn[2] = pSz;
1032
1033    Int_t olttc = oldLttc;
1034    if (oldLttc<=0) {
1035       gGeoManager->FindNode(pSx,pSy,pSz);
1036       olttc = gGeoManager->GetCurrentNodeId()+1;
1037       oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
1038    }
1039
1040    if (gMCGeom->IsDebugging()) {
1041       cout << "g1wr gNstep=" << gNstep << " oldReg="<< oldReg <<" olttc="<< olttc
1042            << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
1043       cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ")  dir: ("
1044            << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl;
1045    }
1046
1047
1048    Int_t ccreg=0,cclat=0;
1049    gMCGeom->GetCurrentRegion(ccreg,cclat);
1050    Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
1051
1052    gMCGeom->SetCurrentRegion(oldReg, olttc);
1053    // Initialize default return values
1054    lttcFlag = 0;
1055    jrLt[lttcFlag] = olttc;
1056    sLt[lttcFlag] = propStep;
1057    jrLt[lttcFlag+1] = -1;
1058    sLt[lttcFlag+1] = 0.;
1059    newReg = oldReg;
1060    newLttc = olttc;
1061    Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE;
1062    Int_t curLttc, curReg;
1063    if (crossedDummy) {
1064    // FLUKA crossed the dummy boundary - update new region/history
1065       retStep = TGeoShape::Tolerance();
1066       saf = 0.;
1067       gMCGeom->GetNextRegion(newReg, newLttc);
1068       gMCGeom->SetMreg(newReg, newLttc);
1069       sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in current region
1070       lttcFlag++;
1071       jrLt[lttcFlag] = newLttc;
1072       sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
1073       jrLt[lttcFlag+1] = -1;
1074       sLt[lttcFlag+1] = 0.; // null step in next region;
1075       if (gMCGeom->IsDebugging()) printf("=> crossed dummy!! newReg=%i newLttc=%i\n", newReg, newLttc);
1076       return;
1077    }
1078
1079    // Reset outside flag
1080    gGeoManager->SetOutside(kFALSE);
1081
1082    curLttc = gGeoManager->GetCurrentNodeId()+1;
1083    curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1084    if (olttc != curLttc) {
1085       // FLUKA crossed the boundary : we trust that the given point is really there,
1086       // so we just update TGeo state
1087       gGeoManager->CdNode(olttc-1);
1088       curLttc = gGeoManager->GetCurrentNodeId()+1;
1089       curReg  = gGeoManager->GetCurrentVolume()->GetNumber();
1090    }
1091    // Now the current TGeo state reflects the FLUKA state
1092
1093    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1094    gGeoManager->SetCurrentDirection(pV);
1095    Double_t pt[3], local[3], ldir[3];
1096    Int_t pid = TRACKR.jtrack;
1097    pt[0] = pSx;
1098    pt[1] = pSy;
1099    pt[2] = pSz;
1100    gGeoManager->MasterToLocal(pt,local);
1101    gGeoManager->MasterToLocalVect(pV,ldir);
1102 /*
1103    Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local);
1104    if (!valid) {
1105       printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid);
1106       printf("local:(%f, %f, %f)  ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]);
1107 //      gGeoManager->FindNode();
1108 //      printf("   -> actual location: %s\n", gGeoManager->GetPath());
1109    }
1110 */
1111    Double_t pstep = propStep;
1112    Double_t snext = propStep;
1113    const Double_t epsil = 0.9999999 * TGeoShape::Tolerance();
1114    // This should never happen !!!
1115    if (pstep<TGeoShape::Tolerance()) {
1116       printf("Proposed step is 0 !!!\n");
1117       pstep = 2.*TGeoShape::Tolerance();
1118    }
1119    if (crossed) {
1120       gGeoManager->FindNextBoundaryAndStep(pstep);
1121       snext = gGeoManager->GetStep();
1122       saf = 0.;
1123       if (snext <= TGeoShape::Tolerance()) {
1124 //         printf("FLUKA crossed boundary but snext=%g\n", snext);
1125          ierr++;
1126          snext = epsil;
1127       } else {
1128          snext += TGeoShape::Tolerance();
1129          ierr = 0;
1130       }
1131    } else {
1132       gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE);
1133       snext = gGeoManager->GetStep();
1134       saf = gGeoManager->GetSafeDistance();
1135       if (snext <= TGeoShape::Tolerance()) {
1136 //         printf("FLUKA put particle on bondary without crossing\n");
1137          ierr++;
1138          snext = epsil;
1139          saf = 0.;
1140       } else {
1141          snext += TGeoShape::Tolerance();
1142          ierr = 0;
1143       }
1144       if (saf<0) saf=0.;
1145       saf -= saf*3.0e-09;
1146    }
1147 //   if (ierr>1) {
1148 //      printf("%d snext=%g\n", ierr, snext);
1149 //   }
1150    if (ierr == 10) {
1151 //      printf("Too many null steps - sending error code -33...\n");
1152       newReg = -33;   // Error code
1153       ierr = 0;
1154       return;
1155    }
1156
1157    PAREM.dist = snext;
1158    NORLAT.distn = snext;
1159    NORLAT.xn[0] += snext*pV[0];
1160    NORLAT.xn[1] += snext*pV[1];
1161    NORLAT.xn[2] += snext*pV[2];
1162    if (!gGeoManager->IsOnBoundary()) {
1163    // Next boundary further than proposed step, which is approved
1164       if (saf>propStep) saf = propStep;
1165       retStep = propStep;
1166       sLt[lttcFlag] = propStep;
1167       return;
1168    }
1169    if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed
1170    gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
1171    newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
1172    newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
1173    if (gMCGeom->IsDebugging()) printf("=> newReg=%i newLttc=%i\n", newReg, newLttc);
1174
1175    // We really crossed the boundary, but is it the same region ?
1176    gMCGeom->SetNextRegion(newReg, newLttc);
1177
1178    if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
1179       // Virtual boundary between replicants
1180       newReg  = gFluka->GetDummyRegion();
1181       newLttc = TFlukaMCGeometry::kLttcVirtual;
1182       if (gMCGeom->IsDebugging()) printf("=> virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1183    }
1184
1185    retStep = snext;
1186    sLt[lttcFlag] = snext;
1187    lttcFlag++;
1188    jrLt[lttcFlag] = newLttc;
1189    sLt[lttcFlag] = snext;
1190    jrLt[lttcFlag+1] = -1;
1191
1192    sLt[lttcFlag+1] = 0.;
1193    gGeoManager->SetOutside(kFALSE);
1194    gGeoManager->CdNode(olttc-1);
1195    if (gMCGeom->IsDebugging()) {
1196       printf("=> snext=%g safe=%g\n", snext, saf);
1197       for (Int_t i=0; i<lttcFlag+1; i++) printf("   jrLt[%i]=%i  sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1198    }
1199 }
1200
1201 //_____________________________________________________________________________
1202 void g1rtwr()
1203 {
1204
1205    if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
1206 }
1207
1208 //_____________________________________________________________________________
1209 void conhwr(Int_t & intHist, Int_t & incrCount)
1210 {
1211    if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR intHist=%d incrCount=%d currentNodeId=%d\n",
1212                                        intHist, incrCount, gGeoManager->GetCurrentNodeId()+1 );
1213 //   if( incrCount != -1 ) {
1214 //      if (intHist==0) gGeoManager->CdTop();
1215 //      else gGeoManager->CdNode(intHist-1);
1216 //   }
1217 //   intHist = gGeoManager->GetCurrentNodeId()+1;
1218 }
1219
1220 //_____________________________________________________________________________
1221 void inihwr(Int_t &intHist)
1222 {
1223    if (gMCGeom->IsDebugging())
1224       printf("========== Inside INIHWR -> reinitializing history: %i \n", intHist);
1225    if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1226    if (intHist<0) {
1227 //      printf("=== wrong history number\n");
1228       return;
1229    }
1230    if (intHist==0) gGeoManager->CdTop();
1231    else gGeoManager->CdNode(intHist-1);
1232    if (gMCGeom->IsDebugging()) {
1233       printf(" --- current path: %s\n", gGeoManager->GetPath());
1234       printf("<= INIHWR\n");
1235    }
1236 }
1237
1238 //_____________________________________________________________________________
1239 void  jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1240              Int_t &flukaReg)
1241 {
1242 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1243 // number of regions (volumes in TGeo)
1244    // build application geometry
1245    if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
1246    flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
1247    if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
1248 }
1249
1250 //_____________________________________________________________________________
1251 void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1252             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1253             Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1254 {
1255    if (gMCGeom->IsDebugging()) {
1256       printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
1257       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1258       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1259    }
1260    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1261 }
1262
1263 //_____________________________________________________________________________
1264 void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1265             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1266             Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1267 {
1268    if (gMCGeom->IsDebugging()) {
1269       printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
1270       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1271       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1272    }
1273    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1274 }
1275
1276 //_____________________________________________________________________________
1277 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1278             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1279             Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1280 {
1281    if (gMCGeom->IsDebugging()) {
1282       printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
1283       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1284       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1285    }
1286    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1287 }
1288
1289 //_____________________________________________________________________________
1290 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1291           Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1292           Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1293 {
1294    if (gMCGeom->IsDebugging()) {
1295       printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
1296       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1297       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1298    }
1299    flagErr = 0;
1300    Double_t epsil = 1000.*TGeoShape::Tolerance();
1301    TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]);
1302    if (gGeoManager->IsOutside()) {
1303       newReg = gMCGeom->NofVolumes()+2;
1304       newLttc = TFlukaMCGeometry::kLttcOutside;
1305       gGeoManager->SetOutside(kFALSE);
1306       if (oldLttc>0 && oldLttc<newLttc) gGeoManager->CdNode(oldLttc-1);
1307       return;
1308    }
1309    gGeoManager->SetOutside(kFALSE);
1310    newReg = node->GetVolume()->GetNumber();
1311    newLttc = gGeoManager->GetCurrentNodeId()+1;
1312    if (oldLttc==TFlukaMCGeometry::kLttcOutside || oldLttc==0) return;
1313
1314    Int_t dummy = gFluka->GetDummyRegion();
1315    if (oldReg==dummy) {
1316       Int_t newreg1, newlttc1;
1317       gMCGeom->GetNextRegion(newreg1, newlttc1);
1318       if (newreg1==newReg && newlttc1==newLttc) {
1319          newReg = dummy;
1320          newLttc = TFlukaMCGeometry::kLttcVirtual;
1321          flagErr = newReg;
1322          if (gMCGeom->IsDebugging()) printf("  virtual boundary (oldReg==dummy) !! newReg=%i newLttc=%i\n", newReg, newLttc);
1323       }
1324       return;
1325    }
1326
1327    if (oldReg==newReg && oldLttc!=newLttc) {
1328       newReg  = dummy;
1329       newLttc = TFlukaMCGeometry::kLttcVirtual;
1330       if (gMCGeom->IsDebugging()) printf("  virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1331    }
1332
1333    if( oldReg!=newReg && oldLttc==newLttc ) {
1334       // this should not happen!! ??? Ernesto
1335 //       cout << "  lkwr      oldReg!=newReg ("<< oldReg <<"!="<< newReg
1336 //            << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
1337       newReg  = dummy;
1338       newLttc = TFlukaMCGeometry::kLttcVirtual;
1339       flagErr = newReg;
1340    }
1341
1342    if (gMCGeom->IsDebugging()) {
1343       printf("  LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
1344    }
1345 }
1346
1347 //_____________________________________________________________________________
1348 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1349             Double_t &pVx, Double_t &pVy, Double_t &pVz,
1350             Double_t *norml, const Int_t &oldReg,
1351             const Int_t &newReg, Int_t &flagErr)
1352 {
1353    if (gMCGeom->IsDebugging()) {
1354       printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1355       printf("                         (%g, %g, %g)\n", NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1356       printf("   oldReg=%i, newReg=%i\n", oldReg,newReg);
1357    }
1358    gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1359    gGeoManager->SetCurrentDirection(pVx, pVy, pVz);
1360    Double_t *dnorm = gGeoManager->FindNormalFast();
1361    flagErr = 0;
1362    if (!dnorm) {
1363       printf("   ERROR: Cannot compute fast normal\n");
1364       flagErr = 1;
1365       norml[0] = -pVx;
1366       norml[1] = -pVy;
1367       norml[2] = -pVz;
1368    } else {
1369       norml[0] = -dnorm[0];
1370       norml[1] = -dnorm[1];
1371       norml[2] = -dnorm[2];
1372    }
1373
1374    if (gMCGeom->IsDebugging()) {
1375       printf("   normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
1376       printf("<= NRMLWR\n");
1377    }
1378
1379 }
1380
1381 //_____________________________________________________________________________
1382 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1383             Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1384 {
1385    if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
1386 }
1387
1388 //_____________________________________________________________________________
1389 Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1390 {
1391 //   from FLUGG:
1392 // Wrapper for saving current navigation history (fCheck=default)
1393 // and returning its pointer. If fCheck=-1 copy of history pointed
1394 // by intHist is made in NavHistWithCount object, and its pointer
1395 // is returned. fCheck=1 and fCheck=2 cases are only in debugging
1396 // version: an array is created by means of FGeometryInit functions
1397 // (but could be a static int * ptrArray = new int[10000] with
1398 // file scope as well) that stores a flag for deleted/undeleted
1399 // histories and at the end of event is checked to verify that
1400 // all saved history objects have been deleted.
1401
1402 // For TGeo, just return the current node ID. No copy need to be made.
1403
1404    if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR check=%d intHist=%d\n", check, intHist);
1405    if (check<0) return intHist;
1406    Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1407    if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
1408    return histInt;
1409 }
1410
1411
1412
1413