TFluka/TFluka.cxx

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