11 #include <simulation/kernel/UserInfo.h>
12 #include <vxd/simulation/SensitiveDetectorBase.h>
13 #include <vxd/dataobjects/VXDElectronDeposit.h>
14 #include <framework/gearbox/Unit.h>
16 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
17 #include <vxd/simulation/SensitiveDetectorDebugHelper.h>
32 const G4Track& track = *
step->GetTrack();
34 const int pdgCode = track.GetDefinition()->GetPDGEncoding();
36 const bool isNeutral = track.GetDefinition()->GetPDGCharge() == 0;
37 const bool isAllowedNeutral = (pdgCode == 22) || (
m_seeNeutrons && (abs(pdgCode) == 2112)) || (abs(pdgCode) == 99666);
39 if (isNeutral && !isAllowedNeutral)
return false;
42 const int trackID = track.GetTrackID();
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;
62 traversal.
setInitial(trackID, pdgCode, isPrimary);
64 if (preStep.GetStepStatus() == fGeomBoundary) traversal.
hasEntered();
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);
71 traversal.
add(postStepPos, postStepMom, electrons,
75 bool isLeaving = (postStep.GetStepStatus() == fGeomBoundary);
77 if (isLeaving) traversal.
hasLeft();
80 if (isLeaving || track.GetTrackStatus() >= fStopAndKill) {
96 SensorTraversal& traversal =
m_tracks.top();
97 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
103 int trueHitIndex = -1;
107 std::vector<std::pair<unsigned int, float>> simhits =
createSimHits();
110 #ifdef VXD_SENSITIVEDETECTOR_DEBUG
111 debug.finishTraversal();
127 SensorTraversal& traversal =
m_tracks.top();
129 std::vector<std::pair<unsigned int, float>> simhits;
137 static std::stack<SensorTraversal::range> stack;
140 stack.push(make_pair(traversal.begin(), traversal.end() - 1));
142 SensorTraversal::iterator firstPoint, finalPoint, splitPoint;
145 while (!stack.empty()) {
147 std::tie(firstPoint, finalPoint) = stack.top();
151 const G4ThreeVector n = (finalPoint->position - firstPoint->position).unit();
153 double maxDistance(0);
154 for (
auto nextPoint = firstPoint + 1; nextPoint != finalPoint; ++nextPoint) {
158 const G4ThreeVector pa = nextPoint->position - firstPoint->position;
159 const double dist = (pa - (pa * n) * n).mag();
161 if (dist > maxDistance) {
162 splitPoint = nextPoint;
170 stack.push(make_pair(splitPoint, finalPoint));
171 stack.push(make_pair(firstPoint, splitPoint));
176 int simHitIndex =
saveSimHit(traversal, std::make_pair(firstPoint, finalPoint));
177 simhits.push_back(std::make_pair(simHitIndex, finalPoint->electrons - firstPoint->electrons));
190 std::vector<unsigned int> electronProfile;
193 SensorTraversal::iterator firstPoint, finalPoint;
195 std::tie(firstPoint, finalPoint) = points;
198 const double electronsOffset = (firstPoint->electrons);
200 const double length = finalPoint->length - firstPoint->length;
202 const double lengthOffset = firstPoint->length;
209 static std::stack <SensorTraversal::range> stack;
215 while (!stack.empty()) {
217 std::tie(firstPoint, finalPoint) = stack.top();
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;
229 const double lengthScale = 1. /
Unit::um * 80;
232 const double slope = segmentElectrons / segmentLength / lengthScale;
238 const double distanceConstant = std::sqrt(slope * slope + 1);
241 double maxDistance(0);
243 SensorTraversal::iterator splitPoint;
246 for (
auto nextPoint = firstPoint + 1; nextPoint != finalPoint; ++nextPoint) {
248 const double x = (nextPoint->length - startLength) * lengthScale;
249 const double dist = fabs(x * slope - nextPoint->electrons + startElectrons) / distanceConstant;
251 if (dist > maxDistance) {
252 splitPoint = nextPoint;
261 stack.push(make_pair(splitPoint, finalPoint));
262 stack.push(make_pair(firstPoint, splitPoint));
266 const double fraction = (finalPoint->length - lengthOffset) / length;
267 const double electrons = (finalPoint->electrons - electronsOffset);
270 return electronProfile;
277 const double midLength = traversal.
getLength() * 0.5;
278 auto after = traversal.begin();
279 while (after->length < midLength) ++after;
283 auto before = after - 1;
285 const double fl = (after->length - midLength) / (after->length - before->length);
286 const double fr = (1 - fl);
288 const double midTime = fl * before->time + fr * after->time;
289 const double midElectrons = fl * before->electrons + fr * after->electrons;
292 const G4ThreeVector& p0 = before->position;
293 const G4ThreeVector& p3 = after->position;
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;
300 const G4ThreeVector midPos = (
302 + 3 * fl * fl * fr * p1
303 + 3 * fl * fr * fr * p2
307 const G4ThreeVector midMom = 1.0 / momentumScale * (
309 + 2 * fl * fr * (p2 - p1)
310 + fr * fr * (p3 - p2)
313 return StepInformation(midPos, midMom, midElectrons, midTime, midLength);