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