Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
//==========================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author : M.Frank
//
//==========================================================================
// Framework include files
#include "DD4hep/DD4hepUnits.h"
#include "DDG4/Geant4SensDetAction.inl"
#include "DDG4/Geant4SteppingAction.h"
#include "DDG4/Geant4TrackingAction.h"
#include "DDG4/Geant4EventAction.h"
#include "G4Event.hh"
#include "G4VSolid.hh"
#include <map>
#include <limits>
#include <sstream>
using namespace std;
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
using namespace Geometry;
/// Geant4 sensitive detector combining all deposits of one G4Track within one sensitive element.
/**
* Geant4SensitiveAction<TrackerWeighted>
*
*
* \author M.Frank
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
/** \addtogroup Geant4SDActionPlugin
*
* @{
* \package Geant4TrackerWeightedAction
*
* \brief Sensitive detector meant for tracking detectors with multiple ways to combine steps
\param integer HitPositionCombination
-# Use energy weights to define the position of the energy deposit
-# Set the hit position between the step pre and post point
-# Set the hit position to the position of the step pre point
-# Set the hit position to the position of the step post point
\param bool CollectSingleDeposits
- If true each step is written out
*
* @}
*/
struct TrackerWeighted {
typedef Geant4HitCollection HitCollection;
/// Enumeration to define the calculation of the position of the energy deposit
enum {
POSITION_WEIGHTED = 1, // Use energy weights to define the position of the energy deposit
POSITION_MIDDLE = 2, // Set the hit position between the step pre and post point
POSITION_PREPOINT = 3, // Set the hit position to the position of the step pre point
POSITION_POSTPOINT = 4 // Set the hit position to the position of the step post point
};
Geant4Tracker::Hit pre, post;
Position mean_pos;
Geant4Sensitive* sensitive = 0;
G4VSensitiveDetector* thisSD = 0;
double distance_to_inside = 0.0;
double distance_to_outside = 0.0;
double mean_time = 0.0;
double step_length = 0.0;
double e_cut = 0.0;
Andre Sailer
committed
int current = -1;
int parent = 0;
int combined = 0;
int hit_position_type = POSITION_MIDDLE;
int hit_flag = 0;
int g4ID = 0;
EInside last_inside = kOutside;
long long int cell = 0;
bool single_deposit_mode = false;
TrackerWeighted() = default;
/// Clear collected information and restart for new hit
TrackerWeighted& clear() {
mean_pos.SetXYZ(0,0,0);
distance_to_inside = 0;
distance_to_outside = 0;
mean_time = 0;
Markus Frank
committed
step_length = 0;
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
post.clear();
pre.clear();
current = -1;
parent = -1;
combined = 0;
cell = 0;
hit_flag = 0;
g4ID = 0;
last_inside = kOutside;
return *this;
}
/// Start a new hit
TrackerWeighted& start(const G4Step* step, const G4StepPoint* point) {
clear();
pre.storePoint(step,point);
pre.truth.deposit = 0.0;
post.truth.deposit = 0.0;
current = pre.truth.trackID;
sensitive->mark(step->GetTrack());
post = pre;
parent = step->GetTrack()->GetParentID();
g4ID = step->GetTrack()->GetTrackID();
return *this;
}
/// Update energy and track information during hit info accumulation
TrackerWeighted& update(const G4Step* step) {
post.storePoint(step,step->GetPostStepPoint());
Markus Frank
committed
Position mean = (post.position+pre.position)*0.5;
double mean_tm = (post.truth.time+pre.truth.time)*0.5;
pre.truth.deposit += post.truth.deposit;
Markus Frank
committed
mean_pos.SetX(mean_pos.x()+mean.x()*post.truth.deposit);
mean_pos.SetY(mean_pos.y()+mean.y()*post.truth.deposit);
mean_pos.SetZ(mean_pos.z()+mean.z()*post.truth.deposit);
mean_time += mean_tm*post.truth.deposit;
step_length += step->GetStepLength();
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
if ( 0 == cell ) {
cell = sensitive->cellID(step);
if ( 0 == cell ) {
cell = sensitive->volumeID(step) ;
sensitive->except("+++ Invalid CELL ID for hit!");
}
}
++combined;
return *this;
}
/// Helper function to decide if the hit has to be extracted and saved in the collection
bool mustSaveTrack(const G4Track* tr) const {
return current > 0 && current != tr->GetTrackID();
}
/// Extract hit information and add the created hit to the collection
void extractHit(EInside ended) {
Geant4HitCollection* collection = sensitive->collection(0);
extractHit(collection, ended);
}
TrackerWeighted& calc_dist_out(const G4VSolid* solid) {
Position v(pre.momentum.unit()), &p=post.position;
double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()),
G4ThreeVector(v.X(),v.Y(),v.Z()));
distance_to_outside = dist;
return *this;
}
TrackerWeighted& calc_dist_in(const G4VSolid* solid) {
Position v(pre.momentum.unit()), &p=pre.position;
double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()),
G4ThreeVector(v.X(),v.Y(),v.Z()));
distance_to_inside = dist;
return *this;
}
void extractHit(Geant4HitCollection* collection, EInside ended) {
double deposit = pre.truth.deposit;
Position pos;
Momentum mom = 0.5 * (pre.momentum + post.momentum);
double time = deposit != 0 ? mean_time / deposit : mean_time;
char dist_in[64], dist_out[64];
switch(hit_position_type) {
case POSITION_WEIGHTED:
pos = deposit != 0 ? mean_pos / deposit : mean_pos;
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
break;
case POSITION_PREPOINT:
pos = pre.position;
break;
case POSITION_POSTPOINT:
pos = post.position;
break;
case POSITION_MIDDLE:
default:
pos = (post.position + pre.position) / 2.0;
break;
}
if ( ended == kSurface || distance_to_outside < numeric_limits<float>::epsilon() )
hit_flag |= Geant4Tracker::Hit::HIT_ENDED_SURFACE;
else if ( ended == kInside )
hit_flag |= Geant4Tracker::Hit::HIT_ENDED_INSIDE;
else if ( ended == kOutside )
hit_flag |= Geant4Tracker::Hit::HIT_ENDED_OUTSIDE;
Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID,
pre.truth.pdgID,
deposit,time);
hit->flag = hit_flag;
hit->position = pos;
hit->momentum = mom;
Markus Frank
committed
hit->length = step_length;
hit->cellID = cell;
hit->g4ID = g4ID;
dist_in[0] = dist_out[0] = 0;
if ( !(hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) )
::snprintf(dist_in,sizeof(dist_in)," [%.2e um]",distance_to_inside/CLHEP::um);
if ( !(hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE) )
::snprintf(dist_out,sizeof(dist_out)," [%.2e um]",distance_to_outside/CLHEP::um);
sensitive->print("+++ G4Track:%5d CREATE hit[%03d]:%3d deps E:"
" %.2e keV Pos:%7.2f %7.2f %7.2f [mm] Start:%s%s%s%s End:%s%s%s%s",
pre.truth.trackID,int(collection->GetSize()),
combined,pre.truth.deposit/CLHEP::keV,
pos.X()/CLHEP::mm,pos.Y()/CLHEP::mm,pos.Z()/CLHEP::mm,
((hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) ? "SURFACE" : ""),
((hit_flag&Geant4Tracker::Hit::HIT_STARTED_OUTSIDE) ? "OUTSIDE" : ""),
((hit_flag&Geant4Tracker::Hit::HIT_STARTED_INSIDE) ? "INSIDE " : ""),
((hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE) ? "SURFACE" : ""),
((hit_flag&Geant4Tracker::Hit::HIT_ENDED_OUTSIDE) ? "OUTSIDE" : ""),
((hit_flag&Geant4Tracker::Hit::HIT_ENDED_INSIDE) ? "INSIDE " : ""),
dist_out);
collection->add(hit);
}
clear();
}
/// Method for generating hit(s) using the information of G4Step object.
G4bool process(const G4Step* step, G4TouchableHistory* ) {
Geant4StepHandler h(step);
// std::cout << " process called - pre pos: " << h.prePos() << " post pos " << h.postPos()
// << " edep: " << h.deposit() << std::endl ;
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
G4VSolid* preSolid = h.solid(h.pre);
G4VSolid* postSolid = h.solid(h.post);
G4ThreeVector local_pre = h.globalToLocalG4(h.prePosG4());
G4ThreeVector local_post = h.globalToLocalG4(h.postPosG4());
EInside pre_inside = preSolid->Inside(local_pre);
EInside post_inside = postSolid->Inside(local_post);
const void* postPV = h.postVolume();
const void* prePV = h.preVolume();
const void* postSD = h.postSD();
const void* preSD = h.preSD();
G4VSolid* solid = (preSD == thisSD) ? preSolid : postSolid;
// 1) Track killed inside SD: trace incomplete. This deposition must be added as well.
if ( current == h.trkID() && !h.trkAlive() ) {
hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
update(step).calc_dist_out(solid).extractHit(post_inside);
return true;
}
// 2) Track leaving SD volume. Sensitive detector changed. Store hit.
else if ( current == h.trkID() && postSD != thisSD ) {
update(step).calc_dist_out(solid).extractHit(kOutside);
return true;
}
// 3) Track leaving SD volume. Store hit.
else if ( current == h.trkID() && postSD == thisSD && post_inside == kSurface ) {
update(step).calc_dist_out(solid).extractHit(kSurface);
return true;
}
// 4) Track leaving SD volume. Store hit.
else if ( current == h.trkID() && postSD == thisSD && post_inside == kOutside ) {
update(step).calc_dist_out(solid).extractHit(post_inside);
return true;
}
// 5) Normal update: either intermediate deposition or track is going to be killed.
else if ( current == h.trkID() && postSD == thisSD && post_inside == kInside ) {
last_inside = post_inside;
update(step).calc_dist_out(solid);
return true;
}
// 6) Track from secondary created in SD volume. Store hit from previous.
// --> New hit will be created, to whom also this deposition belongs
else if ( current != h.trkID() && current >= 0 ) {
extractHit(last_inside);
}
// If the hit got extracted, a new one must be set up
if ( current < 0 ) {
EInside inside = pre_inside;
// Track entering SD volume
if ( preSD != thisSD ) {
start(step, h.post);
inside = post_inside;
sensitive->print("++++++++++ Using POST step volume to start hit -- dubious ?");
}
else {
start(step, h.pre);
}
calc_dist_in(solid);
if ( inside == kSurface )
hit_flag |= Geant4Tracker::Hit::HIT_STARTED_SURFACE;
else if ( inside == kInside )
hit_flag |= Geant4Tracker::Hit::HIT_STARTED_INSIDE;
else if ( inside == kOutside )
hit_flag |= Geant4Tracker::Hit::HIT_STARTED_OUTSIDE;
// New (secondary) track created by some process starting inside the volume
if ( inside == kInside ) {
hit_flag |= Geant4Tracker::Hit::HIT_SECONDARY_TRACK;
}
}
// Update Data
last_inside = post_inside;
update(step);
calc_dist_out(solid);
// Track killed inside SD: trace incomplete. This deposition must be added as well.
if ( !h.trkAlive() ) {
hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
extractHit(post_inside);
}
// Avoid danglich hits if the track leaves the sensitive volume
else if ( post_inside == kSurface ) {
extractHit(post_inside);
}
// Avoid danglich hits if the track leaves the sensitive volume
else if ( thisSD == preSD && (preSD != postSD || prePV != postPV) ) {
extractHit(post_inside);
}
// This should simply not happen!
else if ( thisSD == postSD && (preSD != postSD || prePV != postPV) ) {
sensitive->error("+++++ WRONG!!! Extract. How did we get here?");
extractHit(post_inside);
}
// In single hit mode we MUST write the hit after update
else if ( single_deposit_mode ) {
extractHit(post_inside);
}
return true;
}
/// Post-event action callback
Markus Frank
committed
void endEvent() {
// We need to add the possibly last added hit to the collection here.
// otherwise the last hit would be assigned to the next event and the
// MC truth would be screwed.
//
if ( current > 0 ) {
Geant4HitCollection* coll = sensitive->collection(0);
sensitive->print("++++++++++ Found dangling hit: Is the hit extraction logic correct?");
extractHit(coll,last_inside);
}
}
/// Pre event action callback
Markus Frank
committed
void startEvent() {
thisSD = dynamic_cast<G4VSensitiveDetector*>(&sensitive->detector());
}
};
/// Initialization overload for specialization
template <> void Geant4SensitiveAction<TrackerWeighted>::initialize() {
declareProperty("HitPositionCombination", m_userData.hit_position_type);
declareProperty("CollectSingleDeposits", m_userData.single_deposit_mode);
m_userData.e_cut = m_sensitive.energyCutoff();
m_userData.sensitive = this;
}
Markus Frank
committed
/// G4VSensitiveDetector interface: Method invoked at the begining of each event.
template <> void Geant4SensitiveAction<TrackerWeighted>::begin(G4HCofThisEvent* /* hce */) {
Markus Frank
committed
m_userData.startEvent();
}
/// G4VSensitiveDetector interface: Method invoked at the end of each event.
template <> void Geant4SensitiveAction<TrackerWeighted>::end(G4HCofThisEvent* /* hce */) {
Markus Frank
committed
m_userData.endEvent();
}
/// Define collections created by this sensitivie action object
template <> void Geant4SensitiveAction<TrackerWeighted>::defineCollections() {
m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>();
}
/// Method for generating hit(s) using the information of G4Step object.
template <> void Geant4SensitiveAction<TrackerWeighted>::clear(G4HCofThisEvent* /* hce */) {
m_userData.clear();
}
/// Method for generating hit(s) using the information of G4Step object.
template <> G4bool
Geant4SensitiveAction<TrackerWeighted>::process(G4Step* step, G4TouchableHistory* history) {
return m_userData.process(step, history);
}
typedef Geant4SensitiveAction<TrackerWeighted> Geant4TrackerWeightedAction;
}
}
using namespace DD4hep::Simulation;
#include "DDG4/Factories.h"
DECLARE_GEANT4SENSITIVE(Geant4TrackerWeightedAction)