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