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