TFluka/TFluka.cxx

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