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