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