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