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