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