TFluka/TFlukaMCGeometry.cxx

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