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   Double_t zero, one, two, three;
 568   FILE *pAliceCoreInp, *pAliceFlukaMat, *pAliceInp;
 569
 570   zero  = 0.0;
 571   one   = 1.0;
 572   two   = 2.0;
 573   three = 3.0;
 574
 575   FGeometryInit* geominit = FGeometryInit::GetInstance();
 576   Float_t fLastMaterial = geominit->GetLastMaterialIndex();
 577   printf("   last FLUKA material is %g\n", fLastMaterial);
 578
 579 // construct file names
 580   TString sAliceCoreInp = getenv("ALICE_ROOT");
 581   sAliceCoreInp +="/TFluka/input/";
 582   TString sAliceTmp = "flukaMat.inp";
 583   TString sAliceInp = GetInputFileName();
 584   sAliceCoreInp += GetCoreInputFileName();
 585 /* open files */
 586   if ((pAliceCoreInp = fopen(sAliceCoreInp.Data(),"r")) == NULL) {
 587     printf("\nCannot open file %s\n",sAliceCoreInp.Data());
 588     exit(1);
 589   }
 590   if ((pAliceFlukaMat = fopen(sAliceTmp.Data(),"r")) == NULL) {
 591     printf("\nCannot open file %s\n",sAliceTmp.Data());
 592     exit(1);
 593   }
 594   if ((pAliceInp = fopen(sAliceInp.Data(),"w")) == NULL) {
 595     printf("\nCannot open file %s\n",sAliceInp.Data());
 596     exit(1);
 597   }
 598
 599 // copy core input file
 600   Char_t sLine[255];
 601   Float_t fEventsPerRun;
 602
 603   while ((fgets(sLine,255,pAliceCoreInp)) != NULL) {
 604     if (strncmp(sLine,"GEOEND",6) != 0)
 605       fprintf(pAliceInp,"%s",sLine); // copy until GEOEND card
 606     else {
 607       fprintf(pAliceInp,"GEOEND\n");   // add GEOEND card
 608       goto flukamat;
 609     }
 610   } // end of while until GEOEND card
 611
 612 flukamat:
 613   while ((fgets(sLine,255,pAliceFlukaMat)) != NULL) { // copy flukaMat.inp file
 614     fprintf(pAliceInp,"%s\n",sLine);
 615   }
 616
 617   while ((fgets(sLine,255,pAliceCoreInp)) != NULL) {
 618     if (strncmp(sLine,"START",5) != 0)
 619       fprintf(pAliceInp,"%s\n",sLine);
 620     else {
 621       sscanf(sLine+10,"%10f",&fEventsPerRun);
 622       goto fin;
 623     }
 624   } //end of while until START card
 625
 626 fin:
 627 // in G3 the process control values meaning can be different for
 628 // different processes, but for most of them is:
 629 //   0  process is not activated
 630 //   1  process is activated WITH generation of secondaries
 631 //   2  process is activated WITHOUT generation of secondaries
 632 // if process does not generate secondaries => 1 same as 2
 633 //
 634 // Exceptions:
 635 //   MULS:  also 3
 636 //   LOSS:  also 3, 4
 637 //   RAYL:  only 0,1
 638 //   HADR:  may be > 2
 639 //
 640
 641 // Loop over number of SetProcess calls
 642   fprintf(pAliceInp,"*----------------------------------------------------------------------------- \n");
 643   fprintf(pAliceInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
 644   fprintf(pAliceInp,"*----------------------------------------------------------------------------- \n");
 645   for (i=0; i<iNbOfProc; i++) {
 646
 647     // annihilation
 648     // G3 default value: 1
 649     // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
 650     // Particles: e+
 651     // Physics:   EM
 652     // flag = 0 no annihilation
 653     // flag = 1 annihilation, decays processed
 654     // flag = 2 annihilation, no decay product stored
 655     // gMC ->SetProcess("ANNI",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. ANNH-THR
 656     if (strncmp(&sProcessFlag[i][0],"ANNI",4) == 0) {
 657       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
 658         fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
 659         fprintf(pAliceInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
 660         // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
 661         // zero = not used
 662         // zero = not used
 663         // three = lower bound of the material indices in which the respective thresholds apply
 664         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 665         // one = step length in assigning indices
 666         // "ANNH-THR";
 667         fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,three,fLastMaterial,one);
 668       }
 669       else if (iProcessValue[i] == 0) {
 670         fprintf(pAliceInp,"*\n*No annihilation - no FLUKA card generated\n");
 671         fprintf(pAliceInp,"*Generated from call: SetProcess('ANNI',0)\n");
 672       }
 673       else  {
 674         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
 675         fprintf(pAliceInp,"*No FLUKA card generated\n");
 676       }
 677     }
 678
 679     // bremsstrahlung and pair production are both activated
 680     // G3 default value: 1
 681     // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
 682     //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
 683     //               G4LowEnergyBremstrahlung
 684     // Particles: e-/e+; mu+/mu-
 685     // Physics:   EM
 686     // flag = 0 no bremsstrahlung
 687     // flag = 1 bremsstrahlung, photon processed
 688     // flag = 2 bremsstrahlung, no photon stored
 689     // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
 690                                    // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
 691     // G3 default value: 1
 692     // G4 processes: G4GammaConversion,
 693     //               G4MuPairProduction/G4IMuPairProduction
 694     //               G4LowEnergyGammaConversion
 695     // Particles: gamma, mu
 696     // Physics:   EM
 697     // flag = 0 no delta rays
 698     // flag = 1 delta rays, secondaries processed
 699     // flag = 2 delta rays, no secondaries stored
 700     // gMC ->SetProcess("PAIR",1); // PAIRBREM  1.   0.  0. 3. lastmat
 701                                    // EMFCUT    0.   0. -1. 3. lastmat 0. PHOT-THR
 702     else if ((strncmp(&sProcessFlag[i][0],"PAIR",4) == 0) && (iProcessValue[i] == 1 || iProcessValue[i] == 2)) {
 703       for (j=0; j<iNbOfProc; j++) {
 704         if ((strncmp(&sProcessFlag[j][0],"BREM",4) == 0) && (iProcessValue[j] == 1 || iProcessValue[j] == 2)) {
 705           fprintf(pAliceInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
 706           fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
 707           fprintf(pAliceInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
 708           fprintf(pAliceInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
 709           // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
 710           fprintf(pAliceInp,"PAIRBREM  %10.1f",three);
 711           // direct pair production by muons
 712           // G4 particles: "e-", "e+"
 713           // G3 default value: 0.01 GeV
 714           //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
 715           fCut = 0.0;
 716           for (k=0; k<iNbOfCut; k++) {
 717             if (strncmp(&sCutFlag[k][0],"PPCUTM",6) == 0) fCut = fCutValue[k];
 718           }
 719           fprintf(pAliceInp,"%10.4g",fCut);
 720           // fCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
 721           // muon and hadron bremsstrahlung
 722           // G4 particles: "gamma"
 723           // G3 default value: CUTGAM=0.001 GeV
 724           //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
 725           fCut = 0.0;
 726           for (k=0; k<iNbOfCut; k++) {
 727             if (strncmp(&sCutFlag[k][0],"BCUTM",5) == 0) fCut = fCutValue[k];
 728           }
 729           fprintf(pAliceInp,"%10.4g%10.1f%10.1f\n",fCut,three,fLastMaterial);
 730           // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
 731           // three = lower bound of the material indices in which the respective thresholds apply
 732           // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 733
 734           // for e+ and e-
 735           fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
 736           fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1);\n");
 737           fCut = -1.0;
 738           for (k=0; k<iNbOfCut; k++) {
 739             if (strncmp(&sCutFlag[k][0],"BCUTE",5) == 0) fCut = fCutValue[k];
 740           }
 741           //fCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
 742           // zero = not used
 743           // zero = not used
 744           // three = lower bound of the material indices in which the respective thresholds apply
 745           // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 746           // one = step length in assigning indices
 747           // "ELPO-THR";
 748           fprintf(pAliceInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",fCut,zero,zero,three,fLastMaterial,one);
 749
 750           // for e+ and e-
 751           fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n");
 752           fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1);\n");
 753           fCut = -1.0;
 754           for (j=0; j<iNbOfCut; j++) {
 755             if (strncmp(&sCutFlag[j][0],"CUTGAM",6) == 0) fCut = fCutValue[j];
 756           }
 757           // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
 758           // three = lower bound of the material indices in which the respective thresholds apply
 759           // fLastMaterial =  upper bound of the material indices in which the respective thresholds apply
 760           // one = step length in assigning indices
 761           fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,three,fLastMaterial,one);
 762           goto BOTH;
 763         } // end of if for BREM
 764       } // end of loop for BREM
 765
 766       // only pair production by muons and charged hadrons is activated
 767       fprintf(pAliceInp,"*\n*Pair production by muons and charged hadrons is activated\n");
 768       fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
 769       fprintf(pAliceInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
 770       // direct pair production by muons
 771       // G4 particles: "e-", "e+"
 772       // G3 default value: 0.01 GeV
 773       //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
 774       // one = pair production by muons and charged hadrons is activated
 775       // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
 776       // zero = no explicit bremsstrahlung production is simulated
 777       // three = lower bound of the material indices in which the respective thresholds apply
 778       // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 779       fprintf(pAliceInp,"PAIRBREM  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,three,fLastMaterial);
 780
 781       // for e+ and e-
 782       fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n");
 783       fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
 784       fCut = -1.0;
 785       for (j=0; j<iNbOfCut; j++) {
 786         if (strncmp(&sCutFlag[j][0],"CUTGAM",6) == 0) fCut = fCutValue[j];
 787       }
 788       // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
 789       // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
 790       // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
 791       // three = lower bound of the material indices in which the respective thresholds apply
 792       // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 793       // one = step length in assigning indices
 794       fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,three,fLastMaterial,one);
 795
 796 BOTH:
 797     k = 0;
 798     } // end of if for PAIR
 799
 800
 801
 802     // bremsstrahlung
 803     // G3 default value: 1
 804     // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
 805     //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
 806     //               G4LowEnergyBremstrahlung
 807     // Particles: e-/e+; mu+/mu-
 808     // Physics:   EM
 809     // flag = 0 no bremsstrahlung
 810     // flag = 1 bremsstrahlung, photon processed
 811     // flag = 2 bremsstrahlung, no photon stored
 812     // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
 813                                  // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
 814     else if (strncmp(&sProcessFlag[i][0],"BREM",4) == 0) {
 815       for (j=0; j<iNbOfProc; j++) {
 816         if ((strncmp(&sProcessFlag[j][0],"PAIR",4) == 0) && iProcessValue[j] == 1) goto NOBREM;
 817       }
 818       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
 819         fprintf(pAliceInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
 820         fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
 821         fprintf(pAliceInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
 822         // two = bremsstrahlung by muons and charged hadrons is activated
 823         // zero = no meaning
 824         // muon and hadron bremsstrahlung
 825         // G4 particles: "gamma"
 826         // G3 default value: CUTGAM=0.001 GeV
 827         //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
 828         fCut = 0.0;
 829         for (j=0; j<iNbOfCut; j++) {
 830           if (strncmp(&sCutFlag[j][0],"BCUTM",5) == 0) fCut = fCutValue[j];
 831         }
 832         // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
 833         // three = lower bound of the material indices in which the respective thresholds apply
 834         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 835         fprintf(pAliceInp,"PAIRBREM  %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,fCut,three,fLastMaterial);
 836
 837         // for e+ and e-
 838         fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
 839         fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1);");
 840         // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
 841         // zero = not used
 842         // zero = not used
 843         // three = lower bound of the material indices in which the respective thresholds apply
 844         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 845         // one = step length in assigning indices
 846         //"ELPO-THR";
 847         fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",-one,zero,zero,three,fLastMaterial,one);
 848       }
 849       else if (iProcessValue[i] == 0) {
 850         fprintf(pAliceInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
 851         fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',0)\n");
 852       }
 853       else  {
 854         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
 855         fprintf(pAliceInp,"*No FLUKA card generated\n");
 856       }
 857 NOBREM:
 858       j = 0;
 859     } // end of else if (strncmp(&sProcessFlag[i][0],"BREM",4) == 0)
 860
 861
 862     // Cerenkov photon generation
 863     // G3 default value: 0
 864     // G4 process: G4Cerenkov
 865     //
 866     // Particles: charged
 867     // Physics:   Optical
 868     // flag = 0 no Cerenkov photon generation
 869     // flag = 1 Cerenkov photon generation
 870     // flag = 2 Cerenkov photon generation with primary stopped at each step
 871     //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
 872     else if (strncmp(&sProcessFlag[i][0],"CKOV",4) == 0) {
 873       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
 874         fprintf(pAliceInp,"*\n*Cerenkov photon generation\n");
 875         fprintf(pAliceInp,"*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n");
 876         Double_t emin = 2.07e-9; // minimum Cerenkov photon emission energy (in GeV!). Default: 2.07E-9 GeV (corresponding to 600 nm)
 877         Double_t emax = 4.96e-9; // maximum Cerenkov photon emission energy (in GeV!). Default: 4.96E-9 GeV (corresponding to 250 nm)
 878         fprintf(pAliceInp,"OPT-PROD  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fCERENKOV\n",emin,emax,zero,three,fLastMaterial,one);
 879       }
 880       else if (iProcessValue[i] == 0) {
 881         fprintf(pAliceInp,"*\n*No Cerenkov photon generation\n");
 882         fprintf(pAliceInp,"*Generated from call: SetProcess('CKOV',0)\n");
 883         // zero = not used
 884         // zero = not used
 885         // zero = not used
 886         // three = lower bound of the material indices in which the respective thresholds apply
 887         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 888         // one = step length in assigning indices
 889         //"CERE-OFF";
 890         fprintf(pAliceInp,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,three,fLastMaterial,one);
 891       }
 892       else  {
 893         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
 894         fprintf(pAliceInp,"*No FLUKA card generated\n");
 895       }
 896     } // end of else if (strncmp(&sProcessFlag[i][0],"CKOV",4) == 0)
 897
 898
 899     // Compton scattering
 900     // G3 default value: 1
 901     // G4 processes: G4ComptonScattering,
 902     //               G4LowEnergyCompton,
 903     //               G4PolarizedComptonScattering
 904     // Particles: gamma
 905     // Physics:   EM
 906     // flag = 0 no Compton scattering
 907     // flag = 1 Compton scattering, electron processed
 908     // flag = 2 Compton scattering, no electron stored
 909     // gMC ->SetProcess("COMP",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. PHOT-THR
 910     else if (strncmp(&sProcessFlag[i][0],"COMP",4) == 0) {
 911       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
 912         fprintf(pAliceInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
 913         fprintf(pAliceInp,"*Generated from call: SetProcess('COMP',1);\n");
 914         // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
 915         // zero = not used
 916         // zero = not used
 917         // three = lower bound of the material indices in which the respective thresholds apply
 918         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 919         // one = step length in assigning indices
 920         //"PHOT-THR";
 921         fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",-one,zero,zero,three,fLastMaterial,one);
 922       }
 923       else if (iProcessValue[i] == 0) {
 924         fprintf(pAliceInp,"*\n*No Compton scattering - no FLUKA card generated\n");
 925         fprintf(pAliceInp,"*Generated from call: SetProcess('COMP',0)\n");
 926       }
 927       else  {
 928         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
 929         fprintf(pAliceInp,"*No FLUKA card generated\n");
 930       }
 931     } // end of else if (strncmp(&sProcessFlag[i][0],"COMP",4) == 0)
 932
 933     // decay
 934     // G3 default value: 1
 935     // G4 process: G4Decay
 936     //
 937     // Particles: all which decay is applicable for
 938     // Physics:   General
 939     // flag = 0 no decays
 940     // flag = 1 decays, secondaries processed
 941     // flag = 2 decays, no secondaries stored
 942     //gMC ->SetProcess("DCAY",1); // not available
 943     else if ((strncmp(&sProcessFlag[i][0],"DCAY",4) == 0) && iProcessValue[i] == 1)
 944       cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not avaliable!" << endl;
 945
 946     // delta-ray
 947     // G3 default value: 2
 948     // !! G4 treats delta rays in different way
 949     // G4 processes: G4eIonisation/G4IeIonization,
 950     //               G4MuIonisation/G4IMuIonization,
 951     //               G4hIonisation/G4IhIonisation
 952     // Particles: charged
 953     // Physics:   EM
 954     // flag = 0 no energy loss
 955     // flag = 1 restricted energy loss fluctuations
 956     // flag = 2 complete energy loss fluctuations
 957     // flag = 3 same as 1
 958     // flag = 4 no energy loss fluctuations
 959     // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0.  0. 3. lastmat 0.
 960     else if (strncmp(&sProcessFlag[i][0],"DRAY",4) == 0) {
 961       if (iProcessValue[i] == 0 || iProcessValue[i] == 4) {
 962         fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
 963         fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
 964         fprintf(pAliceInp,"*No delta ray production by muons - threshold set artificially high\n");
 965         Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
 966         // zero = ignored
 967         // zero = ignored
 968         // three = lower bound of the material indices in which the respective thresholds apply
 969         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
 970         // one = step length in assigning indices
 971         fprintf(pAliceInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,three,fLastMaterial,one);
 972       }
 973       else if (iProcessValue[i] == 1 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
 974         fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
 975         fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
 976         fprintf(pAliceInp,"*Delta ray production by muons switched on\n");
 977         fprintf(pAliceInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
 978         fCut = 1.0e+6;
 979         for (j=0; j<iNbOfCut; j++) {
 980           if (strncmp(&sCutFlag[j][0],"DCUTM",5) == 0) fCut = fCutValue[j];
 981         }
 982         // fCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
 983         // zero = ignored
 984         // zero = ignored
 985         // three = lower bound of the material indices in which the respective thresholds apply
 986         // fLastMaterial =  upper bound of the material indices in which the respective thresholds apply
 987         // one = step length in assigning indices
 988         fprintf(pAliceInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCut,zero,zero,three,fLastMaterial,one);
 989       }
 990       else  {
 991         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
 992         fprintf(pAliceInp,"*No FLUKA card generated\n");
 993       }
 994     } // end of else if (strncmp(&sProcessFlag[i][0],"DRAY",4) == 0)
 995
 996     // hadronic process
 997     // G3 default value: 1
 998     // G4 processes: all defined by TG4PhysicsConstructorHadron
 999     //
1000     // Particles: hadrons
1001     // Physics:   Hadron
1002     // flag = 0 no multiple scattering
1003     // flag = 1 hadronic interactions, secondaries processed
1004     // flag = 2 hadronic interactions, no secondaries stored
1005     // gMC ->SetProcess("HADR",1); // ??? hadronic process
1006     //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
1007     else if (strncmp(&sProcessFlag[i][0],"HADR",4) == 0) {
1008       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
1009         fprintf(pAliceInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
1010         fprintf(pAliceInp,"*No FLUKA card generated\n");
1011       }
1012       else if (iProcessValue[i] == 0) {
1013         fprintf(pAliceInp,"*\n*Hadronic interaction is set OFF\n");
1014         fprintf(pAliceInp,"*Generated from call: SetProcess('HADR',0);\n");
1015         // zero = ignored
1016         // three = multiple scattering for hadrons and muons is completely suppressed
1017         // zero = no spin-relativistic corrections
1018         // three = lower bound of the material indices in which the respective thresholds apply
1019         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1020         fprintf(pAliceInp,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,zero,three,fLastMaterial);
1021
1022       }
1023       else  {
1024         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
1025         fprintf(pAliceInp,"*No FLUKA card generated\n");
1026       }
1027     } // end of else if (strncmp(&sProcessFlag[i][0],"HADR",4) == 0)
1028
1029
1030     // energy loss
1031     // G3 default value: 2
1032     // G4 processes: G4eIonisation/G4IeIonization,
1033     //               G4MuIonisation/G4IMuIonization,
1034     //               G4hIonisation/G4IhIonisation
1035     //
1036     // Particles: charged
1037     // Physics:   EM
1038     // flag=0 no energy loss
1039     // flag=1 restricted energy loss fluctuations
1040     // flag=2 complete energy loss fluctuations
1041     // flag=3 same as 1
1042     // flag=4 no energy loss fluctuations
1043     // If the value ILOSS is changed, then (in G3) cross-sections and energy
1044     // loss tables must be recomputed via the command 'PHYSI'
1045     // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1046     else if (strncmp(&sProcessFlag[i][0],"LOSS",4) == 0) {
1047       if (iProcessValue[i] == 2) { // complete energy loss fluctuations
1048         fprintf(pAliceInp,"*\n*Complete energy loss fluctuations do not exist in FLUKA\n");
1049         fprintf(pAliceInp,"*Generated from call: SetProcess('LOSS',2);\n");
1050         fprintf(pAliceInp,"*flag=2=complete energy loss fluctuations\n");
1051         fprintf(pAliceInp,"*No FLUKA card generated\n");
1052       }
1053       else if (iProcessValue[i] == 1 || iProcessValue[i] == 3) { // restricted energy loss fluctuations
1054         fprintf(pAliceInp,"*\n*Restricted energy loss fluctuations\n");
1055         fprintf(pAliceInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
1056         // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1057         // one = restricted energy loss fluctuations (for e+ and e-) switched on
1058         // one = minimal accuracy
1059         // three = lower bound of the material indices in which the respective thresholds apply
1060         // upper bound of the material indices in which the respective thresholds apply
1061         fprintf(pAliceInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,three,fLastMaterial);
1062       }
1063       else if (iProcessValue[i] == 4) { // no energy loss fluctuations
1064         fprintf(pAliceInp,"*\n*No energy loss fluctuations\n");
1065         fprintf(pAliceInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
1066         // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
1067         // - one = restricted energy loss fluctuations (for e+ and e-) switched off
1068         // one = minimal accuracy
1069         // three = lower bound of the material indices in which the respective thresholds apply
1070         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1071         fprintf(pAliceInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,three,fLastMaterial);
1072       }
1073       else  {
1074         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
1075         fprintf(pAliceInp,"*No FLUKA card generated\n");
1076       }
1077     } // end of else if (strncmp(&sProcessFlag[i][0],"LOSS",4) == 0)
1078
1079
1080     // multiple scattering
1081     // G3 default value: 1
1082     // G4 process: G4MultipleScattering/G4IMultipleScattering
1083     //
1084     // Particles: charged
1085     // Physics:   EM
1086     // flag = 0 no multiple scattering
1087     // flag = 1 Moliere or Coulomb scattering
1088     // flag = 2 Moliere or Coulomb scattering
1089     // flag = 3 Gaussian scattering
1090     // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
1091     else if (strncmp(&sProcessFlag[i][0],"MULS",4) == 0) {
1092       if (iProcessValue[i] == 1 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
1093         fprintf(pAliceInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
1094         fprintf(pAliceInp,"*No FLUKA card generated\n");
1095       }
1096       else if (iProcessValue[i] == 0) {
1097         fprintf(pAliceInp,"*\n*Multiple scattering is set OFF\n");
1098         fprintf(pAliceInp,"*Generated from call: SetProcess('MULS',0);\n");
1099         // zero = ignored
1100         // three = multiple scattering for hadrons and muons is completely suppressed
1101         // three = multiple scattering for e+ and e- is completely suppressed
1102         // three = lower bound of the material indices in which the respective thresholds apply
1103         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1104         fprintf(pAliceInp,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,three,fLastMaterial);
1105       }
1106       else  {
1107         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
1108         fprintf(pAliceInp,"*No FLUKA card generated\n");
1109       }
1110     } // end of else if (strncmp(&sProcessFlag[i][0],"MULS",4) == 0)
1111
1112
1113     // muon nuclear interaction
1114     // G3 default value: 0
1115     // G4 processes: G4MuNuclearInteraction,
1116     // G4MuonMinusCaptureAtRest
1117     //
1118     // Particles: mu
1119     // Physics:   Not set
1120     // flag = 0 no muon-nuclear interaction
1121     // flag = 1 nuclear interaction, secondaries processed
1122     // flag = 2 nuclear interaction, secondaries not processed
1123     // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
1124     else if (strncmp(&sProcessFlag[i][0],"MUNU",4) == 0) {
1125       if (iProcessValue[i] == 1) {
1126         fprintf(pAliceInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
1127         fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',1);\n");
1128         // one = full simulation of muon nuclear interactions and production of secondary hadrons
1129         // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1130         // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1131         // three = lower bound of the material indices in which the respective thresholds apply
1132         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1133         fprintf(pAliceInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,three,fLastMaterial);
1134       }
1135       else if (iProcessValue[i] == 2) {
1136         fprintf(pAliceInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
1137         fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',2);\n");
1138         // two = full simulation of muon nuclear interactions and production of secondary hadrons
1139         // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1140         // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1141         // three = lower bound of the material indices in which the respective thresholds apply
1142         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1143         fprintf(pAliceInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,three,fLastMaterial);
1144       }
1145       else if (iProcessValue[i] == 0) {
1146         fprintf(pAliceInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
1147         fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',0)\n");
1148       }
1149       else  {
1150         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
1151         fprintf(pAliceInp,"*No FLUKA card generated\n");
1152       }
1153     } // end of else if (strncmp(&sProcessFlag[i][0],"MUNU",4) == 0)
1154
1155
1156     // photofission
1157     // G3 default value: 0
1158     // G4 process: ??
1159     //
1160     // Particles: gamma
1161     // Physics:   ??
1162     // gMC ->SetProcess("PFIS",0); // PHOTONUC -1.   0.  0. 3. lastmat 0.
1163     // flag = 0 no photon fission
1164     // flag = 1 photon fission, secondaries processed
1165     // flag = 2 photon fission, no secondaries stored
1166     else if (strncmp(&sProcessFlag[i][0],"PFIS",4) == 0) {
1167       if (iProcessValue[i] == 0) {
1168         fprintf(pAliceInp,"*\n*No photonuclear interactions\n");
1169         fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',0);\n");
1170         // - one = no photonuclear interactions
1171         // zero = not used
1172         // zero = not used
1173         // three = lower bound of the material indices in which the respective thresholds apply
1174         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1175         fprintf(pAliceInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,three,fLastMaterial);
1176       }
1177       else if (iProcessValue[i] == 1) {
1178         fprintf(pAliceInp,"*\n*Photon nuclear interactions are activated at all energies\n");
1179         fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',1);\n");
1180         // one = photonuclear interactions are activated at all energies
1181         // zero = not used
1182         // zero = not used
1183         // three = lower bound of the material indices in which the respective thresholds apply
1184         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1185         fprintf(pAliceInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,three,fLastMaterial);
1186       }
1187       else if (iProcessValue[i] == 0) {
1188         fprintf(pAliceInp,"*\n*No photofission - no FLUKA card generated\n");
1189         fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',0)\n");
1190       }
1191       else {
1192         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
1193         fprintf(pAliceInp,"*No FLUKA card generated\n");
1194       }
1195     }
1196
1197
1198     // photo electric effect
1199     // G3 default value: 1
1200     // G4 processes: G4PhotoElectricEffect
1201     //               G4LowEnergyPhotoElectric
1202     // Particles: gamma
1203     // Physics:   EM
1204     // flag = 0 no photo electric effect
1205     // flag = 1 photo electric effect, electron processed
1206     // flag = 2 photo electric effect, no electron stored
1207     // gMC ->SetProcess("PHOT",1); // EMFCUT    0.  -1.  0. 3. lastmat 0. PHOT-THR
1208     else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0) {
1209       if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
1210         fprintf(pAliceInp,"*\n*Photo electric effect is activated\n");
1211         fprintf(pAliceInp,"*Generated from call: SetProcess('PHOT',1);\n");
1212         // zero = ignored
1213         // - one = resets to default=0.
1214         // zero = ignored
1215         // three = lower bound of the material indices in which the respective thresholds apply
1216         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1217         // one = step length in assigning indices
1218         //"PHOT-THR";
1219         fprintf(pAliceInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,-one,zero,three,fLastMaterial,one);
1220       }
1221       else if (iProcessValue[i] == 0) {
1222         fprintf(pAliceInp,"*\n*No photo electric effect - no FLUKA card generated\n");
1223         fprintf(pAliceInp,"*Generated from call: SetProcess('PHOT',0)\n");
1224       }
1225       else {
1226         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
1227         fprintf(pAliceInp,"*No FLUKA card generated\n");
1228       }
1229     } // else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0)
1230
1231
1232     // Rayleigh scattering
1233     // G3 default value: 0
1234     // G4 process: G4OpRayleigh
1235     //
1236     // Particles: optical photon
1237     // Physics:   Optical
1238     // flag = 0 Rayleigh scattering off
1239     // flag = 1 Rayleigh scattering on
1240     //xx gMC ->SetProcess("RAYL",1);
1241     else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0) {
1242       if (iProcessValue[i] == 1) {
1243         fprintf(pAliceInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
1244         fprintf(pAliceInp,"*No FLUKA card generated\n");
1245       }
1246       else if (iProcessValue[i] == 0) {
1247         fprintf(pAliceInp,"*\n*Rayleigh scattering is set OFF\n");
1248         fprintf(pAliceInp,"*Generated from call: SetProcess('RAYL',0);\n");
1249         // - one = no Rayleigh scattering and no binding corrections for Compton
1250         // three = lower bound of the material indices in which the respective thresholds apply
1251         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1252         fprintf(pAliceInp,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n",-one,three,three,fLastMaterial);
1253       }
1254       else  {
1255         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
1256         fprintf(pAliceInp,"*No FLUKA card generated\n");
1257       }
1258     } // end of else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0)
1259
1260
1261     // synchrotron radiation in magnetic field
1262     // G3 default value: 0
1263     // G4 process: G4SynchrotronRadiation
1264     //
1265     // Particles: ??
1266     // Physics:   Not set
1267     // flag = 0 no synchrotron radiation
1268     // flag = 1 synchrotron radiation
1269     //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
1270     else if (strncmp(&sProcessFlag[i][0],"SYNC",4) == 0) {
1271         fprintf(pAliceInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
1272         fprintf(pAliceInp,"*No FLUKA card generated\n");
1273     }
1274
1275
1276     // Automatic calculation of tracking medium parameters
1277     // flag = 0 no automatic calculation
1278     // flag = 1 automatic calculation
1279     //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
1280     else if (strncmp(&sProcessFlag[i][0],"AUTO",4) == 0) {
1281         fprintf(pAliceInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
1282         fprintf(pAliceInp,"*No FLUKA card generated\n");
1283     }
1284
1285
1286     // To control energy loss fluctuation model
1287     // flag = 0 Urban model
1288     // flag = 1 PAI model
1289     // flag = 2 PAI+ASHO model (not active at the moment)
1290     //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1291     else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0) {
1292       if (iProcessValue[i] == 0 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
1293         fprintf(pAliceInp,"*\n*Ionization energy losses calculation is activated\n");
1294         fprintf(pAliceInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
1295         // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1296         // one = restricted energy loss fluctuations (for e+ and e-) switched on
1297         // one = minimal accuracy
1298         // three = lower bound of the material indices in which the respective thresholds apply
1299         // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1300         fprintf(pAliceInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,three,fLastMaterial);
1301       }
1302       else {
1303         fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
1304         fprintf(pAliceInp,"*No FLUKA card generated\n");
1305       }
1306     } // else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0)
1307
1308
1309
1310
1311     else { // processes not yet treated
1312
1313     // light photon absorption (Cerenkov photons)
1314     // it is turned on when Cerenkov process is turned on
1315     // G3 default value: 0
1316     // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1317     //
1318     // Particles: optical photon
1319     // Physics:   Optical
1320     // flag = 0 no absorption of Cerenkov photons
1321     // flag = 1 absorption of Cerenkov photons
1322     // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1323
1324
1325
1326       cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not yet implemented!" << endl;
1327     }
1328   } //end of loop number of SetProcess calls
1329
1330
1331 // Loop over number of SetCut calls
1332   for (Int_t i=0; i<iNbOfCut; i++) {
1333
1334     // cuts used in SetProcess calls
1335     if (strncmp(&sCutFlag[i][0],"BCUTM",5) == 0) continue;
1336     else if (strncmp(&sCutFlag[i][0],"BCUTE",5) == 0) continue;
1337     else if (strncmp(&sCutFlag[i][0],"DCUTM",5) == 0) continue;
1338     else if (strncmp(&sCutFlag[i][0],"PPCUTM",6) == 0) continue;
1339
1340     // gammas
1341     // G4 particles: "gamma"
1342     // G3 default value: 0.001 GeV
1343     //gMC ->SetCut("CUTGAM",cut); // cut for gammas
1344     else if (strncmp(&sCutFlag[i][0],"CUTGAM",6) == 0) {
1345       fprintf(pAliceInp,"*\n*Cut for gamma\n");
1346       fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1347       // -fCutValue[i];
1348       // 7.0 = lower bound of the particle id-numbers to which the cut-off
1349       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f\n",-fCutValue[i],7.0);
1350     }
1351
1352     // electrons
1353     // G4 particles: "e-"
1354     // ?? positrons
1355     // G3 default value: 0.001 GeV
1356     //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1357     else if (strncmp(&sCutFlag[i][0],"CUTELE",6) == 0) {
1358       fprintf(pAliceInp,"*\n*Cut for electrons\n");
1359       fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1360       // -fCutValue[i];
1361       // three = lower bound of the particle id-numbers to which the cut-off
1362       // 4.0 = upper bound of the particle id-numbers to which the cut-off
1363       // one = step length in assigning numbers
1364       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],three,4.0,one);
1365     }
1366
1367     // neutral hadrons
1368     // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1369     // G3 default value: 0.01 GeV
1370     //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1371     else if (strncmp(&sCutFlag[i][0],"CUTNEU",6) == 0) {
1372       fprintf(pAliceInp,"*\n*Cut for neutral hadrons\n");
1373       fprintf(pAliceInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
1374
1375       // 8.0 = Neutron
1376       // 9.0 = Antineutron
1377       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],8.0,9.0);
1378
1379       // 12.0 = Kaon zero long
1380       // 12.0 = Kaon zero long
1381       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],12.0,12.0);
1382
1383       // 17.0 = Lambda, 18.0 = Antilambda
1384       // 19.0 = Kaon zero short
1385       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],17.0,19.0);
1386
1387       // 22.0 = Sigma zero, Pion zero, Kaon zero
1388       // 25.0 = Antikaon zero
1389       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],22.0,25.0);
1390
1391       // 32.0 = Antisigma zero
1392       // 32.0 = Antisigma zero
1393       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],32.0,32.0);
1394
1395       // 34.0 = Xi zero
1396       // 35.0 = AntiXi zero
1397       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],34.0,35.0);
1398
1399       // 47.0 = D zero
1400       // 48.0 = AntiD zero
1401       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],47.0,48.0);
1402
1403       // 53.0 = Xi_c zero
1404       // 53.0 = Xi_c zero
1405       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],53.0,53.0);
1406
1407       // 55.0 = Xi'_c zero
1408       // 56.0 = Omega_c zero
1409       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],55.0,56.0);
1410
1411       // 59.0 = AntiXi_c zero
1412       // 59.0 = AntiXi_c zero
1413       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],59.0,59.0);
1414
1415       // 61.0 = AntiXi'_c zero
1416       // 62.0 = AntiOmega_c zero
1417       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],61.0,62.0);
1418     }
1419
1420     // charged hadrons
1421     // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1422     // G3 default value: 0.01 GeV
1423     //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1424     else if (strncmp(&sCutFlag[i][0],"CUTHAD",6) == 0) {
1425       fprintf(pAliceInp,"*\n*Cut for charged hadrons\n");
1426       fprintf(pAliceInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
1427
1428       // 1.0 = Proton
1429       // 2.0 = Antiproton
1430       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],1.0,2.0);
1431
1432       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
1433       // 16.0 = Negative Kaon
1434       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],13.0,16.0);
1435
1436       // 20.0 = Negative Sigma
1437       // 21.0 = Positive Sigma
1438       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],20.0,21.0);
1439
1440       // 31.0 = Antisigma minus
1441       // 33.0 = Antisigma plus
1442       // 2.0 = step length
1443       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],31.0,33.0,2.0);
1444
1445       // 36.0 = Negative Xi, Positive Xi, Omega minus
1446       // 39.0 = Antiomega
1447       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],36.0,39.0);
1448
1449       // 45.0 = D plus
1450       // 46.0 = D minus
1451       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],45.0,46.0);
1452
1453       // 49.0 = D_s plus, D_s minus, Lambda_c plus
1454       // 52.0 = Xi_c plus
1455       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],49.0,52.0);
1456
1457       // 54.0 = Xi'_c plus
1458       // 60.0 = AntiXi'_c minus
1459       // 6.0 = step length
1460       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],54.0,60.0,6.0);
1461
1462       // 57.0 = Antilambda_c minus
1463       // 58.0 = AntiXi_c minus
1464       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],57.0,58.0);
1465     }
1466
1467     // muons
1468     // G4 particles: "mu+", "mu-"
1469     // G3 default value: 0.01 GeV
1470     //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1471     else if (strncmp(&sCutFlag[i][0],"CUTMUO",6) == 0) {
1472       fprintf(pAliceInp,"*\n*Cut for muons\n");
1473       fprintf(pAliceInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
1474       // 10.0 = Muon+
1475       // 11.0 = Muon-
1476       fprintf(pAliceInp,"PART-THR  %10.4g%10.1f%10.1f\n",-fCutValue[i],10.0,11.0);
1477     }
1478
1479     // delta-rays by electrons
1480     // G4 particles: "e-"
1481     // G3 default value: 10**4 GeV
1482     // gMC ->SetCut("DCUTE",cut);  // cut for deltarays by electrons ???????????????
1483     else if (strncmp(&sCutFlag[i][0],"DCUTE",5) == 0) {
1484       fprintf(pAliceInp,"*\n*Cut for delta rays by electrons ????????????\n");
1485       fprintf(pAliceInp,"*Generated from call: SetCut('DCUTE',cut);\n");
1486       // -fCutValue[i];
1487       // zero = ignored
1488       // zero = ignored
1489       // three = lower bound of the material indices in which the respective thresholds apply
1490       // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
1491       fprintf(pAliceInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f\n",-fCutValue[i],zero,zero,three,fLastMaterial);
1492     }
1493
1494     //
1495     // time of flight cut in seconds
1496     // G4 particles: all
1497     // G3 default value: 0.01 GeV
1498     //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1499     else if (strncmp(&sCutFlag[i][0],"TOFMAX",6) == 0) {
1500       fprintf(pAliceInp,"*\n*Time of flight cuts in seconds\n");
1501       fprintf(pAliceInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
1502       // zero = ignored
1503       // zero = ignored
1504       // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1505       // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1506       fprintf(pAliceInp,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",fCutValue[i]*1.e9,zero,zero,-6.0,64.0);
1507     }
1508
1509     else {
1510       cout << "SetCut for flag=" << &sCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
1511     }
1512   } //end of loop over SetCut calls
1513
1514 // Add START and STOP card
1515    fprintf(pAliceInp,"START     %10.1f\n",fEventsPerRun);
1516    fprintf(pAliceInp,"STOP      \n");
1517
1518 } // end of InitPhysics
1519
1520
1521 void TFluka::SetMaxStep(Double_t)
1522 {
1523 // SetMaxStep is dummy procedure in TFluka !
1524   if (fVerbosityLevel >=3)
1525   cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
1526 }
1527
1528 void TFluka::SetMaxNStep(Int_t)
1529 {
1530 // SetMaxNStep is dummy procedure in TFluka !
1531   if (fVerbosityLevel >=3)
1532   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1533 }
1534
1535 void TFluka::SetUserDecay(Int_t)
1536 {
1537 // SetUserDecay is dummy procedure in TFluka !
1538   if (fVerbosityLevel >=3)
1539   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1540 }
1541
1542 //
1543 // dynamic properties
1544 //
1545 void TFluka::TrackPosition(TLorentzVector& position) const
1546 {
1547 // Return the current position in the master reference frame of the
1548 // track being transported
1549 // TRACKR.atrack = age of the particle
1550 // TRACKR.xtrack = x-position of the last point
1551 // TRACKR.ytrack = y-position of the last point
1552 // TRACKR.ztrack = z-position of the last point
1553   Int_t caller = GetCaller();
1554   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1555     position.SetX(GetXsco());
1556     position.SetY(GetYsco());
1557     position.SetZ(GetZsco());
1558     position.SetT(TRACKR.atrack);
1559   }
1560   else if (caller == 4) { // mgdraw
1561     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1562     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1563     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1564     position.SetT(TRACKR.atrack);
1565   }
1566   else if (caller == 5) { // sodraw
1567     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1568     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1569     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1570     position.SetT(0);
1571   }
1572   else
1573     Warning("TrackPosition","position not available");
1574 }
1575
1576 //
1577 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1578 {
1579 // Return the current position in the master reference frame of the
1580 // track being transported
1581 // TRACKR.atrack = age of the particle
1582 // TRACKR.xtrack = x-position of the last point
1583 // TRACKR.ytrack = y-position of the last point
1584 // TRACKR.ztrack = z-position of the last point
1585   Int_t caller = GetCaller();
1586   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1587     x = GetXsco();
1588     y = GetYsco();
1589     z = GetZsco();
1590   }
1591   else if (caller == 4 || caller == 5) { // mgdraw, sodraw
1592     x = TRACKR.xtrack[TRACKR.ntrack];
1593     y = TRACKR.ytrack[TRACKR.ntrack];
1594     z = TRACKR.ztrack[TRACKR.ntrack];
1595   }
1596   else
1597     Warning("TrackPosition","position not available");
1598 }
1599
1600 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1601 {
1602 // Return the direction and the momentum (GeV/c) of the track
1603 // currently being transported
1604 // TRACKR.ptrack = momentum of the particle (not always defined, if
1605 //               < 0 must be obtained from etrack)
1606 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1607 // TRACKR.etrack = total energy of the particle
1608 // TRACKR.jtrack = identity number of the particle
1609 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1610   Int_t caller = GetCaller();
1611   if (caller != 2) { // not eedraw
1612     if (TRACKR.ptrack >= 0) {
1613       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1614       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1615       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1616       momentum.SetE(TRACKR.etrack);
1617       return;
1618     }
1619     else {
1620       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1621       momentum.SetPx(p*TRACKR.cxtrck);
1622       momentum.SetPy(p*TRACKR.cytrck);
1623       momentum.SetPz(p*TRACKR.cztrck);
1624       momentum.SetE(TRACKR.etrack);
1625       return;
1626     }
1627   }
1628   else
1629     Warning("TrackMomentum","momentum not available");
1630 }
1631
1632 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1633 {
1634 // Return the direction and the momentum (GeV/c) of the track
1635 // currently being transported
1636 // TRACKR.ptrack = momentum of the particle (not always defined, if
1637 //               < 0 must be obtained from etrack)
1638 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1639 // TRACKR.etrack = total energy of the particle
1640 // TRACKR.jtrack = identity number of the particle
1641 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1642   Int_t caller = GetCaller();
1643   if (caller != 2) { // not eedraw
1644     if (TRACKR.ptrack >= 0) {
1645       px = TRACKR.ptrack*TRACKR.cxtrck;
1646       py = TRACKR.ptrack*TRACKR.cytrck;
1647       pz = TRACKR.ptrack*TRACKR.cztrck;
1648       e = TRACKR.etrack;
1649       return;
1650     }
1651     else {
1652       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1653       px = p*TRACKR.cxtrck;
1654       py = p*TRACKR.cytrck;
1655       pz = p*TRACKR.cztrck;
1656       e = TRACKR.etrack;
1657       return;
1658     }
1659   }
1660   else
1661     Warning("TrackMomentum","momentum not available");
1662 }
1663
1664 Double_t TFluka::TrackStep() const
1665 {
1666 // Return the length in centimeters of the current step
1667 // TRACKR.ctrack = total curved path
1668   Int_t caller = GetCaller();
1669   if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
1670     return 0.0;
1671   else if (caller == 4) //mgdraw
1672     return TRACKR.ctrack;
1673   else
1674     return -1.0;
1675 }
1676
1677 Double_t TFluka::TrackLength() const
1678 {
1679 // TRACKR.cmtrck = cumulative curved path since particle birth
1680   Int_t caller = GetCaller();
1681   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1682     return TRACKR.cmtrck;
1683   else
1684     return -1.0;
1685 }
1686
1687 Double_t TFluka::TrackTime() const
1688 {
1689 // Return the current time of flight of the track being transported
1690 // TRACKR.atrack = age of the particle
1691   Int_t caller = GetCaller();
1692   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1693     return TRACKR.atrack;
1694   else
1695     return -1;
1696 }
1697
1698 Double_t TFluka::Edep() const
1699 {
1700 // Energy deposition
1701 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1702 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1703 //    but in the variable "rull" of the procedure "endraw.cxx"
1704 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1705 // -->no energy loss along the track
1706 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1707 // -->energy loss distributed along the track
1708 // TRACKR.dtrack = energy deposition of the jth deposition even
1709
1710   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1711   Int_t caller = GetCaller();
1712   if (caller == 11 || caller==12) return 0.0;
1713   Double_t sum = 0;
1714   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1715     sum +=TRACKR.dtrack[j];
1716   }
1717   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1718     return fRull + sum;
1719   else {
1720     return sum;
1721   }
1722 }
1723
1724 Int_t TFluka::TrackPid() const
1725 {
1726 // Return the id of the particle transported
1727 // TRACKR.jtrack = identity number of the particle
1728   Int_t caller = GetCaller();
1729   if (caller != 2)  // not eedraw
1730     return PDGFromId(TRACKR.jtrack);
1731   else
1732     return -1000;
1733 }
1734
1735 Double_t TFluka::TrackCharge() const
1736 {
1737 // Return charge of the track currently transported
1738 // PAPROP.ichrge = electric charge of the particle
1739 // TRACKR.jtrack = identity number of the particle
1740   Int_t caller = GetCaller();
1741   if (caller != 2)  // not eedraw
1742     return PAPROP.ichrge[TRACKR.jtrack+6];
1743   else
1744     return -1000.0;
1745 }
1746
1747 Double_t TFluka::TrackMass() const
1748 {
1749 // PAPROP.am = particle mass in GeV
1750 // TRACKR.jtrack = identity number of the particle
1751   Int_t caller = GetCaller();
1752   if (caller != 2) {  // not eedraw
1753 //    cout << "JTRACK=" << TRACKR.jtrack << " mass=" << PAPROP.am[TRACKR.jtrack+6] << endl;
1754     return PAPROP.am[TRACKR.jtrack+6];
1755   }
1756   else
1757     return -1000.0;
1758 }
1759
1760 Double_t TFluka::Etot() const
1761 {
1762 // TRACKR.etrack = total energy of the particle
1763   Int_t caller = GetCaller();
1764   if (caller != 2)  // not eedraw
1765     return TRACKR.etrack;
1766   else
1767     return -1000.0;
1768 }
1769
1770 //
1771 // track status
1772 //
1773 Bool_t   TFluka::IsNewTrack() const
1774 {
1775 // Return true for the first call of Stepping()
1776    return fTrackIsNew;
1777 }
1778
1779 Bool_t   TFluka::IsTrackInside() const
1780 {
1781 // True if the track is not at the boundary of the current volume
1782 // In Fluka a step is always inside one kind of material
1783 // If the step would go behind the region of one material,
1784 // it will be shortened to reach only the boundary.
1785 // Therefore IsTrackInside() is always true.
1786   Int_t caller = GetCaller();
1787   if (caller == 11 || caller==12)  // bxdraw
1788     return 0;
1789   else
1790     return 1;
1791 }
1792
1793 Bool_t   TFluka::IsTrackEntering() const
1794 {
1795 // True if this is the first step of the track in the current volume
1796
1797   Int_t caller = GetCaller();
1798   if (caller == 11)  // bxdraw entering
1799     return 1;
1800   else return 0;
1801 }
1802
1803 Bool_t   TFluka::IsTrackExiting() const
1804 {
1805   Int_t caller = GetCaller();
1806   if (caller == 12)  // bxdraw exiting
1807     return 1;
1808   else return 0;
1809 }
1810
1811 Bool_t   TFluka::IsTrackOut() const
1812 {
1813 // True if the track is out of the setup
1814 // means escape
1815 // Icode = 14: escape - call from Kaskad
1816 // Icode = 23: escape - call from Emfsco
1817 // Icode = 32: escape - call from Kasneu
1818 // Icode = 40: escape - call from Kashea
1819 // Icode = 51: escape - call from Kasoph
1820   if (fIcode == 14 ||
1821       fIcode == 23 ||
1822       fIcode == 32 ||
1823       fIcode == 40 ||
1824       fIcode == 51) return 1;
1825   else return 0;
1826 }
1827
1828 Bool_t   TFluka::IsTrackDisappeared() const
1829 {
1830 // means all inelastic interactions and decays
1831 // fIcode from usdraw
1832   if (fIcode == 101 || // inelastic interaction
1833       fIcode == 102 || // particle decay
1834       fIcode == 214 || // in-flight annihilation
1835       fIcode == 215 || // annihilation at rest
1836       fIcode == 217 || // pair production
1837       fIcode == 221) return 1;
1838   else return 0;
1839 }
1840
1841 Bool_t   TFluka::IsTrackStop() const
1842 {
1843 // True if the track energy has fallen below the threshold
1844 // means stopped by signal or below energy threshold
1845 // Icode = 12: stopping particle       - call from Kaskad
1846 // Icode = 15: time kill               - call from Kaskad
1847 // Icode = 21: below threshold, iarg=1 - call from Emfsco
1848 // Icode = 22: below threshold, iarg=2 - call from Emfsco
1849 // Icode = 24: time kill               - call from Emfsco
1850 // Icode = 31: below threshold         - call from Kasneu
1851 // Icode = 33: time kill               - call from Kasneu
1852 // Icode = 41: time kill               - call from Kashea
1853 // Icode = 52: time kill               - call from Kasoph
1854   if (fIcode == 12 ||
1855       fIcode == 15 ||
1856       fIcode == 21 ||
1857       fIcode == 22 ||
1858       fIcode == 24 ||
1859       fIcode == 31 ||
1860       fIcode == 33 ||
1861       fIcode == 41 ||
1862       fIcode == 52) return 1;
1863   else return 0;
1864 }
1865
1866 Bool_t   TFluka::IsTrackAlive() const
1867 {
1868 // means not disappeared or not out
1869   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1870   else return 1;
1871 }
1872
1873 //
1874 // secondaries
1875 //
1876
1877 Int_t TFluka::NSecondaries() const
1878 // Number of secondary particles generated in the current step
1879 // FINUC.np = number of secondaries except light and heavy ions
1880 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
1881 {
1882   Int_t caller = GetCaller();
1883   if (caller == 6)  // valid only after usdraw
1884     return FINUC.np + FHEAVY.npheav;
1885   else
1886     return 0;
1887 } // end of NSecondaries
1888
1889 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1890                 TLorentzVector& position, TLorentzVector& momentum)
1891 {
1892   Int_t caller = GetCaller();
1893   if (caller == 6) {  // valid only after usdraw
1894     if (isec >= 0 && isec < FINUC.np) {
1895       particleId = PDGFromId(FINUC.kpart[isec]);
1896       position.SetX(fXsco);
1897       position.SetY(fYsco);
1898       position.SetZ(fZsco);
1899       position.SetT(TRACKR.atrack);
1900       momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1901       momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1902       momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1903       momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1904     }
1905     else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1906       Int_t jsec = isec - FINUC.np;
1907       particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1908       position.SetX(fXsco);
1909       position.SetY(fYsco);
1910       position.SetZ(fZsco);
1911       position.SetT(TRACKR.atrack);
1912       momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1913       momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1914       momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1915       if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1916         momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1917       else if (FHEAVY.tkheav[jsec] > 6)
1918         momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1919     }
1920     else
1921       Warning("GetSecondary","isec out of range");
1922   }
1923   else
1924     Warning("GetSecondary","no secondaries available");
1925 } // end of GetSecondary
1926
1927 TMCProcess TFluka::ProdProcess(Int_t) const
1928 // Name of the process that has produced the secondary particles
1929 // in the current step
1930 {
1931     const TMCProcess kIpNoProc = kPNoProcess;
1932     const TMCProcess kIpPDecay = kPDecay;
1933     const TMCProcess kIpPPair = kPPair;
1934 //  const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
1935 //  const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
1936     const TMCProcess kIpPCompton = kPCompton;
1937     const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
1938     const TMCProcess kIpPBrem = kPBrem;
1939 //  const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
1940 //  const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
1941     const TMCProcess kIpPDeltaRay = kPDeltaRay;
1942 //  const TMCProcess kIpPMoller = kPMoller;
1943 //  const TMCProcess kIpPBhabha = kPBhabha;
1944     const TMCProcess kIpPAnnihilation = kPAnnihilation;
1945 //  const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
1946 //  const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
1947     const TMCProcess kIpPHadronic = kPHadronic;
1948     const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
1949     const TMCProcess kIpPPhotoFission = kPPhotoFission;
1950     const TMCProcess kIpPRayleigh = kPRayleigh;
1951 //  const TMCProcess kIpPCerenkov = kPCerenkov;
1952 //  const TMCProcess kIpPSynchrotron = kPSynchrotron;
1953
1954     Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
1955     if (fIcode == 102) return kIpPDecay;
1956     else if (fIcode == 104 || fIcode == 217) return kIpPPair;
1957 //  else if (fIcode == 104) return kIpPairFromPhoton;
1958 //  else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
1959     else if (fIcode == 219) return kIpPCompton;
1960     else if (fIcode == 221) return kIpPPhotoelectric;
1961     else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
1962 //  else if (fIcode == 105) return kIpPBremFromHeavy;
1963 //  else if (fIcode == 208) return kPBremFromElectronOrPositron;
1964     else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
1965     else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
1966 //  else if (fIcode == 210) return kIpPMoller;
1967 //  else if (fIcode == 212) return kIpPBhabha;
1968     else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
1969 //  else if (fIcode == 214) return kIpPAnnihilInFlight;
1970 //  else if (fIcode == 215) return kIpPAnnihilAtRest;
1971     else if (fIcode == 101) return kIpPHadronic;
1972     else if (fIcode == 101) {
1973       if (!mugamma) return kIpPHadronic;
1974       else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
1975       else return kIpPMuonNuclear;
1976     }
1977     else if (fIcode == 225) return kIpPRayleigh;
1978 // Fluka codes 100, 300 and 400 still to be investigasted
1979     else return kIpNoProc;
1980 }
1981
1982 //Int_t StepProcesses(TArrayI &proc) const
1983 // Return processes active in the current step
1984 //{
1985 //ck = total energy of the particl ????????????????
1986 //}
1987
1988
1989 Int_t TFluka::VolId2Mate(Int_t id) const
1990 {
1991 //
1992 // Returns the material number for a given volume ID
1993 //
1994     if (fVerbosityLevel >= 3)
1995     printf("VolId2Mate %d %d\n", id, fMediaByRegion[id-1]);
1996     return fMediaByRegion[id-1];
1997 }
1998
1999 const char* TFluka::VolName(Int_t id) const
2000 {
2001 //
2002 // Returns the volume name for a given volume ID
2003 //
2004     FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]);
2005     const char* name = vol->GetName();
2006     if (fVerbosityLevel >= 3)
2007     printf("VolName %d %s \n", id, name);
2008     return name;
2009 }
2010
2011 Int_t TFluka::VolId(const Text_t* volName) const
2012 {
2013 //
2014 // Converts from volume name to volume ID.
2015 // Time consuming. (Only used during set-up)
2016 // Could be replaced by hash-table
2017 //
2018     char tmp[5];
2019     Int_t i =0;
2020     for (i = 0; i < fNVolumes; i++)
2021   {
2022       FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
2023       TString name = vol->GetName();
2024       strcpy(tmp, name.Data());
2025       tmp[4] = '\0';
2026       if (!strcmp(tmp, volName)) break;
2027   }
2028     i++;
2029
2030     return i;
2031 }
2032
2033
2034 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2035 {
2036 //
2037 // Return the logical id and copy number corresponding to the current fluka region
2038 //
2039     int ir = fCurrentFlukaRegion;
2040     int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo);
2041     copyNo++;
2042     if (fVerbosityLevel >= 3)
2043     printf("CurrentVolID: %d %d %d \n", ir, id, copyNo);
2044     return id;
2045 }
2046
2047 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2048 {
2049 //
2050 // Return the logical id and copy number of off'th mother
2051 // corresponding to the current fluka region
2052 //
2053     if (off == 0)
2054         return CurrentVolID(copyNo);
2055
2056     int ir = fCurrentFlukaRegion;
2057     int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo);
2058     copyNo++;
2059     if (fVerbosityLevel >= 3)
2060     printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo);
2061     if (id == -1)
2062         if (fVerbosityLevel >= 0)
2063         printf("CurrentVolOffID: Warning Mother not found !!!\n");
2064     return id;
2065 }
2066
2067
2068 const char* TFluka::CurrentVolName() const
2069 {
2070 //
2071 // Return the current volume name
2072 //
2073     Int_t copy;
2074     Int_t id = TFluka::CurrentVolID(copy);
2075     const char* name = TFluka::VolName(id);
2076     if (fVerbosityLevel >= 3)
2077     printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion,  name);
2078     return name;
2079 }
2080
2081 const char* TFluka::CurrentVolOffName(Int_t off) const
2082 {
2083 //
2084 // Return the volume name of the off'th mother of the current volume
2085 //
2086     Int_t copy;
2087     Int_t id = TFluka::CurrentVolOffID(off, copy);
2088     const char* name = TFluka::VolName(id);
2089     if (fVerbosityLevel >= 3)
2090     printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion,  name);
2091     return name;
2092 }
2093
2094 Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/,
2095                       Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const
2096 {
2097 //
2098 //  Return the current medium number
2099 //
2100     Int_t copy;
2101     Int_t id  =  TFluka::CurrentVolID(copy);
2102     Int_t med =  TFluka::VolId2Mate(id);
2103     if (fVerbosityLevel >= 3)
2104     printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion,  med);
2105     return med;
2106 }
2107
2108 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2109     {
2110 // Transforms a position from the world reference frame
2111 // to the current volume reference frame.
2112 //
2113 //  Geant3 desription:
2114 //  ==================
2115 //       Computes coordinates XD (in DRS)
2116 //       from known coordinates XM in MRS
2117 //       The local reference system can be initialized by
2118 //         - the tracking routines and GMTOD used in GUSTEP
2119 //         - a call to GMEDIA(XM,NUMED)
2120 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2121 //             (inverse routine is GDTOM)
2122 //
2123 //        If IFLAG=1  convert coordinates
2124 //           IFLAG=2  convert direction cosinus
2125 //
2126 // ---
2127         Double_t xmD[3], xdD[3];
2128         xmD[0] = xm[0]; xmD[1] = xm[1]; xmD[2] = xm[2];
2129         (FGeometryInit::GetInstance())->Gmtod(xmD, xdD, iflag);
2130         xd[0] = xdD[0]; xd[1] = xdD[1]; xd[2] = xdD[2];
2131     }
2132
2133
2134 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2135     {
2136 // Transforms a position from the world reference frame
2137 // to the current volume reference frame.
2138 //
2139 //  Geant3 desription:
2140 //  ==================
2141 //       Computes coordinates XD (in DRS)
2142 //       from known coordinates XM in MRS
2143 //       The local reference system can be initialized by
2144 //         - the tracking routines and GMTOD used in GUSTEP
2145 //         - a call to GMEDIA(XM,NUMED)
2146 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2147 //             (inverse routine is GDTOM)
2148 //
2149 //        If IFLAG=1  convert coordinates
2150 //           IFLAG=2  convert direction cosinus
2151 //
2152 // ---
2153         (FGeometryInit::GetInstance())->Gmtod(xm, xd, iflag);
2154     }
2155
2156 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2157     {
2158 // Transforms a position from the current volume reference frame
2159 // to the world reference frame.
2160 //
2161 //  Geant3 desription:
2162 //  ==================
2163 //  Computes coordinates XM (Master Reference System
2164 //  knowing the coordinates XD (Detector Ref System)
2165 //  The local reference system can be initialized by
2166 //    - the tracking routines and GDTOM used in GUSTEP
2167 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2168 //        (inverse routine is GMTOD)
2169 //
2170 //   If IFLAG=1  convert coordinates
2171 //      IFLAG=2  convert direction cosinus
2172 //
2173 // ---
2174         Double_t xmD[3], xdD[3];
2175         xdD[0] = xd[0]; xdD[1] = xd[1]; xdD[2] = xd[2];
2176         (FGeometryInit::GetInstance())->Gdtom(xdD, xmD, iflag);
2177         xm[0] = xmD[0]; xm[1] = xmD[1]; xm[2] = xmD[2];
2178     }
2179 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2180     {
2181 // Transforms a position from the current volume reference frame
2182 // to the world reference frame.
2183 //
2184 //  Geant3 desription:
2185 //  ==================
2186 //  Computes coordinates XM (Master Reference System
2187 //  knowing the coordinates XD (Detector Ref System)
2188 //  The local reference system can be initialized by
2189 //    - the tracking routines and GDTOM used in GUSTEP
2190 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2191 //        (inverse routine is GMTOD)
2192 //
2193 //   If IFLAG=1  convert coordinates
2194 //      IFLAG=2  convert direction cosinus
2195 //
2196 // ---
2197
2198         (FGeometryInit::GetInstance())->Gdtom(xd, xm, iflag);
2199     }
2200
2201 // ===============================================================