TFluka/TFluka.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 //
  19 // Realisation of the TVirtualMC interface for the FLUKA code
  20 // (See official web side http://www.fluka.org/).
  21 //
  22 // This implementation makes use of the TGeo geometry modeller.
  23 // User configuration is via automatic generation of FLUKA input cards.
  24 //
  25 // Authors:
  26 // A. Fasso
  27 // E. Futo
  28 // A. Gheata
  29 // A. Morsch
  30 //
  31
  32 #include <Riostream.h>
  33
  34 #include "TFluka.h"
  35 #include "TCallf77.h"      //For the fortran calls
  36 #include "Fdblprc.h"       //(DBLPRC) fluka common
  37 #include "Fepisor.h"       //(EPISOR) fluka common
  38 #include "Ffinuc.h"        //(FINUC)  fluka common
  39 #include "Fiounit.h"       //(IOUNIT) fluka common
  40 #include "Fpaprop.h"       //(PAPROP) fluka common
  41 #include "Fpart.h"         //(PART)   fluka common
  42 #include "Ftrackr.h"       //(TRACKR) fluka common
  43 #include "Fpaprop.h"       //(PAPROP) fluka common
  44 #include "Ffheavy.h"       //(FHEAVY) fluka common
  45 #include "Fopphst.h"       //(OPPHST) fluka common
  46 #include "Fstack.h"        //(STACK)  fluka common
  47 #include "Fstepsz.h"       //(STEPSZ) fluka common
  48 #include "Fopphst.h"       //(OPPHST) fluka common
  49
  50 #include "TVirtualMC.h"
  51 #include "TMCProcess.h"
  52 #include "TGeoManager.h"
  53 #include "TGeoMaterial.h"
  54 #include "TGeoMedium.h"
  55 #include "TFlukaMCGeometry.h"
  56 #include "TGeoMCGeometry.h"
  57 #include "TFlukaCerenkov.h"
  58 #include "TFlukaConfigOption.h"
  59 #include "TFlukaScoringOption.h"
  60 #include "TLorentzVector.h"
  61 #include "TArrayI.h"
  62
  63 // Fluka methods that may be needed.
  64 #ifndef WIN32
  65 # define flukam  flukam_
  66 # define fluka_openinp fluka_openinp_
  67 # define fluka_openout fluka_openout_
  68 # define fluka_closeinp fluka_closeinp_
  69 # define mcihad mcihad_
  70 # define mpdgha mpdgha_
  71 # define newplo newplo_
  72 #else
  73 # define flukam  FLUKAM
  74 # define fluka_openinp FLUKA_OPENINP
  75 # define fluka_openout FLUKA_OPENOUT
  76 # define fluka_closeinp FLUKA_CLOSEINP
  77 # define mcihad MCIHAD
  78 # define mpdgha MPDGHA
  79 # define newplo NEWPLO
  80 #endif
  81
  82 extern "C"
  83 {
  84   //
  85   // Prototypes for FLUKA functions
  86   //
  87   void type_of_call flukam(const int&);
  88   void type_of_call newplo();
  89   void type_of_call fluka_openinp(const int&, DEFCHARA);
  90   void type_of_call fluka_openout(const int&, DEFCHARA);
  91   void type_of_call fluka_closeinp(const int&);
  92   int  type_of_call mcihad(const int&);
  93   int  type_of_call mpdgha(const int&);
  94 }
  95
  96 //
  97 // Class implementation for ROOT
  98 //
  99 ClassImp(TFluka)
 100
 101 //
 102 //----------------------------------------------------------------------------
 103 // TFluka constructors and destructors.
 104 //______________________________________________________________________________
 105 TFluka::TFluka()
 106   :TVirtualMC(),
 107    fVerbosityLevel(0),
 108    fInputFileName(""),
 109    fUserConfig(0),
 110    fUserScore(0)
 111 {
 112   //
 113   // Default constructor
 114   //
 115    fGeneratePemf = kFALSE;
 116    fNVolumes = 0;
 117    fCurrentFlukaRegion = -1;
 118    fGeom = 0;
 119    fMCGeo = 0;
 120    fMaterials = 0;
 121    fDummyBoundary = 0;
 122    fFieldFlag = 1;
 123    fStopped   = 0;
 124    fStopEvent = 0;
 125    fStopRun   = 0;
 126    fNEvent    = 0;
 127 }
 128
 129 //______________________________________________________________________________
 130 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
 131   :TVirtualMC("TFluka",title, isRootGeometrySupported),
 132    fVerbosityLevel(verbosity),
 133    fInputFileName(""),
 134    fTrackIsEntering(0),
 135    fTrackIsExiting(0),
 136    fTrackIsNew(0),
 137    fUserConfig(new TObjArray(100)),
 138    fUserScore(new TObjArray(100))
 139 {
 140   // create geometry interface
 141    if (fVerbosityLevel >=3)
 142        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 143    SetCoreInputFileName();
 144    SetInputFileName();
 145    SetGeneratePemf(kFALSE);
 146    fNVolumes      = 0;
 147    fCurrentFlukaRegion = -1;
 148    fDummyBoundary = 0;
 149    fFieldFlag = 1;
 150    fGeneratePemf = kFALSE;
 151    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
 152    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
 153    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
 154    fMaterials = 0;
 155    fStopped   = 0;
 156    fStopEvent = 0;
 157    fStopRun   = 0;
 158    fNEvent    = 0;
 159    PrintHeader();
 160 }
 161
 162 //______________________________________________________________________________
 163 TFluka::~TFluka() {
 164 // Destructor
 165     if (fVerbosityLevel >=3)
 166         cout << "<== TFluka::~TFluka() destructor called." << endl;
 167
 168     delete fGeom;
 169     delete fMCGeo;
 170
 171     if (fUserConfig) {
 172         fUserConfig->Delete();
 173         delete fUserConfig;
 174     }
 175
 176     if (fUserScore) {
 177         fUserScore->Delete();
 178         delete fUserScore;
 179     }
 180 }
 181
 182 //
 183 //______________________________________________________________________________
 184 // TFluka control methods
 185 //______________________________________________________________________________
 186 void TFluka::Init() {
 187 //
 188 //  Geometry initialisation
 189 //
 190     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
 191
 192     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
 193     fApplication->ConstructGeometry();
 194     TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
 195     gGeoManager->SetTopVolume(top);
 196     gGeoManager->CloseGeometry("di");
 197     gGeoManager->DefaultColors();  // to be removed
 198
 199     // Now we have TGeo geometry created and we have to patch FlukaVmc.inp
 200     // with the material mapping file FlukaMat.inp
 201
 202     fNVolumes = fGeom->NofVolumes();
 203     fGeom->CreateFlukaMatFile("flukaMat.inp");
 204     if (fVerbosityLevel >=3) {
 205        printf("== Number of volumes: %i\n ==", fNVolumes);
 206        cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
 207     }
 208
 209     fApplication->InitGeometry();
 210
 211 }
 212
 213
 214 //______________________________________________________________________________
 215 void TFluka::FinishGeometry() {
 216 //
 217 // Build-up table with region to medium correspondance
 218 //
 219   if (fVerbosityLevel >=3) {
 220     cout << "==> TFluka::FinishGeometry() called." << endl;
 221     printf("----FinishGeometry - nothing to do with TGeo\n");
 222     cout << "<== TFluka::FinishGeometry() called." << endl;
 223   }
 224 }
 225
 226 //______________________________________________________________________________
 227 void TFluka::BuildPhysics() {
 228 //
 229 //  Prepare FLUKA input files and call FLUKA physics initialisation
 230 //
 231
 232     if (fVerbosityLevel >=3)
 233         cout << "==> TFluka::BuildPhysics() called." << endl;
 234
 235
 236     if (fVerbosityLevel >=3) {
 237         TList *medlist = gGeoManager->GetListOfMedia();
 238         TIter next(medlist);
 239         TGeoMedium*   med = 0x0;
 240         TGeoMaterial* mat = 0x0;
 241         Int_t ic = 0;
 242
 243         while((med = (TGeoMedium*)next()))
 244         {
 245             mat = med->GetMaterial();
 246             printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
 247             ic++;
 248         }
 249     }
 250
 251     //
 252     // At this stage we have the information on materials and cuts available.
 253     // Now create the pemf file
 254
 255     if (fGeneratePemf) fGeom->CreatePemfFile();
 256
 257     //
 258     // Prepare input file with the current physics settings
 259
 260     InitPhysics();
 261
 262     cout << "\t* InitPhysics() - Prepare input file was called" << endl;
 263
 264     if (fVerbosityLevel >=2)
 265         cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
 266              << ") in fluka..." << endl;
 267     GLOBAL.lfdrtr = true;
 268
 269     if (fVerbosityLevel >=2)
 270         cout << "\t* Opening file " << fInputFileName << endl;
 271     const char* fname = fInputFileName;
 272
 273     fluka_openinp(lunin, PASSCHARA(fname));
 274     fluka_openout(11, PASSCHARA("fluka.out"));
 275
 276     if (fVerbosityLevel >=2)
 277         cout << "\t* Calling flukam..." << endl;
 278     flukam(1);
 279
 280     if (fVerbosityLevel >=2)
 281         cout << "\t* Closing file " << fInputFileName << endl;
 282     fluka_closeinp(lunin);
 283
 284     FinishGeometry();
 285
 286     if (fVerbosityLevel >=3)
 287         cout << "<== TFluka::Init() called." << endl;
 288
 289     if (fVerbosityLevel >=3)
 290         cout << "<== TFluka::BuildPhysics() called." << endl;
 291 }
 292
 293 //______________________________________________________________________________
 294 void TFluka::ProcessEvent() {
 295 //
 296 // Process one event
 297 //
 298     if (fStopRun) {
 299         printf("User Run Abortion: No more events handled !\n");
 300         fNEvent += 1;
 301         return;
 302     }
 303
 304     if (fVerbosityLevel >=3)
 305         cout << "==> TFluka::ProcessEvent() called." << endl;
 306     fApplication->GeneratePrimaries();
 307     EPISOR.lsouit = true;
 308     flukam(1);
 309     if (fVerbosityLevel >=3)
 310         cout << "<== TFluka::ProcessEvent() called." << endl;
 311     //
 312     // Increase event number
 313     //
 314     fNEvent += 1;
 315 }
 316
 317 //______________________________________________________________________________
 318 Bool_t TFluka::ProcessRun(Int_t nevent) {
 319 //
 320 // Run steering
 321 //
 322
 323   if (fVerbosityLevel >=3)
 324     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 325          << endl;
 326
 327   if (fVerbosityLevel >=2) {
 328     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 329     cout << "\t* Calling flukam again..." << endl;
 330   }
 331
 332   Int_t todo = TMath::Abs(nevent);
 333   for (Int_t ev = 0; ev < todo; ev++) {
 334       fApplication->BeginEvent();
 335       ProcessEvent();
 336       fApplication->FinishEvent();
 337   }
 338
 339   if (fVerbosityLevel >=3)
 340     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 341          << endl;
 342   // Write fluka specific scoring output
 343   newplo();
 344
 345   return kTRUE;
 346 }
 347
 348 //_____________________________________________________________________________
 349 // methods for building/management of geometry
 350
 351 // functions from GCONS
 352 //____________________________________________________________________________
 353 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 354                     Float_t &dens, Float_t &radl, Float_t &absl,
 355                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 356 //
 357    TGeoMaterial *mat;
 358    TIter next (gGeoManager->GetListOfMaterials());
 359    while ((mat = (TGeoMaterial*)next())) {
 360      if (mat->GetUniqueID() == (UInt_t)imat) break;
 361    }
 362    if (!mat) {
 363       Error("Gfmate", "no material with index %i found", imat);
 364       return;
 365    }
 366    sprintf(name, "%s", mat->GetName());
 367    a = mat->GetA();
 368    z = mat->GetZ();
 369    dens = mat->GetDensity();
 370    radl = mat->GetRadLen();
 371    absl = mat->GetIntLen();
 372 }
 373
 374 //______________________________________________________________________________
 375 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 376                     Double_t &dens, Double_t &radl, Double_t &absl,
 377                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 378 //
 379    TGeoMaterial *mat;
 380    TIter next (gGeoManager->GetListOfMaterials());
 381    while ((mat = (TGeoMaterial*)next())) {
 382      if (mat->GetUniqueID() == (UInt_t)imat) break;
 383    }
 384    if (!mat) {
 385       Error("Gfmate", "no material with index %i found", imat);
 386       return;
 387    }
 388    sprintf(name, "%s", mat->GetName());
 389    a = mat->GetA();
 390    z = mat->GetZ();
 391    dens = mat->GetDensity();
 392    radl = mat->GetRadLen();
 393    absl = mat->GetIntLen();
 394 }
 395
 396 // detector composition
 397 //______________________________________________________________________________
 398 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 399                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 400                       Float_t* buf, Int_t nwbuf) {
 401 //
 402    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
 403    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 404    delete [] dbuf;
 405 }
 406
 407 //______________________________________________________________________________
 408 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 409                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 410                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
 411 //
 412 // Define a material
 413   TGeoMaterial *mat;
 414   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 415   if ((z-Int_t(z)) > 1E-3) {
 416      mat = fGeom->GetMakeWrongMaterial(z);
 417      if (mat) {
 418         mat->SetRadLen(radl,absl);
 419         mat->SetUniqueID(kmat);
 420         return;
 421      }
 422   }
 423   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 424 }
 425
 426 //______________________________________________________________________________
 427 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 428                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 429 //
 430 // Define a material mixture
 431 //
 432   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
 433   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
 434   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
 435
 436   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 437   for (Int_t i=0; i<nlmat; i++) {
 438     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 439   }
 440
 441   delete [] da;
 442   delete [] dz;
 443   delete [] dwmat;
 444 }
 445
 446 //______________________________________________________________________________
 447 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 448                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 449 //
 450   // Defines mixture OR COMPOUND IMAT as composed by
 451   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 452   //
 453   // If NLMAT > 0 then wmat contains the proportion by
 454   // weights of each basic material in the mixture.
 455   //
 456   // If nlmat < 0 then WMAT contains the number of atoms
 457   // of a given kind into the molecule of the COMPOUND
 458   // In this case, WMAT in output is changed to relative
 459   // weigths.
 460   //
 461   Int_t i,j;
 462   if (nlmat < 0) {
 463      nlmat = - nlmat;
 464      Double_t amol = 0;
 465      for (i=0;i<nlmat;i++) {
 466         amol += a[i]*wmat[i];
 467      }
 468      for (i=0;i<nlmat;i++) {
 469         wmat[i] *= a[i]/amol;
 470      }
 471   }
 472   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 473   // Check if we have elements with fractional Z
 474   TGeoMaterial *mat = 0;
 475   TGeoMixture *mix = 0;
 476   Bool_t mixnew = kFALSE;
 477   for (i=0; i<nlmat; i++) {
 478      if (z[i]-Int_t(z[i]) < 1E-3) continue;
 479      // We have found an element with fractional Z -> loop mixtures to look for it
 480      for (j=0; j<kmat; j++) {
 481         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
 482         if (!mat) break;
 483         if (!mat->IsMixture()) continue;
 484         mix = (TGeoMixture*)mat;
 485         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
 486 //        printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
 487         mixnew = kTRUE;
 488         break;
 489      }
 490      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
 491      break;
 492   }
 493   if (mixnew) {
 494      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
 495      Double_t *anew = new Double_t[nlmatnew];
 496      Double_t *znew = new Double_t[nlmatnew];
 497      Double_t *wmatnew = new Double_t[nlmatnew];
 498      Int_t ind=0;
 499      for (j=0; j<nlmat; j++) {
 500         if (j==i) continue;
 501         anew[ind] = a[j];
 502         znew[ind] = z[j];
 503         wmatnew[ind] = wmat[j];
 504         ind++;
 505      }
 506      for (j=0; j<mix->GetNelements(); j++) {
 507         anew[ind] = mix->GetAmixt()[j];
 508         znew[ind] = mix->GetZmixt()[j];
 509         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
 510         ind++;
 511      }
 512      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 513      delete [] anew;
 514      delete [] znew;
 515      delete [] wmatnew;
 516      return;
 517   }
 518   // Now we need to compact identical elements within the mixture
 519   // First check if this happens
 520   mixnew = kFALSE;
 521   for (i=0; i<nlmat-1; i++) {
 522      for (j=i+1; j<nlmat; j++) {
 523         if (z[i] == z[j]) {
 524            mixnew = kTRUE;
 525            break;
 526         }
 527      }
 528      if (mixnew) break;
 529   }
 530   if (mixnew) {
 531      Int_t nlmatnew = 0;
 532      Double_t *anew = new Double_t[nlmat];
 533      Double_t *znew = new Double_t[nlmat];
 534      memset(znew, 0, nlmat*sizeof(Double_t));
 535      Double_t *wmatnew = new Double_t[nlmat];
 536      Bool_t skipi;
 537      for (i=0; i<nlmat; i++) {
 538         skipi = kFALSE;
 539         for (j=0; j<nlmatnew; j++) {
 540            if (z[i] == z[j]) {
 541               wmatnew[j] += wmat[i];
 542               skipi = kTRUE;
 543               break;
 544            }
 545         }
 546         if (skipi) continue;
 547         anew[nlmatnew] = a[i];
 548         znew[nlmatnew] = z[i];
 549         wmatnew[nlmatnew] = wmat[i];
 550         nlmatnew++;
 551      }
 552      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 553      delete [] anew;
 554      delete [] znew;
 555      delete [] wmatnew;
 556      return;
 557    }
 558    gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 559 }
 560
 561 //______________________________________________________________________________
 562 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 563                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 564                     Double_t stemax, Double_t deemax, Double_t epsil,
 565                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 566   // Define a medium
 567   //
 568   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 569   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 570              epsil, stmin, ubuf, nbuf);
 571 }
 572
 573 //______________________________________________________________________________
 574 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 575                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 576                     Double_t stemax, Double_t deemax, Double_t epsil,
 577                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 578   // Define a medium
 579   //
 580   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 581   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 582              epsil, stmin, ubuf, nbuf);
 583 }
 584
 585 //______________________________________________________________________________
 586 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 587                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 588                     Double_t phiZ) {
 589 //
 590   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 591   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 592 }
 593
 594 //______________________________________________________________________________
 595 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 596 //
 597 //
 598 // Check if material is used
 599    if (fVerbosityLevel >= 3)
 600        printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
 601    Int_t* reglist;
 602    Int_t nreg;
 603    reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg);
 604    if (nreg == 0) {
 605        return;
 606    }
 607
 608 //
 609    Bool_t process = kFALSE;
 610    if (strncmp(param, "DCAY",  4) == 0 ||
 611        strncmp(param, "PAIR",  4) == 0 ||
 612        strncmp(param, "COMP",  4) == 0 ||
 613        strncmp(param, "PHOT",  4) == 0 ||
 614        strncmp(param, "PFIS",  4) == 0 ||
 615        strncmp(param, "DRAY",  4) == 0 ||
 616        strncmp(param, "ANNI",  4) == 0 ||
 617        strncmp(param, "BREM",  4) == 0 ||
 618        strncmp(param, "MUNU",  4) == 0 ||
 619        strncmp(param, "CKOV",  4) == 0 ||
 620        strncmp(param, "HADR",  4) == 0 ||
 621        strncmp(param, "LOSS",  4) == 0 ||
 622        strncmp(param, "MULS",  4) == 0 ||
 623        strncmp(param, "RAYL",  4) == 0)
 624    {
 625        process = kTRUE;
 626    }
 627    if (process) {
 628        SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
 629    } else {
 630        SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
 631    }
 632 }
 633
 634 // functions from GGEOM
 635 //_____________________________________________________________________________
 636 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 637 {
 638   // Set visualisation attributes for one volume
 639   char vname[5];
 640   fGeom->Vname(name,vname);
 641   char vatt[5];
 642   fGeom->Vname(att,vatt);
 643   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 644 }
 645
 646 //______________________________________________________________________________
 647 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 648                      Float_t *upar, Int_t np)  {
 649 //
 650     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 651 }
 652
 653 //______________________________________________________________________________
 654 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 655                      Double_t *upar, Int_t np)  {
 656 //
 657     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 658 }
 659
 660 //______________________________________________________________________________
 661 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 662                    Int_t iaxis) {
 663 //
 664     fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
 665 }
 666
 667 //______________________________________________________________________________
 668 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 669                     Int_t iaxis, Double_t c0i, Int_t numed) {
 670 //
 671     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 672 }
 673
 674 //______________________________________________________________________________
 675 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 676                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 677 //
 678     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 679 }
 680
 681 //______________________________________________________________________________
 682 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 683                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 684 //
 685     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 686 }
 687
 688 //______________________________________________________________________________
 689 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 690 //
 691 // Nothing to do with TGeo
 692 }
 693
 694 //______________________________________________________________________________
 695 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 696                    Double_t x, Double_t y, Double_t z, Int_t irot,
 697                    const char *konly) {
 698 //
 699   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
 700 }
 701
 702 //______________________________________________________________________________
 703 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 704                     Double_t x, Double_t y, Double_t z, Int_t irot,
 705                     const char *konly, Float_t *upar, Int_t np)  {
 706   //
 707   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 708 }
 709
 710 //______________________________________________________________________________
 711 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 712                     Double_t x, Double_t y, Double_t z, Int_t irot,
 713                     const char *konly, Double_t *upar, Int_t np)  {
 714   //
 715   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 716 }
 717
 718 //______________________________________________________________________________
 719 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
 720 //
 721 // Nothing to do with TGeo
 722 }
 723
 724 //______________________________________________________________________________
 725 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 726                          Float_t* absco, Float_t* effic, Float_t* rindex) {
 727 //
 728 // Set Cerenkov properties for medium itmed
 729 //
 730 // npckov: number of sampling points
 731 // ppckov: energy values
 732 // absco:  absorption length
 733 // effic:  quantum efficiency
 734 // rindex: refraction index
 735 //
 736 //
 737 //
 738 //  Create object holding Cerenkov properties
 739 //
 740     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 741 //
 742 //  Pass object to medium
 743     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 744     medium->SetCerenkovProperties(cerenkovProperties);
 745 }
 746
 747 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 748                          Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
 749 //
 750 // Set Cerenkov properties for medium itmed
 751 //
 752 // npckov: number of sampling points
 753 // ppckov: energy values
 754 // absco:  absorption length
 755 // effic:  quantum efficiency
 756 // rindex: refraction index
 757 // rfl:    reflectivity for boundary to medium itmed
 758 //
 759 //
 760 //  Create object holding Cerenkov properties
 761 //
 762     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
 763 //
 764 //  Pass object to medium
 765     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 766     medium->SetCerenkovProperties(cerenkovProperties);
 767 }
 768
 769
 770 //______________________________________________________________________________
 771 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
 772                          Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
 773 //
 774 //  Double_t version not implemented
 775 }
 776
 777 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
 778                          Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) {
 779 //
 780 // //  Double_t version not implemented
 781 }
 782
 783 // Euclid
 784 //______________________________________________________________________________
 785 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
 786                           Int_t /*number*/, Int_t /*nlevel*/) {
 787 //
 788 // Not with TGeo
 789    Warning("WriteEuclid", "Not implemented with TGeo");
 790 }
 791
 792
 793
 794 //_____________________________________________________________________________
 795 // methods needed by the stepping
 796 //____________________________________________________________________________
 797
 798 Int_t TFluka::GetMedium() const {
 799 //
 800 //  Get the medium number for the current fluka region
 801 //
 802     return fGeom->GetMedium(); // this I need to check due to remapping !!!
 803 }
 804
 805
 806
 807 //____________________________________________________________________________
 808 // particle table usage
 809 // ID <--> PDG transformations
 810 //_____________________________________________________________________________
 811 Int_t TFluka::IdFromPDG(Int_t pdg) const
 812 {
 813     //
 814     // Return Fluka code from PDG and pseudo ENDF code
 815
 816     // Catch the feedback photons
 817     if (pdg == 50000051) return (-1);
 818     // MCIHAD() goes from pdg to fluka internal.
 819     Int_t intfluka = mcihad(pdg);
 820     // KPTOIP array goes from internal to official
 821     return GetFlukaKPTOIP(intfluka);
 822 }
 823
 824 //______________________________________________________________________________
 825 Int_t TFluka::PDGFromId(Int_t id) const
 826 {
 827   //
 828   // Return PDG code and pseudo ENDF code from Fluka code
 829   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon
 830     Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
 831   // IPTOKP array goes from official to internal
 832
 833     if (id == -1) {
 834 // Cerenkov photon
 835         if (fVerbosityLevel >= 3)
 836             printf("\n PDGFromId: Cerenkov Photon \n");
 837         return  50000050;
 838     }
 839 // Error id
 840     if (id == 0 || id < -6 || id > 250) {
 841         if (fVerbosityLevel >= 3)
 842             printf("PDGFromId: Error id = 0\n");
 843         return -1;
 844     }
 845 // Good id
 846     if (id > 0) {
 847         Int_t intfluka = GetFlukaIPTOKP(id);
 848         if (intfluka == 0) {
 849             if (fVerbosityLevel >= 3)
 850                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
 851             return -1;
 852         } else if (intfluka < 0) {
 853             if (fVerbosityLevel >= 3)
 854                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
 855             return -1;
 856         }
 857         if (fVerbosityLevel >= 3)
 858             printf("mpdgha called with %d %d \n", id, intfluka);
 859         // MPDGHA() goes from fluka internal to pdg.
 860         return mpdgha(intfluka);
 861     } else {
 862         // ions and optical photons
 863         return idSpecial[id + 6];
 864     }
 865 }
 866
 867 void TFluka::StopTrack()
 868 {
 869     // Set stopping conditions
 870     // Works for photons and charged particles
 871     fStopped = kTRUE;
 872 }
 873
 874 //_____________________________________________________________________________
 875 // methods for physics management
 876 //____________________________________________________________________________
 877 //
 878 // set methods
 879 //
 880
 881 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
 882 {
 883 //  Set process user flag for material imat
 884 //
 885 //
 886 //  Update if already in the list
 887 //
 888     TIter next(fUserConfig);
 889     TFlukaConfigOption* proc;
 890     while((proc = (TFlukaConfigOption*)next()))
 891     {
 892         if (proc->Medium() == imed) {
 893             proc->SetProcess(flagName, flagValue);
 894             return;
 895         }
 896     }
 897     proc = new TFlukaConfigOption(imed);
 898     proc->SetProcess(flagName, flagValue);
 899     fUserConfig->Add(proc);
 900 }
 901
 902 //______________________________________________________________________________
 903 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
 904 {
 905 //  Set process user flag
 906 //
 907 //
 908     SetProcess(flagName, flagValue, -1);
 909     return kTRUE;
 910 }
 911
 912 //______________________________________________________________________________
 913 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
 914 {
 915 // Set user cut value for material imed
 916 //
 917     TIter next(fUserConfig);
 918     TFlukaConfigOption* proc;
 919     while((proc = (TFlukaConfigOption*)next()))
 920     {
 921         if (proc->Medium() == imed) {
 922             proc->SetCut(cutName, cutValue);
 923             return;
 924         }
 925     }
 926
 927     proc = new TFlukaConfigOption(imed);
 928     proc->SetCut(cutName, cutValue);
 929     fUserConfig->Add(proc);
 930 }
 931
 932 //______________________________________________________________________________
 933 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
 934 {
 935 // Set user cut value
 936 //
 937 //
 938     SetCut(cutName, cutValue, -1);
 939     return kTRUE;
 940 }
 941
 942
 943 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what)
 944 {
 945 //
 946 // Adds a user scoring option to the list
 947 //
 948     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what);
 949     fUserScore->Add(opt);
 950 }
 951 //______________________________________________________________________________
 952 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3)
 953 {
 954 //
 955 // Adds a user scoring option to the list
 956 //
 957     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3);
 958     fUserScore->Add(opt);
 959 }
 960
 961 //______________________________________________________________________________
 962 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
 963 {
 964   printf("WARNING: Xsec not yet implemented !\n"); return -1.;
 965 }
 966
 967
 968 //______________________________________________________________________________
 969 void TFluka::InitPhysics()
 970 {
 971 //
 972 // Physics initialisation with preparation of FLUKA input cards
 973 //
 974     printf("=>InitPhysics\n");
 975
 976 // Construct file names
 977     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
 978     TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
 979     sFlukaVmcCoreInp +="/TFluka/input/";
 980     TString sFlukaVmcTmp = "flukaMat.inp";
 981     TString sFlukaVmcInp = GetInputFileName();
 982     sFlukaVmcCoreInp += GetCoreInputFileName();
 983
 984 // Open files
 985     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
 986         printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
 987         exit(1);
 988     }
 989     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
 990         printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
 991         exit(1);
 992     }
 993     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
 994         printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
 995         exit(1);
 996     }
 997
 998 // Copy core input file
 999     Char_t sLine[255];
1000     Float_t fEventsPerRun;
1001
1002     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1003         if (strncmp(sLine,"GEOEND",6) != 0)
1004             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1005         else {
1006             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
1007             goto flukamat;
1008         }
1009     } // end of while until GEOEND card
1010
1011
1012  flukamat:
1013     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1014         fprintf(pFlukaVmcInp,"%s\n",sLine);
1015     }
1016
1017     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1018         if (strncmp(sLine,"START",5) != 0)
1019             fprintf(pFlukaVmcInp,"%s\n",sLine);
1020         else {
1021             sscanf(sLine+10,"%10f",&fEventsPerRun);
1022             goto fin;
1023         }
1024     } //end of while until START card
1025
1026  fin:
1027
1028
1029 // Pass information to configuration objects
1030
1031     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1032     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1033
1034     TIter next(fUserConfig);
1035     TFlukaConfigOption* proc;
1036     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1037 //
1038 // Process Fluka specific scoring options
1039 //
1040     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1041     Float_t loginp        = 49.0;
1042     Int_t inp             = 0;
1043     Int_t nscore          = fUserScore->GetEntries();
1044
1045     TFlukaScoringOption *mopo = 0x0;
1046     TFlukaScoringOption *mopi = 0x0;
1047
1048     for (Int_t isc = 0; isc < nscore; isc++)
1049     {
1050         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1051         char*    fileName = mopo->GetFileName();
1052         Int_t    size     = strlen(fileName);
1053         Float_t  lun      = -1.;
1054 //
1055 // Check if new output file has to be opened
1056         for (Int_t isci = 0; isci < isc; isci++) {
1057
1058
1059             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1060             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1061                 //
1062                 // No, the file already exists
1063                 lun = mopi->GetLun();
1064                 mopo->SetLun(lun);
1065                 break;
1066             }
1067         } // inner loop
1068
1069         if (lun == -1.) {
1070             // Open new output file
1071             inp++;
1072             mopo->SetLun(loginp + inp);
1073             mopo->WriteOpenFlukaFile();
1074         }
1075         mopo->WriteFlukaInputCards();
1076     }
1077
1078 // Add START and STOP card
1079     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1080     fprintf(pFlukaVmcInp,"STOP      \n");
1081
1082
1083 // Close files
1084    fclose(pFlukaVmcCoreInp);
1085    fclose(pFlukaVmcFlukaMat);
1086    fclose(pFlukaVmcInp);
1087
1088
1089 //
1090 // Initialisation needed for Cerenkov photon production and transport
1091     TObjArray *matList = GetFlukaMaterials();
1092     Int_t nmaterial =  matList->GetEntriesFast();
1093     fMaterials = new Int_t[nmaterial+3];
1094
1095     for (Int_t im = 0; im < nmaterial; im++)
1096     {
1097         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1098         Int_t idmat = material->GetIndex();
1099         fMaterials[idmat] = im;
1100     }
1101 } // end of InitPhysics
1102
1103
1104 //______________________________________________________________________________
1105 void TFluka::SetMaxStep(Double_t step)
1106 {
1107 // Set the maximum step size
1108     if (step > 1.e4) return;
1109
1110     Int_t mreg, latt;
1111     fGeom->GetCurrentRegion(mreg, latt);
1112     STEPSZ.stepmx[mreg - 1] = step;
1113 }
1114
1115
1116 Double_t TFluka::MaxStep() const
1117 {
1118 // Return the maximum for current medium
1119     Int_t mreg, latt;
1120     fGeom->GetCurrentRegion(mreg, latt);
1121     return (STEPSZ.stepmx[mreg - 1]);
1122 }
1123
1124 //______________________________________________________________________________
1125 void TFluka::SetMaxNStep(Int_t)
1126 {
1127 // SetMaxNStep is dummy procedure in TFluka !
1128   if (fVerbosityLevel >=3)
1129   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1130 }
1131
1132 //______________________________________________________________________________
1133 void TFluka::SetUserDecay(Int_t)
1134 {
1135 // SetUserDecay is dummy procedure in TFluka !
1136   if (fVerbosityLevel >=3)
1137   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1138 }
1139
1140 //
1141 // dynamic properties
1142 //
1143 //______________________________________________________________________________
1144 void TFluka::TrackPosition(TLorentzVector& position) const
1145 {
1146 // Return the current position in the master reference frame of the
1147 // track being transported
1148 // TRACKR.atrack = age of the particle
1149 // TRACKR.xtrack = x-position of the last point
1150 // TRACKR.ytrack = y-position of the last point
1151 // TRACKR.ztrack = z-position of the last point
1152   Int_t caller = GetCaller();
1153   if (caller == 3 || caller == 6 || caller == 11 || caller == 12 || caller == 50) { //bxdraw,endraw,usdraw,ckov
1154     position.SetX(GetXsco());
1155     position.SetY(GetYsco());
1156     position.SetZ(GetZsco());
1157     position.SetT(TRACKR.atrack);
1158   }
1159   else if (caller == 4) { // mgdraw,mgdraw resuming
1160     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1161     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1162     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1163     position.SetT(TRACKR.atrack);
1164   }
1165   else if (caller == 5) { // sodraw
1166     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1167     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1168     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1169     position.SetT(0);
1170   } else if (caller == 40) { // mgdraw resuming transport
1171     position.SetX(TRACKR.spausr[0]);
1172     position.SetY(TRACKR.spausr[1]);
1173     position.SetZ(TRACKR.spausr[2]);
1174     position.SetT(TRACKR.spausr[3]);
1175   }
1176   else
1177     Warning("TrackPosition","position not available");
1178 }
1179
1180 //______________________________________________________________________________
1181 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1182 {
1183 // Return the current position in the master reference frame of the
1184 // track being transported
1185 // TRACKR.atrack = age of the particle
1186 // TRACKR.xtrack = x-position of the last point
1187 // TRACKR.ytrack = y-position of the last point
1188 // TRACKR.ztrack = z-position of the last point
1189   Int_t caller = GetCaller();
1190   if (caller == 3 || caller == 6 || caller == 11 || caller == 12 || caller == 50) { //bxdraw,endraw,usdraw,ckov
1191     x = GetXsco();
1192     y = GetYsco();
1193     z = GetZsco();
1194   }
1195   else if (caller == 4 || caller == 5) { // mgdraw, sodraw, mgdraw resuming
1196     x = TRACKR.xtrack[TRACKR.ntrack];
1197     y = TRACKR.ytrack[TRACKR.ntrack];
1198     z = TRACKR.ztrack[TRACKR.ntrack];
1199   }
1200   else if (caller == 40) { // mgdraw resuming transport
1201     x = TRACKR.spausr[0];
1202     y = TRACKR.spausr[1];
1203     z = TRACKR.spausr[2];
1204   }
1205   else
1206     Warning("TrackPosition","position not available");
1207 }
1208
1209 //______________________________________________________________________________
1210 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1211 {
1212 // Return the direction and the momentum (GeV/c) of the track
1213 // currently being transported
1214 // TRACKR.ptrack = momentum of the particle (not always defined, if
1215 //               < 0 must be obtained from etrack)
1216 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1217 // TRACKR.etrack = total energy of the particle
1218 // TRACKR.jtrack = identity number of the particle
1219 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1220   Int_t caller = GetCaller();
1221   if (caller != 2 && caller != 40) { // not eedraw or mgdraw resuming
1222     if (TRACKR.ptrack >= 0) {
1223       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1224       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1225       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1226       momentum.SetE(TRACKR.etrack);
1227       return;
1228     }
1229     else {
1230       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1231       momentum.SetPx(p*TRACKR.cxtrck);
1232       momentum.SetPy(p*TRACKR.cytrck);
1233       momentum.SetPz(p*TRACKR.cztrck);
1234       momentum.SetE(TRACKR.etrack);
1235       return;
1236     }
1237   } else if  (caller == 40) { // mgdraw resuming transport
1238     momentum.SetPx(TRACKR.spausr[4]);
1239     momentum.SetPy(TRACKR.spausr[5]);
1240     momentum.SetPz(TRACKR.spausr[6]);
1241     momentum.SetE (TRACKR.spausr[7]);
1242     return;
1243   }
1244   else
1245     Warning("TrackMomentum","momentum not available");
1246 }
1247
1248 //______________________________________________________________________________
1249 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1250 {
1251 // Return the direction and the momentum (GeV/c) of the track
1252 // currently being transported
1253 // TRACKR.ptrack = momentum of the particle (not always defined, if
1254 //               < 0 must be obtained from etrack)
1255 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1256 // TRACKR.etrack = total energy of the particle
1257 // TRACKR.jtrack = identity number of the particle
1258 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1259   Int_t caller = GetCaller();
1260   if (caller != 2 && caller != 40) { // not eedraw and not mgdraw resuming
1261     if (TRACKR.ptrack >= 0) {
1262       px = TRACKR.ptrack*TRACKR.cxtrck;
1263       py = TRACKR.ptrack*TRACKR.cytrck;
1264       pz = TRACKR.ptrack*TRACKR.cztrck;
1265       e = TRACKR.etrack;
1266       return;
1267     }
1268     else {
1269       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1270       px = p*TRACKR.cxtrck;
1271       py = p*TRACKR.cytrck;
1272       pz = p*TRACKR.cztrck;
1273       e = TRACKR.etrack;
1274       return;
1275     }
1276   } else if (caller == 40) { // mgdraw resuming transport
1277       px = TRACKR.spausr[4];
1278       py = TRACKR.spausr[5];
1279       pz = TRACKR.spausr[6];
1280       e  = TRACKR.spausr[7];
1281       return;
1282   }
1283   else
1284     Warning("TrackMomentum","momentum not available");
1285 }
1286
1287 //______________________________________________________________________________
1288 Double_t TFluka::TrackStep() const
1289 {
1290 // Return the length in centimeters of the current step
1291 // TRACKR.ctrack = total curved path
1292   Int_t caller = GetCaller();
1293   if (caller == 11 || caller==12 || caller == 3 || caller == 6 || caller == 50 || caller == 40) //bxdraw,endraw,usdraw, ckov
1294     return 0.0;
1295   else if (caller == 4) //mgdraw
1296     return TRACKR.ctrack;
1297   else {
1298     Warning("TrackStep", "track step not available");
1299     return 0.0;
1300   }
1301 }
1302
1303 //______________________________________________________________________________
1304 Double_t TFluka::TrackLength() const
1305 {
1306 // TRACKR.cmtrck = cumulative curved path since particle birth
1307   Int_t caller = GetCaller();
1308   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov
1309     return TRACKR.cmtrck;
1310   else if (caller == 40) // mgdraw resuming transport
1311     return TRACKR.spausr[8];
1312   else {
1313     Warning("TrackLength", "track length not available");
1314     return 0.0;
1315   }
1316 }
1317
1318 //______________________________________________________________________________
1319 Double_t TFluka::TrackTime() const
1320 {
1321 // Return the current time of flight of the track being transported
1322 // TRACKR.atrack = age of the particle
1323   Int_t caller = GetCaller();
1324   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov
1325     return TRACKR.atrack;
1326   else if (caller == 40)
1327     return TRACKR.spausr[3];
1328   else {
1329     Warning("TrackTime", "track time not available");
1330     return 0.0;
1331   }
1332 }
1333
1334 //______________________________________________________________________________
1335 Double_t TFluka::Edep() const
1336 {
1337 // Energy deposition
1338 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1339 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1340 //    but in the variable "rull" of the procedure "endraw.cxx"
1341 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1342 // -->no energy loss along the track
1343 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1344 // -->energy loss distributed along the track
1345 // TRACKR.dtrack = energy deposition of the jth deposition event
1346
1347   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1348   // If coming from usdraw we just signal particle production - no edep
1349   // If just first time after resuming, no edep for the primary
1350   Int_t caller = GetCaller();
1351   if (caller == 11 || caller==12 || caller==6 || caller == 40) return 0.0;
1352   Double_t sum = 0;
1353   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1354       sum +=TRACKR.dtrack[j];
1355   }
1356   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1357       return fRull + sum;
1358   else {
1359       return sum;
1360   }
1361 }
1362
1363 //______________________________________________________________________________
1364 Int_t TFluka::TrackPid() const
1365 {
1366 // Return the id of the particle transported
1367 // TRACKR.jtrack = identity number of the particle
1368   Int_t caller = GetCaller();
1369   if (caller != 2) { // not eedraw
1370       return PDGFromId(TRACKR.jtrack);
1371   }
1372   else
1373     return -1000;
1374 }
1375
1376 //______________________________________________________________________________
1377 Double_t TFluka::TrackCharge() const
1378 {
1379 // Return charge of the track currently transported
1380 // PAPROP.ichrge = electric charge of the particle
1381 // TRACKR.jtrack = identity number of the particle
1382   Int_t caller = GetCaller();
1383   if (caller != 2)  // not eedraw
1384     return PAPROP.ichrge[TRACKR.jtrack+6];
1385   else
1386     return -1000.0;
1387 }
1388
1389 //______________________________________________________________________________
1390 Double_t TFluka::TrackMass() const
1391 {
1392 // PAPROP.am = particle mass in GeV
1393 // TRACKR.jtrack = identity number of the particle
1394   Int_t caller = GetCaller();
1395   if (caller != 2)  // not eedraw
1396     return PAPROP.am[TRACKR.jtrack+6];
1397   else
1398     return -1000.0;
1399 }
1400
1401 //______________________________________________________________________________
1402 Double_t TFluka::Etot() const
1403 {
1404 // TRACKR.etrack = total energy of the particle
1405   Int_t caller = GetCaller();
1406   if (caller != 2)  // not eedraw
1407     return TRACKR.etrack;
1408   else
1409     return -1000.0;
1410 }
1411
1412 //
1413 // track status
1414 //
1415 //______________________________________________________________________________
1416 Bool_t   TFluka::IsNewTrack() const
1417 {
1418 // Return true for the first call of Stepping()
1419    return fTrackIsNew;
1420 }
1421
1422 void     TFluka::SetTrackIsNew(Bool_t flag)
1423 {
1424 // Return true for the first call of Stepping()
1425    fTrackIsNew = flag;
1426
1427 }
1428
1429
1430 //______________________________________________________________________________
1431 Bool_t   TFluka::IsTrackInside() const
1432 {
1433 // True if the track is not at the boundary of the current volume
1434 // In Fluka a step is always inside one kind of material
1435 // If the step would go behind the region of one material,
1436 // it will be shortened to reach only the boundary.
1437 // Therefore IsTrackInside() is always true.
1438   Int_t caller = GetCaller();
1439   if (caller == 11 || caller==12)  // bxdraw
1440     return 0;
1441   else
1442     return 1;
1443 }
1444
1445 //______________________________________________________________________________
1446 Bool_t   TFluka::IsTrackEntering() const
1447 {
1448 // True if this is the first step of the track in the current volume
1449
1450   Int_t caller = GetCaller();
1451   if (caller == 11)  // bxdraw entering
1452     return 1;
1453   else return 0;
1454 }
1455
1456 //______________________________________________________________________________
1457 Bool_t   TFluka::IsTrackExiting() const
1458 {
1459 // True if track is exiting volume
1460 //
1461   Int_t caller = GetCaller();
1462   if (caller == 12)  // bxdraw exiting
1463     return 1;
1464   else return 0;
1465 }
1466
1467 //______________________________________________________________________________
1468 Bool_t   TFluka::IsTrackOut() const
1469 {
1470 // True if the track is out of the setup
1471 // means escape
1472 // Icode = 14: escape - call from Kaskad
1473 // Icode = 23: escape - call from Emfsco
1474 // Icode = 32: escape - call from Kasneu
1475 // Icode = 40: escape - call from Kashea
1476 // Icode = 51: escape - call from Kasoph
1477   if (fIcode == 14 ||
1478       fIcode == 23 ||
1479       fIcode == 32 ||
1480       fIcode == 40 ||
1481       fIcode == 51) return 1;
1482   else return 0;
1483 }
1484
1485 //______________________________________________________________________________
1486 Bool_t   TFluka::IsTrackDisappeared() const
1487 {
1488 // means all inelastic interactions and decays
1489 // fIcode from usdraw
1490   if (fIcode == 101 || // inelastic interaction
1491       fIcode == 102 || // particle decay
1492       fIcode == 103 || // delta ray generation by hadron
1493       fIcode == 104 || // direct pair production
1494       fIcode == 105 || // bremsstrahlung (muon)
1495       fIcode == 208 || // bremsstrahlung (electron)
1496       fIcode == 214 || // in-flight annihilation
1497       fIcode == 215 || // annihilation at rest
1498       fIcode == 217 || // pair production
1499       fIcode == 219 || // Compton scattering
1500       fIcode == 221 || // Photoelectric effect
1501       fIcode == 300 || // hadronic interaction
1502       fIcode == 400    // delta-ray
1503       ) return 1;
1504   else return 0;
1505 }
1506
1507 //______________________________________________________________________________
1508 Bool_t   TFluka::IsTrackStop() const
1509 {
1510 // True if the track energy has fallen below the threshold
1511 // means stopped by signal or below energy threshold
1512 // Icode = 12: stopping particle       - call from Kaskad
1513 // Icode = 15: time kill               - call from Kaskad
1514 // Icode = 21: below threshold, iarg=1 - call from Emfsco
1515 // Icode = 22: below threshold, iarg=2 - call from Emfsco
1516 // Icode = 24: time kill               - call from Emfsco
1517 // Icode = 31: below threshold         - call from Kasneu
1518 // Icode = 33: time kill               - call from Kasneu
1519 // Icode = 41: time kill               - call from Kashea
1520 // Icode = 52: time kill               - call from Kasoph
1521   if (fIcode == 12 ||
1522       fIcode == 15 ||
1523       fIcode == 21 ||
1524       fIcode == 22 ||
1525       fIcode == 24 ||
1526       fIcode == 31 ||
1527       fIcode == 33 ||
1528       fIcode == 41 ||
1529       fIcode == 52) return 1;
1530   else return 0;
1531 }
1532
1533 //______________________________________________________________________________
1534 Bool_t   TFluka::IsTrackAlive() const
1535 {
1536 // means not disappeared or not out
1537   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1538   else return 1;
1539 }
1540
1541 //
1542 // secondaries
1543 //
1544
1545 //______________________________________________________________________________
1546 Int_t TFluka::NSecondaries() const
1547
1548 {
1549 // Number of secondary particles generated in the current step
1550 // FINUC.np = number of secondaries except light and heavy ions
1551 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
1552     Int_t caller = GetCaller();
1553     if (caller == 6)  // valid only after usdraw
1554         return FINUC.np + FHEAVY.npheav;
1555     else if (caller == 50) {
1556         // Cerenkov Photon production
1557         return fNCerenkov;
1558     }
1559     return 0;
1560 } // end of NSecondaries
1561
1562 //______________________________________________________________________________
1563 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1564                 TLorentzVector& position, TLorentzVector& momentum)
1565 {
1566 // Copy particles from secondary stack to vmc stack
1567 //
1568
1569     Int_t caller = GetCaller();
1570     if (caller == 6) {  // valid only after usdraw
1571         if (FINUC.np > 0) {
1572             // Hadronic interaction
1573             if (isec >= 0 && isec < FINUC.np) {
1574                 particleId = PDGFromId(FINUC.kpart[isec]);
1575                 position.SetX(fXsco);
1576                 position.SetY(fYsco);
1577                 position.SetZ(fZsco);
1578                 position.SetT(TRACKR.atrack);
1579                 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1580                 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1581                 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1582                 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1583             }
1584             else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1585                 Int_t jsec = isec - FINUC.np;
1586                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1587                 position.SetX(fXsco);
1588                 position.SetY(fYsco);
1589                 position.SetZ(fZsco);
1590                 position.SetT(TRACKR.atrack);
1591                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1592                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1593                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1594                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1595                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1596                 else if (FHEAVY.tkheav[jsec] > 6)
1597                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1598             }
1599             else
1600                 Warning("GetSecondary","isec out of range");
1601         }
1602     } else if (caller == 50) {
1603         Int_t index = OPPHST.lstopp - isec;
1604         position.SetX(OPPHST.xoptph[index]);
1605         position.SetY(OPPHST.yoptph[index]);
1606         position.SetZ(OPPHST.zoptph[index]);
1607         position.SetT(OPPHST.agopph[index]);
1608         Double_t p = OPPHST.poptph[index];
1609
1610         momentum.SetPx(p * OPPHST.txopph[index]);
1611         momentum.SetPy(p * OPPHST.tyopph[index]);
1612         momentum.SetPz(p * OPPHST.tzopph[index]);
1613         momentum.SetE(p);
1614     }
1615     else
1616         Warning("GetSecondary","no secondaries available");
1617
1618 } // end of GetSecondary
1619
1620
1621 //______________________________________________________________________________
1622 TMCProcess TFluka::ProdProcess(Int_t) const
1623
1624 {
1625 // Name of the process that has produced the secondary particles
1626 // in the current step
1627
1628     Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
1629
1630     if      (fIcode == 102)                  return kPDecay;
1631     else if (fIcode == 104 || fIcode == 217) return kPPair;
1632     else if (fIcode == 219)                  return kPCompton;
1633     else if (fIcode == 221)                  return kPPhotoelectric;
1634     else if (fIcode == 105 || fIcode == 208) return kPBrem;
1635     else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
1636     else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
1637     else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
1638     else if (fIcode == 101)                  return kPHadronic;
1639     else if (fIcode == 101) {
1640         if (!mugamma)                        return kPHadronic;
1641         else if (TRACKR.jtrack == 7)         return kPPhotoFission;
1642         else return kPMuonNuclear;
1643     }
1644     else if (fIcode == 225)                  return kPRayleigh;
1645 // Fluka codes 100, 300 and 400 still to be investigasted
1646     else                                     return kPNoProcess;
1647 }
1648
1649
1650 Int_t TFluka::StepProcesses(TArrayI &proc) const
1651 {
1652   //
1653   // Return processes active in the current step
1654   //
1655     proc.Set(1);
1656     TMCProcess iproc;
1657     switch (fIcode) {
1658     case 15:
1659     case 24:
1660     case 33:
1661     case 41:
1662     case 52:
1663         iproc =  kPTOFlimit;
1664         break;
1665     case 12:
1666     case 14:
1667     case 21:
1668     case 22:
1669     case 23:
1670     case 31:
1671     case 32:
1672     case 40:
1673     case 51:
1674         iproc = kPStop;
1675         break;
1676     case 50:
1677         iproc = kPLightAbsorption;
1678         break;
1679     case 59:
1680         iproc = kPLightRefraction;
1681     case 20:
1682         iproc = kPPhotoelectric;
1683         break;
1684     default:
1685         iproc = ProdProcess(0);
1686     }
1687     proc[0] = iproc;
1688     return 1;
1689 }
1690 //______________________________________________________________________________
1691 Int_t TFluka::VolId2Mate(Int_t id) const
1692 {
1693 //
1694 // Returns the material number for a given volume ID
1695 //
1696    return fMCGeo->VolId2Mate(id);
1697 }
1698
1699 //______________________________________________________________________________
1700 const char* TFluka::VolName(Int_t id) const
1701 {
1702 //
1703 // Returns the volume name for a given volume ID
1704 //
1705    return fMCGeo->VolName(id);
1706 }
1707
1708 //______________________________________________________________________________
1709 Int_t TFluka::VolId(const Text_t* volName) const
1710 {
1711 //
1712 // Converts from volume name to volume ID.
1713 // Time consuming. (Only used during set-up)
1714 // Could be replaced by hash-table
1715 //
1716     char sname[20];
1717     Int_t len;
1718     strncpy(sname, volName, len = strlen(volName));
1719     sname[len] = 0;
1720     while (sname[len - 1] == ' ') sname[--len] = 0;
1721     return fMCGeo->VolId(sname);
1722 }
1723
1724 //______________________________________________________________________________
1725 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1726 {
1727 //
1728 // Return the logical id and copy number corresponding to the current fluka region
1729 //
1730   if (gGeoManager->IsOutside()) return 0;
1731   TGeoNode *node = gGeoManager->GetCurrentNode();
1732   copyNo = node->GetNumber();
1733   Int_t id = node->GetVolume()->GetNumber();
1734   return id;
1735 }
1736
1737 //______________________________________________________________________________
1738 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1739 {
1740 //
1741 // Return the logical id and copy number of off'th mother
1742 // corresponding to the current fluka region
1743 //
1744   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1745   if (off==0) return CurrentVolID(copyNo);
1746   TGeoNode *node = gGeoManager->GetMother(off);
1747   if (!node) return 0;
1748   copyNo = node->GetNumber();
1749   return node->GetVolume()->GetNumber();
1750 }
1751
1752 //______________________________________________________________________________
1753 const char* TFluka::CurrentVolName() const
1754 {
1755 //
1756 // Return the current volume name
1757 //
1758   if (gGeoManager->IsOutside()) return 0;
1759   return gGeoManager->GetCurrentVolume()->GetName();
1760 }
1761
1762 //______________________________________________________________________________
1763 const char* TFluka::CurrentVolOffName(Int_t off) const
1764 {
1765 //
1766 // Return the volume name of the off'th mother of the current volume
1767 //
1768   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1769   if (off==0) return CurrentVolName();
1770   TGeoNode *node = gGeoManager->GetMother(off);
1771   if (!node) return 0;
1772   return node->GetVolume()->GetName();
1773 }
1774
1775 //______________________________________________________________________________
1776 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
1777                       Float_t & dens, Float_t & radl, Float_t & absl) const
1778 {
1779 //
1780 //  Return the current medium number and material properties
1781 //
1782   Int_t copy;
1783   Int_t id  =  TFluka::CurrentVolID(copy);
1784   Int_t med =  TFluka::VolId2Mate(id);
1785   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
1786   TGeoMaterial*   mat = vol->GetMaterial();
1787   a    = mat->GetA();
1788   z    = mat->GetZ();
1789   dens = mat->GetDensity();
1790   radl = mat->GetRadLen();
1791   absl = mat->GetIntLen();
1792
1793   return med;
1794 }
1795
1796 //______________________________________________________________________________
1797 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1798 {
1799 // Transforms a position from the world reference frame
1800 // to the current volume reference frame.
1801 //
1802 //  Geant3 desription:
1803 //  ==================
1804 //       Computes coordinates XD (in DRS)
1805 //       from known coordinates XM in MRS
1806 //       The local reference system can be initialized by
1807 //         - the tracking routines and GMTOD used in GUSTEP
1808 //         - a call to GMEDIA(XM,NUMED)
1809 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1810 //             (inverse routine is GDTOM)
1811 //
1812 //        If IFLAG=1  convert coordinates
1813 //           IFLAG=2  convert direction cosinus
1814 //
1815 // ---
1816    Double_t xmL[3], xdL[3];
1817    Int_t i;
1818    for (i=0;i<3;i++) xmL[i]=xm[i];
1819    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
1820    else            gGeoManager->MasterToLocalVect(xmL,xdL);
1821    for (i=0;i<3;i++) xd[i] = xdL[i];
1822 }
1823
1824 //______________________________________________________________________________
1825 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1826 {
1827    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
1828    else            gGeoManager->MasterToLocalVect(xm,xd);
1829 }
1830
1831 //______________________________________________________________________________
1832 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1833 {
1834 // Transforms a position from the current volume reference frame
1835 // to the world reference frame.
1836 //
1837 //  Geant3 desription:
1838 //  ==================
1839 //  Computes coordinates XM (Master Reference System
1840 //  knowing the coordinates XD (Detector Ref System)
1841 //  The local reference system can be initialized by
1842 //    - the tracking routines and GDTOM used in GUSTEP
1843 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1844 //        (inverse routine is GMTOD)
1845 //
1846 //   If IFLAG=1  convert coordinates
1847 //      IFLAG=2  convert direction cosinus
1848 //
1849 // ---
1850    Double_t xmL[3], xdL[3];
1851    Int_t i;
1852    for (i=0;i<3;i++) xdL[i] = xd[i];
1853    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
1854    else            gGeoManager->LocalToMasterVect(xdL,xmL);
1855    for (i=0;i<3;i++) xm[i]=xmL[i];
1856 }
1857
1858 //______________________________________________________________________________
1859 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1860 {
1861    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
1862    else            gGeoManager->LocalToMasterVect(xd,xm);
1863 }
1864
1865 //______________________________________________________________________________
1866 TObjArray *TFluka::GetFlukaMaterials()
1867 {
1868    return fGeom->GetMatList();
1869 }
1870
1871 //______________________________________________________________________________
1872 void TFluka::SetMreg(Int_t l)
1873 {
1874 // Set current fluka region
1875    fCurrentFlukaRegion = l;
1876    fGeom->SetMreg(l);
1877 }
1878
1879
1880
1881
1882 TString TFluka::ParticleName(Int_t pdg) const
1883 {
1884     // Return particle name for particle with pdg code pdg.
1885     Int_t ifluka = IdFromPDG(pdg);
1886     return TString((CHPPRP.btype[ifluka+6]), 8);
1887 }
1888
1889
1890 Double_t TFluka::ParticleMass(Int_t pdg) const
1891 {
1892     // Return particle mass for particle with pdg code pdg.
1893     Int_t ifluka = IdFromPDG(pdg);
1894     return (PAPROP.am[ifluka+6]);
1895 }
1896
1897 Double_t TFluka::ParticleCharge(Int_t pdg) const
1898 {
1899     // Return particle charge for particle with pdg code pdg.
1900     Int_t ifluka = IdFromPDG(pdg);
1901     return Double_t(PAPROP.ichrge[ifluka+6]);
1902 }
1903
1904 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
1905 {
1906     // Return particle lifetime for particle with pdg code pdg.
1907     Int_t ifluka = IdFromPDG(pdg);
1908     return (PAPROP.thalf[ifluka+6]);
1909 }
1910
1911 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
1912 {
1913     // Retrieve particle properties for particle with pdg code pdg.
1914
1915     strcpy(name, ParticleName(pdg).Data());
1916     type   = ParticleMCType(pdg);
1917     mass   = ParticleMass(pdg);
1918     charge = ParticleCharge(pdg);
1919     tlife  = ParticleLifeTime(pdg);
1920 }
1921
1922 void TFluka::PrintHeader()
1923 {
1924     //
1925     // Print a header
1926     printf("\n");
1927     printf("\n");
1928     printf("------------------------------------------------------------------------------\n");
1929     printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA.           -\n");
1930     printf("- Please see the file fluka.out for FLUKA output and licensing information.  -\n");
1931     printf("------------------------------------------------------------------------------\n");
1932     printf("\n");
1933     printf("\n");
1934 }
1935
1936
1937
1938 #define pushcerenkovphoton pushcerenkovphoton_
1939 #define usersteppingckv    usersteppingckv_
1940
1941
1942 extern "C" {
1943     void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
1944                             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
1945                             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
1946     {
1947         //
1948         // Pushes one cerenkov photon to the stack
1949         //
1950
1951         TFluka* fluka =  (TFluka*) gMC;
1952         TVirtualMCStack* cppstack = fluka->GetStack();
1953         Int_t parent =  TRACKR.ispusr[mkbmx2-1];
1954         cppstack->PushTrack(0, parent, 50000050,
1955                             px, py, pz, e,
1956                             vx, vy, vz, tof,
1957                             polx, poly, polz,
1958                             kPCerenkov, ntr, wgt, 0);
1959     }
1960
1961     void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
1962     {
1963         //
1964         // Calls stepping in order to signal cerenkov production
1965         //
1966         TFluka *fluka = (TFluka*)gMC;
1967         fluka->SetMreg(mreg);
1968         fluka->SetXsco(x);
1969         fluka->SetYsco(y);
1970         fluka->SetZsco(z);
1971         fluka->SetNCerenkov(nphot);
1972         fluka->SetCaller(50);
1973         printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f %s\n", nphot, mreg, x, y, z, fluka->CurrentVolName());
1974         (TVirtualMCApplication::Instance())->Stepping();
1975     }
1976 }
1977