TGeant4/TG4StepManager.cxx

   1 // $Id$
   2 // Category: event
   3 //
   4 // See the class description in the header file.
   5
   6 #include "TG4StepManager.h"
   7 #include "TG4GeometryManager.h"
   8 #include "TG4PhysicsManager.h"
   9 #include "TG4VSensitiveDetector.h"
  10 #include "TG4Globals.h"
  11 #include "TG3Units.h"
  12
  13 #include <G4Step.hh>
  14 #include <G4SteppingManager.hh>
  15 #include <G4UserLimits.hh>
  16 #include <G4ParticleTable.hh>
  17 #include <G4UImanager.hh>
  18 #include <G4AffineTransform.hh>
  19 //#include <G4TransportationManager.hh>
  20 //#include <G4Navigator.hh>
  21
  22 #include <Randomize.hh>
  23
  24 #include <TLorentzVector.h>
  25
  26 TG4StepManager* TG4StepManager::fgInstance = 0;
  27
  28 TG4StepManager::TG4StepManager()
  29   : fStep(0),
  30     fSteppingManager(0)
  31 {
  32 //
  33   if (fgInstance) {
  34     TG4Globals::Exception(
  35       "TG4StepManager: attempt to create two instances of singleton.");
  36   }
  37
  38   fgInstance = this;
  39 }
  40
  41 TG4StepManager::TG4StepManager(const TG4StepManager& right) {
  42 //
  43   TG4Globals::Exception(
  44     "Attempt to copy TG4StepManager singleton.");
  45 }
  46
  47 TG4StepManager::~TG4StepManager() {
  48 //
  49 }
  50
  51 // operators
  52
  53 TG4StepManager& TG4StepManager::operator=(const TG4StepManager& right)
  54 {
  55   // check assignement to self
  56   if (this == &right) return *this;
  57
  58   TG4Globals::Exception(
  59     "Attempt to assign TG4StepManager singleton.");
  60
  61   return *this;
  62 }
  63
  64 // private methods
  65
  66 void TG4StepManager::SetTLorentzVector(G4ThreeVector xyz, G4double t,
  67                                        TLorentzVector& lv) const
  68 {
  69 // Fills TLorentzVector with G4ThreeVector and G4double.
  70 // ---
  71
  72    lv[0] = xyz.x();
  73    lv[1] = xyz.y();
  74    lv[2] = xyz.z();
  75    lv[3] = t;
  76 }
  77
  78 // public methods
  79
  80 void TG4StepManager::StopTrack()
  81 {
  82 // Stops the current track and skips to the next.
  83 // ?? do we want to invoke rest processes?
  84 // ?? do we want to stop secondaries too?
  85 //   possible "stop" track status from G4:
  86 //   fStopButAlive,      // Invoke active rest physics processes and
  87 //                       // and kill the current track afterward
  88 //   fStopAndKill,       // Kill the current track
  89 //   fKillTrackAndSecondaries,
  90 //                       // Kill the current track and also associated
  91 //                       // secondaries.
  92 // ---
  93
  94   if (fStep) {
  95     // check
  96     fStep->GetTrack()->SetTrackStatus(fStopAndKill);
  97     // fStep->GetTrack()->SetTrackStatus(fStopButAlive);
  98     // fStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
  99   }
 100   else {
 101     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 102   }
 103 }
 104
 105 void TG4StepManager::StopEvent()
 106 {
 107 // Aborts the current event processing.
 108 // ---
 109
 110   if (fStep) {
 111     fStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
 112       //StopTrack();   // cannot be used as it keeps secondaries
 113     G4UImanager::GetUIpointer()->ApplyCommand("/event/abort");
 114     G4UImanager::GetUIpointer()->ApplyCommand("/alStacking/clearStack");
 115   }
 116   else {
 117     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 118   }
 119 }
 120
 121 void TG4StepManager::Rndm(Float_t* array, const Int_t size) const
 122 {
 123 // Random numbers array of the specified size.
 124 // ---
 125
 126   G4double* const kpDoubleArray = new G4double[size];
 127   RandFlat::shootArray(size,kpDoubleArray);
 128   for (G4int i=0; i<size; i++) {
 129     array[i] = kpDoubleArray[i];
 130   }
 131   delete [] kpDoubleArray;
 132 }
 133
 134 void TG4StepManager::SetMaxStep(Float_t step)
 135 {
 136 // Maximum step allowed in the current logical volume.
 137 // The maximum value is kept for following tracks - is it ok ??
 138 // ---
 139
 140   if (fStep) {
 141     // check this
 142     G4LogicalVolume* curLogVolume
 143       = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
 144     G4UserLimits* userLimits
 145       = curLogVolume->GetUserLimits();
 146     if (userLimits == 0)
 147     {
 148       userLimits = new G4UserLimits(step);
 149       curLogVolume->SetUserLimits(userLimits);
 150     }
 151     else
 152     { userLimits->SetMaxAllowedStep(step); }
 153   }
 154   else {
 155     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 156   }
 157 }
 158
 159 void TG4StepManager::SetMaxNStep(Int_t maxNofSteps)
 160 {
 161 // Not yet implemented.
 162 // ---
 163
 164   TG4Globals::Warning(
 165     "TG4StepManager::SetMaxNStep(..) is not yet implemented.");
 166 }
 167
 168 void TG4StepManager::SetUserDecay(Int_t pdg)
 169 {
 170 // Not yet implemented.
 171 // ---
 172
 173   TG4Globals::Exception(
 174     "TG4StepManager::SetUserDecay(..) is not yet implemented.");
 175 }
 176
 177 Int_t TG4StepManager::CurrentVolID(Int_t& copyNo) const
 178 {
 179 // Returns the current sensitive detector ID
 180 // and the copy number of the current physical volume.
 181 // ---
 182
 183   if (fStep) {
 184     // check this
 185     G4VPhysicalVolume* physVolume
 186       = fStep->GetPreStepPoint()->GetPhysicalVolume();
 187     copyNo = physVolume->GetCopyNo();
 188
 189     // sensitive detector ID
 190     G4VSensitiveDetector* sd
 191         = physVolume->GetLogicalVolume()->GetSensitiveDetector();
 192     if (sd) {
 193       TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
 194       if (tsd)
 195         return tsd->GetID();
 196       else {
 197         TG4Globals::Exception(
 198           "TG4StepManager::CurrentVol: Unknown sensitive detector type");
 199         return 0;
 200       }
 201     }
 202     else {
 203       G4String text = "TG4StepManager::CurrentVol: \n";
 204       text = text + "    Volume " + physVolume->GetName();
 205       text = text + " has not a sensitive detector.";
 206       TG4Globals::Exception(text);
 207       return 0;
 208     }
 209   }
 210   else {
 211     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 212     return 0;
 213   }
 214 }
 215
 216 Int_t TG4StepManager::CurrentVolOffID(Int_t off, Int_t&  copyNo) const
 217 {
 218 // Returns the off-th mother's of the current volume
 219 // the sensitive detector ID and the copy number.
 220 // ---
 221
 222   if (off == 0) return CurrentVolID(copyNo);
 223
 224   if (fStep) {
 225     // check this !!
 226     // mother of physical volume may not be set ?!!
 227
 228     G4VPhysicalVolume* physVolume
 229       = fStep->GetPreStepPoint()->GetPhysicalVolume();
 230
 231     G4VPhysicalVolume* mother = 0;
 232     Int_t level = off;
 233     while (level > 0) {
 234       mother = physVolume->GetMother();
 235       level--;
 236     }
 237
 238     if (mother) {
 239       copyNo = mother->GetCopyNo();
 240
 241       // sensitive detector ID
 242       G4VSensitiveDetector* sd
 243         = mother->GetLogicalVolume()->GetSensitiveDetector();
 244       if (sd) {
 245         TG4VSensitiveDetector* tsd = dynamic_cast<TG4VSensitiveDetector*>(sd);
 246         if (tsd)
 247           return tsd->GetID();
 248         else {
 249           TG4Globals::Exception(
 250             "TG4StepManager::CurrentVolOff: Unknown sensitive detector type");
 251           return 0;
 252         }
 253       }
 254       else {
 255         G4String text = "TG4StepManager::CurrentVolOff: \n";
 256         text = text + "    Volume " + mother->GetName();
 257         text = text + " has not a sensitive detector.";
 258         TG4Globals::Exception(text);
 259         return 0;
 260       }
 261     }
 262     else {
 263       G4String text = "TG4StepManager::CurrentVolOff: Volume ";
 264       text = text + physVolume->GetName();
 265       text = text + " has not defined mother in required level.";
 266       TG4Globals::Exception(text);
 267       return 0;
 268     }
 269   }
 270   else {
 271     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 272     return 0;
 273   }
 274 }
 275
 276 const char* TG4StepManager::CurrentVolName() const
 277 {
 278 // Returns the current physical volume name.
 279 // ---
 280
 281   if (fStep) {
 282     G4VPhysicalVolume* physVolume
 283       = fStep->GetPreStepPoint()->GetPhysicalVolume();
 284     return physVolume->GetName();
 285   }
 286   else {
 287     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 288     return 0;
 289   }
 290 }
 291
 292 const char* TG4StepManager::CurrentVolOffName(Int_t off) const
 293 {
 294 // Returns the off-th mother's physical volume name.
 295 // ---
 296
 297   if (fStep) {
 298     G4VPhysicalVolume* physVolume
 299       = fStep->GetPreStepPoint()->GetPhysicalVolume();
 300
 301     G4VPhysicalVolume* mother = 0;
 302     Int_t level = off;
 303     while (level > 0) {
 304       mother = physVolume->GetMother();
 305       level--;
 306     }
 307
 308     if (mother) {
 309       return mother->GetName();
 310     }
 311     else {
 312       TG4Globals::Exception("TG4StepManager::CurrentVolOff: wrong usage.");
 313       return 0;
 314     }
 315   }
 316   else {
 317     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 318     return 0;
 319   }
 320 }
 321
 322 Int_t TG4StepManager::CurrentMaterial(Float_t &a, Float_t &z, Float_t &dens,
 323                           Float_t &radl, Float_t &absl) const
 324 {
 325 // Returns the parameters of the current material during transport
 326 // the return value is the number of elements in the mixture.
 327 // ---
 328
 329   if (fStep) {
 330     // check this
 331     G4LogicalVolume* curLogVolume
 332       = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
 333     G4Material* material
 334       = curLogVolume->GetMaterial();
 335
 336     // this if may be redundant - check
 337     if (material)
 338     {
 339       G4int nofElements = material->GetNumberOfElements();
 340       TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
 341       a = pGeometryManager->GetEffA(material);
 342       z = pGeometryManager->GetEffZ(material);
 343
 344       // density
 345       dens = material->GetDensity();
 346       dens /= TG3Units::MassDensity();
 347
 348       // radiation length
 349       radl = material->GetRadlen();
 350       radl /= TG3Units::Length();
 351
 352       absl = 0.;  // this parameter is not defined in Geant4
 353       return nofElements;
 354     }
 355     else {
 356       TG4Globals::Exception(
 357        "TG4StepManager::CurrentMaterial(..): material is not defined.");
 358       return 0;
 359     }
 360   }
 361   else {
 362     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 363     return 0;
 364   }
 365 }
 366
 367 void TG4StepManager::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
 368 {
 369 // Transforms a position from the world reference frame
 370 // to the current volume reference frame.
 371 //
 372 //  Geant3 desription:
 373 //  ==================
 374 //       Computes coordinates XD (in DRS)
 375 //       from known coordinates XM in MRS
 376 //       The local reference system can be initialized by
 377 //         - the tracking routines and GMTOD used in GUSTEP
 378 //         - a call to GMEDIA(XM,NUMED)
 379 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
 380 //             (inverse routine is GDTOM)
 381 //
 382 //        If IFLAG=1  convert coordinates
 383 //           IFLAG=2  convert direction cosinus
 384 //
 385 // ---
 386
 387   if (fStep) {
 388
 389     G4ThreeVector theGlobalPoint(xm[0],xm[1],xm[2]);
 390
 391     //const G4NavigationHistory* history
 392     //   =  fStep->GetPreStepPoint()->GetTouchable()->GetHistory();
 393     G4AffineTransform affineTransform
 394       = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
 395         ->GetTopTransform();
 396
 397     G4ThreeVector theLocalPoint;
 398     if(iflag == 1)
 399          theLocalPoint = affineTransform.TransformPoint(theGlobalPoint);
 400     else if ( iflag == 2)
 401          theLocalPoint = affineTransform.TransformAxis(theGlobalPoint);
 402     else
 403       TG4Globals::Exception(
 404         "TG4StepManager::Gmtod(..,iflag): iflag is not in 1..2");
 405
 406     xd[0] = theLocalPoint.x();
 407     xd[1] = theLocalPoint.y();
 408     xd[2] = theLocalPoint.z();
 409
 410  /*
 411    // does not work ???
 412    G4ThreeVector direction(0,0,0);
 413    G4bool RelativeSearch = true;
 414    G4Navigator* theNavigator =
 415      G4TransportationManager::GetTransportationManager()->
 416      GetNavigatorForTracking();
 417
 418    G4VPhysicalVolume* pPhysVol;
 419    pPhysVol
 420      = LocateGlobalPointAndSetup(theGlobalPoint, &direction, RelativeSearch);
 421    //LocateGlobalPointWithinVolume(theGlobalPoint);
 422
 423    G4AffineTransform at
 424      = theNavigator->GetGlobalToLocalTransform();
 425  */
 426   }
 427   else {
 428     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 429   }
 430 }
 431
 432 void TG4StepManager::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
 433 {
 434 // Transforms a position from the current volume reference frame
 435 // to the world reference frame.
 436 //
 437 //  Geant3 desription:
 438 //  ==================
 439 //  Computes coordinates XM (Master Reference System
 440 //  knowing the coordinates XD (Detector Ref System)
 441 //  The local reference system can be initialized by
 442 //    - the tracking routines and GDTOM used in GUSTEP
 443 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
 444 //        (inverse routine is GMTOD)
 445 //
 446 //   If IFLAG=1  convert coordinates
 447 //      IFLAG=2  convert direction cosinus
 448 //
 449 // ---
 450
 451   if (fStep) {
 452     // check this
 453
 454    G4ThreeVector theLocalPoint(xd[0],xd[1],xd[2]);
 455
 456    G4AffineTransform affineTransform
 457       = fStep->GetPreStepPoint()->GetTouchable()->GetHistory()
 458         ->GetTopTransform().Inverse();
 459
 460    G4ThreeVector theGlobalPoint;
 461    if(iflag == 1)
 462          theGlobalPoint = affineTransform.TransformPoint(theLocalPoint);
 463    else if( iflag == 2)
 464          theGlobalPoint = affineTransform.TransformAxis(theLocalPoint);
 465    else
 466       TG4Globals::Warning(
 467         "TG4StepManager::Gdtom(...,iflag): iflag is not in 1..2");
 468
 469    xm[0] = theGlobalPoint.x();
 470    xm[1] = theGlobalPoint.y();
 471    xm[2] = theGlobalPoint.z();
 472   }
 473   else {
 474     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 475   }
 476 }
 477
 478 Float_t TG4StepManager::MaxStep() const
 479 {
 480 // Returns maximum step allowed in the current logical volume
 481 // by User Limits.
 482 // ---
 483
 484   if (fStep) {
 485     // check this
 486     G4LogicalVolume* curLogVolume
 487       = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
 488     G4UserLimits* userLimits
 489       = curLogVolume->GetUserLimits();
 490     G4double maxStep;
 491     if (userLimits == 0)
 492     {
 493       G4String text = "User Limits are not defined for log volume ";
 494       text = text + curLogVolume->GetName();
 495       TG4Globals::Warning(text);
 496       return DBL_MAX;
 497     }
 498     else
 499     {
 500       const G4Track& trackRef = *(fStep->GetTrack());
 501       maxStep = userLimits->GetMaxAllowedStep(trackRef);
 502       maxStep /= TG3Units::Length();
 503       return maxStep;
 504     }
 505   }
 506   else {
 507     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 508     return DBL_MAX;
 509   }
 510 }
 511
 512 Int_t TG4StepManager::GetMaxNStep() const
 513 {
 514 // Not yet implemented.
 515 // ---
 516
 517   TG4Globals::Warning(
 518     "Method GetMaxNStep is not yet implemented in TG4StepManager.");
 519   return 0;
 520 }
 521
 522 void TG4StepManager::TrackPosition(TLorentzVector& position) const
 523 {
 524 // Current particle position (in the world reference frame)
 525 // and the local time since the current track is created
 526 // (position of the PostStepPoint).
 527 // ---
 528
 529   if (fStep) {
 530     // position
 531     G4ThreeVector positionVector
 532        = fStep->GetPostStepPoint()->GetPosition();
 533     positionVector *= 1./(TG3Units::Length());
 534
 535     // local time
 536     G4double time
 537       = fStep->GetTrack()->GetLocalTime();
 538     time /= TG3Units::Time();
 539
 540     SetTLorentzVector(positionVector, time, position);
 541   }
 542   else {
 543     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 544   }
 545 }
 546
 547 Int_t TG4StepManager::GetMedium() const
 548 {
 549 // Returns the second index of the current material (corresponding to
 550 // G3 tracking medium index).
 551 // ---
 552
 553   if (fStep) {
 554     // current material
 555     G4Material* curMaterial
 556       = fStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();
 557
 558     // medium index
 559     TG4GeometryManager* pGeometryManager = TG4GeometryManager::Instance();
 560     return pGeometryManager->GetMediumId(curMaterial);
 561   }
 562   else {
 563     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 564     return 0;
 565   }
 566 }
 567
 568 void TG4StepManager::TrackMomentum(TLorentzVector& momentum) const
 569 {
 570 // Current particle "momentum" (px, py, pz, Etot).
 571 // ---
 572
 573   if (fStep) {
 574     G4ThreeVector momentumVector
 575       = fStep->GetTrack()->GetMomentum();
 576
 577     G4double energy
 578       = fStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
 579     energy /= TG3Units::Energy();
 580
 581     SetTLorentzVector(momentumVector, energy, momentum);
 582   }
 583   else {
 584     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 585   }
 586 }
 587
 588 void TG4StepManager::TrackVertexPosition(TLorentzVector& position) const
 589 {
 590 // The vertex particle position (in the world reference frame)
 591 // and the local time since the current track is created.
 592 // ---
 593
 594   if (fStep) {
 595     // position
 596     G4ThreeVector positionVector
 597        = fStep->GetTrack()->GetVertexPosition();
 598     positionVector *= 1./(TG3Units::Length());
 599
 600     // local time
 601     // to be checked
 602     G4double time
 603       = fStep->GetTrack()->GetLocalTime();
 604     time /= TG3Units::Time();
 605
 606     SetTLorentzVector(positionVector, time, position);
 607   }
 608   else {
 609     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 610   }
 611 }
 612
 613 void TG4StepManager::TrackVertexMomentum(TLorentzVector& momentum) const
 614 {
 615 // The vertex particle "momentum" (px, py, pz, Ekin)
 616 // to do: change Ekin -> Etot
 617 // ---
 618
 619   if (fStep) {
 620     G4ThreeVector momentumVector
 621       = fStep->GetTrack()->GetVertexMomentumDirection();
 622
 623     G4double energy
 624       = fStep->GetTrack()->GetVertexKineticEnergy();
 625     energy /= TG3Units::Energy();
 626
 627     SetTLorentzVector(momentumVector, energy, momentum);
 628   }
 629   else {
 630     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 631   }
 632 }
 633
 634 Float_t TG4StepManager::TrackStep() const
 635 {
 636 // Returns the current step length.
 637 // ---
 638
 639   if (fStep) {
 640     // check
 641     G4double length = fStep->GetStepLength();
 642     length /= TG3Units::Length();
 643     return length;
 644   }
 645   else {
 646     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 647     return 0.;
 648   }
 649 }
 650
 651 Float_t TG4StepManager::TrackLength() const
 652 {
 653 // Returns the length of the current track from its origin.
 654 // ---
 655
 656   if (fStep) {
 657     G4double length
 658       = fStep->GetTrack()->GetTrackLength();
 659     length /= TG3Units::Length();
 660     return length;
 661   }
 662   else {
 663     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 664     return 0;
 665   }
 666 }
 667
 668 Float_t TG4StepManager::TrackTime() const
 669 {
 670 // Returns the local time since the current track is created.
 671 // Comment:
 672 // in Geant4: there is also defined proper time as
 673 // the proper time of the dynamical particle of the current track.
 674 // ---
 675
 676   if (fStep) {
 677     G4double time
 678       = fStep->GetTrack()->GetLocalTime();
 679     time /= TG3Units::Time();
 680     return time;
 681   }
 682   else {
 683     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 684     return 0;
 685   }
 686 }
 687
 688 Float_t TG4StepManager::Edep() const
 689 {
 690 // Returns total energy deposit in this step.
 691 // ---
 692
 693   if (fStep) {
 694     // G4double deltaEnergy = fStep->GetDeltaEnergy();
 695     G4double deltaEnergy = fStep->GetTotalEnergyDeposit();
 696     deltaEnergy /= TG3Units::Energy();
 697     return deltaEnergy;
 698   }
 699   else {
 700     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 701     return 0.;
 702   }
 703 }
 704
 705 Int_t TG4StepManager::TrackPid() const
 706 {
 707 // Returns the current particle PDG encoding.
 708 // ---
 709
 710   if (fStep) {
 711     G4ParticleDefinition* particle
 712       = fStep->GetTrack()->GetDynamicParticle()->GetDefinition();
 713
 714     // ask TG4PhysicsManager to get PDG encoding
 715     // (in order to get PDG from extended TDatabasePDG
 716     // in case the standard PDG code is not defined)
 717     TG4PhysicsManager* pPhysicsManager = TG4PhysicsManager::Instance();
 718     G4int pdgEncoding = pPhysicsManager->GetPDGEncodingFast(particle);
 719
 720     return pdgEncoding;
 721   }
 722   else {
 723     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 724     return 0;
 725   }
 726 }
 727
 728 Float_t TG4StepManager::TrackCharge() const
 729 {
 730 // Returns the current particle charge.
 731 // ---
 732
 733   if (fStep) {
 734     G4double charge
 735       = fStep->GetTrack()->GetDynamicParticle()->GetDefinition()
 736         ->GetPDGCharge();
 737     charge /= TG3Units::Charge();
 738     return charge;
 739   }
 740   else {
 741     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 742     return 0;
 743   }
 744 }
 745
 746 Float_t TG4StepManager::TrackMass() const
 747 {
 748 // Returns current particle rest mass.
 749 // ---
 750
 751   if (fStep) {
 752     G4double mass
 753       = fStep->GetTrack()->GetDynamicParticle()->GetDefinition()
 754         ->GetPDGMass();
 755     mass /= TG3Units::Mass();
 756     return mass;
 757   }
 758   else {
 759     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 760     return 0;
 761   }
 762 }
 763
 764 Float_t TG4StepManager::Etot() const
 765 {
 766 // Returns total energy of the current particle.
 767 // ---
 768
 769   if (fStep) {
 770     G4double energy
 771       = fStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
 772     energy /= TG3Units::Energy();
 773     return energy;
 774   }
 775   else {
 776     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 777     return 0.;
 778   }
 779 }
 780
 781 Bool_t TG4StepManager::IsTrackInside() const
 782 {
 783 // Returns true if particle does not cross geometrical boundary
 784 // at both pre-step and post-step points.
 785 // ---
 786
 787   if ( !(IsTrackEntering()) && !(IsTrackExiting()) )
 788   { return true; }
 789   else
 790   { return false; }
 791 }
 792
 793 Bool_t TG4StepManager::IsTrackEntering() const
 794 {
 795 // Returns true if particle cross a geometrical boundary
 796 // at pre-step point.
 797 // ---
 798
 799   if (fStep) {
 800     G4StepStatus status
 801        = fStep->GetPreStepPoint()->GetStepStatus();
 802     if (status == fGeomBoundary)
 803     { return true; }
 804     else
 805     { return false; }
 806   }
 807   else {
 808     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 809     return false;
 810   }
 811 }
 812
 813 Bool_t TG4StepManager::IsTrackExiting() const
 814 {
 815 // Returns true if particle cross a geometrical boundary
 816 // at post-step point.
 817 // ---
 818
 819   if (fStep) {
 820     // check
 821     G4StepStatus status
 822        = fStep->GetPostStepPoint()->GetStepStatus();
 823     if (status == fGeomBoundary)
 824     { return true; }
 825     else
 826     { return false; }
 827   }
 828   else {
 829     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 830     return false;
 831   }
 832 }
 833
 834 Bool_t TG4StepManager::IsTrackOut() const
 835 {
 836 // Returns true if particle cross the Hall frame
 837 // at post-step point.
 838 // ---
 839
 840   if (fStep) {
 841     // check
 842     G4String volName
 843        = fStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
 844     if (volName == "pHallFrame")
 845     { return true; }
 846     else
 847     { return false; }
 848   }
 849   else {
 850     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 851     return false;
 852   }
 853 }
 854
 855 Bool_t TG4StepManager::IsTrackStop() const
 856 {
 857 // Returns true if particle has stopped
 858 // or has been killed, suspended or postponed to next event.
 859 //
 860 // Possible track status from G4:
 861 //   fAlive,             // Continue the tracking
 862 //   fStopButAlive,      // Invoke active rest physics processes and
 863 //                            // and kill the current track afterward
 864 //   fStopAndKill,       // Kill the current track
 865 //   fKillTrackAndSecondaries,
 866 //                       // Kill the current track and also associated
 867 //                       // secondaries.
 868 //   fSuspend,           // Suspend the current track
 869 //   fPostponeToNextEvent
 870 //                       // Postpones the tracking of thecurrent track
 871 //                       // to the next event.
 872 // ---
 873
 874   if (fStep) {
 875     // check
 876     G4TrackStatus status
 877        = fStep->GetTrack()->GetTrackStatus();
 878     if ((status == fStopAndKill) ||
 879         (status == fKillTrackAndSecondaries) ||
 880         (status == fSuspend) ||
 881         (status == fPostponeToNextEvent))
 882     { return true; }
 883     else
 884     { return false; }
 885   }
 886   else {
 887     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 888     return false;
 889   }
 890 }
 891
 892 Bool_t TG4StepManager::IsTrackDisappeared() const
 893 {
 894 // Returns true if particle has disappeared
 895 // (due to any physical process)
 896 // or has been killed, suspended or postponed to next event.
 897 // ---
 898
 899   if (fStep) {
 900     // check
 901     G4TrackStatus status
 902        = fStep->GetTrack()->GetTrackStatus();
 903     if ((status == fStopButAlive) ||
 904         (status == fKillTrackAndSecondaries) ||
 905         (status == fSuspend) ||
 906         (status == fPostponeToNextEvent))
 907     { return true; }
 908     else
 909     { return false; }
 910   }
 911   else {
 912     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 913     return false;
 914   }
 915 }
 916
 917 Bool_t TG4StepManager::IsTrackAlive() const
 918 {
 919 // Returns true if particle continues tracking.
 920 // ---
 921
 922   if (fStep) {
 923     G4TrackStatus status
 924        = fStep->GetTrack()->GetTrackStatus();
 925     if (status == fAlive)
 926     { return true; }
 927     else
 928     { return false; }
 929   }
 930   else {
 931     TG4Globals::Exception("TG4StepManager: Step is not defined.");
 932     return false;
 933   }
 934 }
 935
 936 Bool_t TG4StepManager::IsNewTrack() const
 937 {
 938 // Returns true when track performs the first step.
 939 // ---
 940
 941   if (fStep->GetTrack()->GetCurrentStepNumber() == 1)
 942     return true;
 943   else
 944     return false;
 945 }
 946
 947 Int_t TG4StepManager::NSecondaries() const
 948 {
 949 // Returns the number of secondary particles generated
 950 // in the current step.
 951 // ---
 952
 953   if (fSteppingManager)
 954   {
 955     G4int nofSecondaries = 0;
 956     nofSecondaries += fSteppingManager->GetfN2ndariesAtRestDoIt();
 957     nofSecondaries += fSteppingManager->GetfN2ndariesAlongStepDoIt();
 958     nofSecondaries += fSteppingManager->GetfN2ndariesPostStepDoIt();
 959
 960     return nofSecondaries;
 961   }
 962   else {
 963     TG4Globals::Exception("TG4StepManager: SteppingManager is not defined.");
 964     return 0;
 965   }
 966 }
 967
 968 void TG4StepManager::GetSecondary(Int_t index, Int_t& particleId,
 969                           TLorentzVector& position, TLorentzVector& momentum)
 970 {
 971 // Fills the parameters (particle encoding, position, momentum)
 972 // of the generated secondary particle which is specified by index.
 973 // !! Check if indexing of secondaries is same !!
 974 // ---
 975
 976   if (fSteppingManager)
 977   {
 978     G4int nofSecondaries = NSecondaries();
 979     G4TrackVector* secondaryTracks = fSteppingManager->GetSecondary();
 980
 981     if (secondaryTracks)
 982     {
 983       if (index < nofSecondaries)
 984       {
 985         // the index of the first secondary of this step
 986         G4int startIndex
 987           = secondaryTracks->entries() - nofSecondaries;
 988                // (the secondaryTracks vector contains secondaries
 989                // produced by the track at previous steps, too)
 990         G4Track* track
 991           = (*secondaryTracks)[startIndex + index];
 992
 993         // particle encoding
 994         particleId
 995           = track->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
 996
 997         // position & time
 998         G4ThreeVector positionVector = track->GetPosition();
 999         positionVector *= 1./(TG3Units::Length());
1000         G4double time = track->GetLocalTime();
1001         time /= TG3Units::Time();
1002         SetTLorentzVector(positionVector, time, position);
1003
1004         // momentum & energy
1005         G4ThreeVector momentumVector = track->GetMomentum();
1006         G4double energy = track->GetDynamicParticle()->GetTotalEnergy();
1007         energy /= TG3Units::Energy();
1008         SetTLorentzVector(momentumVector, energy, momentum);
1009       }
1010       else {
1011         TG4Globals::Exception(
1012           "TG4StepManager::GetSecondary(): wrong secondary track index.");
1013       }
1014     }
1015     else {
1016       TG4Globals::Exception(
1017         "TG4StepManager::GetSecondary(): secondary tracks vector is empty");
1018     }
1019   }
1020   else {
1021     TG4Globals::Exception("TG4StepManager: SteppingManager is not defined.");
1022   }
1023 }
1024
1025 const char* TG4StepManager::ProdProcess() const
1026 {
1027 // Returns the name of the process that defined current step
1028 // (and may produce the secondary particles).
1029 // ---
1030
1031   if (fStep)
1032   {
1033     const G4VProcess* curProcess
1034       = fStep->GetPostStepPoint()->GetProcessDefinedStep();
1035
1036     G4String g4Name = curProcess->GetProcessName();
1037     return g4Name;
1038   }
1039   else {
1040     TG4Globals::Exception("TG4StepManager: Step is not defined.");
1041     return "";
1042   }
1043 }