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