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