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 #include <TList.h>
  34
  35 #include "TFluka.h"
  36 #include "TFlukaIon.h"
  37 #include "TFlukaCodes.h"
  38 #include "TCallf77.h"      //For the fortran calls
  39 #include "Fdblprc.h"       //(DBLPRC) fluka common
  40 #include "Fsourcm.h"       //(SOURCM) fluka common
  41 #include "Fgenstk.h"       //(GENSTK)  fluka common
  42 #include "Fiounit.h"       //(IOUNIT) fluka common
  43 #include "Fpaprop.h"       //(PAPROP) fluka common
  44 #include "Fpart.h"         //(PART)   fluka common
  45 #include "Ftrackr.h"       //(TRACKR) fluka common
  46 #include "Fpaprop.h"       //(PAPROP) fluka common
  47 #include "Ffheavy.h"       //(FHEAVY) fluka common
  48 #include "Fopphst.h"       //(OPPHST) fluka common
  49 #include "Fflkstk.h"       //(FLKSTK) fluka common
  50 #include "Fstepsz.h"       //(STEPSZ) fluka common
  51 #include "Fopphst.h"       //(OPPHST) fluka common
  52 #include "Fltclcm.h"       //(LTCLCM) fluka common
  53 #include "Falldlt.h"       //(ALLDLT) fluka common
  54
  55 #include "TVirtualMC.h"
  56 #include "TMCProcess.h"
  57 #include "TGeoManager.h"
  58 #include "TGeoMaterial.h"
  59 #include "TGeoMedium.h"
  60 #include "TFlukaMCGeometry.h"
  61 #include "TGeoMCGeometry.h"
  62 #include "TFlukaCerenkov.h"
  63 #include "TFlukaConfigOption.h"
  64 #include "TFlukaScoringOption.h"
  65 #include "TLorentzVector.h"
  66 #include "TArrayI.h"
  67 #include "TArrayD.h"
  68 #include "TDatabasePDG.h"
  69 #include "TStopwatch.h"
  70
  71
  72 // Fluka methods that may be needed.
  73 #ifndef WIN32
  74 # define flukam  flukam_
  75 # define fluka_openinp fluka_openinp_
  76 # define fluka_openout fluka_openout_
  77 # define fluka_closeinp fluka_closeinp_
  78 # define mcihad mcihad_
  79 # define mpdgha mpdgha_
  80 # define newplo newplo_
  81 # define genout genout_
  82 # define flkend flkend_
  83 #else
  84 # define flukam  FLUKAM
  85 # define fluka_openinp FLUKA_OPENINP
  86 # define fluka_openout FLUKA_OPENOUT
  87 # define fluka_closeinp FLUKA_CLOSEINP
  88 # define mcihad MCIHAD
  89 # define mpdgha MPDGHA
  90 # define newplo NEWPLO
  91 # define genout GENOUT
  92 # define flkend FLKEND
  93 #endif
  94
  95 extern "C"
  96 {
  97   //
  98   // Prototypes for FLUKA functions
  99   //
 100   void type_of_call flukam(const int&);
 101   void type_of_call newplo();
 102   void type_of_call genout();
 103   void type_of_call flkend();
 104   void type_of_call fluka_openinp(const int&, DEFCHARA);
 105   void type_of_call fluka_openout(const int&, DEFCHARA);
 106   void type_of_call fluka_closeinp(const int&);
 107   int  type_of_call mcihad(const int&);
 108   int  type_of_call mpdgha(const int&);
 109 }
 110
 111 //
 112 // Class implementation for ROOT
 113 //
 114 ClassImp(TFluka)
 115
 116 //
 117 //----------------------------------------------------------------------------
 118 // TFluka constructors and destructors.
 119 //______________________________________________________________________________
 120 TFluka::TFluka()
 121   :TVirtualMC(),
 122    fVerbosityLevel(0),
 123    fNEvent(0),
 124    fInputFileName(""),
 125    fCoreInputFileName(""),
 126    fCaller(kNoCaller),
 127    fIcode(kNoProcess),
 128    fNewReg(-1),
 129    fRull(0),
 130    fXsco(0),
 131    fYsco(0),
 132    fZsco(0),
 133    fPItime(0),
 134    fPIlength(0),
 135    fNPI(0),
 136    fTrackIsEntering(kFALSE),
 137    fTrackIsExiting(kFALSE),
 138    fTrackIsNew(kFALSE),
 139    fFieldFlag(kTRUE),
 140    fDummyBoundary(kFALSE),
 141    fStopped(kFALSE),
 142    fStopEvent(kFALSE),
 143    fStopRun(kFALSE),
 144    fPrimaryElectronIndex(-1),
 145    fLowEnergyNeutronTransport(kFALSE),
 146    fMaterials(0),
 147    fNVolumes(0),
 148    fCurrentFlukaRegion(-1),
 149    fNCerenkov(0),
 150    fGeom(0),
 151    fMCGeo(0),
 152    fUserConfig(0),
 153    fUserScore(0),
 154    fUserIon(0)
 155 {
 156   //
 157   // Default constructor
 158   //
 159     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
 160 }
 161
 162 //______________________________________________________________________________
 163 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
 164   :TVirtualMC("TFluka",title, isRootGeometrySupported),
 165    fVerbosityLevel(verbosity),
 166    fNEvent(0),
 167    fInputFileName(""),
 168    fCoreInputFileName(""),
 169    fCaller(kNoCaller),
 170    fIcode(kNoProcess),
 171    fNewReg(-1),
 172    fRull(0),
 173    fXsco(0),
 174    fYsco(0),
 175    fZsco(0),
 176    fPItime(0),
 177    fPIlength(0),
 178    fNPI(0),
 179    fTrackIsEntering(kFALSE),
 180    fTrackIsExiting(kFALSE),
 181    fTrackIsNew(kFALSE),
 182    fFieldFlag(kTRUE),
 183    fDummyBoundary(kFALSE),
 184    fStopped(kFALSE),
 185    fStopEvent(kFALSE),
 186    fStopRun(kFALSE),
 187    fPrimaryElectronIndex(-1),
 188    fLowEnergyNeutronTransport(kFALSE),
 189    fMaterials(0),
 190    fNVolumes(0),
 191    fCurrentFlukaRegion(-1),
 192    fNCerenkov(0),
 193    fGeom(0),
 194    fMCGeo(0),
 195    fUserConfig(new TObjArray(100)),
 196    fUserScore(new TObjArray(100)),
 197    fUserIon(0)
 198 {
 199   // create geometry interface
 200     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
 201
 202    if (fVerbosityLevel >=3)
 203        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 204    SetCoreInputFileName();
 205    SetInputFileName();
 206    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
 207    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
 208    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
 209    PrintHeader();
 210 }
 211
 212 //______________________________________________________________________________
 213 TFluka::~TFluka()
 214 {
 215     // Destructor
 216     if (fVerbosityLevel >=3)
 217         cout << "<== TFluka::~TFluka() destructor called." << endl;
 218     if (fMaterials) delete [] fMaterials;
 219
 220 //    delete fGeom;
 221 //    delete fMCGeo;
 222
 223     if (fUserConfig) {
 224         fUserConfig->Delete();
 225         delete fUserConfig;
 226     }
 227
 228     if (fUserScore) {
 229         fUserScore->Delete();
 230         delete fUserScore;
 231     }
 232 }
 233
 234 //
 235 //______________________________________________________________________________
 236 // TFluka control methods
 237 //______________________________________________________________________________
 238 void TFluka::Init() {
 239 //
 240 //  Geometry initialisation
 241 //
 242     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
 243
 244     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
 245     fApplication->ConstructGeometry();
 246     if (!gGeoManager->IsClosed()) {
 247        TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
 248        gGeoManager->SetTopVolume(top);
 249        gGeoManager->CloseGeometry("di");
 250     } else {
 251        TGeoNodeCache *cache = gGeoManager->GetCache();
 252        if (!cache->HasIdArray()) {
 253           Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
 254           cache->BuildIdArray();
 255        }
 256     }
 257     fNVolumes = fGeom->NofVolumes();
 258     fGeom->CreateFlukaMatFile("flukaMat.inp");
 259     if (fVerbosityLevel >=3) {
 260        printf("== Number of volumes: %i\n ==", fNVolumes);
 261        cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
 262     }
 263
 264     fApplication->InitGeometry();
 265     fApplication->ConstructOpGeometry();
 266     //
 267     // Add ions to PDG Data base
 268     //
 269      AddParticlesToPdgDataBase();
 270     //
 271 }
 272
 273
 274 //______________________________________________________________________________
 275 void TFluka::FinishGeometry() {
 276 //
 277 // Build-up table with region to medium correspondance
 278 //
 279   if (fVerbosityLevel >=3) {
 280     cout << "==> TFluka::FinishGeometry() called." << endl;
 281     printf("----FinishGeometry - applying misalignment if any\n");
 282     cout << "<== TFluka::FinishGeometry() called." << endl;
 283   }
 284   TVirtualMCApplication::Instance()->MisalignGeometry();
 285 }
 286
 287 //______________________________________________________________________________
 288 void TFluka::BuildPhysics() {
 289 //
 290 //  Prepare FLUKA input files and call FLUKA physics initialisation
 291 //
 292
 293     if (fVerbosityLevel >=3)
 294         cout << "==> TFluka::BuildPhysics() called." << endl;
 295
 296
 297     if (fVerbosityLevel >=3) {
 298         TList *medlist = gGeoManager->GetListOfMedia();
 299         TIter next(medlist);
 300         TGeoMedium*   med = 0x0;
 301         TGeoMaterial* mat = 0x0;
 302         Int_t ic = 0;
 303
 304         while((med = (TGeoMedium*)next()))
 305         {
 306             mat = med->GetMaterial();
 307             printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
 308             ic++;
 309         }
 310     }
 311
 312
 313     // Prepare input file with the current physics settings
 314
 315     InitPhysics();
 316 //  Open fortran files
 317     const char* fname = fInputFileName;
 318     fluka_openinp(lunin, PASSCHARA(fname));
 319     fluka_openout(11, PASSCHARA("fluka.out"));
 320 //  Read input cards
 321     cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
 322     TStopwatch timer;
 323     timer.Start();
 324     GLOBAL.lfdrtr = true;
 325     flukam(1);
 326     cout << "<== TFluka::BuildPhysics() Read input cards End"
 327          << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
 328 //  Close input file
 329     fluka_closeinp(lunin);
 330 //  Finish geometry
 331     FinishGeometry();
 332 }
 333
 334 //______________________________________________________________________________
 335 void TFluka::ProcessEvent() {
 336 //
 337 // Process one event
 338 //
 339     if (fStopRun) {
 340         Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
 341         fNEvent += 1;
 342         return;
 343     }
 344
 345     if (fVerbosityLevel >=3)
 346         cout << "==> TFluka::ProcessEvent() called." << endl;
 347     fApplication->GeneratePrimaries();
 348     SOURCM.lsouit = true;
 349     flukam(1);
 350     if (fVerbosityLevel >=3)
 351         cout << "<== TFluka::ProcessEvent() called." << endl;
 352     //
 353     // Increase event number
 354     //
 355     fNEvent += 1;
 356 }
 357
 358 //______________________________________________________________________________
 359 Bool_t TFluka::ProcessRun(Int_t nevent) {
 360 //
 361 // Run steering
 362 //
 363
 364   if (fVerbosityLevel >=3)
 365     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 366          << endl;
 367
 368   if (fVerbosityLevel >=2) {
 369     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 370     cout << "\t* Calling flukam again..." << endl;
 371   }
 372
 373   Int_t todo = TMath::Abs(nevent);
 374   for (Int_t ev = 0; ev < todo; ev++) {
 375       TStopwatch timer;
 376       timer.Start();
 377       fApplication->BeginEvent();
 378       ProcessEvent();
 379       fApplication->FinishEvent();
 380       cout << "Event: "<< ev
 381            << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
 382   }
 383
 384   if (fVerbosityLevel >=3)
 385     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 386          << endl;
 387
 388   // Write fluka specific scoring output
 389   genout();
 390   newplo();
 391   flkend();
 392
 393   return kTRUE;
 394 }
 395
 396 //_____________________________________________________________________________
 397 // methods for building/management of geometry
 398
 399 // functions from GCONS
 400 //____________________________________________________________________________
 401 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 402                     Float_t &dens, Float_t &radl, Float_t &absl,
 403                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 404 //
 405    TGeoMaterial *mat;
 406    TIter next (gGeoManager->GetListOfMaterials());
 407    while ((mat = (TGeoMaterial*)next())) {
 408      if (mat->GetUniqueID() == (UInt_t)imat) break;
 409    }
 410    if (!mat) {
 411       Error("Gfmate", "no material with index %i found", imat);
 412       return;
 413    }
 414    sprintf(name, "%s", mat->GetName());
 415    a = mat->GetA();
 416    z = mat->GetZ();
 417    dens = mat->GetDensity();
 418    radl = mat->GetRadLen();
 419    absl = mat->GetIntLen();
 420 }
 421
 422 //______________________________________________________________________________
 423 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 424                     Double_t &dens, Double_t &radl, Double_t &absl,
 425                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 426 //
 427    TGeoMaterial *mat;
 428    TIter next (gGeoManager->GetListOfMaterials());
 429    while ((mat = (TGeoMaterial*)next())) {
 430      if (mat->GetUniqueID() == (UInt_t)imat) break;
 431    }
 432    if (!mat) {
 433       Error("Gfmate", "no material with index %i found", imat);
 434       return;
 435    }
 436    sprintf(name, "%s", mat->GetName());
 437    a = mat->GetA();
 438    z = mat->GetZ();
 439    dens = mat->GetDensity();
 440    radl = mat->GetRadLen();
 441    absl = mat->GetIntLen();
 442 }
 443
 444 // detector composition
 445 //______________________________________________________________________________
 446 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 447                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 448                       Float_t* buf, Int_t nwbuf) {
 449 //
 450    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
 451    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 452    delete [] dbuf;
 453 }
 454
 455 //______________________________________________________________________________
 456 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 457                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 458                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
 459 //
 460 // Define a material
 461   TGeoMaterial *mat;
 462   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 463   if ((z-Int_t(z)) > 1E-3) {
 464      mat = fGeom->GetMakeWrongMaterial(z);
 465      if (mat) {
 466         mat->SetRadLen(radl,absl);
 467         mat->SetUniqueID(kmat);
 468         return;
 469      }
 470   }
 471   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 472 }
 473
 474 //______________________________________________________________________________
 475 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 476                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 477 //
 478 // Define a material mixture
 479 //
 480   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
 481   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
 482   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
 483
 484   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 485   for (Int_t i=0; i<nlmat; i++) {
 486     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 487   }
 488
 489   delete [] da;
 490   delete [] dz;
 491   delete [] dwmat;
 492 }
 493
 494 //______________________________________________________________________________
 495 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 496                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 497 //
 498   // Defines mixture OR COMPOUND IMAT as composed by
 499   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 500   //
 501   // If NLMAT > 0 then wmat contains the proportion by
 502   // weights of each basic material in the mixture.
 503   //
 504   // If nlmat < 0 then WMAT contains the number of atoms
 505   // of a given kind into the molecule of the COMPOUND
 506   // In this case, WMAT in output is changed to relative
 507   // weigths.
 508   //
 509   Int_t i,j;
 510   if (nlmat < 0) {
 511      nlmat = - nlmat;
 512      Double_t amol = 0;
 513      for (i=0;i<nlmat;i++) {
 514         amol += a[i]*wmat[i];
 515      }
 516      for (i=0;i<nlmat;i++) {
 517         wmat[i] *= a[i]/amol;
 518      }
 519   }
 520   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 521   // Check if we have elements with fractional Z
 522   TGeoMaterial *mat = 0;
 523   TGeoMixture *mix = 0;
 524   Bool_t mixnew = kFALSE;
 525   for (i=0; i<nlmat; i++) {
 526      if (z[i]-Int_t(z[i]) < 1E-3) continue;
 527      // We have found an element with fractional Z -> loop mixtures to look for it
 528      for (j=0; j<kmat; j++) {
 529         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
 530         if (!mat) break;
 531         if (!mat->IsMixture()) continue;
 532         mix = (TGeoMixture*)mat;
 533         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
 534         mixnew = kTRUE;
 535         break;
 536      }
 537      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
 538      break;
 539   }
 540   if (mixnew) {
 541      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
 542      Double_t *anew = new Double_t[nlmatnew];
 543      Double_t *znew = new Double_t[nlmatnew];
 544      Double_t *wmatnew = new Double_t[nlmatnew];
 545      Int_t ind=0;
 546      for (j=0; j<nlmat; j++) {
 547         if (j==i) continue;
 548         anew[ind] = a[j];
 549         znew[ind] = z[j];
 550         wmatnew[ind] = wmat[j];
 551         ind++;
 552      }
 553      for (j=0; j<mix->GetNelements(); j++) {
 554         anew[ind] = mix->GetAmixt()[j];
 555         znew[ind] = mix->GetZmixt()[j];
 556         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
 557         ind++;
 558      }
 559      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 560      delete [] anew;
 561      delete [] znew;
 562      delete [] wmatnew;
 563      return;
 564   }
 565   // Now we need to compact identical elements within the mixture
 566   // First check if this happens
 567   mixnew = kFALSE;
 568   for (i=0; i<nlmat-1; i++) {
 569      for (j=i+1; j<nlmat; j++) {
 570         if (z[i] == z[j]) {
 571            mixnew = kTRUE;
 572            break;
 573         }
 574      }
 575      if (mixnew) break;
 576   }
 577   if (mixnew) {
 578      Int_t nlmatnew = 0;
 579      Double_t *anew = new Double_t[nlmat];
 580      Double_t *znew = new Double_t[nlmat];
 581      memset(znew, 0, nlmat*sizeof(Double_t));
 582      Double_t *wmatnew = new Double_t[nlmat];
 583      Bool_t skipi;
 584      for (i=0; i<nlmat; i++) {
 585         skipi = kFALSE;
 586         for (j=0; j<nlmatnew; j++) {
 587            if (z[i] == z[j]) {
 588               wmatnew[j] += wmat[i];
 589               skipi = kTRUE;
 590               break;
 591            }
 592         }
 593         if (skipi) continue;
 594         anew[nlmatnew] = a[i];
 595         znew[nlmatnew] = z[i];
 596         wmatnew[nlmatnew] = wmat[i];
 597         nlmatnew++;
 598      }
 599      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 600      delete [] anew;
 601      delete [] znew;
 602      delete [] wmatnew;
 603      return;
 604    }
 605   gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 606 }
 607
 608 //______________________________________________________________________________
 609 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 610                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 611                     Double_t stemax, Double_t deemax, Double_t epsil,
 612                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 613   // Define a medium
 614   //
 615   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 616   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 617              epsil, stmin, ubuf, nbuf);
 618 }
 619
 620 //______________________________________________________________________________
 621 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 622                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 623                     Double_t stemax, Double_t deemax, Double_t epsil,
 624                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 625   // Define a medium
 626   //
 627   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 628   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 629              epsil, stmin, ubuf, nbuf);
 630 }
 631
 632 //______________________________________________________________________________
 633 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 634                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 635                     Double_t phiZ) {
 636 //
 637   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 638   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 639 }
 640
 641 //______________________________________________________________________________
 642 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 643 //
 644 //
 645 //
 646    Bool_t process = kFALSE;
 647    Bool_t modelp  = kFALSE;
 648
 649    if (strncmp(param, "DCAY",  4) == 0 ||
 650        strncmp(param, "PAIR",  4) == 0 ||
 651        strncmp(param, "COMP",  4) == 0 ||
 652        strncmp(param, "PHOT",  4) == 0 ||
 653        strncmp(param, "PFIS",  4) == 0 ||
 654        strncmp(param, "DRAY",  4) == 0 ||
 655        strncmp(param, "ANNI",  4) == 0 ||
 656        strncmp(param, "BREM",  4) == 0 ||
 657        strncmp(param, "MUNU",  4) == 0 ||
 658        strncmp(param, "CKOV",  4) == 0 ||
 659        strncmp(param, "HADR",  4) == 0 ||
 660        strncmp(param, "LOSS",  4) == 0 ||
 661        strncmp(param, "MULS",  4) == 0 ||
 662        strncmp(param, "RAYL",  4) == 0 ||
 663        strncmp(param, "STRA",  4) == 0)
 664    {
 665        process = kTRUE;
 666    }
 667
 668    if (strncmp(param, "PRIMIO_N",  8) == 0 ||
 669        strncmp(param, "PRIMIO_E",  8) == 0)
 670    {
 671        modelp = kTRUE;
 672    }
 673
 674    if (process) {
 675        // Process switch
 676        SetProcess(param, Int_t (parval), itmed);
 677    } else if (modelp) {
 678        // Model parameters
 679        SetModelParameter(param, parval, itmed);
 680    } else {
 681        // Cuts
 682        SetCut(param, parval, itmed);
 683    }
 684
 685
 686 }
 687
 688 // functions from GGEOM
 689 //_____________________________________________________________________________
 690 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 691 {
 692   // Set visualisation attributes for one volume
 693   char vname[5];
 694   fGeom->Vname(name,vname);
 695   char vatt[5];
 696   fGeom->Vname(att,vatt);
 697   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 698 }
 699
 700 //______________________________________________________________________________
 701 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 702                      Float_t *upar, Int_t np)  {
 703 //
 704     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 705 }
 706
 707 //______________________________________________________________________________
 708 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 709                      Double_t *upar, Int_t np)  {
 710 //
 711     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 712 }
 713
 714 //______________________________________________________________________________
 715 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 716                    Int_t iaxis) {
 717 //
 718     fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
 719 }
 720
 721 //______________________________________________________________________________
 722 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 723                     Int_t iaxis, Double_t c0i, Int_t numed) {
 724 //
 725     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 726 }
 727
 728 //______________________________________________________________________________
 729 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 730                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 731 //
 732     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 733 }
 734
 735 //______________________________________________________________________________
 736 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 737                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 738 //
 739     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 740 }
 741
 742 //______________________________________________________________________________
 743 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 744 //
 745 // Nothing to do with TGeo
 746 }
 747
 748 //______________________________________________________________________________
 749 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 750                    Double_t x, Double_t y, Double_t z, Int_t irot,
 751                    const char *konly) {
 752 //
 753   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
 754 }
 755
 756 //______________________________________________________________________________
 757 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 758                     Double_t x, Double_t y, Double_t z, Int_t irot,
 759                     const char *konly, Float_t *upar, Int_t np)  {
 760   //
 761   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 762 }
 763
 764 //______________________________________________________________________________
 765 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 766                     Double_t x, Double_t y, Double_t z, Int_t irot,
 767                     const char *konly, Double_t *upar, Int_t np)  {
 768   //
 769   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 770 }
 771
 772 //______________________________________________________________________________
 773 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
 774 //
 775 // Nothing to do with TGeo
 776 }
 777
 778 //______________________________________________________________________
 779 Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
 780 {
 781     // Returns the Transformation matrix between the volume specified
 782     // by the path volumePath and the Top or mater volume. The format
 783     // of the path volumePath is as follows (assuming ALIC is the Top volume)
 784     // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
 785     // or master volume which has only 1 instance of. Of all of the daughter
 786     // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
 787     // the daughter volume of DDIP is S05I copy #2 and so on.
 788     // Inputs:
 789     //   TString& volumePath  The volume path to the specific volume
 790     //                        for which you want the matrix. Volume name
 791     //                        hierarchy is separated by "/" while the
 792     //                        copy number is appended using a "_".
 793     // Outputs:
 794     //  TGeoHMatrix &mat      A matrix with its values set to those
 795     //                        appropriate to the Local to Master transformation
 796     // Return:
 797     //   A logical value if kFALSE then an error occurred and no change to
 798     //   mat was made.
 799
 800    // We have to preserve the modeler state
 801    return fMCGeo->GetTransformation(volumePath, mat);
 802 }
 803
 804 //______________________________________________________________________
 805 Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
 806                         TArrayD &par)
 807 {
 808     // Returns the shape and its parameters for the volume specified
 809     // by volumeName.
 810     // Inputs:
 811     //   TString& volumeName  The volume name
 812     // Outputs:
 813     //   TString &shapeType   Shape type
 814     //   TArrayD &par         A TArrayD of parameters with all of the
 815     //                        parameters of the specified shape.
 816     // Return:
 817     //   A logical indicating whether there was an error in getting this
 818     //   information
 819    return fMCGeo->GetShape(volumePath, shapeType, par);
 820 }
 821
 822 //______________________________________________________________________
 823 Bool_t TFluka::GetMaterial(const TString &volumeName,
 824                             TString &name,Int_t &imat,
 825                             Double_t &a,Double_t &z,Double_t &dens,
 826                             Double_t &radl,Double_t &inter,TArrayD &par)
 827 {
 828     // Returns the Material and its parameters for the volume specified
 829     // by volumeName.
 830     // Note, Geant3 stores and uses mixtures as an element with an effective
 831     // Z and A. Consequently, if the parameter Z is not integer, then
 832     // this material represents some sort of mixture.
 833     // Inputs:
 834     //   TString& volumeName  The volume name
 835     // Outputs:
 836     //   TSrting   &name       Material name
 837     //   Int_t     &imat       Material index number
 838     //   Double_t  &a          Average Atomic mass of material
 839     //   Double_t  &z          Average Atomic number of material
 840     //   Double_t  &dens       Density of material [g/cm^3]
 841     //   Double_t  &radl       Average radiation length of material [cm]
 842     //   Double_t  &inter      Average interaction length of material [cm]
 843     //   TArrayD   &par        A TArrayD of user defined parameters.
 844     // Return:
 845     //   kTRUE if no errors
 846    return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
 847 }
 848
 849 //______________________________________________________________________
 850 Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
 851                          Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
 852                          Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
 853                          Double_t &deemax,Double_t &epsil, Double_t &stmin,
 854                          TArrayD &par)
 855 {
 856     // Returns the Medium and its parameters for the volume specified
 857     // by volumeName.
 858     // Inputs:
 859     //   TString& volumeName  The volume name.
 860     // Outputs:
 861     //   TString  &name       Medium name
 862     //   Int_t    &nmat       Material number defined for this medium
 863     //   Int_t    &imed       The medium index number
 864     //   Int_t    &isvol      volume number defined for this medium
 865     //   Int_t    &iflield    Magnetic field flag
 866     //   Double_t &fieldm     Magnetic field strength
 867     //   Double_t &tmaxfd     Maximum angle of deflection per step
 868     //   Double_t &stemax     Maximum step size
 869     //   Double_t &deemax     Maximum fraction of energy allowed to be lost
 870     //                        to continuous process.
 871     //   Double_t &epsil      Boundary crossing precision
 872     //   Double_t &stmin      Minimum step size allowed
 873     //   TArrayD  &par        A TArrayD of user parameters with all of the
 874     //                        parameters of the specified medium.
 875     // Return:
 876     //   kTRUE if there where no errors
 877    return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
 878 }
 879
 880 //______________________________________________________________________________
 881 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 882                          Float_t* absco, Float_t* effic, Float_t* rindex) {
 883 //
 884 // Set Cerenkov properties for medium itmed
 885 //
 886 // npckov: number of sampling points
 887 // ppckov: energy values
 888 // absco:  absorption length
 889 // effic:  quantum efficiency
 890 // rindex: refraction index
 891 //
 892 //
 893 //
 894 //  Create object holding Cerenkov properties
 895 //
 896
 897     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 898 //
 899 //  Pass object to medium
 900     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 901     medium->SetCerenkovProperties(cerenkovProperties);
 902 }
 903
 904 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 905                          Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
 906 //
 907 // Set Cerenkov properties for medium itmed
 908 //
 909 // npckov: number of sampling points
 910 // ppckov: energy values
 911 // absco:  absorption length
 912 // effic:  quantum efficiency
 913 // rindex: refraction index
 914 // rfl:    reflectivity for boundary to medium itmed
 915 //
 916 //
 917 //  Create object holding Cerenkov properties
 918 //
 919     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
 920 //
 921 //  Pass object to medium
 922     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 923     medium->SetCerenkovProperties(cerenkovProperties);
 924 }
 925
 926
 927 //______________________________________________________________________________
 928 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
 929                          Double_t *absco, Double_t *effic, Double_t *rindex) {
 930 //
 931 // Set Cerenkov properties for medium itmed
 932 //
 933 // npckov: number of sampling points
 934 // ppckov: energy values
 935 // absco:  absorption length
 936 // effic:  quantum efficiency
 937 // rindex: refraction index
 938 //
 939
 940 //
 941 //  Double_t version
 942   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
 943   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
 944   Float_t* feffic  = CreateFloatArray(effic,  npckov);
 945   Float_t* frindex = CreateFloatArray(rindex, npckov);
 946
 947   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
 948
 949   delete [] fppckov;
 950   delete [] fabsco;
 951   delete [] feffic;
 952   delete [] frindex;
 953 }
 954
 955 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
 956                          Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
 957 //
 958 // Set Cerenkov properties for medium itmed
 959 //
 960 // npckov: number of sampling points
 961 // ppckov: energy values
 962 // absco:  absorption length
 963 // effic:  quantum efficiency
 964 // rindex: refraction index
 965 // rfl:    reflectivity for boundary to medium itmed
 966 //
 967
 968 //
 969 // //  Double_t version
 970   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
 971   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
 972   Float_t* feffic  = CreateFloatArray(effic,  npckov);
 973   Float_t* frindex = CreateFloatArray(rindex, npckov);
 974   Float_t* frfl    = CreateFloatArray(rfl,    npckov);
 975
 976   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
 977
 978   delete [] fppckov;
 979   delete [] fabsco;
 980   delete [] feffic;
 981   delete [] frindex;
 982   delete [] frfl;
 983 }
 984
 985 // Euclid
 986 //______________________________________________________________________________
 987 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
 988                           Int_t /*number*/, Int_t /*nlevel*/) {
 989 //
 990 // Not with TGeo
 991    Warning("WriteEuclid", "Not implemented !");
 992 }
 993
 994
 995
 996 //_____________________________________________________________________________
 997 // methods needed by the stepping
 998 //____________________________________________________________________________
 999
1000 Int_t TFluka::GetMedium() const {
1001 //
1002 //  Get the medium number for the current fluka region
1003 //
1004     if (gGeoManager->IsOutside()) {
1005         return (-1);
1006     } else {
1007         return (fGeom->GetMedium()); // this I need to check due to remapping !!!
1008     }
1009 }
1010
1011 //____________________________________________________________________________
1012 Int_t TFluka::GetDummyRegion() const
1013 {
1014 // Returns index of the dummy region.
1015    return fGeom->GetDummyRegion();
1016 }
1017
1018 //____________________________________________________________________________
1019 Int_t TFluka::GetDummyLattice() const
1020 {
1021 // Returns index of the dummy lattice.
1022    return fGeom->GetDummyLattice();
1023 }
1024
1025 //____________________________________________________________________________
1026 // particle table usage
1027 // ID <--> PDG transformations
1028 //_____________________________________________________________________________
1029 Int_t TFluka::IdFromPDG(Int_t pdg) const
1030 {
1031
1032     //
1033     // Return Fluka code from PDG and pseudo ENDF code
1034     Int_t idSpecial[4] = {TFlukaIon::GetIonPdg(2,4),
1035                           TFlukaIon::GetIonPdg(2,3),
1036                           TFlukaIon::GetIonPdg(1,3),
1037                           TFlukaIon::GetIonPdg(1,2)};
1038     // Catch the feedback photons
1039     if (pdg == 50000051) return (kFLUKAoptical);
1040     // Ion as primary
1041     for (Int_t i = 0; i < 4; i++) {
1042         if (pdg == idSpecial[i]) return (i + kFLUKAcodemin);
1043     }
1044
1045     if ((!fUserIon && pdg == TFlukaIon::GetIonPdg(6,12)) ||
1046         ( fUserIon && pdg == fUserIon->GetPdgCode()))
1047         return (-2);
1048
1049     // MCIHAD() goes from pdg to fluka internal.
1050     Int_t intfluka = mcihad(pdg);
1051     // KPTOIP array goes from internal to official
1052     return GetFlukaKPTOIP(intfluka);
1053 }
1054
1055 //______________________________________________________________________________
1056 Int_t TFluka::PDGFromId(Int_t id) const
1057 {
1058   //
1059   // Return PDG code and pseudo ENDF code from Fluka code
1060     Int_t idSpecial[6] = {TFlukaIon::GetIonPdg(2,4), // alpha
1061                           TFlukaIon::GetIonPdg(2,3), // He3
1062                           TFlukaIon::GetIonPdg(1,3), // triton
1063                           TFlukaIon::GetIonPdg(1,2), // deuteron
1064                           TFlukaIon::GetIonPdg(0,0), // gen. ion
1065                           50000050};
1066   // IPTOKP array goes from official to internal
1067
1068     if (id == kFLUKAoptical) {
1069 // Cerenkov photon
1070 //        if (fVerbosityLevel >= 3)
1071 //            printf("\n PDGFromId: Cerenkov Photon \n");
1072         return  50000050;
1073     }
1074 // Error id
1075     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
1076         if (fVerbosityLevel >= 3)
1077             printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
1078         return -1;
1079     }
1080 // Good id
1081     if (id > 0) {
1082         Int_t intfluka = GetFlukaIPTOKP(id);
1083         if (intfluka == 0) {
1084             if (fVerbosityLevel >= 3)
1085                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
1086             return -1;
1087         } else if (intfluka < 0) {
1088             if (fVerbosityLevel >= 3)
1089                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
1090             return -1;
1091         }
1092 //        if (fVerbosityLevel >= 3)
1093 //            printf("mpdgha called with %d %d \n", id, intfluka);
1094         return mpdgha(intfluka);
1095     } else {
1096         // ions and optical photons
1097         return idSpecial[id - kFLUKAcodemin];
1098     }
1099 }
1100
1101 void TFluka::StopTrack()
1102 {
1103     // Set stopping conditions
1104     // Works for photons and charged particles
1105     fStopped = kTRUE;
1106 }
1107
1108 //_____________________________________________________________________________
1109 // methods for physics management
1110 //____________________________________________________________________________
1111 //
1112 // set methods
1113 //
1114
1115 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
1116 {
1117 //  Set process user flag for material imat
1118 //
1119 //
1120 //  Update if already in the list
1121 //
1122     TIter next(fUserConfig);
1123     TFlukaConfigOption* proc;
1124     while((proc = (TFlukaConfigOption*)next()))
1125     {
1126         if (proc->Medium() == imed) {
1127             proc->SetProcess(flagName, flagValue);
1128             return;
1129         }
1130     }
1131     proc = new TFlukaConfigOption(imed);
1132     proc->SetProcess(flagName, flagValue);
1133     fUserConfig->Add(proc);
1134 }
1135
1136 //______________________________________________________________________________
1137 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
1138 {
1139 //  Set process user flag
1140 //
1141 //
1142     SetProcess(flagName, flagValue, -1);
1143     return kTRUE;
1144 }
1145
1146 //______________________________________________________________________________
1147 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1148 {
1149 // Set user cut value for material imed
1150 //
1151     TIter next(fUserConfig);
1152     TFlukaConfigOption* proc;
1153     while((proc = (TFlukaConfigOption*)next()))
1154     {
1155         if (proc->Medium() == imed) {
1156             proc->SetCut(cutName, cutValue);
1157             return;
1158         }
1159     }
1160
1161     proc = new TFlukaConfigOption(imed);
1162     proc->SetCut(cutName, cutValue);
1163     fUserConfig->Add(proc);
1164 }
1165
1166
1167 //______________________________________________________________________________
1168 void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t imed)
1169 {
1170 // Set model parameter for material imed
1171 //
1172     TIter next(fUserConfig);
1173     TFlukaConfigOption* proc;
1174     while((proc = (TFlukaConfigOption*)next()))
1175     {
1176         if (proc->Medium() == imed) {
1177             proc->SetModelParameter(parName, parValue);
1178             return;
1179         }
1180     }
1181
1182     proc = new TFlukaConfigOption(imed);
1183     proc->SetModelParameter(parName, parValue);
1184     fUserConfig->Add(proc);
1185 }
1186
1187 //______________________________________________________________________________
1188 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1189 {
1190 // Set user cut value
1191 //
1192 //
1193     SetCut(cutName, cutValue, -1);
1194     return kTRUE;
1195 }
1196
1197
1198 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
1199 {
1200 //
1201 // Adds a user scoring option to the list
1202 //
1203     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
1204     fUserScore->Add(opt);
1205 }
1206 //______________________________________________________________________________
1207 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what,
1208                             const char* det1, const char* det2, const char* det3)
1209 {
1210 //
1211 // Adds a user scoring option to the list
1212 //
1213     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
1214     fUserScore->Add(opt);
1215 }
1216
1217 //______________________________________________________________________________
1218 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1219 {
1220   Warning("Xsec", "Not yet implemented.!\n"); return -1.;
1221 }
1222
1223
1224 //______________________________________________________________________________
1225 void TFluka::InitPhysics()
1226 {
1227 //
1228 // Physics initialisation with preparation of FLUKA input cards
1229 //
1230 // Construct file names
1231     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
1232     TString sFlukaVmcTmp = "flukaMat.inp";
1233     TString sFlukaVmcInp = GetInputFileName();
1234     TString sFlukaVmcCoreInp = GetCoreInputFileName();
1235
1236 // Open files
1237     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
1238         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
1239         exit(1);
1240     }
1241     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
1242         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
1243         exit(1);
1244     }
1245     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
1246         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
1247         exit(1);
1248     }
1249
1250 // Copy core input file
1251     Char_t sLine[255];
1252     Float_t fEventsPerRun;
1253
1254     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1255         if (strncmp(sLine,"GEOEND",6) != 0)
1256             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1257         else {
1258             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
1259             goto flukamat;
1260         }
1261     } // end of while until GEOEND card
1262
1263
1264  flukamat:
1265     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1266         fprintf(pFlukaVmcInp,"%s\n",sLine);
1267     }
1268
1269     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1270         if (strncmp(sLine,"START",5) != 0)
1271             fprintf(pFlukaVmcInp,"%s\n",sLine);
1272         else {
1273             sscanf(sLine+10,"%10f",&fEventsPerRun);
1274             goto fin;
1275         }
1276     } //end of while until START card
1277
1278  fin:
1279
1280
1281 // Pass information to configuration objects
1282
1283     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1284     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1285
1286     TIter next(fUserConfig);
1287     TFlukaConfigOption* proc;
1288     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1289 //
1290 // Process Fluka specific scoring options
1291 //
1292     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1293     Float_t loginp        = -49.0;
1294     Int_t inp             = 0;
1295     Int_t nscore          = fUserScore->GetEntries();
1296
1297     TFlukaScoringOption *mopo = 0;
1298     TFlukaScoringOption *mopi = 0;
1299
1300     for (Int_t isc = 0; isc < nscore; isc++)
1301     {
1302         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1303         char*    fileName = mopo->GetFileName();
1304         Int_t    size     = strlen(fileName);
1305         Float_t  lun      = -1.;
1306 //
1307 // Check if new output file has to be opened
1308         for (Int_t isci = 0; isci < isc; isci++) {
1309
1310
1311             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1312             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1313                 //
1314                 // No, the file already exists
1315                 lun = mopi->GetLun();
1316                 mopo->SetLun(lun);
1317                 break;
1318             }
1319         } // inner loop
1320
1321         if (lun == -1.) {
1322             // Open new output file
1323             inp++;
1324             mopo->SetLun(loginp + inp);
1325             mopo->WriteOpenFlukaFile();
1326         }
1327         mopo->WriteFlukaInputCards();
1328     }
1329
1330 // Add RANDOMIZ card
1331     fprintf(pFlukaVmcInp,"RANDOMIZ  %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
1332 // User defined ion
1333     if (fUserIon) fUserIon->WriteUserInputCard(pFlukaVmcInp);
1334 // Add START and STOP card
1335     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1336     fprintf(pFlukaVmcInp,"STOP      \n");
1337
1338
1339 // Close files
1340    fclose(pFlukaVmcCoreInp);
1341    fclose(pFlukaVmcFlukaMat);
1342    fclose(pFlukaVmcInp);
1343
1344
1345 //
1346 // Initialisation needed for Cerenkov photon production and transport
1347     TObjArray *matList = GetFlukaMaterials();
1348     Int_t nmaterial =  matList->GetEntriesFast();
1349     fMaterials = new Int_t[nmaterial+25];
1350
1351     for (Int_t im = 0; im < nmaterial; im++)
1352     {
1353         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1354         Int_t idmat = material->GetIndex();
1355         fMaterials[idmat] = im;
1356     }
1357 } // end of InitPhysics
1358
1359
1360 //______________________________________________________________________________
1361 void TFluka::SetMaxStep(Double_t step)
1362 {
1363 // Set the maximum step size
1364 //    if (step > 1.e4) return;
1365
1366 //    Int_t mreg=0, latt=0;
1367 //    fGeom->GetCurrentRegion(mreg, latt);
1368
1369
1370     Int_t mreg = fGeom->GetCurrentRegion();
1371     STEPSZ.stepmx[mreg - 1] = step;
1372 }
1373
1374
1375 Double_t TFluka::MaxStep() const
1376 {
1377 // Return the maximum for current medium
1378     Int_t mreg, latt;
1379     fGeom->GetCurrentRegion(mreg, latt);
1380     return (STEPSZ.stepmx[mreg - 1]);
1381 }
1382
1383 //______________________________________________________________________________
1384 void TFluka::SetMaxNStep(Int_t)
1385 {
1386 // SetMaxNStep is dummy procedure in TFluka !
1387   if (fVerbosityLevel >=3)
1388   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1389 }
1390
1391 //______________________________________________________________________________
1392 void TFluka::SetUserDecay(Int_t)
1393 {
1394 // SetUserDecay is dummy procedure in TFluka !
1395   if (fVerbosityLevel >=3)
1396   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1397 }
1398
1399 //
1400 // dynamic properties
1401 //
1402 //______________________________________________________________________________
1403 void TFluka::TrackPosition(TLorentzVector& position) const
1404 {
1405 // Return the current position in the master reference frame of the
1406 // track being transported
1407 // TRACKR.atrack = age of the particle
1408 // TRACKR.xtrack = x-position of the last point
1409 // TRACKR.ytrack = y-position of the last point
1410 // TRACKR.ztrack = z-position of the last point
1411   FlukaCallerCode_t caller = GetCaller();
1412   if (caller == kENDRAW    || caller == kUSDRAW ||
1413       caller == kBXExiting || caller == kBXEntering ||
1414       caller == kUSTCKV) {
1415       position.SetX(GetXsco());
1416       position.SetY(GetYsco());
1417       position.SetZ(GetZsco());
1418       position.SetT(TRACKR.atrack);
1419   }
1420   else if (caller == kMGDRAW) {
1421       Int_t i = -1;
1422       if ((i = fPrimaryElectronIndex) > -1) {
1423           // Primary Electron Ionisation
1424           Double_t x, y, z, t;
1425           GetPrimaryElectronPosition(i, x, y, z, t);
1426           position.SetX(x);
1427           position.SetY(y);
1428           position.SetZ(z);
1429           position.SetT(t);
1430       } else {
1431           position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1432           position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1433           position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1434           position.SetT(TRACKR.atrack);
1435       }
1436   }
1437   else if (caller == kSODRAW) {
1438       Int_t ist = FLKSTK.npflka;
1439       position.SetX(FLKSTK.xflk[ist]);
1440       position.SetY(FLKSTK.yflk[ist]);
1441       position.SetZ(FLKSTK.zflk[ist]);
1442       position.SetT(FLKSTK.agestk[ist]);
1443   } else if (caller == kMGResumedTrack) {
1444       position.SetX(TRACKR.spausr[0]);
1445       position.SetY(TRACKR.spausr[1]);
1446       position.SetZ(TRACKR.spausr[2]);
1447       position.SetT(TRACKR.spausr[3]);
1448   }
1449   else
1450       Warning("TrackPosition","position not available");
1451 }
1452
1453 //______________________________________________________________________________
1454 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1455 {
1456 // Return the current position in the master reference frame of the
1457 // track being transported
1458 // TRACKR.atrack = age of the particle
1459 // TRACKR.xtrack = x-position of the last point
1460 // TRACKR.ytrack = y-position of the last point
1461 // TRACKR.ztrack = z-position of the last point
1462   FlukaCallerCode_t caller = GetCaller();
1463   if (caller == kENDRAW    || caller == kUSDRAW ||
1464       caller == kBXExiting || caller == kBXEntering ||
1465       caller == kUSTCKV) {
1466       x = GetXsco();
1467       y = GetYsco();
1468       z = GetZsco();
1469   }
1470   else if (caller == kMGDRAW) {
1471       Int_t i = -1;
1472       if ((i = fPrimaryElectronIndex) > -1) {
1473           Double_t t;
1474           GetPrimaryElectronPosition(i, x, y, z, t);
1475       } else {
1476           x = TRACKR.xtrack[TRACKR.ntrack];
1477           y = TRACKR.ytrack[TRACKR.ntrack];
1478           z = TRACKR.ztrack[TRACKR.ntrack];
1479       }
1480   }
1481   else if (caller == kSODRAW) {
1482       Int_t ist = FLKSTK.npflka;
1483       x = FLKSTK.xflk[ist];
1484       y = FLKSTK.yflk[ist];
1485       z = FLKSTK.zflk[ist];
1486   }
1487   else if (caller == kMGResumedTrack) {
1488       x = TRACKR.spausr[0];
1489       y = TRACKR.spausr[1];
1490       z = TRACKR.spausr[2];
1491   }
1492   else
1493       Warning("TrackPosition","position not available");
1494 }
1495
1496 //______________________________________________________________________________
1497 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1498 {
1499 // Return the direction and the momentum (GeV/c) of the track
1500 // currently being transported
1501 // TRACKR.ptrack = momentum of the particle (not always defined, if
1502 //               < 0 must be obtained from etrack)
1503 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1504 // TRACKR.etrack = total energy of the particle
1505 // TRACKR.jtrack = identity number of the particle
1506 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1507   FlukaCallerCode_t  caller = GetCaller();
1508   FlukaProcessCode_t icode  = GetIcode();
1509
1510   if (caller  != kEEDRAW         &&
1511       caller  != kMGResumedTrack &&
1512       caller  != kSODRAW         &&
1513       caller  != kUSDRAW         &&
1514       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1515       if (TRACKR.ptrack >= 0) {
1516           momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1517           momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1518           momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1519           momentum.SetE(TRACKR.etrack);
1520           return;
1521       }
1522       else {
1523           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1524           momentum.SetPx(p*TRACKR.cxtrck);
1525           momentum.SetPy(p*TRACKR.cytrck);
1526           momentum.SetPz(p*TRACKR.cztrck);
1527           momentum.SetE(TRACKR.etrack);
1528           return;
1529       }
1530   } else if  (caller == kMGResumedTrack) {
1531       momentum.SetPx(TRACKR.spausr[4]);
1532       momentum.SetPy(TRACKR.spausr[5]);
1533       momentum.SetPz(TRACKR.spausr[6]);
1534       momentum.SetE (TRACKR.spausr[7]);
1535       return;
1536   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1537       momentum.SetPx(0.);
1538       momentum.SetPy(0.);
1539       momentum.SetPz(0.);
1540       momentum.SetE(TrackMass());
1541
1542   } else if (caller == kSODRAW) {
1543       Int_t ist  = FLKSTK.npflka;
1544       Double_t p = FLKSTK.pmoflk[ist];
1545       Int_t ifl  = FLKSTK.iloflk[ist];
1546       Double_t m = PAPROP.am[ifl + 6];
1547       Double_t e = TMath::Sqrt(p * p + m * m);
1548       momentum.SetPx(p * FLKSTK.txflk[ist]);
1549       momentum.SetPy(p * FLKSTK.tyflk[ist]);
1550       momentum.SetPz(p * FLKSTK.tzflk[ist]);
1551       momentum.SetE(e);
1552   } else if (caller == kUSDRAW) {
1553       if (icode == kEMFSCObrems  ||
1554           icode == kEMFSCOmoller ||
1555           icode == kEMFSCObhabha ||
1556           icode == kEMFSCOcompton )
1557       {
1558           momentum.SetPx(fPint[0]);
1559           momentum.SetPy(fPint[1]);
1560           momentum.SetPz(fPint[2]);
1561           momentum.SetE(fPint[3]);
1562       } else if (icode == kKASKADdray  ||
1563                  icode == kKASKADbrems ||
1564                  icode == kKASKADpair) {
1565           momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
1566           momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
1567           momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
1568           momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
1569       } else {
1570           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack
1571                             - ParticleMassFPC(TRACKR.jtrack)
1572                             * ParticleMassFPC(TRACKR.jtrack));
1573           momentum.SetPx(p*TRACKR.cxtrck);
1574           momentum.SetPy(p*TRACKR.cytrck);
1575           momentum.SetPz(p*TRACKR.cztrck);
1576           momentum.SetE(TRACKR.etrack);
1577       }
1578   }
1579   else
1580     Warning("TrackMomentum","momentum not available");
1581 }
1582
1583 //______________________________________________________________________________
1584 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1585 {
1586 // Return the direction and the momentum (GeV/c) of the track
1587 // currently being transported
1588 // TRACKR.ptrack = momentum of the particle (not always defined, if
1589 //               < 0 must be obtained from etrack)
1590 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1591 // TRACKR.etrack = total energy of the particle
1592 // TRACKR.jtrack = identity number of the particle
1593 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1594   FlukaCallerCode_t   caller = GetCaller();
1595   FlukaProcessCode_t  icode  = GetIcode();
1596   if (caller != kEEDRAW         &&
1597       caller != kMGResumedTrack &&
1598       caller != kSODRAW         &&
1599       caller != kUSDRAW         &&
1600       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1601     if (TRACKR.ptrack >= 0) {
1602       px = TRACKR.ptrack*TRACKR.cxtrck;
1603       py = TRACKR.ptrack*TRACKR.cytrck;
1604       pz = TRACKR.ptrack*TRACKR.cztrck;
1605       e  = TRACKR.etrack;
1606       return;
1607     }
1608     else {
1609       Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1610       px = p*TRACKR.cxtrck;
1611       py = p*TRACKR.cytrck;
1612       pz = p*TRACKR.cztrck;
1613       e  = TRACKR.etrack;
1614       return;
1615     }
1616   } else if (caller == kMGResumedTrack) {
1617       px = TRACKR.spausr[4];
1618       py = TRACKR.spausr[5];
1619       pz = TRACKR.spausr[6];
1620       e  = TRACKR.spausr[7];
1621       return;
1622   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1623       px = 0.;
1624       py = 0.;
1625       pz = 0.;
1626       e  = TrackMass();
1627   } else if (caller == kSODRAW) {
1628       Int_t ist  = FLKSTK.npflka;
1629       Double_t p = FLKSTK.pmoflk[ist];
1630       Int_t ifl  = FLKSTK.iloflk[ist];
1631       Double_t m = PAPROP.am[ifl + 6];
1632                e = TMath::Sqrt(p * p + m * m);
1633       px = p * FLKSTK.txflk[ist];
1634       py = p * FLKSTK.tyflk[ist];
1635       pz = p * FLKSTK.tzflk[ist];
1636   } else if (caller == kUSDRAW) {
1637       if (icode == kEMFSCObrems  ||
1638           icode == kEMFSCOmoller ||
1639           icode == kEMFSCObhabha ||
1640           icode == kEMFSCOcompton )
1641       {
1642           px = fPint[0];
1643           py = fPint[1];
1644           pz = fPint[2];
1645           e  = fPint[3];
1646       } else if (icode == kKASKADdray  ||
1647                  icode == kKASKADbrems ||
1648                  icode == kKASKADpair) {
1649           px = GENSTK.plr[0] * GENSTK.cxr[0];
1650           py = GENSTK.plr[0] * GENSTK.cyr[0];
1651           pz = GENSTK.plr[0] * GENSTK.czr[0];
1652           e  = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
1653       } else {
1654           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1655           px = p*TRACKR.cxtrck;
1656           py = p*TRACKR.cytrck;
1657           pz = p*TRACKR.cztrck;
1658           e  = TRACKR.etrack;
1659       }
1660   }
1661   else
1662       Warning("TrackMomentum","momentum not available");
1663 }
1664
1665 //______________________________________________________________________________
1666 Double_t TFluka::TrackStep() const
1667 {
1668 // Return the length in centimeters of the current step
1669 // TRACKR.ctrack = total curved path
1670     FlukaCallerCode_t caller = GetCaller();
1671     if (caller == kMGDRAW) {
1672         Int_t i;
1673         if ((i = fPrimaryElectronIndex) > -1) {
1674             if (i > 0) {
1675                 return (fPIlength[i] - fPIlength[i-1]);
1676             } else {
1677                 Double_t s (TRACKR.ctrack - (fPIlength[fNPI - 1] - fPIlength[0]));
1678                 return s;
1679             }
1680         } else {
1681             return TRACKR.ctrack;
1682         }
1683     } else if (caller == kBXEntering || caller == kBXExiting ||
1684                caller == kENDRAW     || caller == kUSDRAW ||
1685                caller == kUSTCKV     || caller == kMGResumedTrack ||
1686                caller == kSODRAW)
1687     {
1688         return 0.0;
1689     } else {
1690         Warning("TrackStep", "track step not available");
1691         return 0.0;
1692     }
1693 }
1694
1695 //______________________________________________________________________________
1696 Double_t TFluka::TrackLength() const
1697 {
1698 // TRACKR.cmtrck = cumulative curved path since particle birth
1699   FlukaCallerCode_t caller = GetCaller();
1700   if (caller == kMGDRAW) {
1701       Int_t i;
1702       if ((i = fPrimaryElectronIndex) > -1) {
1703           return fPIlength[i];
1704       } else {
1705           return TRACKR.cmtrck;
1706       }
1707
1708   } else if (caller == kBXEntering || caller == kBXExiting ||
1709              caller == kENDRAW || caller == kUSDRAW || caller == kUSTCKV)
1710       return TRACKR.cmtrck;
1711   else if (caller == kMGResumedTrack)
1712       return TRACKR.spausr[8];
1713   else if (caller == kSODRAW)
1714       return 0.0;
1715   else {
1716       Warning("TrackLength", "track length not available for caller %5d \n", caller);
1717       return 0.0;
1718   }
1719 }
1720
1721
1722 //______________________________________________________________________________
1723 Double_t TFluka::TrackTime() const
1724 {
1725 // Return the current time of flight of the track being transported
1726 // TRACKR.atrack = age of the particle
1727   FlukaCallerCode_t caller = GetCaller();
1728   if (caller == kMGDRAW) {
1729       Int_t i;
1730       if ((i = fPrimaryElectronIndex) > -1) {
1731           Double_t t = fPItime[i];
1732           return t;
1733       } else {
1734           return TRACKR.atrack;
1735       }
1736   } else if (caller == kBXEntering || caller == kBXExiting ||
1737              caller == kENDRAW     || caller == kUSDRAW    ||
1738              caller == kUSTCKV)
1739     return TRACKR.atrack;
1740   else if (caller == kMGResumedTrack)
1741     return TRACKR.spausr[3];
1742   else if (caller == kSODRAW) {
1743       return (FLKSTK.agestk[FLKSTK.npflka]);
1744   }
1745   else {
1746     Warning("TrackTime", "track time not available");
1747     return 0.0;
1748   }
1749 }
1750
1751 //______________________________________________________________________________
1752 Double_t TFluka::Edep() const
1753 {
1754 // Energy deposition
1755 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1756 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1757 //    but in the variable "rull" of the procedure "endraw.cxx"
1758 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1759 // -->no energy loss along the track
1760 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1761 // -->energy loss distributed along the track
1762 // TRACKR.dtrack = energy deposition of the jth deposition event
1763
1764   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1765   // If coming from usdraw we just signal particle production - no edep
1766   // If just first time after resuming, no edep for the primary
1767   FlukaCallerCode_t caller = GetCaller();
1768
1769   if (caller == kBXExiting || caller == kBXEntering ||
1770       caller == kUSDRAW    || caller == kMGResumedTrack ||
1771       caller == kSODRAW)
1772       return 0.0;
1773   Double_t sum = 0;
1774   Int_t i = -1;
1775
1776   // Material with primary ionisation activated but number of primary electrons nprim = 0
1777   if (fPrimaryElectronIndex == -2) return 0.0;
1778   // nprim > 0
1779   if ((i = fPrimaryElectronIndex) > -1) {
1780       // Primary ionisation
1781       sum = GetPrimaryElectronKineticEnergy(i);
1782       if (sum > 100.) {
1783           printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
1784       }
1785       return sum;
1786   } else {
1787       // Normal ionisation
1788       if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
1789
1790       for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1791           sum +=TRACKR.dtrack[j];
1792       }
1793       if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1794           return fRull + sum;
1795       else {
1796           return sum;
1797       }
1798   }
1799 }
1800
1801 //______________________________________________________________________________
1802 Int_t TFluka::CorrectFlukaId() const
1803 {
1804    // since we don't put photons and e- created bellow transport cut on the vmc stack
1805    // and there is a call to endraw for energy deposition for each of them
1806    // and they have the track number of their parent, but different identity (pdg)
1807    // so we want to assign also their parent identity.
1808
1809    if( (IsTrackStop())
1810         && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
1811         && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
1812       if (fVerbosityLevel >=3)
1813          cout << "CorrectFlukaId() for icode=" << GetIcode()
1814                << " track=" << TRACKR.ispusr[mkbmx2 - 1]
1815                << " current PDG=" << PDGFromId(TRACKR.jtrack)
1816                << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
1817       return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
1818    }
1819    if (TRACKR.jtrack <= 64){
1820        return TRACKR.jtrack;
1821    } else {
1822        return TRACKR.j0trck;
1823    }
1824 }
1825
1826
1827 //______________________________________________________________________________
1828 Int_t TFluka::TrackPid() const
1829 {
1830 // Return the id of the particle transported
1831 // TRACKR.jtrack = identity number of the particle
1832   FlukaCallerCode_t caller = GetCaller();
1833   if (caller != kEEDRAW && caller != kSODRAW) {
1834      return PDGFromId( CorrectFlukaId() );
1835   }
1836   else if (caller == kSODRAW) {
1837       return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1838   }
1839   else
1840     return -1000;
1841 }
1842
1843 //______________________________________________________________________________
1844 Double_t TFluka::TrackCharge() const
1845 {
1846 // Return charge of the track currently transported
1847 // PAPROP.ichrge = electric charge of the particle
1848 // TRACKR.jtrack = identity number of the particle
1849
1850   FlukaCallerCode_t caller = GetCaller();
1851   if (caller != kEEDRAW && caller != kSODRAW)
1852      return PAPROP.ichrge[CorrectFlukaId() + 6];
1853   else if (caller == kSODRAW) {
1854       Int_t ifl =  PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1855       return PAPROP.ichrge[ifl + 6];
1856   }
1857   else
1858     return -1000.0;
1859 }
1860
1861 //______________________________________________________________________________
1862 Double_t TFluka::TrackMass() const
1863 {
1864 // PAPROP.am = particle mass in GeV
1865 // TRACKR.jtrack = identity number of the particle
1866   FlukaCallerCode_t caller = GetCaller();
1867   if (caller != kEEDRAW && caller != kSODRAW)
1868      return PAPROP.am[CorrectFlukaId()+6];
1869   else if (caller == kSODRAW) {
1870       Int_t ifl =  FLKSTK.iloflk[FLKSTK.npflka];
1871       return PAPROP.am[ifl + 6];
1872   }
1873   else
1874     return -1000.0;
1875 }
1876
1877 //______________________________________________________________________________
1878 Double_t TFluka::Etot() const
1879 {
1880 // TRACKR.etrack = total energy of the particle
1881   FlukaCallerCode_t  caller = GetCaller();
1882   FlukaProcessCode_t icode  = GetIcode();
1883   if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
1884   {
1885       return TRACKR.etrack;
1886   } else if (caller == kUSDRAW) {
1887       if (icode == kEMFSCObrems  ||
1888           icode == kEMFSCOmoller ||
1889           icode == kEMFSCObhabha ||
1890           icode == kEMFSCOcompton ) {
1891           return  fPint[3];
1892       }
1893       else if (icode == kKASKADdray  ||
1894                icode == kKASKADbrems ||
1895                icode == kKASKADpair) {
1896           return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
1897       } else {
1898           return TRACKR.etrack;
1899       }
1900
1901   }
1902   else if (caller == kSODRAW) {
1903       Int_t ist  = FLKSTK.npflka;
1904       Double_t p = FLKSTK.pmoflk[ist];
1905       Int_t ifl  = FLKSTK.iloflk[ist];
1906       Double_t m = PAPROP.am[ifl + 6];
1907       Double_t e = TMath::Sqrt(p * p + m * m);
1908       return e;
1909   }
1910   printf("Etot %5d %5d \n", caller, icode);
1911
1912   return -1000.0;
1913 }
1914
1915 //
1916 // track status
1917 //
1918 //______________________________________________________________________________
1919 Bool_t   TFluka::IsNewTrack() const
1920 {
1921 // Return true for the first call of Stepping()
1922    return fTrackIsNew;
1923 }
1924
1925 void     TFluka::SetTrackIsNew(Bool_t flag)
1926 {
1927 // Return true for the first call of Stepping()
1928    fTrackIsNew = flag;
1929
1930 }
1931
1932
1933 //______________________________________________________________________________
1934 Bool_t   TFluka::IsTrackInside() const
1935 {
1936 // True if the track is not at the boundary of the current volume
1937 // In Fluka a step is always inside one kind of material
1938 // If the step would go behind the region of one material,
1939 // it will be shortened to reach only the boundary.
1940 // Therefore IsTrackInside() is always true.
1941   FlukaCallerCode_t caller = GetCaller();
1942   if (caller == kBXEntering || caller == kBXExiting)
1943     return 0;
1944   else
1945     return 1;
1946 }
1947
1948 //______________________________________________________________________________
1949 Bool_t   TFluka::IsTrackEntering() const
1950 {
1951 // True if this is the first step of the track in the current volume
1952
1953   FlukaCallerCode_t caller = GetCaller();
1954   if (caller == kBXEntering)
1955     return 1;
1956   else return 0;
1957 }
1958
1959 //______________________________________________________________________________
1960 Bool_t   TFluka::IsTrackExiting() const
1961 {
1962 // True if track is exiting volume
1963 //
1964   FlukaCallerCode_t caller = GetCaller();
1965   if (caller == kBXExiting)
1966     return 1;
1967   else return 0;
1968 }
1969
1970 //______________________________________________________________________________
1971 Bool_t   TFluka::IsTrackOut() const
1972 {
1973 // True if the track is out of the setup
1974 // means escape
1975   FlukaProcessCode_t icode = GetIcode();
1976
1977   if (icode == kKASKADescape ||
1978       icode == kEMFSCOescape ||
1979       icode == kKASNEUescape ||
1980       icode == kKASHEAescape ||
1981       icode == kKASOPHescape)
1982        return 1;
1983   else return 0;
1984 }
1985
1986 //______________________________________________________________________________
1987 Bool_t   TFluka::IsTrackDisappeared() const
1988 {
1989 // All inelastic interactions and decays
1990 // fIcode from usdraw
1991   FlukaProcessCode_t icode = GetIcode();
1992   if (icode == kKASKADinelint    || // inelastic interaction
1993       icode == kKASKADdecay      || // particle decay
1994       icode == kKASKADdray       || // delta ray generation by hadron
1995       icode == kKASKADpair       || // direct pair production
1996       icode == kKASKADbrems      || // bremsstrahlung (muon)
1997       icode == kEMFSCObrems      || // bremsstrahlung (electron)
1998       icode == kEMFSCOmoller     || // Moller scattering
1999       icode == kEMFSCObhabha     || // Bhaba scattering
2000       icode == kEMFSCOanniflight || // in-flight annihilation
2001       icode == kEMFSCOannirest   || // annihilation at rest
2002       icode == kEMFSCOpair       || // pair production
2003       icode == kEMFSCOcompton    || // Compton scattering
2004       icode == kEMFSCOphotoel    || // Photoelectric effect
2005       icode == kKASNEUhadronic   || // hadronic interaction
2006       icode == kKASHEAdray          // delta-ray
2007       ) return 1;
2008   else return 0;
2009 }
2010
2011 //______________________________________________________________________________
2012 Bool_t   TFluka::IsTrackStop() const
2013 {
2014 // True if the track energy has fallen below the threshold
2015 // means stopped by signal or below energy threshold
2016   FlukaProcessCode_t icode = GetIcode();
2017   if (icode == kKASKADstopping  || // stopping particle
2018       icode == kKASKADtimekill  || // time kill
2019       icode == kEMFSCOstopping1 || // below user-defined cut-off
2020       icode == kEMFSCOstopping2 || // below user cut-off
2021       icode == kEMFSCOtimekill  || // time kill
2022       icode == kKASNEUstopping  || // neutron below threshold
2023       icode == kKASNEUtimekill  || // time kill
2024       icode == kKASHEAtimekill  || // time kill
2025       icode == kKASOPHtimekill) return 1; // time kill
2026   else return 0;
2027 }
2028
2029 //______________________________________________________________________________
2030 Bool_t   TFluka::IsTrackAlive() const
2031 {
2032 // Means not disappeared or not out
2033     FlukaProcessCode_t icode = GetIcode();
2034
2035     if (IsTrackOut()               ||
2036         IsTrackStop()              ||
2037         icode == kKASKADinelint    || // inelastic interaction
2038         icode == kKASKADdecay      || // particle decay
2039         icode == kEMFSCOanniflight || // in-flight annihilation
2040         icode == kEMFSCOannirest   || // annihilation at rest
2041         icode == kEMFSCOpair       || // pair production
2042         icode == kEMFSCOphotoel    || // Photoelectric effect
2043         icode == kKASNEUhadronic      // hadronic interaction
2044         )
2045     {
2046         // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
2047         return 0;
2048     } else {
2049         return 1;
2050     }
2051 }
2052
2053 //
2054 // secondaries
2055 //
2056
2057 //______________________________________________________________________________
2058 Int_t TFluka::NSecondaries() const
2059
2060 {
2061 // Number of secondary particles generated in the current step
2062 // GENSTK.np = number of secondaries except light and heavy ions
2063 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2064     FlukaCallerCode_t caller = GetCaller();
2065     if (caller == kUSDRAW)  // valid only after usdraw
2066         return GENSTK.np + FHEAVY.npheav;
2067     else if (caller == kUSTCKV) {
2068         // Cerenkov Photon production
2069         return fNCerenkov;
2070     }
2071     return 0;
2072 } // end of NSecondaries
2073
2074 //______________________________________________________________________________
2075 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2076                 TLorentzVector& position, TLorentzVector& momentum)
2077 {
2078 // Copy particles from secondary stack to vmc stack
2079 //
2080
2081     FlukaCallerCode_t caller = GetCaller();
2082     if (caller == kUSDRAW) {  // valid only after usdraw
2083         if (GENSTK.np > 0) {
2084             // Hadronic interaction
2085             if (isec >= 0 && isec < GENSTK.np) {
2086                 particleId = PDGFromId(GENSTK.kpart[isec]);
2087                 position.SetX(fXsco);
2088                 position.SetY(fYsco);
2089                 position.SetZ(fZsco);
2090                 position.SetT(TRACKR.atrack);
2091                 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
2092                 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
2093                 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
2094                 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
2095             }
2096             else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
2097                 Int_t jsec = isec - GENSTK.np;
2098                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2099                 position.SetX(fXsco);
2100                 position.SetY(fYsco);
2101                 position.SetZ(fZsco);
2102                 position.SetT(TRACKR.atrack);
2103                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2104                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2105                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2106                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2107                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2108                 else if (FHEAVY.tkheav[jsec] > 6)
2109                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2110             }
2111             else
2112                 Warning("GetSecondary","isec out of range");
2113         }
2114     } else if (caller == kUSTCKV) {
2115         Int_t index = OPPHST.lstopp - isec;
2116         position.SetX(OPPHST.xoptph[index]);
2117         position.SetY(OPPHST.yoptph[index]);
2118         position.SetZ(OPPHST.zoptph[index]);
2119         position.SetT(OPPHST.agopph[index]);
2120         Double_t p = OPPHST.poptph[index];
2121
2122         momentum.SetPx(p * OPPHST.txopph[index]);
2123         momentum.SetPy(p * OPPHST.tyopph[index]);
2124         momentum.SetPz(p * OPPHST.tzopph[index]);
2125         momentum.SetE(p);
2126     }
2127     else
2128         Warning("GetSecondary","no secondaries available");
2129
2130 } // end of GetSecondary
2131
2132
2133 //______________________________________________________________________________
2134 TMCProcess TFluka::ProdProcess(Int_t) const
2135
2136 {
2137 // Name of the process that has produced the secondary particles
2138 // in the current step
2139
2140     Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton ||
2141                      TRACKR.jtrack == kFLUKAmuplus ||
2142                      TRACKR.jtrack == kFLUKAmuminus);
2143     FlukaProcessCode_t icode = GetIcode();
2144
2145     if      (icode == kKASKADdecay)                                   return kPDecay;
2146     else if (icode == kKASKADpair || icode == kEMFSCOpair)            return kPPair;
2147     else if (icode == kEMFSCOcompton)                                 return kPCompton;
2148     else if (icode == kEMFSCOphotoel)                                 return kPPhotoelectric;
2149     else if (icode == kKASKADbrems      || icode == kEMFSCObrems)     return kPBrem;
2150     else if (icode == kKASKADdray       || icode == kKASHEAdray)      return kPDeltaRay;
2151     else if (icode == kEMFSCOmoller     || icode == kEMFSCObhabha)    return kPDeltaRay;
2152     else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest)  return kPAnnihilation;
2153     else if (icode == kKASKADinelint) {
2154         if (!mugamma)                                                 return kPHadronic;
2155         else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
2156         else                                                          return kPMuonNuclear;
2157     }
2158     else if (icode == kEMFSCOrayleigh)                                return kPRayleigh;
2159 // Fluka codes 100, 300 and 400 still to be investigasted
2160     else                                                              return kPNoProcess;
2161 }
2162
2163
2164 Int_t TFluka::StepProcesses(TArrayI &proc) const
2165 {
2166   //
2167   // Return processes active in the current step
2168   //
2169     FlukaProcessCode_t icode   = GetIcode();
2170     FlukaCallerCode_t  caller  = GetCaller();
2171     proc.Set(1);
2172     TMCProcess iproc;
2173     if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
2174         iproc = kPTransportation;
2175     }
2176     else if (caller == kUSTCKV) {
2177         iproc = kPCerenkov;
2178     } else {
2179         switch (icode) {
2180         case kEMFSCO:
2181             if (Edep() > 0.) {
2182                 iproc = kPEnergyLoss;
2183             } else {
2184                 iproc = kPTransportation;
2185             }
2186             break;
2187         case kKASKAD:
2188             if (Edep() > 0.) {
2189                 iproc = kPEnergyLoss;
2190             } else {
2191                 iproc = kPTransportation;
2192             }
2193             break;
2194         case kKASHEA:
2195         case kKASNEU:
2196         case kKASOPH:
2197         case kKASKADescape:
2198         case kEMFSCOescape:
2199         case kKASNEUescape:
2200         case kKASHEAescape:
2201         case kKASOPHescape:
2202             iproc = kPTransportation;
2203             break;
2204         case kKASKADtimekill:
2205         case kEMFSCOtimekill:
2206         case kKASNEUtimekill:
2207         case kKASHEAtimekill:
2208         case kKASOPHtimekill:
2209             iproc =  kPTOFlimit;
2210             break;
2211         case kKASKADstopping:
2212         case kEMFSCOstopping1:
2213         case kEMFSCOstopping2:
2214         case kKASNEUstopping:
2215             iproc = kPStop;
2216             break;
2217         case kKASKADinelint:
2218         case kKASNEUhadronic:
2219             iproc = kPHadronic;
2220             break;
2221         case kKASKADinelarecoil:
2222             iproc = kPHadronic;
2223             break;
2224         case kKASKADnelint:
2225             iproc = kPHElastic;
2226             break;
2227         case kKASOPHabsorption:
2228             iproc = kPLightAbsorption;
2229             break;
2230         case kKASOPHrefraction:
2231             iproc = kPLightRefraction;
2232             break;
2233         case kEMFSCOlocaldep :
2234             iproc = kPPhotoelectric;
2235             break;
2236         default:
2237             iproc = ProdProcess(0);
2238         }
2239     }
2240
2241     proc[0] = iproc;
2242     return 1;
2243 }
2244 //______________________________________________________________________________
2245 Int_t TFluka::VolId2Mate(Int_t id) const
2246 {
2247 //
2248 // Returns the material number for a given volume ID
2249 //
2250    return fMCGeo->VolId2Mate(id);
2251 }
2252
2253 //______________________________________________________________________________
2254 const char* TFluka::VolName(Int_t id) const
2255 {
2256 //
2257 // Returns the volume name for a given volume ID
2258 //
2259    return fMCGeo->VolName(id);
2260 }
2261
2262 Int_t TFluka::MediumId(const Text_t* mediumName) const
2263 {
2264     //
2265     // Return the unique medium id for medium with name mediumName
2266     TList *medlist = gGeoManager->GetListOfMedia();
2267     TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
2268     if (med) {
2269         return (med->GetId());
2270     } else {
2271         return (-1);
2272     }
2273 }
2274
2275 //______________________________________________________________________________
2276 Int_t TFluka::VolId(const Text_t* volName) const
2277 {
2278 //
2279 // Converts from volume name to volume ID.
2280 // Time consuming. (Only used during set-up)
2281 // Could be replaced by hash-table
2282 //
2283     char sname[20];
2284     Int_t len;
2285     strncpy(sname, volName, len = strlen(volName));
2286     sname[len] = 0;
2287     while (sname[len - 1] == ' ') sname[--len] = 0;
2288     return fMCGeo->VolId(sname);
2289 }
2290
2291 //______________________________________________________________________________
2292 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2293 {
2294 //
2295 // Return the logical id and copy number corresponding to the current fluka region
2296 //
2297   if (gGeoManager->IsOutside()) return 0;
2298   TGeoNode *node = gGeoManager->GetCurrentNode();
2299   copyNo = node->GetNumber();
2300   Int_t id = node->GetVolume()->GetNumber();
2301   return id;
2302 }
2303
2304 //______________________________________________________________________________
2305 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2306 {
2307 //
2308 // Return the logical id and copy number of off'th mother
2309 // corresponding to the current fluka region
2310 //
2311   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2312   if (off==0) return CurrentVolID(copyNo);
2313   TGeoNode *node = gGeoManager->GetMother(off);
2314   if (!node) return 0;
2315   copyNo = node->GetNumber();
2316   return node->GetVolume()->GetNumber();
2317 }
2318
2319 //______________________________________________________________________________
2320 const char* TFluka::CurrentVolName() const
2321 {
2322 //
2323 // Return the current volume name
2324 //
2325   if (gGeoManager->IsOutside()) return "OutOfWorld";
2326   return gGeoManager->GetCurrentVolume()->GetName();
2327 }
2328
2329 //______________________________________________________________________________
2330 const char* TFluka::CurrentVolOffName(Int_t off) const
2331 {
2332 //
2333 // Return the volume name of the off'th mother of the current volume
2334 //
2335   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2336   if (off==0) return CurrentVolName();
2337   TGeoNode *node = gGeoManager->GetMother(off);
2338   if (!node) return 0;
2339   return node->GetVolume()->GetName();
2340 }
2341
2342 const char* TFluka::CurrentVolPath() {
2343   // Return the current volume path
2344   return gGeoManager->GetPath();
2345 }
2346 //______________________________________________________________________________
2347 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
2348                       Float_t & dens, Float_t & radl, Float_t & absl) const
2349 {
2350 //
2351 //  Return the current medium number and material properties
2352 //
2353   Int_t copy;
2354   Int_t id  =  TFluka::CurrentVolID(copy);
2355   Int_t med =  TFluka::VolId2Mate(id);
2356   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
2357   TGeoMaterial*   mat = vol->GetMaterial();
2358   a    = mat->GetA();
2359   z    = mat->GetZ();
2360   dens = mat->GetDensity();
2361   radl = mat->GetRadLen();
2362   absl = mat->GetIntLen();
2363
2364   return med;
2365 }
2366
2367 //______________________________________________________________________________
2368 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2369 {
2370 // Transforms a position from the world reference frame
2371 // to the current volume reference frame.
2372 //
2373 //  Geant3 desription:
2374 //  ==================
2375 //       Computes coordinates XD (in DRS)
2376 //       from known coordinates XM in MRS
2377 //       The local reference system can be initialized by
2378 //         - the tracking routines and GMTOD used in GUSTEP
2379 //         - a call to GMEDIA(XM,NUMED)
2380 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2381 //             (inverse routine is GDTOM)
2382 //
2383 //        If IFLAG=1  convert coordinates
2384 //           IFLAG=2  convert direction cosinus
2385 //
2386 // ---
2387    Double_t xmL[3], xdL[3];
2388    Int_t i;
2389    for (i=0;i<3;i++) xmL[i]=xm[i];
2390    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2391    else            gGeoManager->MasterToLocalVect(xmL,xdL);
2392    for (i=0;i<3;i++) xd[i] = xdL[i];
2393 }
2394
2395 //______________________________________________________________________________
2396 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2397 {
2398 //
2399 // See Gmtod(Float_t*, Float_t*, Int_t)
2400 //
2401    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2402    else            gGeoManager->MasterToLocalVect(xm,xd);
2403 }
2404
2405 //______________________________________________________________________________
2406 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2407 {
2408 // Transforms a position from the current volume reference frame
2409 // to the world reference frame.
2410 //
2411 //  Geant3 desription:
2412 //  ==================
2413 //  Computes coordinates XM (Master Reference System
2414 //  knowing the coordinates XD (Detector Ref System)
2415 //  The local reference system can be initialized by
2416 //    - the tracking routines and GDTOM used in GUSTEP
2417 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2418 //        (inverse routine is GMTOD)
2419 //
2420 //   If IFLAG=1  convert coordinates
2421 //      IFLAG=2  convert direction cosinus
2422 //
2423 // ---
2424    Double_t xmL[3], xdL[3];
2425    Int_t i;
2426    for (i=0;i<3;i++) xdL[i] = xd[i];
2427    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2428    else            gGeoManager->LocalToMasterVect(xdL,xmL);
2429    for (i=0;i<3;i++) xm[i]=xmL[i];
2430 }
2431
2432 //______________________________________________________________________________
2433 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2434 {
2435 //
2436 // See Gdtom(Float_t*, Float_t*, Int_t)
2437 //
2438    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2439    else            gGeoManager->LocalToMasterVect(xd,xm);
2440 }
2441
2442 //______________________________________________________________________________
2443 TObjArray *TFluka::GetFlukaMaterials()
2444 {
2445 //
2446 // Get array of Fluka materials
2447    return fGeom->GetMatList();
2448 }
2449
2450 //______________________________________________________________________________
2451 void TFluka::SetMreg(Int_t l, Int_t lttc)
2452 {
2453 // Set current fluka region
2454    fCurrentFlukaRegion = l;
2455    fGeom->SetMreg(l,lttc);
2456 }
2457
2458
2459
2460
2461 //______________________________________________________________________________
2462 TString TFluka::ParticleName(Int_t pdg) const
2463 {
2464     // Return particle name for particle with pdg code pdg.
2465     Int_t ifluka = IdFromPDG(pdg);
2466     return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
2467 }
2468
2469
2470 //______________________________________________________________________________
2471 Double_t TFluka::ParticleMass(Int_t pdg) const
2472 {
2473     // Return particle mass for particle with pdg code pdg.
2474     Int_t ifluka = IdFromPDG(pdg);
2475     return (PAPROP.am[ifluka - kFLUKAcodemin]);
2476 }
2477
2478 //______________________________________________________________________________
2479 Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2480 {
2481     // Return particle mass for particle with Fluka particle code fpc
2482     return (PAPROP.am[fpc - kFLUKAcodemin]);
2483 }
2484
2485 //______________________________________________________________________________
2486 Double_t TFluka::ParticleCharge(Int_t pdg) const
2487 {
2488     // Return particle charge for particle with pdg code pdg.
2489     Int_t ifluka = IdFromPDG(pdg);
2490     return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
2491 }
2492
2493 //______________________________________________________________________________
2494 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2495 {
2496     // Return particle lifetime for particle with pdg code pdg.
2497     Int_t ifluka = IdFromPDG(pdg);
2498     return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
2499 }
2500
2501 //______________________________________________________________________________
2502 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2503 {
2504     // Retrieve particle properties for particle with pdg code pdg.
2505
2506     strcpy(name, ParticleName(pdg).Data());
2507     type   = ParticleMCType(pdg);
2508     mass   = ParticleMass(pdg);
2509     charge = ParticleCharge(pdg);
2510     tlife  = ParticleLifeTime(pdg);
2511 }
2512
2513 //______________________________________________________________________________
2514 void TFluka::PrintHeader()
2515 {
2516     //
2517     // Print a header
2518     printf("\n");
2519     printf("\n");
2520     printf("------------------------------------------------------------------------------\n");
2521     printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA.           -\n");
2522     printf("- Please see the file fluka.out for FLUKA output and licensing information.  -\n");
2523     printf("------------------------------------------------------------------------------\n");
2524     printf("\n");
2525     printf("\n");
2526 }
2527
2528
2529 #define pshckp pshckp_
2530 #define ustckv ustckv_
2531
2532
2533 extern "C" {
2534   void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2535               Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2536               Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2537   {
2538     //
2539     // Pushes one cerenkov photon to the stack
2540     //
2541
2542     TFluka* fluka =  (TFluka*) gMC;
2543     TVirtualMCStack* cppstack = fluka->GetStack();
2544     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2545     cppstack->PushTrack(0, parent, 50000050,
2546                         px, py, pz, e,
2547                         vx, vy, vz, tof,
2548                         polx, poly, polz,
2549                         kPCerenkov, ntr, wgt, 0);
2550     if (fluka->GetVerbosityLevel() >= 3)
2551             printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
2552   }
2553
2554     void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
2555     {
2556         //
2557         // Calls stepping in order to signal cerenkov production
2558         //
2559         TFluka *fluka = (TFluka*)gMC;
2560         fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
2561         fluka->SetXsco(x);
2562         fluka->SetYsco(y);
2563         fluka->SetZsco(z);
2564         fluka->SetNCerenkov(nphot);
2565         fluka->SetCaller(kUSTCKV);
2566         if (fluka->GetVerbosityLevel() >= 3)
2567             printf("ustckv: %10d mreg=%d lattc=%d  newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
2568                     nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
2569
2570     // check region lattice consistency (debug Ernesto)
2571     // *****************************************************
2572    Int_t nodeId;
2573    Int_t volId = fluka->CurrentVolID(nodeId);
2574    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
2575
2576    if( mreg != volId  && !gGeoManager->IsOutside() ) {
2577        cout << "  ustckv:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
2578             << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
2579             << "               fluka   mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
2580             << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
2581             << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
2582             << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
2583             << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
2584             << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
2585         if( TRACKR.lt1trk == crtlttc ) cout << "   *************************************************************" << endl;
2586     }
2587     // *****************************************************
2588
2589
2590
2591         (TVirtualMCApplication::Instance())->Stepping();
2592     }
2593 }
2594
2595 //______________________________________________________________________________
2596 void TFluka::AddParticlesToPdgDataBase() const
2597 {
2598
2599 //
2600 // Add particles to the PDG data base
2601
2602     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
2603
2604     const Double_t kAu2Gev   = 0.9314943228;
2605     const Double_t khSlash   = 1.0545726663e-27;
2606     const Double_t kErg2Gev  = 1/1.6021773349e-3;
2607     const Double_t khShGev   = khSlash*kErg2Gev;
2608     const Double_t kYear2Sec = 3600*24*365.25;
2609 //
2610 // Ions
2611 //
2612   pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
2613                      0,3,"Ion",TFlukaIon::GetIonPdg(1,2));
2614   pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
2615                      khShGev/(12.33*kYear2Sec),3,"Ion",TFlukaIon::GetIonPdg(1,3));
2616   pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
2617                      khShGev/(12.33*kYear2Sec),6,"Ion",TFlukaIon::GetIonPdg(2,4));
2618   pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
2619                      0,6,"Ion",TFlukaIon::GetIonPdg(2,3));
2620 //
2621 // Default user ion
2622   TFlukaIon::AddIon(12, 6);
2623
2624 //
2625 //
2626 //
2627 // Special particles
2628 //
2629   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
2630                      0,0,"Special",GetSpecialPdg(50));
2631   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
2632                      0,0,"Special",GetSpecialPdg(51));
2633 }
2634
2635
2636 //
2637 // Info about primary ionization electrons
2638 //
2639
2640 //______________________________________________________________________________
2641 Int_t TFluka::GetNPrimaryElectrons()
2642 {
2643     // Get number of primary electrons
2644     return ALLDLT.nalldl;
2645 }
2646
2647 //______________________________________________________________________________
2648 Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
2649 {
2650     // Returns kinetic energy of primary electron i
2651
2652     Double_t ekin = -1.;
2653
2654     if (i >= 0 && i < ALLDLT.nalldl) {
2655         ekin =  ALLDLT.talldl[i];
2656     } else {
2657         Warning("GetPrimaryElectronKineticEnergy",
2658                 "Primary electron index out of range %d %d \n",
2659                 i, ALLDLT.nalldl);
2660     }
2661     return ekin;
2662 }
2663
2664 void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const
2665 {
2666     // Returns position  of primary electron i
2667         if (i >= 0 && i < ALLDLT.nalldl) {
2668             x = ALLDLT.xalldl[i];
2669             y = ALLDLT.yalldl[i];
2670             z = ALLDLT.zalldl[i];
2671             t = ALLDLT.talldl[i];
2672             return;
2673         } else {
2674             Warning("GetPrimaryElectronPosition",
2675                     "Primary electron index out of range %d %d \n",
2676                     i, ALLDLT.nalldl);
2677             return;
2678         }
2679         return;
2680 }
2681
2682
2683 //__________________________________________________________________
2684 Int_t TFluka::GetSpecialPdg(Int_t number) const
2685 {
2686 // Numbering for special particles
2687
2688   return 50000000 + number;
2689 }
2690
2691
2692 void  TFluka::PrimaryIonisationStepping(Int_t nprim)
2693 {
2694 // Call Stepping for primary ionisation electrons
2695 // Protection against nprim > mxalld
2696 // Multiple steps for nprim > 0
2697     Int_t i;
2698     fNPI = nprim;
2699     if (nprim > 0) {
2700         CalcPrimaryIonisationTime();
2701         for (i = 0; i < nprim; i++) {
2702             SetCurrentPrimaryElectronIndex(i);
2703             (TVirtualMCApplication::Instance())->Stepping();
2704             if (i == 0) SetTrackIsNew(kFALSE);
2705         }
2706     } else {
2707         // No primary electron ionisation
2708         // Call Stepping anyway but flag nprim = 0 as index = -2
2709         SetCurrentPrimaryElectronIndex(-2);
2710         (TVirtualMCApplication::Instance())->Stepping();
2711     }
2712     // Reset the index
2713     SetCurrentPrimaryElectronIndex(-1);
2714 }
2715
2716 //______________________________________________________________________
2717 Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
2718 {
2719 // Converts Double_t* array to Float_t*,
2720 // !! The new array has to be deleted by user.
2721 // ---
2722
2723   Float_t* floatArray;
2724   if (size>0) {
2725     floatArray = new Float_t[size];
2726     for (Int_t i=0; i<size; i++)
2727       if (array[i] >= FLT_MAX )
2728         floatArray[i] = FLT_MAX/100.;
2729       else
2730         floatArray[i] = array[i];
2731   }
2732   else {
2733     //floatArray = 0;
2734     floatArray = new Float_t[1];
2735   }
2736   return floatArray;
2737 }
2738
2739 void TFluka::CalcPrimaryIonisationTime()
2740 {
2741     // Calculates the primary ionisation time
2742     if (fPItime) delete [] fPItime;
2743     fPItime = new Double_t[fNPI];
2744     if (fPIlength) delete [] fPIlength;
2745     fPIlength = new Double_t[fNPI];
2746     //
2747     Double_t px, py, pz, e, t;
2748     TrackMomentum(px, py, pz, e);
2749     Double_t p    = TMath::Sqrt(px * px + py * py + pz * pz);
2750     Double_t beta = p / e;
2751     Double_t x0, y0, z0;
2752     fPItime[fNPI -1]   = TRACKR.atrack;
2753     fPIlength[fNPI -1] = TRACKR.cmtrck;
2754     GetPrimaryElectronPosition(fNPI - 1, x0, y0, z0, t);
2755     if (fNPI > 1) {
2756         for (Int_t i = fNPI - 2; i > -1; i--) {
2757             Double_t x, y, z, t;
2758             GetPrimaryElectronPosition(i, x, y, z, t);
2759             Double_t ds = TMath::Sqrt((x-x0) * (x-x0) + (y-y0) * (y-y0) + (z-z0) * (z-z0));
2760             fPItime[i]   = fPItime[i+1]   - ds / (beta * 2.99792458e10);
2761             fPIlength[i] = fPIlength[i+1] - ds;
2762             x0 = x; y0 = y; z0 = z;
2763         }
2764     }
2765
2766 }
2767
2768 Bool_t TFluka::DefineIon(const char* name , Int_t z, Int_t a, Int_t q, Double_t exE, Double_t mass)
2769 {
2770     // User defined ion that can be used as a primary
2771     if (fUserIon) {
2772         Warning("DefineIon", "Only one user ion can be defined !");
2773         return kFALSE;
2774     } else {
2775         fUserIon = new TFlukaIon(name, z, a, q, exE, mass);
2776         return kTRUE;
2777     }
2778 }