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    fProcesses(0),
 107    fCuts(0),
 108    fUserScore(0)
 109 {
 110   //
 111   // Default constructor
 112   //
 113    fGeneratePemf = kFALSE;
 114    fNVolumes = 0;
 115    fCurrentFlukaRegion = -1;
 116    fGeom = 0;
 117    fMCGeo = 0;
 118    fMaterials = 0;
 119    fDummyBoundary = 0;
 120    fFieldFlag = 1;
 121    fStopped   = 0;
 122    fStopEvent = 0;
 123    fStopRun   = 0;
 124    fNEvent    = 0;
 125 }
 126
 127 //______________________________________________________________________________
 128 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
 129   :TVirtualMC("TFluka",title, isRootGeometrySupported),
 130    fVerbosityLevel(verbosity),
 131    fInputFileName(""),
 132    fTrackIsEntering(0),
 133    fTrackIsExiting(0),
 134    fTrackIsNew(0),
 135    fProcesses(new TObjArray(100)),
 136    fCuts(new TObjArray(100)),
 137    fUserScore(new TObjArray(100))
 138 {
 139   // create geometry interface
 140    if (fVerbosityLevel >=3)
 141        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 142    SetCoreInputFileName();
 143    SetInputFileName();
 144    SetGeneratePemf(kFALSE);
 145    fNVolumes      = 0;
 146    fCurrentFlukaRegion = -1;
 147    fDummyBoundary = 0;
 148    fFieldFlag = 1;
 149    fGeneratePemf = kFALSE;
 150    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
 151    fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
 152    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
 153    fMaterials = 0;
 154    fStopped   = 0;
 155    fStopEvent = 0;
 156    fStopRun   = 0;
 157    fNEvent    = 0;
 158 }
 159
 160 //______________________________________________________________________________
 161 TFluka::~TFluka() {
 162 // Destructor
 163     if (fVerbosityLevel >=3)
 164         cout << "<== TFluka::~TFluka() destructor called." << endl;
 165
 166     delete fGeom;
 167     delete fMCGeo;
 168
 169     if (fCuts) {
 170         fCuts->Delete();
 171         delete fCuts;
 172     }
 173
 174     if (fProcesses) {
 175         fProcesses->Delete();
 176         delete fProcesses;
 177     }
 178
 179
 180 }
 181
 182 //
 183 //______________________________________________________________________________
 184 // TFluka control methods
 185 //______________________________________________________________________________
 186 void TFluka::Init() {
 187 //
 188 //  Geometry initialisation
 189 //
 190     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
 191
 192     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
 193     fApplication->ConstructGeometry();
 194     TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
 195     gGeoManager->SetTopVolume(top);
 196     gGeoManager->CloseGeometry("di");
 197     gGeoManager->DefaultColors();  // to be removed
 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     // now we have TGeo geometry created and we have to patch FlukaVmc.inp
 205     // with the material mapping file FlukaMat.inp
 206 }
 207
 208
 209 //______________________________________________________________________________
 210 void TFluka::FinishGeometry() {
 211 //
 212 // Build-up table with region to medium correspondance
 213 //
 214   if (fVerbosityLevel >=3) {
 215     cout << "==> TFluka::FinishGeometry() called." << endl;
 216     printf("----FinishGeometry - nothing to do with TGeo\n");
 217     cout << "<== TFluka::FinishGeometry() called." << endl;
 218   }
 219 }
 220
 221 //______________________________________________________________________________
 222 void TFluka::BuildPhysics() {
 223 //
 224 //  Prepare FLUKA input files and call FLUKA physics initialisation
 225 //
 226
 227     if (fVerbosityLevel >=3)
 228         cout << "==> TFluka::BuildPhysics() called." << endl;
 229 // Prepare input file with the current physics settings
 230     InitPhysics();
 231     cout << "\t* InitPhysics() - Prepare input file was called" << endl;
 232
 233     if (fVerbosityLevel >=2)
 234         cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
 235              << ") in fluka..." << endl;
 236     GLOBAL.lfdrtr = true;
 237
 238     if (fVerbosityLevel >=2)
 239         cout << "\t* Opening file " << fInputFileName << endl;
 240     const char* fname = fInputFileName;
 241     fluka_openinp(lunin, PASSCHARA(fname));
 242
 243     if (fVerbosityLevel >=2)
 244         cout << "\t* Calling flukam..." << endl;
 245     flukam(1);
 246
 247     if (fVerbosityLevel >=2)
 248         cout << "\t* Closing file " << fInputFileName << endl;
 249     fluka_closeinp(lunin);
 250
 251     FinishGeometry();
 252
 253     if (fVerbosityLevel >=3)
 254         cout << "<== TFluka::Init() called." << endl;
 255
 256
 257     if (fVerbosityLevel >=3)
 258         cout << "<== TFluka::BuildPhysics() called." << endl;
 259 }
 260
 261 //______________________________________________________________________________
 262 void TFluka::ProcessEvent() {
 263 //
 264 // Process one event
 265 //
 266     if (fStopRun) {
 267         printf("User Run Abortion: No more events handled !\n");
 268         fNEvent += 1;
 269         return;
 270     }
 271
 272     if (fVerbosityLevel >=3)
 273         cout << "==> TFluka::ProcessEvent() called." << endl;
 274     fApplication->GeneratePrimaries();
 275     EPISOR.lsouit = true;
 276     flukam(1);
 277     if (fVerbosityLevel >=3)
 278         cout << "<== TFluka::ProcessEvent() called." << endl;
 279     //
 280     // Increase event number
 281     //
 282     fNEvent += 1;
 283 }
 284
 285 //______________________________________________________________________________
 286 Bool_t TFluka::ProcessRun(Int_t nevent) {
 287 //
 288 // Run steering
 289 //
 290
 291   if (fVerbosityLevel >=3)
 292     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 293          << endl;
 294
 295   if (fVerbosityLevel >=2) {
 296     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 297     cout << "\t* Calling flukam again..." << endl;
 298   }
 299
 300   fApplication->InitGeometry();
 301   Int_t todo = TMath::Abs(nevent);
 302   for (Int_t ev = 0; ev < todo; ev++) {
 303       fApplication->BeginEvent();
 304       ProcessEvent();
 305       fApplication->FinishEvent();
 306   }
 307
 308   if (fVerbosityLevel >=3)
 309     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 310          << endl;
 311   // Write fluka specific scoring output
 312   newplo();
 313
 314   return kTRUE;
 315 }
 316
 317 //_____________________________________________________________________________
 318 // methods for building/management of geometry
 319
 320 // functions from GCONS
 321 //____________________________________________________________________________
 322 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 323                     Float_t &dens, Float_t &radl, Float_t &absl,
 324                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 325 //
 326    TGeoMaterial *mat;
 327    TIter next (gGeoManager->GetListOfMaterials());
 328    while ((mat = (TGeoMaterial*)next())) {
 329      if (mat->GetUniqueID() == (UInt_t)imat) break;
 330    }
 331    if (!mat) {
 332       Error("Gfmate", "no material with index %i found", imat);
 333       return;
 334    }
 335    sprintf(name, "%s", mat->GetName());
 336    a = mat->GetA();
 337    z = mat->GetZ();
 338    dens = mat->GetDensity();
 339    radl = mat->GetRadLen();
 340    absl = mat->GetIntLen();
 341 }
 342
 343 //______________________________________________________________________________
 344 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 345                     Double_t &dens, Double_t &radl, Double_t &absl,
 346                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 347 //
 348    TGeoMaterial *mat;
 349    TIter next (gGeoManager->GetListOfMaterials());
 350    while ((mat = (TGeoMaterial*)next())) {
 351      if (mat->GetUniqueID() == (UInt_t)imat) break;
 352    }
 353    if (!mat) {
 354       Error("Gfmate", "no material with index %i found", imat);
 355       return;
 356    }
 357    sprintf(name, "%s", mat->GetName());
 358    a = mat->GetA();
 359    z = mat->GetZ();
 360    dens = mat->GetDensity();
 361    radl = mat->GetRadLen();
 362    absl = mat->GetIntLen();
 363 }
 364
 365 // detector composition
 366 //______________________________________________________________________________
 367 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 368                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 369                       Float_t* buf, Int_t nwbuf) {
 370 //
 371    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
 372    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 373    delete [] dbuf;
 374 }
 375
 376 //______________________________________________________________________________
 377 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 378                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 379                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
 380 //
 381   TGeoMaterial *mat;
 382   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 383   if ((z-Int_t(z)) > 1E-3) {
 384      mat = fGeom->GetMakeWrongMaterial(z);
 385      if (mat) {
 386         mat->SetRadLen(radl,absl);
 387         mat->SetUniqueID(kmat);
 388         return;
 389      }
 390   }
 391   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 392 }
 393
 394 //______________________________________________________________________________
 395 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 396                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 397 //
 398   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
 399   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
 400   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
 401
 402   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 403   for (Int_t i=0; i<nlmat; i++) {
 404     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 405   }
 406
 407   delete [] da;
 408   delete [] dz;
 409   delete [] dwmat;
 410 }
 411
 412 //______________________________________________________________________________
 413 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 414                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 415 //
 416   // Defines mixture OR COMPOUND IMAT as composed by
 417   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 418   //
 419   // If NLMAT > 0 then wmat contains the proportion by
 420   // weights of each basic material in the mixture.
 421   //
 422   // If nlmat < 0 then WMAT contains the number of atoms
 423   // of a given kind into the molecule of the COMPOUND
 424   // In this case, WMAT in output is changed to relative
 425   // weigths.
 426   //
 427   Int_t i,j;
 428   if (nlmat < 0) {
 429      nlmat = - nlmat;
 430      Double_t amol = 0;
 431      for (i=0;i<nlmat;i++) {
 432         amol += a[i]*wmat[i];
 433      }
 434      for (i=0;i<nlmat;i++) {
 435         wmat[i] *= a[i]/amol;
 436      }
 437   }
 438   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 439   // Check if we have elements with fractional Z
 440   TGeoMaterial *mat = 0;
 441   TGeoMixture *mix = 0;
 442   Bool_t mixnew = kFALSE;
 443   for (i=0; i<nlmat; i++) {
 444      if (z[i]-Int_t(z[i]) < 1E-3) continue;
 445      // We have found an element with fractional Z -> loop mixtures to look for it
 446      for (j=0; j<kmat; j++) {
 447         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
 448         if (!mat) break;
 449         if (!mat->IsMixture()) continue;
 450         mix = (TGeoMixture*)mat;
 451         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
 452 //        printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
 453         mixnew = kTRUE;
 454         break;
 455      }
 456      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
 457      break;
 458   }
 459   if (mixnew) {
 460      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
 461      Double_t *anew = new Double_t[nlmatnew];
 462      Double_t *znew = new Double_t[nlmatnew];
 463      Double_t *wmatnew = new Double_t[nlmatnew];
 464      Int_t ind=0;
 465      for (j=0; j<nlmat; j++) {
 466         if (j==i) continue;
 467         anew[ind] = a[j];
 468         znew[ind] = z[j];
 469         wmatnew[ind] = wmat[j];
 470         ind++;
 471      }
 472      for (j=0; j<mix->GetNelements(); j++) {
 473         anew[ind] = mix->GetAmixt()[j];
 474         znew[ind] = mix->GetZmixt()[j];
 475         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
 476         ind++;
 477      }
 478      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 479      delete [] anew;
 480      delete [] znew;
 481      delete [] wmatnew;
 482      return;
 483   }
 484   // Now we need to compact identical elements within the mixture
 485   // First check if this happens
 486   mixnew = kFALSE;
 487   for (i=0; i<nlmat-1; i++) {
 488      for (j=i+1; j<nlmat; j++) {
 489         if (z[i] == z[j]) {
 490            mixnew = kTRUE;
 491            break;
 492         }
 493      }
 494      if (mixnew) break;
 495   }
 496   if (mixnew) {
 497      Int_t nlmatnew = 0;
 498      Double_t *anew = new Double_t[nlmat];
 499      Double_t *znew = new Double_t[nlmat];
 500      memset(znew, 0, nlmat*sizeof(Double_t));
 501      Double_t *wmatnew = new Double_t[nlmat];
 502      Bool_t skipi;
 503      for (i=0; i<nlmat; i++) {
 504         skipi = kFALSE;
 505         for (j=0; j<nlmatnew; j++) {
 506            if (z[i] == z[j]) {
 507               wmatnew[j] += wmat[i];
 508               skipi = kTRUE;
 509               break;
 510            }
 511         }
 512         if (skipi) continue;
 513         anew[nlmatnew] = a[i];
 514         znew[nlmatnew] = z[i];
 515         wmatnew[nlmatnew] = wmat[i];
 516         nlmatnew++;
 517      }
 518      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 519      delete [] anew;
 520      delete [] znew;
 521      delete [] wmatnew;
 522      return;
 523    }
 524    gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 525 }
 526
 527 //______________________________________________________________________________
 528 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 529                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 530                     Double_t stemax, Double_t deemax, Double_t epsil,
 531                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 532   // Define a medium
 533   //
 534   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 535   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 536              epsil, stmin, ubuf, nbuf);
 537 }
 538
 539 //______________________________________________________________________________
 540 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 541                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 542                     Double_t stemax, Double_t deemax, Double_t epsil,
 543                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 544   // Define a medium
 545   //
 546   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 547   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 548              epsil, stmin, ubuf, nbuf);
 549 }
 550
 551 //______________________________________________________________________________
 552 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 553                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 554                     Double_t phiZ) {
 555 //
 556   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 557   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 558 }
 559
 560 //______________________________________________________________________________
 561 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 562 //
 563 //
 564 // Check if material is used
 565    if (fVerbosityLevel >=3)
 566        printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
 567    Int_t* reglist;
 568    Int_t nreg;
 569    reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg);
 570    if (nreg == 0) {
 571        return;
 572    }
 573
 574 //
 575    Bool_t process = kFALSE;
 576    if (strncmp(param, "DCAY",  4) == 0 ||
 577        strncmp(param, "PAIR",  4) == 0 ||
 578        strncmp(param, "COMP",  4) == 0 ||
 579        strncmp(param, "PHOT",  4) == 0 ||
 580        strncmp(param, "PFIS",  4) == 0 ||
 581        strncmp(param, "DRAY",  4) == 0 ||
 582        strncmp(param, "ANNI",  4) == 0 ||
 583        strncmp(param, "BREM",  4) == 0 ||
 584        strncmp(param, "MUNU",  4) == 0 ||
 585        strncmp(param, "CKOV",  4) == 0 ||
 586        strncmp(param, "HADR",  4) == 0 ||
 587        strncmp(param, "LOSS",  4) == 0 ||
 588        strncmp(param, "MULS",  4) == 0 ||
 589        strncmp(param, "RAYL",  4) == 0)
 590    {
 591        process = kTRUE;
 592    }
 593    if (process) {
 594        SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
 595    } else {
 596        SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
 597    }
 598 }
 599
 600 // functions from GGEOM
 601 //_____________________________________________________________________________
 602 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 603 {
 604   // Set visualisation attributes for one volume
 605   char vname[5];
 606   fGeom->Vname(name,vname);
 607   char vatt[5];
 608   fGeom->Vname(att,vatt);
 609   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 610 }
 611
 612 //______________________________________________________________________________
 613 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 614                      Float_t *upar, Int_t np)  {
 615 //
 616     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 617 }
 618
 619 //______________________________________________________________________________
 620 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 621                      Double_t *upar, Int_t np)  {
 622 //
 623     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 624 }
 625
 626 //______________________________________________________________________________
 627 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 628                    Int_t iaxis) {
 629 //
 630     fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
 631 }
 632
 633 //______________________________________________________________________________
 634 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 635                     Int_t iaxis, Double_t c0i, Int_t numed) {
 636 //
 637     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 638 }
 639
 640 //______________________________________________________________________________
 641 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 642                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 643 //
 644     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 645 }
 646
 647 //______________________________________________________________________________
 648 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 649                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 650 //
 651     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 652 }
 653
 654 //______________________________________________________________________________
 655 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 656 //
 657 // Nothing to do with TGeo
 658 }
 659
 660 //______________________________________________________________________________
 661 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 662                    Double_t x, Double_t y, Double_t z, Int_t irot,
 663                    const char *konly) {
 664 //
 665   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
 666 }
 667
 668 //______________________________________________________________________________
 669 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 670                     Double_t x, Double_t y, Double_t z, Int_t irot,
 671                     const char *konly, Float_t *upar, Int_t np)  {
 672   //
 673   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 674 }
 675
 676 //______________________________________________________________________________
 677 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 678                     Double_t x, Double_t y, Double_t z, Int_t irot,
 679                     const char *konly, Double_t *upar, Int_t np)  {
 680   //
 681   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 682 }
 683
 684 //______________________________________________________________________________
 685 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
 686 //
 687 // Nothing to do with TGeo
 688 }
 689
 690 //______________________________________________________________________________
 691 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 692                          Float_t* absco, Float_t* effic, Float_t* rindex) {
 693 //
 694 // Set Cerenkov properties for medium itmed
 695 //
 696 // npckov: number of sampling points
 697 // ppckov: energy values
 698 // absco:  absorption length
 699 // effic:  quantum efficiency
 700 // rindex: refraction index
 701 //
 702 //
 703 //
 704 //  Create object holding Cerenkov properties
 705 //
 706     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 707 //
 708 //  Pass object to medium
 709     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 710     medium->SetCerenkovProperties(cerenkovProperties);
 711 }
 712
 713 //______________________________________________________________________________
 714 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
 715                          Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
 716 //
 717 // Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
 718    Warning("SetCerenkov", "Not implemented with TGeo");
 719 }
 720
 721 // Euclid
 722 //______________________________________________________________________________
 723 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
 724                           Int_t /*number*/, Int_t /*nlevel*/) {
 725 //
 726 // Not with TGeo
 727    Warning("WriteEuclid", "Not implemented with TGeo");
 728 }
 729
 730
 731
 732 //_____________________________________________________________________________
 733 // methods needed by the stepping
 734 //____________________________________________________________________________
 735
 736 Int_t TFluka::GetMedium() const {
 737 //
 738 //  Get the medium number for the current fluka region
 739 //
 740     return fGeom->GetMedium(); // this I need to check due to remapping !!!
 741 }
 742
 743
 744
 745 //____________________________________________________________________________
 746 // particle table usage
 747 // ID <--> PDG transformations
 748 //_____________________________________________________________________________
 749 Int_t TFluka::IdFromPDG(Int_t pdg) const
 750 {
 751     //
 752     // Return Fluka code from PDG and pseudo ENDF code
 753
 754     // Catch the feedback photons
 755     if (pdg == 50000051) return (-1);
 756     // MCIHAD() goes from pdg to fluka internal.
 757     Int_t intfluka = mcihad(pdg);
 758     // KPTOIP array goes from internal to official
 759     return GetFlukaKPTOIP(intfluka);
 760 }
 761
 762 //______________________________________________________________________________
 763 Int_t TFluka::PDGFromId(Int_t id) const
 764 {
 765   //
 766   // Return PDG code and pseudo ENDF code from Fluka code
 767   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon
 768     Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
 769   // IPTOKP array goes from official to internal
 770
 771     if (id == -1) {
 772 // Cerenkov photon
 773         if (fVerbosityLevel >= 3)
 774             printf("\n PDGFromId: Cerenkov Photon \n");
 775         return  50000050;
 776     }
 777 // Error id
 778     if (id == 0 || id < -6 || id > 250) {
 779         if (fVerbosityLevel >= 3)
 780             printf("PDGFromId: Error id = 0\n");
 781         return -1;
 782     }
 783 // Good id
 784     if (id > 0) {
 785         Int_t intfluka = GetFlukaIPTOKP(id);
 786         if (intfluka == 0) {
 787             if (fVerbosityLevel >= 3)
 788                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
 789             return -1;
 790         } else if (intfluka < 0) {
 791             if (fVerbosityLevel >= 3)
 792                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
 793             return -1;
 794         }
 795         if (fVerbosityLevel >= 3)
 796             printf("mpdgha called with %d %d \n", id, intfluka);
 797         // MPDGHA() goes from fluka internal to pdg.
 798         return mpdgha(intfluka);
 799     } else {
 800         // ions and optical photons
 801         return idSpecial[id + 6];
 802     }
 803 }
 804
 805 void TFluka::StopTrack()
 806 {
 807     // Set stopping conditions
 808     // Works for photons and charged particles
 809     fStopped = kTRUE;
 810 }
 811
 812 //_____________________________________________________________________________
 813 // methods for physics management
 814 //____________________________________________________________________________
 815 //
 816 // set methods
 817 //
 818
 819 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
 820 {
 821 //  Set process user flag for material imat
 822 //
 823     TFlukaConfigOption* proc = new TFlukaConfigOption(flagName, flagValue, imed);
 824     fProcesses->Add(proc);
 825 }
 826
 827 //______________________________________________________________________________
 828 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
 829 {
 830 //  Set process user flag
 831 //
 832 //
 833 //  Update if already in the list
 834 //
 835
 836     TIter next(fProcesses);
 837     TFlukaConfigOption* proc;
 838     while((proc = (TFlukaConfigOption*)next()))
 839     {
 840         if (strcmp(proc->GetName(), flagName) == 0) {
 841             proc->SetFlag(flagValue);
 842             proc->SetMedium(-1);
 843             return kTRUE;
 844          }
 845     }
 846 //
 847 // If not create a new process
 848 //
 849
 850     proc = new TFlukaConfigOption(flagName, flagValue);
 851     fProcesses->Add(proc);
 852
 853     return kTRUE;
 854 }
 855
 856 //______________________________________________________________________________
 857 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
 858 {
 859 // Set user cut value for material imed
 860 //
 861     TFlukaConfigOption* cut = new TFlukaConfigOption(cutName, cutValue, imed);
 862     fCuts->Add(cut);
 863 }
 864
 865 //______________________________________________________________________________
 866 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
 867 {
 868 // Set user cut value
 869 //
 870 //
 871 //  Update if already in the list
 872 //
 873
 874     TIter next(fCuts);
 875     TFlukaConfigOption* cut;
 876     while((cut = (TFlukaConfigOption*)next()))
 877     {
 878         if (strcmp(cut->GetName(), cutName) == 0) {
 879             cut->SetCut(cutValue);
 880             return kTRUE;
 881          }
 882     }
 883 //
 884 // If not create a new process
 885 //
 886
 887     cut = new TFlukaConfigOption(cutName, cutValue);
 888     fCuts->Add(cut);
 889
 890     return kTRUE;
 891 }
 892
 893 void TFluka::SetUserScoring(const char* option, Int_t npar, Float_t what[12])
 894 {
 895 //
 896 // Ads a user scoring option to th list
 897 //
 898     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npar, what);
 899     fUserScore->Add(opt);
 900 }
 901
 902
 903 //______________________________________________________________________________
 904 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
 905 {
 906   printf("WARNING: Xsec not yet implemented !\n"); return -1.;
 907 }
 908
 909
 910 //______________________________________________________________________________
 911 void TFluka::InitPhysics()
 912 {
 913 //
 914 // Physics initialisation with preparation of FLUKA input cards
 915 //
 916   printf("=>InitPhysics\n");
 917   Int_t j, k;
 918   Double_t theCut;
 919
 920   FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
 921
 922   Double_t zero  = 0.0;
 923   Double_t one   = 1.0;
 924   Double_t two   = 2.0;
 925   Double_t three = 3.0;
 926
 927   Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
 928   if (fVerbosityLevel >= 3) printf("   last FLUKA material is %g\n", fLastMaterial);
 929
 930   // Prepare  Cerenkov
 931   TObjArray *matList = GetFlukaMaterials();
 932   Int_t nmaterial =  matList->GetEntriesFast();
 933   fMaterials = new Int_t[nmaterial+3];
 934
 935 // construct file names
 936
 937   TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
 938   sFlukaVmcCoreInp +="/TFluka/input/";
 939   TString sFlukaVmcTmp = "flukaMat.inp";
 940   TString sFlukaVmcInp = GetInputFileName();
 941   sFlukaVmcCoreInp += GetCoreInputFileName();
 942
 943 // open files
 944
 945   if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
 946       printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
 947       exit(1);
 948   }
 949   if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
 950       printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
 951       exit(1);
 952   }
 953   if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
 954       printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
 955       exit(1);
 956   }
 957
 958 // copy core input file
 959   Char_t sLine[255];
 960   Float_t fEventsPerRun;
 961
 962   while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
 963       if (strncmp(sLine,"GEOEND",6) != 0)
 964           fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
 965       else {
 966           fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
 967           goto flukamat;
 968       }
 969   } // end of while until GEOEND card
 970
 971
 972  flukamat:
 973   while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
 974       fprintf(pFlukaVmcInp,"%s\n",sLine);
 975   }
 976
 977   while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
 978       if (strncmp(sLine,"START",5) != 0)
 979           fprintf(pFlukaVmcInp,"%s\n",sLine);
 980       else {
 981           sscanf(sLine+10,"%10f",&fEventsPerRun);
 982       goto fin;
 983       }
 984   } //end of while until START card
 985
 986 fin:
 987 // in G3 the process control values meaning can be different for
 988 // different processes, but for most of them is:
 989 //   0  process is not activated
 990 //   1  process is activated WITH generation of secondaries
 991 //   2  process is activated WITHOUT generation of secondaries
 992 // if process does not generate secondaries => 1 same as 2
 993 //
 994 // Exceptions:
 995 //   MULS:  also 3
 996 //   LOSS:  also 3, 4
 997 //   RAYL:  only 0,1
 998 //   HADR:  may be > 2
 999 //
1000
1001 // Loop over number of SetProcess calls
1002   fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
1003   fprintf(pFlukaVmcInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
1004   fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
1005
1006 // Outer loop over processes
1007   TIter next(fProcesses);
1008   TFlukaConfigOption *proc;
1009 // Inner loop over processes
1010   TIter nextp(fProcesses);
1011   TFlukaConfigOption *procp;
1012 // Loop over cuts
1013   TIter nextc(fCuts);
1014   TFlukaConfigOption *cut = 0x0;
1015
1016   while((proc = (TFlukaConfigOption*)next())) {
1017       Float_t matMin = three;
1018       Float_t matMax = fLastMaterial;
1019       Bool_t  global = kTRUE;
1020       if (proc->Medium() != -1) {
1021           Int_t mat;
1022           if ((mat = proc->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
1023           matMin = Float_t(mat);
1024           matMax = matMin;
1025           global = kFALSE;
1026
1027           fprintf(pFlukaVmcInp,"*\n*Material specific process setting for #%8d \n", mat);
1028       }
1029
1030     // annihilation
1031     // G3 default value: 1
1032     // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
1033     // Particles: e+
1034     // Physics:   EM
1035     // flag = 0 no annihilation
1036     // flag = 1 annihilation, decays processed
1037     // flag = 2 annihilation, no decay product stored
1038     // gMC ->SetProcess("ANNI",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. ANNH-THR
1039       if (strncmp(proc->GetName(),"ANNI",4) == 0) {
1040           if (proc->Flag() == 1 || proc->Flag() == 2) {
1041               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
1042               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
1043               // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
1044               // zero = not used
1045               // zero = not used
1046               // matMin = lower bound of the material indices in which the respective thresholds apply
1047               // matMax = upper bound of the material indices in which the respective thresholds apply
1048               // one = step length in assigning indices
1049               // "ANNH-THR";
1050               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one);
1051           }
1052           else if (proc->Flag() == 0) {
1053               fprintf(pFlukaVmcInp,"*\n*No annihilation - no FLUKA card generated\n");
1054               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',0)\n");
1055           }
1056           else  {
1057               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
1058               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1059           }
1060       }
1061
1062     // bremsstrahlung and pair production are both activated
1063     // G3 default value: 1
1064     // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
1065     //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
1066     //               G4LowEnergyBremstrahlung
1067     // Particles: e-/e+; mu+/mu-
1068     // Physics:   EM
1069     // flag = 0 no bremsstrahlung
1070     // flag = 1 bremsstrahlung, photon processed
1071     // flag = 2 bremsstrahlung, no photon stored
1072     // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
1073                                  // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
1074     // G3 default value: 1
1075     // G4 processes: G4GammaConversion,
1076     //               G4MuPairProduction/G4IMuPairProduction
1077     //               G4LowEnergyGammaConversion
1078     // Particles: gamma, mu
1079     // Physics:   EM
1080     // flag = 0 no delta rays
1081     // flag = 1 delta rays, secondaries processed
1082     // flag = 2 delta rays, no secondaries stored
1083
1084     else if ((strncmp(proc->GetName(),"PAIR",4) == 0) && (proc->Flag() == 1 || proc->Flag() == 2)) {
1085
1086         nextp.Reset();
1087
1088         while ((procp = (TFlukaConfigOption*)nextp())) {
1089             if ((strncmp(procp->GetName(),"BREM",4) == 0) &&
1090                 (proc->Flag() == 1) &&
1091                 (procp->Medium() == proc->Medium())) {
1092                 fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
1093                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
1094                 fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
1095                 fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
1096                 // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
1097                 fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f",three);
1098                 // direct pair production by muons
1099                 // G4 particles: "e-", "e+"
1100                 // G3 default value: 0.01 GeV
1101                 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
1102                 theCut = 0.0;
1103                 nextc.Reset();
1104                 while ((cut = (TFlukaConfigOption*)nextc())) {
1105                     if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 &&
1106                         (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1107                 }
1108                 fprintf(pFlukaVmcInp,"%10.4g",theCut);
1109                 // theCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
1110                 // muon and hadron bremsstrahlung
1111                 // G4 particles: "gamma"
1112                 // G3 default value: CUTGAM=0.001 GeV
1113                 //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
1114                 theCut = 0.0;
1115                 nextc.Reset();
1116                 while ((cut = (TFlukaConfigOption*)nextc())) {
1117                     if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
1118                         (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1119                 }
1120                 fprintf(pFlukaVmcInp,"%10.4g%10.1f%10.1f\n",theCut,matMin,matMax);
1121                 // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production
1122                 // matMin = lower bound of the material indices in which the respective thresholds apply
1123                 // matMax = upper bound of the material indices in which the respective thresholds apply
1124
1125                 // for e+ and e-
1126                 fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
1127                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);\n");
1128                 theCut = -1.0;
1129                 nextc.Reset();
1130                 while ((cut = (TFlukaConfigOption*)nextc())) {
1131                     if (strncmp(cut->GetName(), "BCUTE", 5) == 0 &&
1132                         (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1133                 }
1134                 //theCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
1135                 // zero = not used
1136                 // zero = not used
1137                 // matMin = lower bound of the material indices in which the respective thresholds apply
1138                 // matMax = upper bound of the material indices in which the respective thresholds apply
1139                 // one = step length in assigning indices
1140                 // "ELPO-THR";
1141                 fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",theCut,zero,zero,matMin,matMax,one);
1142
1143                 // for gamma -> e+ and e-
1144                 fprintf(pFlukaVmcInp,"*\n*Pair production by photons is activated\n");
1145                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1);\n");
1146                 theCut = -1.0;
1147                 nextc.Reset();
1148                 while ((cut = (TFlukaConfigOption*)nextc())) {
1149                     if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
1150                         (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1151                 }
1152                 // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
1153                 // matMin = lower bound of the material indices in which the respective thresholds apply
1154                 // matMax =  upper bound of the material indices in which the respective thresholds apply
1155                 // one = step length in assigning indices
1156                 fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one);
1157                 goto BOTH;
1158             } // end of if for BREM
1159         } // end of loop for BREM
1160
1161         // only pair production by muons and charged hadrons is activated
1162         fprintf(pFlukaVmcInp,"*\n*Pair production by muons and charged hadrons is activated\n");
1163         fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
1164         fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
1165         // direct pair production by muons
1166         // G4 particles: "e-", "e+"
1167         // G3 default value: 0.01 GeV
1168         //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
1169         // one = pair production by muons and charged hadrons is activated
1170         // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
1171         // zero = no explicit bremsstrahlung production is simulated
1172         // matMin = lower bound of the material indices in which the respective thresholds apply
1173         // matMax = upper bound of the material indices in which the respective thresholds apply
1174         fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1175
1176         // for gamma -> e+ and e-
1177         fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
1178         fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
1179         theCut = -1.0;
1180         nextc.Reset();
1181         while ((cut = (TFlukaConfigOption*)nextc())) {
1182             if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
1183                 (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1184         }
1185         // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
1186         // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
1187         // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
1188         // matMin = lower bound of the material indices in which the respective thresholds apply
1189         // matMax = upper bound of the material indices in which the respective thresholds apply
1190         // one = step length in assigning indices
1191         fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one);
1192
1193     BOTH:
1194         k = 0;
1195     } // end of if for PAIR
1196
1197
1198
1199       // bremsstrahlung
1200       // G3 default value: 1
1201       // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
1202       //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
1203       //               G4LowEnergyBremstrahlung
1204       // Particles: e-/e+; mu+/mu-
1205       // Physics:   EM
1206       // flag = 0 no bremsstrahlung
1207       // flag = 1 bremsstrahlung, photon processed
1208       // flag = 2 bremsstrahlung, no photon stored
1209       // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
1210       // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
1211       else if (strncmp(proc->GetName(),"BREM",4) == 0) {
1212           nextp.Reset();
1213           while((procp = (TFlukaConfigOption*)nextp())) {
1214               if ((strncmp(procp->GetName(),"PAIR",4) == 0) &&
1215                   procp->Flag() == 1 &&
1216                   (procp->Medium() == proc->Medium())) goto NOBREM;
1217           }
1218           if (proc->Flag() == 1) {
1219               fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
1220               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
1221               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
1222               // two = bremsstrahlung by muons and charged hadrons is activated
1223               // zero = no meaning
1224               // muon and hadron bremsstrahlung
1225               // G4 particles: "gamma"
1226               // G3 default value: CUTGAM=0.001 GeV
1227               //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
1228               theCut = 0.0;
1229               nextc.Reset();
1230               while ((cut = (TFlukaConfigOption*)nextc())) {
1231                   if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
1232                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1233               }
1234               // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production
1235               // matMin = lower bound of the material indices in which the respective thresholds apply
1236               // matMax = upper bound of the material indices in which the respective thresholds apply
1237               fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,theCut,matMin,matMax);
1238
1239               // for e+ and e-
1240               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
1241               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);");
1242               // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
1243               // zero = not used
1244               // zero = not used
1245               // matMin = lower bound of the material indices in which the respective thresholds apply
1246               // matMax = upper bound of the material indices in which the respective thresholds apply
1247               // one = step length in assigning indices
1248               //"ELPO-THR";
1249               theCut = -1.0;
1250               nextc.Reset();
1251               while ((cut = (TFlukaConfigOption*)nextc())) {
1252                   if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
1253                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1254               }
1255               fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n", theCut,zero,zero,matMin,matMax,one);
1256           }
1257           else if (proc->Flag() == 0) {
1258               fprintf(pFlukaVmcInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
1259               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',0)\n");
1260           }
1261           else  {
1262               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
1263               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1264           }
1265       NOBREM:
1266           j = 0;
1267       } // end of else if (strncmp(proc->GetName(),"BREM",4) == 0)
1268
1269       // Cerenkov photon generation
1270       // G3 default value: 0
1271       // G4 process: G4Cerenkov
1272       //
1273       // Particles: charged
1274       // Physics:   Optical
1275       // flag = 0 no Cerenkov photon generation
1276       // flag = 1 Cerenkov photon generation
1277       // flag = 2 Cerenkov photon generation with primary stopped at each step
1278       //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
1279
1280       else if (strncmp(proc->GetName(),"CKOV",4) == 0) {
1281           if ((proc->Flag() == 1 || proc->Flag() == 2) && global) {
1282               // Write comments
1283               fprintf(pFlukaVmcInp, "* \n");
1284               fprintf(pFlukaVmcInp, "*Cerenkov photon generation\n");
1285               fprintf(pFlukaVmcInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n");
1286               // Loop over media
1287               for (Int_t im = 0; im < nmaterial; im++)
1288               {
1289                   TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1290                   Int_t idmat = material->GetIndex();
1291
1292                   if (!global && idmat != proc->Medium()) continue;
1293
1294                   fMaterials[idmat] = im;
1295                   // Skip media with no Cerenkov properties
1296                   TFlukaCerenkov* cerenkovProp;
1297                   if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
1298                   //
1299                   // This medium has Cerenkov properties
1300                   //
1301                   //
1302                   // Write OPT-PROD card for each medium
1303                   Float_t  emin  = cerenkovProp->GetMinimumEnergy();
1304                   Float_t  emax  = cerenkovProp->GetMaximumEnergy();
1305                   fprintf(pFlukaVmcInp, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
1306                           Float_t(idmat), Float_t(idmat), 0.);
1307                   //
1308                   // Write OPT-PROP card for each medium
1309                   // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
1310                   //
1311                   fprintf(pFlukaVmcInp, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
1312                           cerenkovProp->GetMinimumWavelength(),
1313                           cerenkovProp->GetMaximumWavelength(),
1314                           cerenkovProp->GetMaximumWavelength(),
1315                           Float_t(idmat), Float_t(idmat), 0.0);
1316
1317                   if (cerenkovProp->IsMetal()) {
1318                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n",
1319                               -100., -100., -100.,
1320                               Float_t(idmat), Float_t(idmat), 0.0);
1321                   } else {
1322                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
1323                               -100., -100., -100.,
1324                               Float_t(idmat), Float_t(idmat), 0.0);
1325                   }
1326
1327
1328                   for (Int_t j = 0; j < 3; j++) {
1329                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n",
1330                               -100., -100., -100.,
1331                               Float_t(idmat), Float_t(idmat), 0.0);
1332                   }
1333                   // Photon detection efficiency user defined
1334
1335                   if (cerenkovProp->IsSensitive())
1336                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n",
1337                               -100., -100., -100.,
1338                               Float_t(idmat), Float_t(idmat), 0.0);
1339
1340               } // materials
1341           } else if (proc->Flag() == 0) {
1342               fprintf(pFlukaVmcInp,"*\n*No Cerenkov photon generation\n");
1343               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('CKOV',0)\n");
1344               // zero = not used
1345               // zero = not used
1346               // zero = not used
1347               // matMin = lower bound of the material indices in which the respective thresholds apply
1348               // matMax = upper bound of the material indices in which the respective thresholds apply
1349               // one = step length in assigning indices
1350               //"CERE-OFF";
1351               fprintf(pFlukaVmcInp,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one);
1352           }
1353           else  {
1354               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
1355               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1356           }
1357       } // end of else if (strncmp(proc->GetName(),"CKOV",4) == 0)
1358
1359       // Compton scattering
1360       // G3 default value: 1
1361       // G4 processes: G4ComptonScattering,
1362       //               G4LowEnergyCompton,
1363       //               G4PolarizedComptonScattering
1364       // Particles: gamma
1365       // Physics:   EM
1366       // flag = 0 no Compton scattering
1367       // flag = 1 Compton scattering, electron processed
1368       // flag = 2 Compton scattering, no electron stored
1369       // gMC ->SetProcess("COMP",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. PHOT-THR
1370       else if (strncmp(proc->GetName(),"COMP",4) == 0) {
1371           if (proc->Flag() == 1 || proc->Flag() == 2) {
1372               fprintf(pFlukaVmcInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
1373               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',1);\n");
1374               // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
1375               // zero = not used
1376               // zero = not used
1377               // matMin = lower bound of the material indices in which the respective thresholds apply
1378               // matMax = upper bound of the material indices in which the respective thresholds apply
1379               // one = step length in assigning indices
1380               //"PHOT-THR";
1381               theCut = -1.0;
1382               nextc.Reset();
1383               while ((cut = (TFlukaConfigOption*)nextc())) {
1384                   if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
1385                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1386               }
1387               fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",theCut,zero,zero,matMin,matMax,one);
1388           }
1389           else if (proc->Flag() == 0) {
1390               fprintf(pFlukaVmcInp,"*\n*No Compton scattering - no FLUKA card generated\n");
1391               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',0)\n");
1392           }
1393           else  {
1394               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
1395               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1396           }
1397       } // end of else if (strncmp(proc->GetName(),"COMP",4) == 0)
1398
1399       // decay
1400       // G3 default value: 1
1401       // G4 process: G4Decay
1402       //
1403       // Particles: all which decay is applicable for
1404       // Physics:   General
1405       // flag = 0 no decays
1406       // flag = 1 decays, secondaries processed
1407       // flag = 2 decays, no secondaries stored
1408       //gMC ->SetProcess("DCAY",0); // not available
1409       else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 0)
1410           cout << "SetProcess for flag =" << proc->GetName() << " value=" << proc->Flag() << " not avaliable!" << endl;
1411       else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 1) {
1412           // Nothing to do decays are switched on by default
1413       }
1414
1415
1416       // delta-ray
1417       // G3 default value: 2
1418       // !! G4 treats delta rays in different way
1419       // G4 processes: G4eIonisation/G4IeIonization,
1420       //               G4MuIonisation/G4IMuIonization,
1421       //               G4hIonisation/G4IhIonisation
1422       // Particles: charged
1423       // Physics:   EM
1424       // flag = 0 no energy loss
1425       // flag = 1 restricted energy loss fluctuations
1426       // flag = 2 complete energy loss fluctuations
1427       // flag = 3 same as 1
1428       // flag = 4 no energy loss fluctuations
1429       // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0.  0. 3. lastmat 0.
1430       else if (strncmp(proc->GetName(),"DRAY",4) == 0) {
1431           if (proc->Flag() == 0 || proc->Flag() == 4) {
1432               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1433               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
1434               fprintf(pFlukaVmcInp,"*No delta ray production by muons - threshold set artificially high\n");
1435               Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1436               // zero = ignored
1437               // zero = ignored
1438               // matMin = lower bound of the material indices in which the respective thresholds apply
1439               // matMax = upper bound of the material indices in which the respective thresholds apply
1440               // one = step length in assigning indices
1441               fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one);
1442           }
1443           else if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
1444               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1445               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
1446               fprintf(pFlukaVmcInp,"*Delta ray production by muons switched on\n");
1447               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
1448               theCut = 1.0e+6;
1449               nextc.Reset();
1450               //
1451               // Check cut one delta-rays from electrons
1452               //
1453               while ((cut = (TFlukaConfigOption*)nextc())) {
1454                   if (strncmp(cut->GetName(), "DCUTM", 5) == 0 &&
1455                       cut->Medium() == proc->Medium()) theCut = cut->Cut();
1456               }
1457               // theCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1458               // zero = ignored
1459               // zero = ignored
1460               // matMin = lower bound of the material indices in which the respective thresholds apply
1461               // matMax =  upper bound of the material indices in which the respective thresholds apply
1462               // one = step length in assigning indices
1463               fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",theCut,zero,zero,matMin,matMax,one);
1464           }
1465           else  {
1466               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
1467               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1468           }
1469       } // end of else if (strncmp(proc->GetName(),"DRAY",4) == 0)
1470
1471       // hadronic process
1472       // G3 default value: 1
1473       // G4 processes: all defined by TG4PhysicsConstructorHadron
1474       //
1475       // Particles: hadrons
1476       // Physics:   Hadron
1477       // flag = 0 no multiple scattering
1478       // flag = 1 hadronic interactions, secondaries processed
1479       // flag = 2 hadronic interactions, no secondaries stored
1480       // gMC ->SetProcess("HADR",1); // ??? hadronic process
1481       //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
1482       else if (strncmp(proc->GetName(),"HADR",4) == 0) {
1483           if (proc->Flag() == 1 || proc->Flag() == 2) {
1484               fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
1485               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1486           }
1487           else if (proc->Flag() == 0) {
1488               fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is set OFF\n");
1489               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('HADR',0);\n");
1490               fprintf(pFlukaVmcInp,"*Switching off hadronic interactions not foreseen in FLUKA\n");
1491               fprintf(pFlukaVmcInp,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero);
1492           }
1493           else  {
1494               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
1495               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1496           }
1497       } // end of else if (strncmp(proc->GetName(),"HADR",4) == 0)
1498
1499
1500       // energy loss
1501       // G3 default value: 2
1502       // G4 processes: G4eIonisation/G4IeIonization,
1503       //               G4MuIonisation/G4IMuIonization,
1504       //               G4hIonisation/G4IhIonisation
1505       //
1506       // Particles: charged
1507       // Physics:   EM
1508       // flag=0 no energy loss
1509       // flag=1 restricted energy loss fluctuations
1510       // flag=2 complete energy loss fluctuations
1511       // flag=3 same as 1
1512       // flag=4 no energy loss fluctuations
1513       // If the value ILOSS is changed, then (in G3) cross-sections and energy
1514       // loss tables must be recomputed via the command 'PHYSI'
1515       // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1516       else if (strncmp(proc->GetName(),"LOSS",4) == 0) {
1517           if (proc->Flag() > 0 || proc->Flag() < 4) { // restricted energy loss fluctuations
1518               fprintf(pFlukaVmcInp,"*\n*Restricted energy loss fluctuations\n");
1519               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
1520               // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1521               // one = restricted energy loss fluctuations (for e+ and e-) switched on
1522               // one = minimal accuracy
1523               // matMin = lower bound of the material indices in which the respective thresholds apply
1524               // upper bound of the material indices in which the respective thresholds apply
1525               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one, one, 4., matMin, matMax);
1526           }
1527           else if (proc->Flag() == 4) { // no energy loss fluctuations
1528               fprintf(pFlukaVmcInp,"*\n*No energy loss fluctuations\n");
1529               fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
1530               // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
1531               // - one = restricted energy loss fluctuations (for e+ and e-) switched off
1532               // one = minimal accuracy
1533               // matMin = lower bound of the material indices in which the respective thresholds apply
1534               // matMax = upper bound of the material indices in which the respective thresholds apply
1535               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax);
1536           }
1537           else  {
1538               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
1539               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1540           }
1541       } // end of else if (strncmp(proc->GetName(),"LOSS",4) == 0)
1542
1543
1544       // multiple scattering
1545       // G3 default value: 1
1546       // G4 process: G4MultipleScattering/G4IMultipleScattering
1547       //
1548       // Particles: charged
1549       // Physics:   EM
1550       // flag = 0 no multiple scattering
1551       // flag = 1 Moliere or Coulomb scattering
1552       // flag = 2 Moliere or Coulomb scattering
1553       // flag = 3 Gaussian scattering
1554       // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
1555       else if (strncmp(proc->GetName(),"MULS",4) == 0) {
1556           if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
1557               fprintf(pFlukaVmcInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
1558               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1559           }
1560           else if (proc->Flag() == 0) {
1561               fprintf(pFlukaVmcInp,"*\n*Multiple scattering is set OFF\n");
1562               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MULS',0);\n");
1563               // zero = ignored
1564               // three = multiple scattering for hadrons and muons is completely suppressed
1565               // three = multiple scattering for e+ and e- is completely suppressed
1566               // matMin = lower bound of the material indices in which the respective thresholds apply
1567               // matMax = upper bound of the material indices in which the respective thresholds apply
1568               fprintf(pFlukaVmcInp,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax);
1569           }
1570           else  {
1571               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
1572               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1573           }
1574       } // end of else if (strncmp(proc->GetName(),"MULS",4) == 0)
1575
1576
1577       // muon nuclear interaction
1578       // G3 default value: 0
1579       // G4 processes: G4MuNuclearInteraction,
1580       // G4MuonMinusCaptureAtRest
1581       //
1582       // Particles: mu
1583       // Physics:   Not set
1584       // flag = 0 no muon-nuclear interaction
1585       // flag = 1 nuclear interaction, secondaries processed
1586       // flag = 2 nuclear interaction, secondaries not processed
1587       // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
1588       else if (strncmp(proc->GetName(),"MUNU",4) == 0) {
1589           if (proc->Flag() == 1) {
1590               fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
1591               fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('MUNU',1);\n");
1592               // one = full simulation of muon nuclear interactions and production of secondary hadrons
1593               // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1594               // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1595               // matMin = lower bound of the material indices in which the respective thresholds apply
1596               // matMax = upper bound of the material indices in which the respective thresholds apply
1597               fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1598           }
1599           else if (proc->Flag() == 2) {
1600               fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
1601               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',2);\n");
1602               // two = full simulation of muon nuclear interactions and production of secondary hadrons
1603               // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1604               // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1605               // matMin = lower bound of the material indices in which the respective thresholds apply
1606               // matMax = upper bound of the material indices in which the respective thresholds apply
1607               fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax);
1608           }
1609           else if (proc->Flag() == 0) {
1610               fprintf(pFlukaVmcInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
1611               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',0)\n");
1612           }
1613           else  {
1614               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
1615               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1616           }
1617       } // end of else if (strncmp(proc->GetName(),"MUNU",4) == 0)
1618
1619
1620       // photofission
1621       // G3 default value: 0
1622       // G4 process: ??
1623       //
1624       // Particles: gamma
1625       // Physics:   ??
1626       // gMC ->SetProcess("PFIS",0); // PHOTONUC -1.   0.  0. 3. lastmat 0.
1627       // flag = 0 no photon fission
1628       // flag = 1 photon fission, secondaries processed
1629       // flag = 2 photon fission, no secondaries stored
1630       else if (strncmp(proc->GetName(),"PFIS",4) == 0) {
1631           if (proc->Flag() == 0) {
1632               fprintf(pFlukaVmcInp,"*\n*No photonuclear interactions\n");
1633               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0);\n");
1634               // - one = no photonuclear interactions
1635               // zero = not used
1636               // zero = not used
1637               // matMin = lower bound of the material indices in which the respective thresholds apply
1638               // matMax = upper bound of the material indices in which the respective thresholds apply
1639               fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax);
1640           }
1641           else if (proc->Flag() == 1) {
1642               fprintf(pFlukaVmcInp,"*\n*Photon nuclear interactions are activated at all energies\n");
1643               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',1);\n");
1644               // one = photonuclear interactions are activated at all energies
1645               // zero = not used
1646               // zero = not used
1647               // matMin = lower bound of the material indices in which the respective thresholds apply
1648               // matMax = upper bound of the material indices in which the respective thresholds apply
1649               fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1650           }
1651           else if (proc->Flag() == 0) {
1652               fprintf(pFlukaVmcInp,"*\n*No photofission - no FLUKA card generated\n");
1653               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0)\n");
1654           }
1655           else {
1656               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
1657               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1658           }
1659       }
1660
1661
1662       // photo electric effect
1663       // G3 default value: 1
1664       // G4 processes: G4PhotoElectricEffect
1665       //               G4LowEnergyPhotoElectric
1666       // Particles: gamma
1667       // Physics:   EM
1668       // flag = 0 no photo electric effect
1669       // flag = 1 photo electric effect, electron processed
1670       // flag = 2 photo electric effect, no electron stored
1671       // gMC ->SetProcess("PHOT",1); // EMFCUT    0.  -1.  0. 3. lastmat 0. PHOT-THR
1672       else if (strncmp(proc->GetName(),"PHOT",4) == 0) {
1673           if (proc->Flag() == 1) {
1674               fprintf(pFlukaVmcInp,"*\n*Photo electric effect is activated\n");
1675               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',1);\n");
1676               // zero = ignored
1677               // - one = resets to default=0.
1678               // zero = ignored
1679               // matMin = lower bound of the material indices in which the respective thresholds apply
1680               // matMax = upper bound of the material indices in which the respective thresholds apply
1681               // one = step length in assigning indices
1682               //"PHOT-THR";
1683               theCut = -1.0;
1684               nextc.Reset();
1685               while ((cut = (TFlukaConfigOption*)nextc())) {
1686                   if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
1687                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
1688               }
1689               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,theCut,zero,matMin,matMax,one);
1690           }
1691           else if (proc->Flag() == 0) {
1692               fprintf(pFlukaVmcInp,"*\n*No photo electric effect - no FLUKA card generated\n");
1693               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',0)\n");
1694           }
1695           else {
1696               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
1697               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1698           }
1699       } // else if (strncmp(proc->GetName(),"PHOT",4) == 0)
1700
1701
1702       // Rayleigh scattering
1703       // G3 default value: 0
1704       // G4 process: G4OpRayleigh
1705       //
1706       // Particles: optical photon
1707       // Physics:   Optical
1708       // flag = 0 Rayleigh scattering off
1709       // flag = 1 Rayleigh scattering on
1710       //xx gMC ->SetProcess("RAYL",1);
1711       else if (strncmp(proc->GetName(),"RAYL",4) == 0) {
1712           if (proc->Flag() == 1) {
1713               fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
1714               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1715           }
1716           else if (proc->Flag() == 0) {
1717               fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is set OFF\n");
1718               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('RAYL',0);\n");
1719               // - one = no Rayleigh scattering and no binding corrections for Compton
1720               // matMin = lower bound of the material indices in which the respective thresholds apply
1721               // matMax = upper bound of the material indices in which the respective thresholds apply
1722               fprintf(pFlukaVmcInp,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax);
1723           }
1724           else  {
1725               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
1726               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1727           }
1728       } // end of else if (strncmp(proc->GetName(),"RAYL",4) == 0)
1729
1730
1731       // synchrotron radiation in magnetic field
1732       // G3 default value: 0
1733       // G4 process: G4SynchrotronRadiation
1734       //
1735       // Particles: ??
1736       // Physics:   Not set
1737       // flag = 0 no synchrotron radiation
1738       // flag = 1 synchrotron radiation
1739       //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
1740       else if (strncmp(proc->GetName(),"SYNC",4) == 0) {
1741           fprintf(pFlukaVmcInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
1742           fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1743       }
1744
1745
1746       // Automatic calculation of tracking medium parameters
1747       // flag = 0 no automatic calculation
1748       // flag = 1 automatic calculation
1749       //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
1750       else if (strncmp(proc->GetName(),"AUTO",4) == 0) {
1751           fprintf(pFlukaVmcInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
1752           fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1753       }
1754
1755
1756       // To control energy loss fluctuation model
1757       // flag = 0 Urban model
1758       // flag = 1 PAI model
1759       // flag = 2 PAI+ASHO model (not active at the moment)
1760       //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1761       else if (strncmp(proc->GetName(),"STRA",4) == 0) {
1762           if (proc->Flag() == 0 || proc->Flag() == 2 || proc->Flag() == 3) {
1763               fprintf(pFlukaVmcInp,"*\n*Ionization energy losses calculation is activated\n");
1764               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
1765               // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1766               // one = restricted energy loss fluctuations (for e+ and e-) switched on
1767               // one = minimal accuracy
1768               // matMin = lower bound of the material indices in which the respective thresholds apply
1769               // matMax = upper bound of the material indices in which the respective thresholds apply
1770               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
1771           }
1772           else {
1773               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
1774               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1775           }
1776       } // else if (strncmp(proc->GetName(),"STRA",4) == 0)
1777
1778
1779
1780
1781       else { // processes not yet treated
1782
1783           // light photon absorption (Cerenkov photons)
1784           // it is turned on when Cerenkov process is turned on
1785           // G3 default value: 0
1786           // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1787           //
1788           // Particles: optical photon
1789           // Physics:   Optical
1790           // flag = 0 no absorption of Cerenkov photons
1791           // flag = 1 absorption of Cerenkov photons
1792           // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1793
1794
1795
1796           cout << "SetProcess for flag=" << proc->GetName() << " value=" << proc->Flag() << " not yet implemented!" << endl;
1797       }
1798   } //end of loop number of SetProcess calls
1799
1800
1801 // Loop over number of SetCut calls
1802
1803   nextc.Reset();
1804   while ((cut = (TFlukaConfigOption*)nextc())) {
1805       Float_t matMin = three;
1806       Float_t matMax = fLastMaterial;
1807       Bool_t global  = kTRUE;
1808       if (cut->Medium() != -1) {
1809           Int_t mat;
1810           if ((mat = cut->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
1811           matMin = Float_t(mat);
1812           matMax = matMin;
1813           global = kFALSE;
1814           TGeoMaterial*    material =  (TGeoMaterial*) (GetFlukaMaterials())->At(GetMaterialIndex(mat));
1815           fprintf(pFlukaVmcInp,"*\n*Material specific cut setting for #%8d %s %s %13.3e\n",
1816                   mat, material->GetName(), cut->GetName(), cut->Cut());
1817
1818       }
1819
1820       // cuts handled in SetProcess calls
1821       if      (strncmp(cut->GetName(),"BCUTM",5)  == 0) continue;
1822       else if (strncmp(cut->GetName(),"BCUTE",5)  == 0) continue;
1823       else if (strncmp(cut->GetName(),"DCUTM",5)  == 0) continue;
1824       else if (strncmp(cut->GetName(),"PPCUTM",6) == 0) continue;
1825
1826       // gammas
1827       // G4 particles: "gamma"
1828       // G3 default value: 0.001 GeV
1829       // gMC ->SetCut("CUTGAM",cut); // cut for gammas
1830
1831       else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && global) {
1832           fprintf(pFlukaVmcInp,"*\n*Cut for gamma\n");
1833           fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1834           fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
1835                   zero, cut->Cut(), zero, zero, Float_t(fGeom->NofVolumes()), one);
1836       }
1837       else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) {
1838           // loop over materials for EMFCUT FLUKA cards
1839           for (j=0; j < matMax-matMin+1; j++) {
1840               Int_t nreg, imat, *reglist;
1841               Float_t ireg;
1842               imat = (Int_t) matMin + j;
1843               reglist = fGeom->GetMaterialList(imat, nreg);
1844               // loop over regions of a given material
1845               for (Int_t k = 0; k < nreg; k++) {
1846                   ireg = reglist[k];
1847                   fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, cut->Cut(), zero, ireg, ireg, one);
1848               }
1849           }
1850       } // end of else if for gamma
1851
1852
1853       // electrons
1854       // G4 particles: "e-"
1855       // ?? positrons
1856       // G3 default value: 0.001 GeV
1857       //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1858       else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && global) {
1859           fprintf(pFlukaVmcInp,"*\n*Cut for electrons\n");
1860           fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1861           fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
1862                   -cut->Cut(), zero, zero, zero, Float_t(fGeom->NofVolumes()), one);
1863       }
1864       else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) {
1865           // loop over materials for EMFCUT FLUKA cards
1866           for (j=0; j < matMax-matMin+1; j++) {
1867               Int_t nreg, imat, *reglist;
1868               Float_t ireg;
1869               imat = (Int_t) matMin + j;
1870               reglist = fGeom->GetMaterialList(imat, nreg);
1871               // loop over regions of a given material
1872               for (k=0; k<nreg; k++) {
1873                   ireg = reglist[k];
1874                   fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -cut->Cut(), zero, zero, ireg, ireg, one);
1875               }
1876           }
1877       } // end of else if for electrons
1878
1879
1880       // neutral hadrons
1881       // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1882       // G3 default value: 0.01 GeV
1883       //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1884       else if (strncmp(cut->GetName(),"CUTNEU",6) == 0 && global) {
1885         fprintf(pFlukaVmcInp,"*\n*Cut for neutral hadrons\n");
1886         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
1887
1888         // 8.0 = Neutron
1889         // 9.0 = Antineutron
1890         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),8.0,9.0);
1891
1892         // 12.0 = Kaon zero long
1893         // 12.0 = Kaon zero long
1894         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),12.0,12.0);
1895
1896         // 17.0 = Lambda, 18.0 = Antilambda
1897         // 19.0 = Kaon zero short
1898         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),17.0,19.0);
1899
1900         // 22.0 = Sigma zero, Pion zero, Kaon zero
1901         // 25.0 = Antikaon zero
1902         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0);
1903
1904         // 32.0 = Antisigma zero
1905         // 32.0 = Antisigma zero
1906         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0);
1907
1908         // 34.0 = Xi zero
1909         // 35.0 = AntiXi zero
1910         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0);
1911
1912         // 47.0 = D zero
1913         // 48.0 = AntiD zero
1914         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),47.0,48.0);
1915
1916         // 53.0 = Xi_c zero
1917         // 53.0 = Xi_c zero
1918         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),53.0,53.0);
1919
1920         // 55.0 = Xi'_c zero
1921         // 56.0 = Omega_c zero
1922         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),55.0,56.0);
1923
1924         // 59.0 = AntiXi_c zero
1925         // 59.0 = AntiXi_c zero
1926         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),59.0,59.0);
1927
1928         // 61.0 = AntiXi'_c zero
1929         // 62.0 = AntiOmega_c zero
1930         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),61.0,62.0);
1931       }
1932
1933       // charged hadrons
1934       // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1935       // G3 default value: 0.01 GeV
1936       //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1937       else if (strncmp(cut->GetName(),"CUTHAD",6) == 0 && global) {
1938         fprintf(pFlukaVmcInp,"*\n*Cut for charged hadrons\n");
1939         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
1940
1941         // 1.0 = Proton
1942         // 2.0 = Antiproton
1943         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),1.0,2.0);
1944
1945         // 13.0 = Positive Pion, Negative Pion, Positive Kaon
1946         // 16.0 = Negative Kaon
1947         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0);
1948
1949         // 20.0 = Negative Sigma
1950         // 21.0 = Positive Sigma
1951         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),20.0,21.0);
1952
1953         // 31.0 = Antisigma minus
1954         // 33.0 = Antisigma plus
1955         // 2.0 = step length
1956         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),31.0,33.0,2.0);
1957
1958         // 36.0 = Negative Xi, Positive Xi, Omega minus
1959         // 39.0 = Antiomega
1960         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0);
1961
1962         // 45.0 = D plus
1963         // 46.0 = D minus
1964         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),45.0,46.0);
1965
1966         // 49.0 = D_s plus, D_s minus, Lambda_c plus
1967         // 52.0 = Xi_c plus
1968         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),49.0,52.0);
1969
1970         // 54.0 = Xi'_c plus
1971         // 60.0 = AntiXi'_c minus
1972         // 6.0 = step length
1973         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),54.0,60.0,6.0);
1974
1975         // 57.0 = Antilambda_c minus
1976         // 58.0 = AntiXi_c minus
1977         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),57.0,58.0);
1978       }
1979
1980       // muons
1981       // G4 particles: "mu+", "mu-"
1982       // G3 default value: 0.01 GeV
1983       //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1984       else if (strncmp(cut->GetName(),"CUTMUO",6)== 0 && global) {
1985         fprintf(pFlukaVmcInp,"*\n*Cut for muons\n");
1986         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
1987         // 10.0 = Muon+
1988         // 11.0 = Muon-
1989         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),10.0,11.0);
1990       }
1991
1992       //
1993       // time of flight cut in seconds
1994       // G4 particles: all
1995       // G3 default value: 0.01 GeV
1996       //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1997       else if (strncmp(cut->GetName(),"TOFMAX",6) == 0) {
1998         fprintf(pFlukaVmcInp,"*\n*Time of flight cuts in seconds\n");
1999         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
2000         // zero = ignored
2001         // zero = ignored
2002         // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
2003         // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
2004         fprintf(pFlukaVmcInp,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut->Cut()*1.e9,zero,zero,-6.0,64.0);
2005       }
2006
2007       else if (global){
2008         cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " not yet implemented!" << endl;
2009       }
2010       else {
2011         cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " (material specific) not yet implemented!" << endl;
2012       }
2013
2014   } //end of loop over SetCut calls
2015
2016 // Add START and STOP card
2017   fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
2018   fprintf(pFlukaVmcInp,"STOP      \n");
2019
2020
2021 // Close files
2022
2023    fclose(pFlukaVmcCoreInp);
2024    fclose(pFlukaVmcFlukaMat);
2025    fclose(pFlukaVmcInp);
2026
2027 } // end of InitPhysics
2028
2029
2030 //______________________________________________________________________________
2031 void TFluka::SetMaxStep(Double_t step)
2032 {
2033 // Set the maximum step size
2034     if (step > 1.e4) return;
2035
2036     Int_t mreg, latt;
2037     fGeom->GetCurrentRegion(mreg, latt);
2038     STEPSZ.stepmx[mreg - 1] = step;
2039 }
2040
2041
2042 Double_t TFluka::MaxStep() const
2043 {
2044 // Return the maximum for current medium
2045     Int_t mreg, latt;
2046     fGeom->GetCurrentRegion(mreg, latt);
2047     return (STEPSZ.stepmx[mreg - 1]);
2048 }
2049
2050 //______________________________________________________________________________
2051 void TFluka::SetMaxNStep(Int_t)
2052 {
2053 // SetMaxNStep is dummy procedure in TFluka !
2054   if (fVerbosityLevel >=3)
2055   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
2056 }
2057
2058 //______________________________________________________________________________
2059 void TFluka::SetUserDecay(Int_t)
2060 {
2061 // SetUserDecay is dummy procedure in TFluka !
2062   if (fVerbosityLevel >=3)
2063   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
2064 }
2065
2066 //
2067 // dynamic properties
2068 //
2069 //______________________________________________________________________________
2070 void TFluka::TrackPosition(TLorentzVector& position) const
2071 {
2072 // Return the current position in the master reference frame of the
2073 // track being transported
2074 // TRACKR.atrack = age of the particle
2075 // TRACKR.xtrack = x-position of the last point
2076 // TRACKR.ytrack = y-position of the last point
2077 // TRACKR.ztrack = z-position of the last point
2078   Int_t caller = GetCaller();
2079   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
2080     position.SetX(GetXsco());
2081     position.SetY(GetYsco());
2082     position.SetZ(GetZsco());
2083     position.SetT(TRACKR.atrack);
2084   }
2085   else if (caller == 4) { // mgdraw
2086     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
2087     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
2088     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
2089     position.SetT(TRACKR.atrack);
2090   }
2091   else if (caller == 5) { // sodraw
2092     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
2093     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
2094     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
2095     position.SetT(0);
2096   }
2097   else
2098     Warning("TrackPosition","position not available");
2099 }
2100
2101 //______________________________________________________________________________
2102 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
2103 {
2104 // Return the current position in the master reference frame of the
2105 // track being transported
2106 // TRACKR.atrack = age of the particle
2107 // TRACKR.xtrack = x-position of the last point
2108 // TRACKR.ytrack = y-position of the last point
2109 // TRACKR.ztrack = z-position of the last point
2110   Int_t caller = GetCaller();
2111   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
2112     x = GetXsco();
2113     y = GetYsco();
2114     z = GetZsco();
2115   }
2116   else if (caller == 4 || caller == 5) { // mgdraw, sodraw
2117     x = TRACKR.xtrack[TRACKR.ntrack];
2118     y = TRACKR.ytrack[TRACKR.ntrack];
2119     z = TRACKR.ztrack[TRACKR.ntrack];
2120   }
2121   else
2122     Warning("TrackPosition","position not available");
2123 }
2124
2125 //______________________________________________________________________________
2126 void TFluka::TrackMomentum(TLorentzVector& momentum) const
2127 {
2128 // Return the direction and the momentum (GeV/c) of the track
2129 // currently being transported
2130 // TRACKR.ptrack = momentum of the particle (not always defined, if
2131 //               < 0 must be obtained from etrack)
2132 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2133 // TRACKR.etrack = total energy of the particle
2134 // TRACKR.jtrack = identity number of the particle
2135 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2136   Int_t caller = GetCaller();
2137   if (caller != 2) { // not eedraw
2138     if (TRACKR.ptrack >= 0) {
2139       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
2140       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
2141       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
2142       momentum.SetE(TRACKR.etrack);
2143       return;
2144     }
2145     else {
2146       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2147       momentum.SetPx(p*TRACKR.cxtrck);
2148       momentum.SetPy(p*TRACKR.cytrck);
2149       momentum.SetPz(p*TRACKR.cztrck);
2150       momentum.SetE(TRACKR.etrack);
2151       return;
2152     }
2153   }
2154   else
2155     Warning("TrackMomentum","momentum not available");
2156 }
2157
2158 //______________________________________________________________________________
2159 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
2160 {
2161 // Return the direction and the momentum (GeV/c) of the track
2162 // currently being transported
2163 // TRACKR.ptrack = momentum of the particle (not always defined, if
2164 //               < 0 must be obtained from etrack)
2165 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2166 // TRACKR.etrack = total energy of the particle
2167 // TRACKR.jtrack = identity number of the particle
2168 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2169   Int_t caller = GetCaller();
2170   if (caller != 2) { // not eedraw
2171     if (TRACKR.ptrack >= 0) {
2172       px = TRACKR.ptrack*TRACKR.cxtrck;
2173       py = TRACKR.ptrack*TRACKR.cytrck;
2174       pz = TRACKR.ptrack*TRACKR.cztrck;
2175       e = TRACKR.etrack;
2176       return;
2177     }
2178     else {
2179       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2180       px = p*TRACKR.cxtrck;
2181       py = p*TRACKR.cytrck;
2182       pz = p*TRACKR.cztrck;
2183       e = TRACKR.etrack;
2184       return;
2185     }
2186   }
2187   else
2188     Warning("TrackMomentum","momentum not available");
2189 }
2190
2191 //______________________________________________________________________________
2192 Double_t TFluka::TrackStep() const
2193 {
2194 // Return the length in centimeters of the current step
2195 // TRACKR.ctrack = total curved path
2196   Int_t caller = GetCaller();
2197   if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
2198     return 0.0;
2199   else if (caller == 4) //mgdraw
2200     return TRACKR.ctrack;
2201   else
2202     return -1.0;
2203 }
2204
2205 //______________________________________________________________________________
2206 Double_t TFluka::TrackLength() const
2207 {
2208 // TRACKR.cmtrck = cumulative curved path since particle birth
2209   Int_t caller = GetCaller();
2210   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2211     return TRACKR.cmtrck;
2212   else
2213     return -1.0;
2214 }
2215
2216 //______________________________________________________________________________
2217 Double_t TFluka::TrackTime() const
2218 {
2219 // Return the current time of flight of the track being transported
2220 // TRACKR.atrack = age of the particle
2221   Int_t caller = GetCaller();
2222   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2223     return TRACKR.atrack;
2224   else
2225     return -1;
2226 }
2227
2228 //______________________________________________________________________________
2229 Double_t TFluka::Edep() const
2230 {
2231 // Energy deposition
2232 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
2233 // -->local energy deposition (the value and the point are not recorded in TRACKR)
2234 //    but in the variable "rull" of the procedure "endraw.cxx"
2235 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
2236 // -->no energy loss along the track
2237 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
2238 // -->energy loss distributed along the track
2239 // TRACKR.dtrack = energy deposition of the jth deposition event
2240
2241   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
2242   Int_t caller = GetCaller();
2243   if (caller == 11 || caller==12) return 0.0;
2244   Double_t sum = 0;
2245   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
2246     sum +=TRACKR.dtrack[j];
2247   }
2248   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
2249     return fRull + sum;
2250   else {
2251     return sum;
2252   }
2253 }
2254
2255 //______________________________________________________________________________
2256 Int_t TFluka::TrackPid() const
2257 {
2258 // Return the id of the particle transported
2259 // TRACKR.jtrack = identity number of the particle
2260   Int_t caller = GetCaller();
2261   if (caller != 2) { // not eedraw
2262       return PDGFromId(TRACKR.jtrack);
2263   }
2264   else
2265     return -1000;
2266 }
2267
2268 //______________________________________________________________________________
2269 Double_t TFluka::TrackCharge() const
2270 {
2271 // Return charge of the track currently transported
2272 // PAPROP.ichrge = electric charge of the particle
2273 // TRACKR.jtrack = identity number of the particle
2274   Int_t caller = GetCaller();
2275   if (caller != 2)  // not eedraw
2276     return PAPROP.ichrge[TRACKR.jtrack+6];
2277   else
2278     return -1000.0;
2279 }
2280
2281 //______________________________________________________________________________
2282 Double_t TFluka::TrackMass() const
2283 {
2284 // PAPROP.am = particle mass in GeV
2285 // TRACKR.jtrack = identity number of the particle
2286   Int_t caller = GetCaller();
2287   if (caller != 2)  // not eedraw
2288     return PAPROP.am[TRACKR.jtrack+6];
2289   else
2290     return -1000.0;
2291 }
2292
2293 //______________________________________________________________________________
2294 Double_t TFluka::Etot() const
2295 {
2296 // TRACKR.etrack = total energy of the particle
2297   Int_t caller = GetCaller();
2298   if (caller != 2)  // not eedraw
2299     return TRACKR.etrack;
2300   else
2301     return -1000.0;
2302 }
2303
2304 //
2305 // track status
2306 //
2307 //______________________________________________________________________________
2308 Bool_t   TFluka::IsNewTrack() const
2309 {
2310 // Return true for the first call of Stepping()
2311    return fTrackIsNew;
2312 }
2313
2314 void     TFluka::SetTrackIsNew(Bool_t flag)
2315 {
2316 // Return true for the first call of Stepping()
2317    fTrackIsNew = flag;
2318
2319 }
2320
2321
2322 //______________________________________________________________________________
2323 Bool_t   TFluka::IsTrackInside() const
2324 {
2325 // True if the track is not at the boundary of the current volume
2326 // In Fluka a step is always inside one kind of material
2327 // If the step would go behind the region of one material,
2328 // it will be shortened to reach only the boundary.
2329 // Therefore IsTrackInside() is always true.
2330   Int_t caller = GetCaller();
2331   if (caller == 11 || caller==12)  // bxdraw
2332     return 0;
2333   else
2334     return 1;
2335 }
2336
2337 //______________________________________________________________________________
2338 Bool_t   TFluka::IsTrackEntering() const
2339 {
2340 // True if this is the first step of the track in the current volume
2341
2342   Int_t caller = GetCaller();
2343   if (caller == 11)  // bxdraw entering
2344     return 1;
2345   else return 0;
2346 }
2347
2348 //______________________________________________________________________________
2349 Bool_t   TFluka::IsTrackExiting() const
2350 {
2351 // True if track is exiting volume
2352 //
2353   Int_t caller = GetCaller();
2354   if (caller == 12)  // bxdraw exiting
2355     return 1;
2356   else return 0;
2357 }
2358
2359 //______________________________________________________________________________
2360 Bool_t   TFluka::IsTrackOut() const
2361 {
2362 // True if the track is out of the setup
2363 // means escape
2364 // Icode = 14: escape - call from Kaskad
2365 // Icode = 23: escape - call from Emfsco
2366 // Icode = 32: escape - call from Kasneu
2367 // Icode = 40: escape - call from Kashea
2368 // Icode = 51: escape - call from Kasoph
2369   if (fIcode == 14 ||
2370       fIcode == 23 ||
2371       fIcode == 32 ||
2372       fIcode == 40 ||
2373       fIcode == 51) return 1;
2374   else return 0;
2375 }
2376
2377 //______________________________________________________________________________
2378 Bool_t   TFluka::IsTrackDisappeared() const
2379 {
2380 // means all inelastic interactions and decays
2381 // fIcode from usdraw
2382   if (fIcode == 101 || // inelastic interaction
2383       fIcode == 102 || // particle decay
2384       fIcode == 103 || // delta ray generation by hadron
2385       fIcode == 104 || // direct pair production
2386       fIcode == 105 || // bremsstrahlung (muon)
2387       fIcode == 208 || // bremsstrahlung (electron)
2388       fIcode == 214 || // in-flight annihilation
2389       fIcode == 215 || // annihilation at rest
2390       fIcode == 217 || // pair production
2391       fIcode == 219 || // Compton scattering
2392       fIcode == 221 || // Photoelectric effect
2393       fIcode == 300 || // hadronic interaction
2394       fIcode == 400    // delta-ray
2395       ) return 1;
2396   else return 0;
2397 }
2398
2399 //______________________________________________________________________________
2400 Bool_t   TFluka::IsTrackStop() const
2401 {
2402 // True if the track energy has fallen below the threshold
2403 // means stopped by signal or below energy threshold
2404 // Icode = 12: stopping particle       - call from Kaskad
2405 // Icode = 15: time kill               - call from Kaskad
2406 // Icode = 21: below threshold, iarg=1 - call from Emfsco
2407 // Icode = 22: below threshold, iarg=2 - call from Emfsco
2408 // Icode = 24: time kill               - call from Emfsco
2409 // Icode = 31: below threshold         - call from Kasneu
2410 // Icode = 33: time kill               - call from Kasneu
2411 // Icode = 41: time kill               - call from Kashea
2412 // Icode = 52: time kill               - call from Kasoph
2413   if (fIcode == 12 ||
2414       fIcode == 15 ||
2415       fIcode == 21 ||
2416       fIcode == 22 ||
2417       fIcode == 24 ||
2418       fIcode == 31 ||
2419       fIcode == 33 ||
2420       fIcode == 41 ||
2421       fIcode == 52) return 1;
2422   else return 0;
2423 }
2424
2425 //______________________________________________________________________________
2426 Bool_t   TFluka::IsTrackAlive() const
2427 {
2428 // means not disappeared or not out
2429   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
2430   else return 1;
2431 }
2432
2433 //
2434 // secondaries
2435 //
2436
2437 //______________________________________________________________________________
2438 Int_t TFluka::NSecondaries() const
2439
2440 {
2441 // Number of secondary particles generated in the current step
2442 // FINUC.np = number of secondaries except light and heavy ions
2443 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2444     Int_t caller = GetCaller();
2445     if (caller == 6)  // valid only after usdraw
2446         return FINUC.np + FHEAVY.npheav;
2447     else if (caller == 50) {
2448         // Cerenkov Photon production
2449         return fNCerenkov;
2450     }
2451     return 0;
2452 } // end of NSecondaries
2453
2454 //______________________________________________________________________________
2455 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2456                 TLorentzVector& position, TLorentzVector& momentum)
2457 {
2458 // Copy particles from secondary stack to vmc stack
2459 //
2460
2461     Int_t caller = GetCaller();
2462     if (caller == 6) {  // valid only after usdraw
2463         if (FINUC.np > 0) {
2464             // Hadronic interaction
2465             if (isec >= 0 && isec < FINUC.np) {
2466                 particleId = PDGFromId(FINUC.kpart[isec]);
2467                 position.SetX(fXsco);
2468                 position.SetY(fYsco);
2469                 position.SetZ(fZsco);
2470                 position.SetT(TRACKR.atrack);
2471                 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
2472                 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
2473                 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
2474                 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
2475             }
2476             else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
2477                 Int_t jsec = isec - FINUC.np;
2478                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2479                 position.SetX(fXsco);
2480                 position.SetY(fYsco);
2481                 position.SetZ(fZsco);
2482                 position.SetT(TRACKR.atrack);
2483                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2484                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2485                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2486                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2487                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2488                 else if (FHEAVY.tkheav[jsec] > 6)
2489                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2490             }
2491             else
2492                 Warning("GetSecondary","isec out of range");
2493         }
2494     } else if (caller == 50) {
2495         Int_t index = OPPHST.lstopp - isec;
2496         position.SetX(OPPHST.xoptph[index]);
2497         position.SetY(OPPHST.yoptph[index]);
2498         position.SetZ(OPPHST.zoptph[index]);
2499         position.SetT(OPPHST.agopph[index]);
2500         Double_t p = OPPHST.poptph[index];
2501
2502         momentum.SetPx(p * OPPHST.txopph[index]);
2503         momentum.SetPy(p * OPPHST.tyopph[index]);
2504         momentum.SetPz(p * OPPHST.tzopph[index]);
2505         momentum.SetE(p);
2506     }
2507     else
2508         Warning("GetSecondary","no secondaries available");
2509
2510 } // end of GetSecondary
2511
2512
2513 //______________________________________________________________________________
2514 TMCProcess TFluka::ProdProcess(Int_t) const
2515
2516 {
2517 // Name of the process that has produced the secondary particles
2518 // in the current step
2519
2520     Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
2521
2522     if      (fIcode == 102)                  return kPDecay;
2523     else if (fIcode == 104 || fIcode == 217) return kPPair;
2524     else if (fIcode == 219)                  return kPCompton;
2525     else if (fIcode == 221)                  return kPPhotoelectric;
2526     else if (fIcode == 105 || fIcode == 208) return kPBrem;
2527     else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
2528     else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
2529     else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
2530     else if (fIcode == 101)                  return kPHadronic;
2531     else if (fIcode == 101) {
2532         if (!mugamma)                        return kPHadronic;
2533         else if (TRACKR.jtrack == 7)         return kPPhotoFission;
2534         else return kPMuonNuclear;
2535     }
2536     else if (fIcode == 225)                  return kPRayleigh;
2537 // Fluka codes 100, 300 and 400 still to be investigasted
2538     else                                     return kPNoProcess;
2539 }
2540
2541
2542 Int_t TFluka::StepProcesses(TArrayI &proc) const
2543 {
2544   //
2545   // Return processes active in the current step
2546   //
2547     proc.Set(1);
2548     TMCProcess iproc;
2549     switch (fIcode) {
2550     case 15:
2551     case 24:
2552     case 33:
2553     case 41:
2554     case 52:
2555         iproc =  kPTOFlimit;
2556         break;
2557     case 12:
2558     case 14:
2559     case 21:
2560     case 22:
2561     case 23:
2562     case 31:
2563     case 32:
2564     case 40:
2565     case 51:
2566         iproc = kPStop;
2567         break;
2568     case 50:
2569         iproc = kPLightAbsorption;
2570         break;
2571     case 59:
2572         iproc = kPLightRefraction;
2573     case 20:
2574         iproc = kPPhotoelectric;
2575         break;
2576     default:
2577         iproc = ProdProcess(0);
2578     }
2579     proc[0] = iproc;
2580     return 1;
2581 }
2582 //______________________________________________________________________________
2583 Int_t TFluka::VolId2Mate(Int_t id) const
2584 {
2585 //
2586 // Returns the material number for a given volume ID
2587 //
2588    return fMCGeo->VolId2Mate(id);
2589 }
2590
2591 //______________________________________________________________________________
2592 const char* TFluka::VolName(Int_t id) const
2593 {
2594 //
2595 // Returns the volume name for a given volume ID
2596 //
2597    return fMCGeo->VolName(id);
2598 }
2599
2600 //______________________________________________________________________________
2601 Int_t TFluka::VolId(const Text_t* volName) const
2602 {
2603 //
2604 // Converts from volume name to volume ID.
2605 // Time consuming. (Only used during set-up)
2606 // Could be replaced by hash-table
2607 //
2608    return fMCGeo->VolId(volName);
2609 }
2610
2611 //______________________________________________________________________________
2612 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2613 {
2614 //
2615 // Return the logical id and copy number corresponding to the current fluka region
2616 //
2617   if (gGeoManager->IsOutside()) return 0;
2618   TGeoNode *node = gGeoManager->GetCurrentNode();
2619   copyNo = node->GetNumber();
2620   Int_t id = node->GetVolume()->GetNumber();
2621   return id;
2622 }
2623
2624 //______________________________________________________________________________
2625 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2626 {
2627 //
2628 // Return the logical id and copy number of off'th mother
2629 // corresponding to the current fluka region
2630 //
2631   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2632   if (off==0) return CurrentVolID(copyNo);
2633   TGeoNode *node = gGeoManager->GetMother(off);
2634   if (!node) return 0;
2635   copyNo = node->GetNumber();
2636   return node->GetVolume()->GetNumber();
2637 }
2638
2639 //______________________________________________________________________________
2640 const char* TFluka::CurrentVolName() const
2641 {
2642 //
2643 // Return the current volume name
2644 //
2645   if (gGeoManager->IsOutside()) return 0;
2646   return gGeoManager->GetCurrentVolume()->GetName();
2647 }
2648
2649 //______________________________________________________________________________
2650 const char* TFluka::CurrentVolOffName(Int_t off) const
2651 {
2652 //
2653 // Return the volume name of the off'th mother of the current volume
2654 //
2655   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2656   if (off==0) return CurrentVolName();
2657   TGeoNode *node = gGeoManager->GetMother(off);
2658   if (!node) return 0;
2659   return node->GetVolume()->GetName();
2660 }
2661
2662 //______________________________________________________________________________
2663 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
2664                       Float_t & dens, Float_t & radl, Float_t & absl) const
2665 {
2666 //
2667 //  Return the current medium number and material properties
2668 //
2669   Int_t copy;
2670   Int_t id  =  TFluka::CurrentVolID(copy);
2671   Int_t med =  TFluka::VolId2Mate(id);
2672   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
2673   TGeoMaterial*   mat = vol->GetMaterial();
2674   a    = mat->GetA();
2675   z    = mat->GetZ();
2676   dens = mat->GetDensity();
2677   radl = mat->GetRadLen();
2678   absl = mat->GetIntLen();
2679
2680   return med;
2681 }
2682
2683 //______________________________________________________________________________
2684 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2685 {
2686 // Transforms a position from the world reference frame
2687 // to the current volume reference frame.
2688 //
2689 //  Geant3 desription:
2690 //  ==================
2691 //       Computes coordinates XD (in DRS)
2692 //       from known coordinates XM in MRS
2693 //       The local reference system can be initialized by
2694 //         - the tracking routines and GMTOD used in GUSTEP
2695 //         - a call to GMEDIA(XM,NUMED)
2696 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2697 //             (inverse routine is GDTOM)
2698 //
2699 //        If IFLAG=1  convert coordinates
2700 //           IFLAG=2  convert direction cosinus
2701 //
2702 // ---
2703    Double_t xmL[3], xdL[3];
2704    Int_t i;
2705    for (i=0;i<3;i++) xmL[i]=xm[i];
2706    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2707    else            gGeoManager->MasterToLocalVect(xmL,xdL);
2708    for (i=0;i<3;i++) xd[i] = xdL[i];
2709 }
2710
2711 //______________________________________________________________________________
2712 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2713 {
2714    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2715    else            gGeoManager->MasterToLocalVect(xm,xd);
2716 }
2717
2718 //______________________________________________________________________________
2719 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2720 {
2721 // Transforms a position from the current volume reference frame
2722 // to the world reference frame.
2723 //
2724 //  Geant3 desription:
2725 //  ==================
2726 //  Computes coordinates XM (Master Reference System
2727 //  knowing the coordinates XD (Detector Ref System)
2728 //  The local reference system can be initialized by
2729 //    - the tracking routines and GDTOM used in GUSTEP
2730 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2731 //        (inverse routine is GMTOD)
2732 //
2733 //   If IFLAG=1  convert coordinates
2734 //      IFLAG=2  convert direction cosinus
2735 //
2736 // ---
2737    Double_t xmL[3], xdL[3];
2738    Int_t i;
2739    for (i=0;i<3;i++) xdL[i] = xd[i];
2740    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2741    else            gGeoManager->LocalToMasterVect(xdL,xmL);
2742    for (i=0;i<3;i++) xm[i]=xmL[i];
2743 }
2744
2745 //______________________________________________________________________________
2746 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2747 {
2748    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2749    else            gGeoManager->LocalToMasterVect(xd,xm);
2750 }
2751
2752 //______________________________________________________________________________
2753 TObjArray *TFluka::GetFlukaMaterials()
2754 {
2755    return fGeom->GetMatList();
2756 }
2757
2758 //______________________________________________________________________________
2759 void TFluka::SetMreg(Int_t l)
2760 {
2761 // Set current fluka region
2762    fCurrentFlukaRegion = l;
2763    fGeom->SetMreg(l);
2764 }
2765
2766
2767
2768
2769 TString TFluka::ParticleName(Int_t pdg) const
2770 {
2771     // Return particle name for particle with pdg code pdg.
2772     Int_t ifluka = IdFromPDG(pdg);
2773     return TString((CHPPRP.btype[ifluka+6]), 8);
2774 }
2775
2776
2777 Double_t TFluka::ParticleMass(Int_t pdg) const
2778 {
2779     // Return particle mass for particle with pdg code pdg.
2780     Int_t ifluka = IdFromPDG(pdg);
2781     return (PAPROP.am[ifluka+6]);
2782 }
2783
2784 Double_t TFluka::ParticleCharge(Int_t pdg) const
2785 {
2786     // Return particle charge for particle with pdg code pdg.
2787     Int_t ifluka = IdFromPDG(pdg);
2788     return Double_t(PAPROP.ichrge[ifluka+6]);
2789 }
2790
2791 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2792 {
2793     // Return particle lifetime for particle with pdg code pdg.
2794     Int_t ifluka = IdFromPDG(pdg);
2795     return (PAPROP.thalf[ifluka+6]);
2796 }
2797
2798 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2799 {
2800     // Retrieve particle properties for particle with pdg code pdg.
2801
2802     strcpy(name, ParticleName(pdg).Data());
2803     type   = ParticleMCType(pdg);
2804     mass   = ParticleMass(pdg);
2805     charge = ParticleCharge(pdg);
2806     tlife  = ParticleLifeTime(pdg);
2807 }
2808
2809
2810
2811 #define pushcerenkovphoton pushcerenkovphoton_
2812 #define usersteppingckv    usersteppingckv_
2813
2814
2815 extern "C" {
2816     void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2817                             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2818                             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2819     {
2820         //
2821         // Pushes one cerenkov photon to the stack
2822         //
2823
2824         TFluka* fluka =  (TFluka*) gMC;
2825         TVirtualMCStack* cppstack = fluka->GetStack();
2826         Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2827         cppstack->PushTrack(0, parent, 50000050,
2828                             px, py, pz, e,
2829                             vx, vy, vz, tof,
2830                             polx, poly, polz,
2831                             kPCerenkov, ntr, wgt, 0);
2832     }
2833
2834     void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
2835     {
2836         //
2837         // Calls stepping in order to signal cerenkov production
2838         //
2839         TFluka *fluka = (TFluka*)gMC;
2840         fluka->SetMreg(mreg);
2841         fluka->SetXsco(x);
2842         fluka->SetYsco(y);
2843         fluka->SetZsco(z);
2844         fluka->SetNCerenkov(nphot);
2845         fluka->SetCaller(50);
2846         printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f\n", nphot, mreg, x, y, z);
2847         (TVirtualMCApplication::Instance())->Stepping();
2848     }
2849 }
2850