TFluka/TFluka.cxx

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