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