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