TFluka/TFluka.cxx

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