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