TFluka/TFlukaMCGeometry.cxx

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