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