Belle II Software  release-05-01-25
SensitiveDetector.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2015 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter, Igal Jaegle *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <beast/fangs/simulation/SensitiveDetector.h>
12 #include <framework/gearbox/Unit.h>
13 #include <G4Track.hh>
14 #include <G4Step.hh>
15 
16 #include <array>
17 
18 namespace Belle2 {
24  namespace fangs {
25  namespace {
26  std::array<float, 3> vecToFloat(const G4ThreeVector& vec)
27  {
28  return std::array<float, 3> {{(float)vec.x(), (float)vec.y(), (float)vec.z()}};
29  }
30  }
31 
33  Simulation::SensitiveDetectorBase("FANGSSensitiveDetector", Const::invalidDetector)
34  {
35  //Register all collections we want to modify and require those we want to use
36  m_mcParticles.registerInDataStore();
37  m_simHits.registerInDataStore();
38  m_relMCSimHit.registerInDataStore();
39 
40  //Register the Relation so that the TrackIDs get replaced by the actual
41  //MCParticle indices after simulating the events. This is needed as
42  //secondary particles might not be stored so everything relating to those
43  //particles will be attributed to the last saved mother particle
45  }
46 
47  bool SensitiveDetector::step(G4Step* step, G4TouchableHistory*)
48  {
49  //Get scintilator ladder and sensor number
50  const G4TouchableHistory* touchable = dynamic_cast<const G4TouchableHistory*>(step->GetPreStepPoint()->GetTouchable());
51  int ladderID = touchable->GetVolume(2)->GetCopyNo();
52  int sensorID = touchable->GetVolume(1)->GetCopyNo();
53 
54  //Get Track information
55  const G4Track& track = *step->GetTrack();
56  const int trackID = track.GetTrackID();
57  const int pdgCode = step->GetTrack()->GetDefinition()->GetPDGEncoding();
58  // deposited Energy in Geant4 Units
59  double depEnergy = step->GetTotalEnergyDeposit();
60  // convert into Belle2 Units
61  depEnergy *= Unit::GeV / CLHEP::GeV;
62 
63  // get pre and post step point
64  const G4StepPoint& preStep = *step->GetPreStepPoint();
65  const G4StepPoint& postStep = *step->GetPostStepPoint();
66 
67  // check if this is the same sensor traversal, otherwise add one to the stack
68  if (m_tracks.empty() || (!m_tracks.top().check(trackID, ladderID, sensorID))) {
69  m_tracks.push(SensorTraversal());
70  }
71  // get the top of the stack
72  SensorTraversal& traversal = m_tracks.top();
73 
74  //If new track, remember values
75  if (traversal.getTrackID() == 0) {
76  bool isPrimary = Simulation::TrackInfo::getInfo(track).hasStatus(MCParticle::c_PrimaryParticle);
77  //Add start position
78  const G4ThreeVector preStepPos = preStep.GetPosition() / CLHEP::mm * Unit::mm;
79  const G4ThreeVector preStepMom = preStep.GetMomentum() / CLHEP::MeV * Unit::MeV;
80  //const G4AffineTransform& localToGlobalTransform = preStep.GetTouchableHandle()->GetHistory()->GetTopTransform().Inverse();
81  const G4AffineTransform& localToGlobalTransform = preStep.GetTouchableHandle()->GetHistory()->GetTopTransform();
82  const G4ThreeVector localpreStepPos = localToGlobalTransform.TransformPoint(preStep.GetPosition()) / CLHEP::mm * Unit::mm;
83  const double time = preStep.GetGlobalTime() / CLHEP::ns * Unit::ns;
84  traversal.setInitial(trackID, ladderID, sensorID, pdgCode, isPrimary, preStepPos, localpreStepPos, preStepMom, time);
85  //Remember if the track came from the outside
86  if (preStep.GetStepStatus() == fGeomBoundary) traversal.hasEntered();
87  }
88  //Add current step
89  const G4ThreeVector postStepPos = postStep.GetPosition() / CLHEP::mm * Unit::mm;
90  const double length = step->GetStepLength() / CLHEP::cm * Unit::cm;
91  traversal.add(postStepPos, depEnergy, length);
92 
93  //check if we are leaving the volume
94  bool isLeaving = (postStep.GetStepStatus() == fGeomBoundary);
95  if (isLeaving) traversal.hasLeft();
96 
97  //if we are leaving or the track is stopped, finish it
98  if (isLeaving || track.GetTrackStatus() >= fStopAndKill) {
99  bool contained = traversal.isContained();
100  bool saved = finishTrack();
101  //we mark all particles as important which created a hit and entered or left the volume
102  if (saved && !contained) {
103  Simulation::TrackInfo::getInfo(track).setIgnore(false);
104  }
105  return saved;
106  }
107  // Track not finished, return false right now
108  return false;
109  }
110 
111  bool SensitiveDetector::finishTrack()
112  {
113  SensorTraversal& traversal = m_tracks.top();
114  // ignore everything below 1eV
115  bool save = traversal.getDepEnergy() > Unit::eV;
116  if (save) {
117  auto momEntry = vecToFloat(traversal.getEntryMomentum());
118  auto posEntry = vecToFloat(traversal.getEntryPosition());
119  auto localposEntry = vecToFloat(traversal.getLocalEntryPosition());
120  auto posExit = vecToFloat(traversal.getExitPosition());
121  int hitIndex = m_simHits.getEntries();
122  m_simHits.appendNew(
123  traversal.getTrackID(),
124  traversal.getLadderID(), traversal.getSensorID(),
125  traversal.getPDGCode(), traversal.getEntryTime(),
126  traversal.getDepEnergy(),
127  traversal.getLength(), posEntry.data(),
128  localposEntry.data(),
129  posExit.data(), momEntry.data()
130  );
131  m_relMCSimHit.add(traversal.getTrackID(), hitIndex, traversal.getDepEnergy());
132  }
133  //No we just need to take care of the stack of traversals
134  if (m_tracks.size() == 1) {
135  //reuse traversal to keep memory if this is the first one
136  traversal.reset();
137  } else {
138  //this should only happen if the parent track got suspended. As this
139  //rarely happens in PXD we do not care for re-usability here
140  m_tracks.pop();
141  }
142  return save;
143  }
144  } //claw namespace
146 } //Belle2 namespace
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::Unit::ns
static const double ns
Standard of [time].
Definition: Unit.h:58
Belle2::fangs::SensorTraversal::hasLeft
void hasLeft()
indicate that the track left the current volume
Definition: SensorTraversal.h:82
Belle2::RelationArray::c_negativeWeight
@ c_negativeWeight
Flip the sign of the weight to become negative if the original element got re-attributed.
Definition: RelationArray.h:89
Belle2::Simulation::UserInfo::getInfo
static Payload getInfo(Carrier &obj)
Static function to just return UserInformation attached to the obj of type Carrier.
Definition: UserInfo.h:110
Belle2::fangs::SensorTraversal::getTrackID
int getTrackID() const
get Geant4 trackID
Definition: SensorTraversal.h:50
Belle2::fangs::SensorTraversal::getDepEnergy
double getDepEnergy() const
get total energy deposition
Definition: SensorTraversal.h:58
Belle2::Simulation::SensitiveDetectorBase::registerMCParticleRelation
static void registerMCParticleRelation(const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
Register an relation involving MCParticles.
Definition: SensitiveDetectorBase.cc:24
Belle2::fangs::SensorTraversal::isContained
bool isContained() const
return whether the track was contained in the volume so far
Definition: SensorTraversal.h:75
Belle2::fangs::SensorTraversal::hasEntered
void hasEntered()
indicate that the track originated outisde the current volume
Definition: SensorTraversal.h:80
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::Unit::eV
static const double eV
[electronvolt]
Definition: Unit.h:122
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::fangs::SensorTraversal::setInitial
void setInitial(int trackID, int ladderID, int sensorID, int pdgCode, bool primary, const G4ThreeVector &position, const G4ThreeVector &localposition, const G4ThreeVector &momentum, double time)
set initial values for a new track
Definition: SensorTraversal.h:85
Belle2::fangs::SensorTraversal
Class to keep track of the traversal of the sensitive volume for one track.
Definition: SensorTraversal.h:35
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::dosi::SensitiveDetector::step
bool step(G4Step *step, G4TouchableHistory *) override
Step processing method.
Definition: SensitiveDetector.cc:80
Belle2::Unit::GeV
static const double GeV
Standard of [energy, momentum, mass].
Definition: Unit.h:61
Belle2::fangs::SensitiveDetector::SensitiveDetector
SensitiveDetector()
Constructor.
Definition: SensitiveDetector.cc:40
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::fangs::SensorTraversal::add
void add(const G4ThreeVector &position, double depEnergy, double length)
add a new step
Definition: SensorTraversal.h:42