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