TFluka/TFluka.cxx

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