TGeant4/TG4StepManager.cxx

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