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