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