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