TFluka/TFluka.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.9  2002/12/06 12:41:29  morsch
  19 Mess from last merge cleaned up.
  20
  21 Revision 1.8  2002/12/06 12:28:44  morsch
  22 Region to media mapping corrected and improved.
  23
  24 Revision 1.7  2002/12/06 12:21:32  morsch
  25 User stepping methods added (E. Futo)
  26
  27 Revision 1.6  2002/11/21 18:40:06  iglez2
  28 Media handling added
  29
  30 Revision 1.5  2002/11/07 17:59:10  iglez2
  31 Included the geometry through geant4_vmc/FLUGG
  32
  33 Revision 1.4  2002/11/04 16:00:46  iglez2
  34 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
  35
  36 Revision 1.3  2002/10/22 15:12:14  alibrary
  37 Introducing Riostream.h
  38
  39 Revision 1.2  2002/10/14 14:57:40  hristov
  40 Merging the VirtualMC branch to the main development branch (HEAD)
  41
  42 Revision 1.1.2.8  2002/10/08 16:33:17  iglez2
  43 LSOUIT is set to true before the second call to flukam.
  44
  45 Revision 1.1.2.7  2002/10/08 09:30:37  iglez2
  46 Solved stupid missing ;
  47
  48 Revision 1.1.2.6  2002/10/07 13:40:22  iglez2
  49 First implementations of the PDG <--> Fluka Id conversion routines
  50
  51 Revision 1.1.2.5  2002/09/26 16:26:03  iglez2
  52 Added verbosity
  53 Call to gAlice->Generator()->Generate()
  54
  55 Revision 1.1.2.4  2002/09/26 13:22:23  iglez2
  56 Naive implementation of ProcessRun and ProcessEvent
  57 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
  58
  59 Revision 1.1.2.3  2002/09/20 15:35:51  iglez2
  60 Modification of LFDRTR. Value is passed to FLUKA !!!
  61
  62 Revision 1.1.2.2  2002/09/18 14:34:44  iglez2
  63 Revised version with all pure virtual methods implemented
  64
  65 Revision 1.1.2.1  2002/07/24 08:49:41  alibrary
  66 Adding TFluka to VirtualMC
  67
  68 Revision 1.1  2002/07/05 13:10:07  morsch
  69 First commit of Fluka interface.
  70
  71 */
  72
  73 #include <Riostream.h>
  74
  75 #include "TClonesArray.h"
  76 #include "TFluka.h"
  77 #include "TCallf77.h"      //For the fortran calls
  78 #include "Fdblprc.h"       //(DBLPRC) fluka common
  79 #include "Fepisor.h"       //(EPISOR) fluka common
  80 #include "Ffinuc.h"        //(FINUC) fluka common
  81 #include "Fiounit.h"       //(IOUNIT) fluka common
  82 #include "Fpaprop.h"       //(PAPROP) fluka common
  83 #include "Fpart.h"         //(PART)   fluka common
  84 #include "Ftrackr.h"       //(TRACKR) fluka common
  85 #include "Fpaprop.h"       //(PAPROP) fluka common
  86 #include "Ffheavy.h"       //(FHEAVY) fluka common
  87
  88 #include "TVirtualMC.h"
  89 #include "TG4GeometryManager.h" //For the geometry management
  90 #include "TG4DetConstruction.h" //For the detector construction
  91
  92 #include "FGeometryInit.hh"
  93 #include "TLorentzVector.h"
  94 #include "FlukaVolume.h"
  95
  96 // Fluka methods that may be needed.
  97 #ifndef WIN32
  98 # define flukam  flukam_
  99 # define fluka_openinp fluka_openinp_
 100 # define fluka_closeinp fluka_closeinp_
 101 # define mcihad mcihad_
 102 # define mpdgha mpdgha_
 103 #else
 104 # define flukam  FLUKAM
 105 # define fluka_openinp FLUKA_OPENINP
 106 # define fluka_closeinp FLUKA_CLOSEINP
 107 # define mcihad MCIHAD
 108 # define mpdgha MPDGHA
 109 #endif
 110
 111 extern "C"
 112 {
 113   //
 114   // Prototypes for FLUKA functions
 115   //
 116   void type_of_call flukam(const int&);
 117   void type_of_call fluka_openinp(const int&, DEFCHARA);
 118   void type_of_call fluka_closeinp(const int&);
 119   int  type_of_call mcihad(const int&);
 120   int  type_of_call mpdgha(const int&);
 121 }
 122
 123 //
 124 // Class implementation for ROOT
 125 //
 126 ClassImp(TFluka)
 127
 128 //
 129 //----------------------------------------------------------------------------
 130 // TFluka constructors and destructors.
 131 //____________________________________________________________________________
 132 TFluka::TFluka()
 133   :TVirtualMC(),
 134    fVerbosityLevel(0),
 135    fInputFileName(""),
 136    fDetector(0),
 137    fCurrentFlukaRegion(-1)
 138 {
 139   //
 140   // Default constructor
 141   //
 142 }
 143
 144 TFluka::TFluka(const char *title, Int_t verbosity)
 145   :TVirtualMC("TFluka",title),
 146    fVerbosityLevel(verbosity),
 147    fInputFileName(""),
 148    fDetector(0),
 149    fCurrentFlukaRegion(-1)
 150 {
 151   if (fVerbosityLevel >=3)
 152     cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
 153
 154
 155   // create geometry manager
 156   if (fVerbosityLevel >=2)
 157     cout << "\t* Creating G4 Geometry manager..." << endl;
 158   fGeometryManager = new TG4GeometryManager();
 159   if (fVerbosityLevel >=2)
 160     cout << "\t* Creating G4 Detector..." << endl;
 161   fDetector = new TG4DetConstruction();
 162   FGeometryInit* geominit = FGeometryInit::GetInstance();
 163   if (geominit)
 164     geominit->setDetConstruction(fDetector);
 165   else {
 166     cerr << "ERROR: Could not create FGeometryInit!" << endl;
 167     cerr << "       Exiting!!!" << endl;
 168     abort();
 169   }
 170
 171   if (fVerbosityLevel >=3)
 172     cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 173
 174   fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
 175   fNVolumes      = 0;
 176   fMediaByRegion = 0;
 177 }
 178
 179 TFluka::~TFluka() {
 180   if (fVerbosityLevel >=3)
 181     cout << "==> TFluka::~TFluka() destructor called." << endl;
 182
 183   delete fGeometryManager;
 184   fVolumeMediaMap->Delete();
 185   delete  fVolumeMediaMap;
 186
 187
 188   if (fVerbosityLevel >=3)
 189     cout << "<== TFluka::~TFluka() destructor called." << endl;
 190 }
 191
 192 //
 193 //_____________________________________________________________________________
 194 // TFluka control methods
 195 //____________________________________________________________________________
 196 void TFluka::Init() {
 197   if (fVerbosityLevel >=3)
 198     cout << "==> TFluka::Init() called." << endl;
 199
 200   if (fVerbosityLevel >=2)
 201     cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
 202          << ") in fluka..." << endl;
 203   GLOBAL.lfdrtr = true;
 204
 205   if (fVerbosityLevel >=2)
 206     cout << "\t* Opening file " << fInputFileName << endl;
 207   const char* fname = fInputFileName;
 208   fluka_openinp(lunin, PASSCHARA(fname));
 209
 210   if (fVerbosityLevel >=2)
 211     cout << "\t* Calling flukam..." << endl;
 212   flukam(1);
 213
 214   if (fVerbosityLevel >=2)
 215     cout << "\t* Closing file " << fInputFileName << endl;
 216   fluka_closeinp(lunin);
 217
 218   if (fVerbosityLevel >=3)
 219     cout << "<== TFluka::Init() called." << endl;
 220
 221   FinishGeometry();
 222
 223 }
 224
 225 void TFluka::FinishGeometry() {
 226 //
 227 // Build-up table with region to medium correspondance
 228 //
 229     char tmp[5];
 230
 231   if (fVerbosityLevel >=3)
 232     cout << "==> TFluka::FinishGeometry() called." << endl;
 233
 234 //  fGeometryManager->Ggclos();
 235
 236   FGeometryInit* flugg = FGeometryInit::GetInstance();
 237
 238   fMediaByRegion = new Int_t[fNVolumes+2];
 239   for (Int_t i = 0; i < fNVolumes; i++)
 240   {
 241       FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
 242       TString volName = vol->GetName();
 243       Int_t   media   = vol->GetMedium();
 244       printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
 245       strcpy(tmp, volName.Data());
 246       tmp[4] = '\0';
 247       flugg->SetMediumFromName(tmp, media);
 248   }
 249
 250   flugg->BuildMediaMap();
 251
 252   if (fVerbosityLevel >=3)
 253     cout << "<== TFluka::FinishGeometry() called." << endl;
 254 }
 255
 256 void TFluka::BuildPhysics() {
 257   if (fVerbosityLevel >=3)
 258     cout << "==> TFluka::BuildPhysics() called." << endl;
 259
 260
 261   if (fVerbosityLevel >=3)
 262     cout << "<== TFluka::BuildPhysics() called." << endl;
 263 }
 264
 265 void TFluka::ProcessEvent() {
 266   if (fVerbosityLevel >=3)
 267     cout << "==> TFluka::ProcessEvent() called." << endl;
 268
 269   if (fVerbosityLevel >=3)
 270     cout << "<== TFluka::ProcessEvent() called." << endl;
 271 }
 272
 273
 274 void TFluka::ProcessRun(Int_t nevent) {
 275   if (fVerbosityLevel >=3)
 276     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 277          << endl;
 278
 279   if (fVerbosityLevel >=2) {
 280     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 281     cout << "\t* Calling flukam again..." << endl;
 282   }
 283   fApplication->GeneratePrimaries();
 284   EPISOR.lsouit = true;
 285   flukam(1);
 286
 287   if (fVerbosityLevel >=3)
 288     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 289          << endl;
 290 }
 291
 292 //_____________________________________________________________________________
 293 // methods for building/management of geometry
 294 //____________________________________________________________________________
 295 // functions from GCONS
 296 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 297                     Float_t &dens, Float_t &radl, Float_t &absl,
 298                     Float_t* ubuf, Int_t& nbuf) {
 299 //
 300   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
 301 }
 302
 303 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 304                     Double_t &dens, Double_t &radl, Double_t &absl,
 305                     Double_t* ubuf, Int_t& nbuf) {
 306 //
 307   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
 308 }
 309
 310 // detector composition
 311 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 312                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 313                       Float_t* buf, Int_t nwbuf) {
 314 //
 315   fGeometryManager
 316     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
 317 }
 318 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 319                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 320                       Double_t* buf, Int_t nwbuf) {
 321 //
 322   fGeometryManager
 323     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
 324 }
 325
 326 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 327                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 328 //
 329    fGeometryManager
 330      ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
 331 }
 332 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 333                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 334 //
 335    fGeometryManager
 336      ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
 337 }
 338
 339 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 340                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 341                     Double_t stemax, Double_t deemax, Double_t epsil,
 342                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 343   //
 344   fGeometryManager
 345     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 346              epsil, stmin, ubuf, nbuf);
 347 }
 348 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 349                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 350                     Double_t stemax, Double_t deemax, Double_t epsil,
 351                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 352   //
 353   fGeometryManager
 354     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 355              epsil, stmin, ubuf, nbuf);
 356 }
 357
 358 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 359                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 360                     Double_t phiZ) {
 361 //
 362   fGeometryManager
 363     ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 364 }
 365
 366 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
 367 //
 368   fGeometryManager->Gstpar(itmed, param, parval);
 369 }
 370
 371 // functions from GGEOM
 372 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 373                      Float_t *upar, Int_t np)  {
 374 //
 375 //  fVolumeMediaMap[TString(name)] = nmed;
 376     TClonesArray &lvols = *fVolumeMediaMap;
 377     new(lvols[fNVolumes++])
 378         FlukaVolume(name, nmed);
 379     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
 380 }
 381 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 382                      Double_t *upar, Int_t np)  {
 383 //
 384     TClonesArray &lvols = *fVolumeMediaMap;
 385     new(lvols[fNVolumes++])
 386         FlukaVolume(name, nmed);
 387
 388     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
 389 }
 390
 391 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 392                    Int_t iaxis) {
 393 //
 394     fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
 395 }
 396
 397 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 398                     Int_t iaxis, Double_t c0i, Int_t numed) {
 399 //
 400     TClonesArray &lvols = *fVolumeMediaMap;
 401     new(lvols[fNVolumes++])
 402         FlukaVolume(name, numed);
 403     fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 404 }
 405
 406 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 407                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 408 //
 409     TClonesArray &lvols = *fVolumeMediaMap;
 410     new(lvols[fNVolumes++])
 411         FlukaVolume(name, numed);
 412     fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 413 }
 414
 415 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 416                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 417 //
 418     TClonesArray &lvols = *fVolumeMediaMap;
 419     new(lvols[fNVolumes++])
 420         FlukaVolume(name, numed);
 421     fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 422 }
 423
 424 void TFluka::Gsord(const char *name, Int_t iax) {
 425 //
 426   fGeometryManager->Gsord(name, iax);
 427 }
 428
 429 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 430                    Double_t x, Double_t y, Double_t z, Int_t irot,
 431                    const char *konly) {
 432 //
 433   fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly);
 434 }
 435
 436 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 437                     Double_t x, Double_t y, Double_t z, Int_t irot,
 438                     const char *konly, Float_t *upar, Int_t np)  {
 439   //
 440   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 441 }
 442 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 443                     Double_t x, Double_t y, Double_t z, Int_t irot,
 444                     const char *konly, Double_t *upar, Int_t np)  {
 445   //
 446   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 447 }
 448
 449 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
 450 //
 451   fGeometryManager->Gsbool(onlyVolName, manyVolName);
 452 }
 453
 454 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
 455                          Float_t *absco, Float_t *effic, Float_t *rindex) {
 456 //
 457   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
 458 }
 459 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
 460                          Double_t *absco, Double_t *effic, Double_t *rindex) {
 461 //
 462   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
 463 }
 464
 465 // Euclid
 466 void TFluka::WriteEuclid(const char* fileName, const char* topVol,
 467                           Int_t number, Int_t nlevel) {
 468 //
 469   fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel);
 470 }
 471
 472
 473
 474 //_____________________________________________________________________________
 475 // methods needed by the stepping
 476 //____________________________________________________________________________
 477
 478 Int_t TFluka::GetMedium() const {
 479     FGeometryInit* flugg = FGeometryInit::GetInstance();
 480     return flugg->GetMedium(fCurrentFlukaRegion);
 481 }
 482
 483
 484
 485 //____________________________________________________________________________
 486 // ID <--> PDG transformations
 487 //_____________________________________________________________________________
 488 Int_t TFluka::IdFromPDG(Int_t pdg) const
 489 {
 490   //
 491   // Return Fluka code from PDG and pseudo ENDF code
 492
 493   // MCIHAD() goes from pdg to fluka internal.
 494   Int_t intfluka = mcihad(pdg);
 495   // KPTOIP array goes from internal to official
 496   return GetFlukaKPTOIP(intfluka);
 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   Int_t intfluka = GetFlukaIPTOKP(id);
 506   //MPKDHA() goes from internal to PDG
 507   return mpdgha(intfluka);
 508 }
 509
 510 //_____________________________________________________________________________
 511 // methods for step management
 512 //____________________________________________________________________________
 513 //
 514 // set methods
 515 //
 516 void TFluka::SetMaxStep(Double_t)
 517 {
 518 // SetMaxStep is dummy procedure in TFluka !
 519   cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
 520 }
 521
 522 void TFluka::SetMaxNStep(Int_t)
 523 {
 524 // SetMaxNStep is dummy procedure in TFluka !
 525   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
 526 }
 527
 528 void TFluka::SetUserDecay(Int_t)
 529 {
 530 // SetUserDecay is dummy procedure in TFluka !
 531   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
 532 }
 533
 534 //
 535 // dynamic properties
 536 //
 537 void TFluka::TrackPosition(TLorentzVector& position) const
 538 {
 539 // Return the current position in the master reference frame of the
 540 // track being transported
 541 // TRACKR.atrack = age of the particle
 542 // TRACKR.xtrack = x-position of the last point
 543 // TRACKR.ytrack = y-position of the last point
 544 // TRACKR.ztrack = z-position of the last point
 545   position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
 546   position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
 547   position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
 548   position.SetT(TRACKR.atrack);
 549 }
 550
 551 void TFluka::TrackMomentum(TLorentzVector& momentum) const
 552 {
 553 // Return the direction and the momentum (GeV/c) of the track
 554 // currently being transported
 555 // TRACKR.ptrack = momentum of the particle (not always defined, if
 556 //               < 0 must be obtained from etrack)
 557 // TRACKR.cx,y,ztrck = direction cosines of the current particle
 558 // TRACKR.etrack = total energy of the particle
 559 // TRACKR.jtrack = identity number of the particle
 560 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
 561   if (TRACKR.ptrack >= 0) {
 562     momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
 563     momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
 564     momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
 565     momentum.SetE(TRACKR.etrack);
 566     return;
 567   }
 568   else {
 569     Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
 570     momentum.SetPx(p*TRACKR.cxtrck);
 571     momentum.SetPy(p*TRACKR.cytrck);
 572     momentum.SetPz(p*TRACKR.cztrck);
 573     momentum.SetE(TRACKR.etrack);
 574     return;
 575   }
 576 }
 577
 578 Double_t TFluka::TrackStep() const
 579 {
 580 // Return the length in centimeters of the current step
 581 // TRACKR.ctrack = total curved path
 582     return TRACKR.ctrack;
 583 }
 584
 585 Double_t TFluka::TrackLength() const
 586 {
 587 // Still wrong !!!
 588 // This is the sum of substeps !!!
 589 // TRACKR.ctrack = total curved path of the current step
 590 // Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
 591 // The sum of all step length starting from the beginning of the track
 592 // for the time being returns only the length in centimeters of the current step
 593     Double_t sum = 0;
 594     for ( Int_t j=0;j<TRACKR.ntrack;j++) {
 595       sum +=TRACKR.ttrack[j];
 596     }
 597     return sum;
 598 }
 599
 600 Double_t TFluka::TrackTime() const
 601 {
 602 // Return the current time of flight of the track being transported
 603 // TRACKR.atrack = age of the particle
 604   return TRACKR.atrack;
 605 }
 606
 607 Double_t TFluka::Edep() const
 608 {
 609 // Energy deposition
 610 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
 611 // -->local energy deposition (the value and the point are not recorded in TRACKR)
 612 //    but in the variable "rull" of the procedure "endraw.cxx"
 613 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
 614 // -->no energy loss along the track
 615 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
 616 // -->energy loss distributed along the track
 617 // TRACKR.dtrack = energy deposition of the jth deposition even
 618   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
 619     return fRull;
 620   else {
 621     Double_t sum = 0;
 622     for ( Int_t j=0;j<TRACKR.mtrack;j++) {
 623       sum +=TRACKR.dtrack[j];
 624     }
 625     return sum;
 626   }
 627 }
 628
 629 Int_t TFluka::TrackPid() const
 630 {
 631 // Return the id of the particle transported
 632 // TRACKR.jtrack = identity number of the particle
 633   return PDGFromId(TRACKR.jtrack);
 634 }
 635
 636 Double_t TFluka::TrackCharge() const
 637 {
 638 // Return charge of the track currently transported
 639 // PAPROP.ichrge = electric charge of the particle
 640 // TRACKR.jtrack = identity number of the particle
 641   return PAPROP.ichrge[TRACKR.jtrack+6];
 642 }
 643
 644 Double_t TFluka::TrackMass() const
 645 {
 646 // PAPROP.am = particle mass in GeV
 647 // TRACKR.jtrack = identity number of the particle
 648   return PAPROP.am[TRACKR.jtrack+6];
 649 }
 650
 651 Double_t TFluka::Etot() const
 652 {
 653 // TRACKR.etrack = total energy of the particle
 654   return TRACKR.etrack;
 655 }
 656
 657 //
 658 // track status
 659 //
 660 Bool_t   TFluka::IsNewTrack() const
 661 {
 662 // ???????????????,
 663 // True if the track is not at the boundary of the current volume
 664 // Not true in some cases in bxdraw - to be solved
 665   return 1;
 666 }
 667
 668 Bool_t   TFluka::IsTrackInside() const
 669 {
 670 // True if the track is not at the boundary of the current volume
 671 // In Fluka a step is always inside one kind of material
 672 // If the step would go behind the region of one material,
 673 // it will be shortened to reach only the boundary.
 674 // Therefore IsTrackInside() is always true.
 675 // Not true in some cases in bxdraw - to be solved
 676   return 1;
 677 }
 678
 679 Bool_t   TFluka::IsTrackEntering() const
 680 {
 681 // True if this is the first step of the track in the current volume
 682 // Boundary- (X) crossing
 683 // Icode = 19: boundary crossing - call from Kaskad
 684 // Icode = 29: boundary crossing - call from Emfsco
 685 // Icode = 39: boundary crossing - call from Kasneu
 686 // Icode = 49: boundary crossing - call from Kashea
 687 // Icode = 59: boundary crossing - call from Kasoph
 688   if (fIcode == 19 ||
 689       fIcode == 29 ||
 690       fIcode == 39 ||
 691       fIcode == 49 ||
 692       fIcode == 59) return 1;
 693   else return 0;
 694 }
 695
 696 Bool_t   TFluka::IsTrackExiting() const
 697 {
 698 // True if this is the last step of the track in the current volume
 699 // Boundary- (X) crossing
 700 // Icode = 19: boundary crossing - call from Kaskad
 701 // Icode = 29: boundary crossing - call from Emfsco
 702 // Icode = 39: boundary crossing - call from Kasneu
 703 // Icode = 49: boundary crossing - call from Kashea
 704 // Icode = 59: boundary crossing - call from Kasoph
 705   if (fIcode == 19 ||
 706       fIcode == 29 ||
 707       fIcode == 39 ||
 708       fIcode == 49 ||
 709       fIcode == 59) return 1;
 710   else return 0;
 711 }
 712
 713 Bool_t   TFluka::IsTrackOut() const
 714 {
 715 // True if the track is out of the setup
 716 // means escape
 717 // Icode = 14: escape - call from Kaskad
 718 // Icode = 23: escape - call from Emfsco
 719 // Icode = 32: escape - call from Kasneu
 720 // Icode = 40: escape - call from Kashea
 721 // Icode = 51: escape - call from Kasoph
 722   if (fIcode == 14 ||
 723       fIcode == 23 ||
 724       fIcode == 32 ||
 725       fIcode == 40 ||
 726       fIcode == 51) return 1;
 727   else return 0;
 728 }
 729
 730 Bool_t   TFluka::IsTrackDisappeared() const
 731 {
 732 // means all inelastic interactions and decays
 733 // fIcode from usdraw
 734   if (fIcode == 101 || // inelastic interaction
 735       fIcode == 102 || // particle decay
 736       fIcode == 214 || // in-flight annihilation
 737       fIcode == 215 || // annihilation at rest
 738       fIcode == 217 || // pair production
 739       fIcode == 221) return 1;
 740   else return 0;
 741 }
 742
 743 Bool_t   TFluka::IsTrackStop() const
 744 {
 745 // True if the track energy has fallen below the threshold
 746 // means stopped by signal or below energy threshold
 747 // Icode = 12: stopping particle       - call from Kaskad
 748 // Icode = 15: time kill               - call from Kaskad
 749 // Icode = 21: below threshold, iarg=1 - call from Emfsco
 750 // Icode = 22: below threshold, iarg=2 - call from Emfsco
 751 // Icode = 24: time kill               - call from Emfsco
 752 // Icode = 31: below threshold         - call from Kasneu
 753 // Icode = 33: time kill               - call from Kasneu
 754 // Icode = 41: time kill               - call from Kashea
 755 // Icode = 52: time kill               - call from Kasoph
 756   if (fIcode == 12 ||
 757       fIcode == 15 ||
 758       fIcode == 21 ||
 759       fIcode == 22 ||
 760       fIcode == 24 ||
 761       fIcode == 31 ||
 762       fIcode == 33 ||
 763       fIcode == 41 ||
 764       fIcode == 52) return 1;
 765   else return 0;
 766 }
 767
 768 Bool_t   TFluka::IsTrackAlive() const
 769 {
 770 // means not disappeared or not out
 771   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
 772   else return 1;
 773 }
 774
 775 //
 776 // secondaries
 777 //
 778
 779 Int_t TFluka::NSecondaries() const
 780 // Number of secondary particles generated in the current step
 781 // FINUC.np = number of secondaries except light and heavy ions
 782 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
 783 {
 784   return FINUC.np + FHEAVY.npheav;
 785 }
 786
 787 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
 788                 TLorentzVector& position, TLorentzVector& momentum)
 789 {
 790   if (isec >= 0 && isec < FINUC.np) {
 791           // more fine condition depending on icode
 792           // icode = 100 ?
 793           // icode = 101 OK
 794           // icode = 102 OK
 795           // icode = 103 ?
 796           // icode = 104 ?
 797           // icode = 105 ?
 798           // icode = 208 ?
 799           // icode = 210 ?
 800           // icode = 212 ?
 801           // icode = 214 OK
 802           // icode = 215 OK
 803           // icode = 219 ?
 804           // icode = 221 OK
 805           // icode = 225 ?
 806           // icode = 300 ?
 807           // icode = 400 ?
 808
 809     particleId = PDGFromId(FINUC.kpart[isec]);
 810     position.SetX(fXsco);
 811     position.SetY(fYsco);
 812     position.SetZ(fZsco);
 813     position.SetT(TRACKR.atrack);
 814 //  position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
 815     momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
 816     momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
 817     momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
 818     momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
 819   }
 820   if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
 821     Int_t jsec = isec - FINUC.np;
 822     particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
 823     position.SetX(fXsco);
 824     position.SetY(fYsco);
 825     position.SetZ(fZsco);
 826     position.SetT(TRACKR.atrack);
 827 //  position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
 828     momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
 829     momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
 830     momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
 831     if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
 832       momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
 833     else if (FHEAVY.tkheav[jsec] > 6)
 834       momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
 835     }
 836 }
 837
 838 TMCProcess TFluka::ProdProcess(Int_t isec) const
 839 // Name of the process that has produced the secondary particles
 840 // in the current step
 841 {
 842   const TMCProcess kIpNoProc = kPNoProcess;
 843   const TMCProcess kIpPDecay = kPDecay;
 844   const TMCProcess kIpPPair = kPPair;
 845 //const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
 846 //const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
 847   const TMCProcess kIpPCompton = kPCompton;
 848   const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
 849   const TMCProcess kIpPBrem = kPBrem;
 850 //const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
 851 //const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
 852   const TMCProcess kIpPDeltaRay = kPDeltaRay;
 853 //const TMCProcess kIpPMoller = kPMoller;
 854 //const TMCProcess kIpPBhabha = kPBhabha;
 855   const TMCProcess kIpPAnnihilation = kPAnnihilation;
 856 //const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
 857 //const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
 858   const TMCProcess kIpPHadronic = kPHadronic;
 859   const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
 860   const TMCProcess kIpPPhotoFission = kPPhotoFission;
 861   const TMCProcess kIpPRayleigh = kPRayleigh;
 862   const TMCProcess kIpPCerenkov = kPCerenkov;
 863   const TMCProcess kIpPSynchrotron = kPSynchrotron;
 864
 865   Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
 866   if (fIcode == 102) return kIpPDecay;
 867   else if (fIcode == 104 || fIcode == 217) return kIpPPair;
 868 //else if (fIcode == 104) return kIpPairFromPhoton;
 869 //else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
 870   else if (fIcode == 219) return kIpPCompton;
 871   else if (fIcode == 221) return kIpPPhotoelectric;
 872   else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
 873 //else if (fIcode == 105) return kIpPBremFromHeavy;
 874 //else if (fIcode == 208) return kPBremFromElectronOrPositron;
 875   else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
 876   else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
 877 //else if (fIcode == 210) return kIpPMoller;
 878 //else if (fIcode == 212) return kIpPBhabha;
 879   else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
 880 //else if (fIcode == 214) return kIpPAnnihilInFlight;
 881 //else if (fIcode == 215) return kIpPAnnihilAtRest;
 882   else if (fIcode == 101) return kIpPHadronic;
 883   else if (fIcode == 101) {
 884     if (!mugamma) return kIpPHadronic;
 885     else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
 886     else return kIpPMuonNuclear;
 887   }
 888   else if (fIcode == 225) return kIpPRayleigh;
 889 // Fluka codes 100, 300 and 400 still to be investigasted
 890   else return kIpNoProc;
 891 }
 892
 893 //Int_t StepProcesses(TArrayI &proc) const
 894 // Return processes active in the current step
 895 //{
 896 //ck = total energy of the particl ????????????????
 897 //}
 898
 899
 900 // ===============================================================
 901 void TFluka::FutoTest()
 902 {
 903   Int_t icode, mreg, newreg, particleId;
 904 //  Int_t medium;
 905   Double_t rull, xsco, ysco, zsco;
 906   TLorentzVector position, momentum;
 907   icode = GetIcode();
 908   if (icode == 0) {
 909     cout << " icode=" << icode << endl;
 910     /*
 911     cout << "TLorentzVector positionX=" << position.X()
 912        << "positionY=" << position.Y()
 913        << "positionZ=" << position.Z()
 914        << "timeT=" << position.T() << endl;
 915     cout << "TLorentzVector momentumX=" << momentum.X()
 916        << "momentumY=" << momentum.Y()
 917        << "momentumZ=" << momentum.Z()
 918        << "energyE=" << momentum.E() << endl;
 919     cout << "TrackPid=" << TrackPid() << endl;
 920     */
 921   }
 922
 923   else if (icode > 0 && icode <= 5) {
 924 // mgdraw
 925     mreg = GetMreg();
 926 //    medium = GetMedium();
 927     cout << " icode=" << icode
 928          << " mreg=" << mreg
 929 //       << " medium=" << medium
 930          << endl;
 931   TrackPosition(position);
 932   TrackMomentum(momentum);
 933   cout << "TLorentzVector positionX=" << position.X()
 934        << "positionY=" << position.Y()
 935        << "positionZ=" << position.Z()
 936        << "timeT=" << position.T() << endl;
 937   cout << "TLorentzVector momentumX=" << momentum.X()
 938        << "momentumY=" << momentum.Y()
 939        << "momentumZ=" << momentum.Z()
 940        << "energyE=" << momentum.E() << endl;
 941   cout << "TrackStep=" << TrackStep() << endl;
 942   cout << "TrackLength=" << TrackLength() << endl;
 943   cout << "TrackTime=" << TrackTime() << endl;
 944   cout << "Edep=" << Edep() << endl;
 945   cout << "TrackPid=" << TrackPid() << endl;
 946   cout << "TrackCharge=" << TrackCharge() << endl;
 947   cout << "TrackMass=" << TrackMass() << endl;
 948   cout << "Etot=" << Etot() << endl;
 949   cout << "IsNewTrack=" << IsNewTrack() << endl;
 950   cout << "IsTrackInside=" << IsTrackInside() << endl;
 951   cout << "IsTrackEntering=" << IsTrackEntering() << endl;
 952   cout << "IsTrackExiting=" << IsTrackExiting() << endl;
 953   cout << "IsTrackOut=" << IsTrackOut() << endl;
 954   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
 955   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
 956   }
 957
 958   else if((icode >= 10 && icode <= 15) ||
 959           (icode >= 20 && icode <= 24) ||
 960           (icode >= 30 && icode <= 33) ||
 961           (icode >= 40 && icode <= 41) ||
 962           (icode >= 50 && icode <= 52)) {
 963 // endraw
 964     mreg = GetMreg();
 965 //    medium = GetMedium();
 966     rull = GetRull();
 967     xsco = GetXsco();
 968     ysco = GetYsco();
 969     zsco = GetZsco();
 970     cout << " icode=" << icode
 971          << " mreg=" << mreg
 972 //       << " medium=" << medium
 973          << " rull=" << rull
 974          << " xsco=" << xsco
 975          << " ysco=" << ysco
 976          << " zsco=" << zsco << endl;
 977   TrackPosition(position);
 978   TrackMomentum(momentum);
 979   cout << "Edep=" << Edep() << endl;
 980   cout << "Etot=" << Etot() << endl;
 981   cout << "TrackPid=" << TrackPid() << endl;
 982   cout << "TrackCharge=" << TrackCharge() << endl;
 983   cout << "TrackMass=" << TrackMass() << endl;
 984   cout << "IsTrackOut=" << IsTrackOut() << endl;
 985   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
 986   cout << "IsTrackStop=" << IsTrackStop() << endl;
 987   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
 988   }
 989
 990   else if((icode >= 100 && icode <= 105) ||
 991            (icode == 208) ||
 992            (icode == 210) ||
 993            (icode == 212) ||
 994            (icode >= 214 && icode <= 215) ||
 995            (icode == 217) ||
 996            (icode == 219) ||
 997            (icode == 221) ||
 998            (icode == 225) ||
 999            (icode == 300) ||
1000            (icode == 400)) {
1001 // usdraw
1002     mreg = GetMreg();
1003 //    medium = GetMedium();
1004     xsco = GetXsco();
1005     ysco = GetYsco();
1006     zsco = GetZsco();
1007     cout << " icode=" << icode
1008          << " mreg=" << mreg
1009 //       << " medium=" << medium
1010          << " xsco=" << xsco
1011          << " ysco=" << ysco
1012          << " zsco=" << zsco << endl;
1013     cout << "TrackPid=" << TrackPid() << endl;
1014     cout << "NSecondaries=" << NSecondaries() << endl;
1015     for (Int_t isec=0; isec< NSecondaries(); isec++) {
1016       TFluka::GetSecondary(isec, particleId, position, momentum);
1017       cout << "TLorentzVector positionX=" << position.X()
1018            << "positionY=" << position.Y()
1019            << "positionZ=" << position.Z()
1020            << "timeT=" << position.T() << endl;
1021       cout << "TLorentzVector momentumX=" << momentum.X()
1022            << "momentumY=" << momentum.Y()
1023            << "momentumZ=" << momentum.Z()
1024            << "energyE=" << momentum.E() << endl;
1025       cout << "TrackPid=" << particleId << endl;
1026
1027     }
1028   }
1029
1030   else if((icode == 19) ||
1031           (icode == 29) ||
1032           (icode == 39) ||
1033           (icode == 49) ||
1034           (icode == 59)) {
1035     mreg = GetMreg();
1036 //    medium = GetMedium();
1037     newreg = GetNewreg();
1038     xsco = GetXsco();
1039     ysco = GetYsco();
1040     zsco = GetZsco();
1041     cout << " icode=" << icode
1042          << " mreg=" << mreg
1043 //       << " medium=" << medium
1044          << " newreg=" << newreg
1045          << " xsco=" << xsco
1046          << " ysco=" << ysco
1047          << " zsco=" << zsco << endl;
1048   }
1049 //
1050 // ====================================================================
1051 //
1052
1053
1054
1055 } // end of FutoTest