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