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