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