Belle II Software  release-05-01-25
SensitiveDetectorBase.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2014 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <simulation/kernel/UserInfo.h>
12 #include <vxd/simulation/SensitiveDetectorBase.h>
13 #include <vxd/dataobjects/VXDElectronDeposit.h>
14 #include <framework/gearbox/Unit.h>
15 
16 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
17 #include <vxd/simulation/SensitiveDetectorDebugHelper.h>
18 #endif
19 
20 #include <G4Step.hh>
21 
22 namespace Belle2 {
27  namespace VXD {
28 
29  bool SensitiveDetectorBase::step(G4Step* step, G4TouchableHistory*)
30  {
31  // Get track
32  const G4Track& track = *step->GetTrack();
33  // Get particle PDG code
34  const int pdgCode = track.GetDefinition()->GetPDGEncoding();
35  // Get particle charge (only keep charged tracks, photons and magnetic monopoles)
36  const bool isNeutral = track.GetDefinition()->GetPDGCharge() == 0;
37  const bool isAllowedNeutral = (pdgCode == 22) || (m_seeNeutrons && (abs(pdgCode) == 2112)) || (abs(pdgCode) == 99666);
38  //Not interested in neutral particles except photons, magnetic monopoles and maybe neutrons
39  if (isNeutral && !isAllowedNeutral) return false;
40 
41  // Get track ID
42  const int trackID = track.GetTrackID();
43  //Get deposited energy
44  const double electrons = step->GetTotalEnergyDeposit() * Unit::MeV / Const::ehEnergy;
45  // Get step information
46  const G4StepPoint& postStep = *step->GetPostStepPoint();
47  const G4StepPoint& preStep = *step->GetPreStepPoint();
48  const G4AffineTransform& topTransform = preStep.GetTouchableHandle()->GetHistory()->GetTopTransform();
49  const G4ThreeVector postStepPos = topTransform.TransformPoint(postStep.GetPosition()) * Unit::mm;
50  const G4ThreeVector postStepMom = topTransform.TransformAxis(postStep.GetMomentum()) * Unit::MeV;
51 
52  //Get the current track info
53  //if none present or not matching current trackID create one first
54  //if trackID is 0 we can reuse the existing top
55  if (m_tracks.empty() || (m_tracks.top().getTrackID() > 0 && m_tracks.top().getTrackID() != trackID)) {
56  m_tracks.push(SensorTraversal());
57  }
58  SensorTraversal& traversal = m_tracks.top();
59  if (traversal.getTrackID() == 0) {
60  bool isPrimary = Simulation::TrackInfo::getInfo(track).hasStatus(MCParticle::c_PrimaryParticle);
61  //If new track, remember values
62  traversal.setInitial(trackID, pdgCode, isPrimary);
63  //Remember if the track came from the outside
64  if (preStep.GetStepStatus() == fGeomBoundary) traversal.hasEntered();
65  //Add start position
66  const G4ThreeVector preStepPos = topTransform.TransformPoint(preStep.GetPosition()) * Unit::mm;
67  const G4ThreeVector preStepMom = topTransform.TransformAxis(preStep.GetMomentum()) * Unit::MeV;
68  traversal.add(preStepPos, preStepMom, 0, preStep.GetGlobalTime() * Unit::ns, 0);
69  }
70  //Add new track
71  traversal.add(postStepPos, postStepMom, electrons,
72  postStep.GetGlobalTime() * Unit::ns,
73  step->GetStepLength() * Unit::mm);
74  //check if we are leaving the volume
75  bool isLeaving = (postStep.GetStepStatus() == fGeomBoundary);
76  //remember that in the track info
77  if (isLeaving) traversal.hasLeft();
78  //If this step is to the boundary or track gets killed, save hits and
79  //return whether or not we saved something
80  if (isLeaving || track.GetTrackStatus() >= fStopAndKill) {
81  bool saved = finishTrack();
82  //If we saved at least one simhit and the track is not contained inside
83  //the sensor volume: keep MCParticle
84  if (saved && !traversal.isContained()) {
85  Simulation::TrackInfo::getInfo(track).setIgnore(false);
86  }
87  return saved;
88  }
89  //Track not finished so we do not save anything, let's return false for
90  //now
91  return false;
92  }
93 
95  {
96  SensorTraversal& traversal = m_tracks.top();
97 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
98  SensitiveDetectorDebugHelper& debug = SensitiveDetectorDebugHelper::getInstance();
99  debug.startTraversal(getSensorID(), traversal);
100 #endif
101  bool save = traversal.getElectrons() >= m_minimumElectrons || m_seeNeutrons;
102  if (save) {
103  int trueHitIndex = -1;
104  if (!m_onlyPrimaryTrueHits || traversal.isPrimary()) {
105  trueHitIndex = saveTrueHit(traversal);
106  }
107  std::vector<std::pair<unsigned int, float>> simhits = createSimHits();
108  saveRelations(traversal, trueHitIndex, simhits);
109  }
110 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
111  debug.finishTraversal();
112 #endif
113  //No we just need to take care of the stack of traversals
114  if (m_tracks.size() == 1) {
115  //reuse traversal to keep memory if this is the first one
116  traversal.reset();
117  } else {
118  //this should only happen if the parent track got suspended. As this
119  //rarely happens in PXD we do not care for re-usability here
120  m_tracks.pop();
121  }
122  return save;
123  }
124 
125  std::vector<std::pair<unsigned int, float>> SensitiveDetectorBase::createSimHits()
126  {
127  SensorTraversal& traversal = m_tracks.top();
128  //List of created simhit indices to be able to create relations
129  std::vector<std::pair<unsigned int, float>> simhits;
130 
131  //We need to check how close the steps would be to a linear approximation
132  //of the traversal and split the track if they are too far away. To do
133  //this we check the whole track and split it at the point with the
134  //largest distance recursively until the largest distance is inside the
135  //tolerance. For that we need a stack of segments and check each segment
136  //in turn until they fulfill the criteria.
137  static std::stack<SensorTraversal::range> stack;
138  //Lets push the full segment, first step point to last step point, on the
139  //stack. We use inclusive range, so the second iterator is still included
140  stack.push(make_pair(traversal.begin(), traversal.end() - 1));
141  //Iterators needed for checking
142  SensorTraversal::iterator firstPoint, finalPoint, splitPoint;
143 
144  //Check all segments ...
145  while (!stack.empty()) {
146  //Get first and last point
147  std::tie(firstPoint, finalPoint) = stack.top();
148  //Remove segment from stack
149  stack.pop();
150  //Direction of the segment
151  const G4ThreeVector n = (finalPoint->position - firstPoint->position).unit();
152  //find largest distance to segment by looping over all intermediate points
153  double maxDistance(0);
154  for (auto nextPoint = firstPoint + 1; nextPoint != finalPoint; ++nextPoint) {
155  //Calculate distances between point p and line segment,
156  //x = a + t*n, where a is a point on the line and n is the unit vector
157  //pointing in line direction. distance = ||(p-a) - ((p-a)*n)*n||
158  const G4ThreeVector pa = nextPoint->position - firstPoint->position;
159  const double dist = (pa - (pa * n) * n).mag();
160  //Update splitpoint if distance is larger then previously known
161  if (dist > maxDistance) {
162  splitPoint = nextPoint;
163  maxDistance = dist;
164  }
165  }
166  //If we are above the tolerance we split the track
167  if (maxDistance > m_distanceTolerance) {
168  //If we split in this order, all created simhits will be in correct
169  //order. That is not a requirement but a nice side effect.
170  stack.push(make_pair(splitPoint, finalPoint));
171  stack.push(make_pair(firstPoint, splitPoint));
172  continue;
173  }
174  //Otherwise we create a SimHit
175  //FIXME: check for m_minimumElectrons?
176  int simHitIndex = saveSimHit(traversal, std::make_pair(firstPoint, finalPoint));
177  simhits.push_back(std::make_pair(simHitIndex, finalPoint->electrons - firstPoint->electrons));
178  }
179  return simhits;
180  }
181 
182  std::vector<unsigned int> SensitiveDetectorBase::simplifyEnergyDeposit(const SensorTraversal::range& points)
183  {
184  //At the end we want a list of points indication how many electrons where
185  //deposited so far along the track, e.g. [(0.5,100), (1.0, 2000)] would
186  //indicate that after half the step length 100 electrons were deposited
187  //and at the end of the step we have 2000 electrons. To save some memory
188  //we encode this information as unsigned int using ElectronDeposit later
189  //on.
190  std::vector<unsigned int> electronProfile;
191 
192  //We need an iterator to the first and last point in a segment
193  SensorTraversal::iterator firstPoint, finalPoint;
194  //for now let's extract the full segment
195  std::tie(firstPoint, finalPoint) = points;
196  //If this is not the first SimHit for a sensor traversal we need to
197  //subtract the electrons already taken care of
198  const double electronsOffset = (firstPoint->electrons);
199  //We also need the length of the step
200  const double length = finalPoint->length - firstPoint->length;
201  //And the start length
202  const double lengthOffset = firstPoint->length;
203 
204  //We need to keep track of sub segments and where they should insert a
205  //midpoint if required. So we need a tuple of three iterators: insert
206  //position, first point and last point of the segment. We store those
207  //in a stack which we declare static to avoid some relocations if
208  //possible.
209  static std::stack <SensorTraversal::range> stack;
210  //And we push the full segment on the stack
211  stack.push(points);
212 
213  //Now we check all segments we encounter until none exceed the maximum
214  //tolerance anymore
215  while (!stack.empty()) {
216  //Get next segment to be checked
217  std::tie(firstPoint, finalPoint) = stack.top();
218  //And remove it from the stack
219  stack.pop();
220 
221  //Some variables we need
222  const double startElectrons = firstPoint->electrons;
223  const double startLength = firstPoint->length;
224  const double segmentLength = finalPoint->length - startLength;
225  const double segmentElectrons = finalPoint->electrons - startElectrons;
226  //We want to give the tolerance in electrons so we need to convert the
227  //step length into electrons. We do this by using the average number of
228  //electrons created per micrometer for a minimum ionizing particle.
229  const double lengthScale = 1. / Unit::um * 80;
230  //we need the slope for the linear approximation:
231  //electrons= slope*(length-startLength)*lengthScale+startElectrons;
232  const double slope = segmentElectrons / segmentLength / lengthScale;
233  //Distance between point x0,y0 and line a*x+b*y+c=0 is given as
234  //abs(a*x0+b*y0+c)/sqrt(a^2+b^2). In our case:
235  //x0=(length-startLength)*lengthScale, y0=electrons
236  //a=slope, b=-1, c=startElectrons.
237  //Nominator is independent of the actual point so we calculate it now
238  const double distanceConstant = std::sqrt(slope * slope + 1);
239 
240  //Variable to remember maximum distance from linear approximation
241  double maxDistance(0);
242  //and also the point with the largest distance
243  SensorTraversal::iterator splitPoint;
244 
245  //No go through all intermediate points
246  for (auto nextPoint = firstPoint + 1; nextPoint != finalPoint; ++nextPoint) {
247  //and check their distance from the linear approximation
248  const double x = (nextPoint->length - startLength) * lengthScale;
249  const double dist = fabs(x * slope - nextPoint->electrons + startElectrons) / distanceConstant;
250  //and remember if it is the largest distance
251  if (dist > maxDistance) {
252  splitPoint = nextPoint;
253  maxDistance = dist;
254  }
255  }
256  //if the largest distance is above the tolerance we split the segment
257  if (maxDistance > m_electronTolerance) {
258  //And add the two sub segments to the stack in the correcto order to
259  //ensure sorted processing: We always add the segments at the front
260  //last so they will be processed first
261  stack.push(make_pair(splitPoint, finalPoint));
262  stack.push(make_pair(firstPoint, splitPoint));
263  continue;
264  }
265  //otherwise we add the endpoint of the current segment to the list of points.
266  const double fraction = (finalPoint->length - lengthOffset) / length;
267  const double electrons = (finalPoint->electrons - electronsOffset);
268  electronProfile.push_back(VXDElectronDeposit(fraction, electrons));
269  }
270  return electronProfile;
271  }
272 
274  {
275  //We want the middle of the track so we need to find the two steps
276  //which enclose this
277  const double midLength = traversal.getLength() * 0.5;
278  auto after = traversal.begin();
279  while (after->length < midLength) ++after;
280  //Now we have an iterator after the half length. We know that the first
281  //step contains length 0 so we can savely subtract one to get the one
282  //before
283  auto before = after - 1;
284  //the midpoint is in between so calculate the fractions from both sides
285  const double fl = (after->length - midLength) / (after->length - before->length);
286  const double fr = (1 - fl);
287  //we calculate the time and electrons using weighted average
288  const double midTime = fl * before->time + fr * after->time;
289  const double midElectrons = fl * before->electrons + fr * after->electrons;
290  //now we use 3rd order bezier curve to approximate mid position and momentum.
291  //The two base points are easy ...
292  const G4ThreeVector& p0 = before->position;
293  const G4ThreeVector& p3 = after->position;
294  //And the two control points are the base points +- the appropriately
295  //scaled momenta
296  const double momentumScale = (p3 - p0).mag() / before->momentum.mag() / 3;
297  const G4ThreeVector p1 = p0 + momentumScale * before->momentum;
298  const G4ThreeVector p2 = p3 - momentumScale * after->momentum;
299  // The curve is B(t) = (1-t)^3*P0 + 3*(1-t)^2*t*P1 + 3*(1-t)*t^2*P2 + t^3*P3
300  const G4ThreeVector midPos = (
301  fl * fl * fl * p0
302  + 3 * fl * fl * fr * p1
303  + 3 * fl * fr * fr * p2
304  + fr * fr * fr * p3
305  );
306  // The derivative is dB(t)/dt = 3*[(1-t)^2*(P1-P0)+2*(1-t)*t*(P2-P1)+t^2*(P3-P2)]
307  const G4ThreeVector midMom = 1.0 / momentumScale * (
308  fl * fl * (p1 - p0)
309  + 2 * fl * fr * (p2 - p1)
310  + fr * fr * (p3 - p2)
311  );
312  //Okay, we now have everything
313  return StepInformation(midPos, midMom, midElectrons, midTime, midLength);
314  }
315 
316  } //VXD namespace
318 } //Belle2 namespace
Belle2::Unit::ns
static const double ns
Standard of [time].
Definition: Unit.h:58
Belle2::VXD::SensitiveDetectorBase::m_minimumElectrons
float m_minimumElectrons
minimum number of electrons a track must deposit for SimHit/TrueHits to be created
Definition: SensitiveDetectorBase.h:184
Belle2::SensorTraversal
Class to keep track of the traversal of the sensitive volume for one track.
Definition: SensorTraversal.h:54
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::VXD::SensitiveDetectorBase::createSimHits
std::vector< std::pair< unsigned int, float > > createSimHits()
Determine which SimHits to create.
Definition: SensitiveDetectorBase.cc:133
Belle2::VXD::SensitiveDetectorBase::finishTrack
bool finishTrack()
Process a track once all steps are known.
Definition: SensitiveDetectorBase.cc:102
Belle2::VXD::SensitiveDetectorBase::getSensorID
VxdID getSensorID() const
Return the VxdID belonging to this sensitive detector.
Definition: SensitiveDetectorBase.h:90
Belle2::SensorTraversal::getLength
double getLength() const
get flight length so far
Definition: SensorTraversal.h:81
Belle2::SensorTraversal::add
void add(const G4ThreeVector &position, const G4ThreeVector &momentum, double electrons, double time, double length)
add a new step
Definition: SensorTraversal.h:67
Belle2::SensorTraversal::hasLeft
void hasLeft()
indicate that the track left the current volume
Definition: SensorTraversal.h:90
Belle2::VXD::SensitiveDetectorBase::findMidPoint
StepInformation findMidPoint(const SensorTraversal &traversal)
Find the mid-point of the track traversal.
Definition: SensitiveDetectorBase.cc:281
Belle2::VXDElectronDeposit
Packed class to represent energy deposit along a path in electrons.
Definition: VXDElectronDeposit.h:34
Belle2::SensorTraversal::setInitial
void setInitial(int trackID, int pdgCode, bool primary)
set initial values for a new track
Definition: SensorTraversal.h:93
Belle2::VXD::SensitiveDetectorBase::m_seeNeutrons
bool m_seeNeutrons
also create SimHit/TrueHit objects for neutrons (or charged particles which deposit less than m_minim...
Definition: SensitiveDetectorBase.h:187
Belle2::StepInformation
Simple struct to keep information about steps in the sensitive detector.
Definition: SensorTraversal.h:27
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::SensorTraversal::getTrackID
int getTrackID() const
get Geant4 trackID
Definition: SensorTraversal.h:75
Belle2::VXD::SensitiveDetectorBase::m_electronTolerance
float m_electronTolerance
maximum relative difference between electron density of two steps where they can be considered simila...
Definition: SensitiveDetectorBase.h:181
Belle2::VXD::SensitiveDetectorBase::m_onlyPrimaryTrueHits
bool m_onlyPrimaryTrueHits
only create TrueHits for primary particles
Definition: SensitiveDetectorBase.h:189
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SensorTraversal::range
std::pair< iterator, iterator > range
Iterator pair for a set of points.
Definition: SensorTraversal.h:57
Belle2::VXD::SensitiveDetectorBase::step
bool step(G4Step *step, G4TouchableHistory *) override
Process a single Geant4 Step.
Definition: SensitiveDetectorBase.cc:37
Belle2::VXD::SensitiveDetectorDebugHelper::getInstance
static SensitiveDetectorDebugHelper & getInstance()
Singleton class: get instance.
Definition: SensitiveDetectorDebugHelper.cc:30
Belle2::VXD::SensitiveDetectorBase::saveTrueHit
virtual int saveTrueHit(const SensorTraversal &traversal)=0
Save the actual TrueHit for this sensor traversal.
Belle2::VXD::SensitiveDetectorBase::simplifyEnergyDeposit
std::vector< unsigned int > simplifyEnergyDeposit(const SensorTraversal::range &points)
Simplify the energy deposition profile using Douglas-Peuker-Algorithm We normally force a Geant4 step...
Definition: SensitiveDetectorBase.cc:190
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::VXD::SensitiveDetectorBase::m_tracks
std::stack< SensorTraversal > m_tracks
stack of SensorTraversal information for all tracks not finished so far
Definition: SensitiveDetectorBase.h:175
Belle2::VXD::SensitiveDetectorBase::saveSimHit
virtual int saveSimHit(const SensorTraversal &traversal, const SensorTraversal::range &points)=0
Save a SimHit for this track including the given points.
Belle2::SensorTraversal::hasEntered
void hasEntered()
indicate that the track originated outisde the current volume
Definition: SensorTraversal.h:88
Belle2::SensorTraversal::isContained
bool isContained() const
return whether the track was contained in the volume so far
Definition: SensorTraversal.h:83
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::VXD::SensitiveDetectorBase::saveRelations
virtual void saveRelations(const SensorTraversal &traversal, int trueHitIndex, std::vector< std::pair< unsigned int, float >> simHitIndices)=0
Save the relations between MCParticle, TrueHit and SimHits.
Belle2::VXD::SensitiveDetectorBase::m_distanceTolerance
float m_distanceTolerance
maximum distance between step point and linear interpolation of sensor traversal before a new simhit ...
Definition: SensitiveDetectorBase.h:178
Belle2::Const::ehEnergy
static const double ehEnergy
Energy needed to create an electron-hole pair in Si at std.
Definition: Const.h:570
Belle2::Unit::um
static const double um
[micrometers]
Definition: Unit.h:81