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