DDG4/plugins/Geant4TrackerWeightedSD.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/DD4hepUnits.h>
0016 #include <DDG4/Geant4SensDetAction.inl>
0017 #include <DDG4/Geant4SteppingAction.h>
0018 #include <DDG4/Geant4TrackingAction.h>
0019 #include <DDG4/Geant4EventAction.h>
0020 #include <G4Event.hh>
0021 #include <G4VSolid.hh>
0022
0023 #include <map>
0024 #include <limits>
0025 #include <sstream>
0026
0027 /// Namespace for the AIDA detector description toolkit
0028 namespace dd4hep {
0029
0030   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0031   namespace sim   {
0032
0033     using namespace detail;
0034
0035     /**
0036      *  \addtogroup Geant4SDActionPlugin
0037      *
0038      *  @{
0039      *  \package Geant4TrackerWeightedAction
0040      *
0041      *  \brief Sensitive detector meant for tracking detectors with multiple ways to combine steps
0042      *
0043      *
0044      *  \param integer HitPositionCombination
0045      *   -# Use energy weights to define the position of the energy deposit
0046      *   -# Set the hit position between the step pre and post point
0047      *   -# Set the hit position to the position of the step pre point
0048      *   -# Set the hit position to the position of the step post point
0049      *
0050      *  \param bool CollectSingleDeposits
0051      *   - If true each step is written out
0052      *
0053      * @}
0054      */
0055     /// Geant4 sensitive detector combining all deposits of one G4Track within one sensitive element.
0056     /**
0057      *  Geant4SensitiveAction<TrackerWeighted>
0058      *
0059      *
0060      *  \author  M.Frank
0061      *  \version 1.0
0062      *  \ingroup DD4HEP_SIMULATION
0063      */
0064     struct TrackerWeighted {
0065
0066       /// Enumeration to define the calculation of the position of the energy deposit
0067       enum  {
0068         POSITION_WEIGHTED  = 1, // Use energy weights to define the position of the energy deposit
0069         POSITION_MIDDLE    = 2, // Set the hit position between the step pre and post point
0070         POSITION_PREPOINT  = 3, // Set the hit position to the position of the step pre point
0071         POSITION_POSTPOINT = 4  // Set the hit position to the position of the step post point
0072       };
0073
0074       Geant4Tracker::Hit    pre, post;
0075       Position              mean_pos;
0076       Geant4Sensitive*      sensitive            = 0;
0077       G4VSensitiveDetector* thisSD               = 0;
0078       G4VPhysicalVolume*    thisPV               = 0;
0079       double                distance_to_inside   = 0.0;
0080       double                distance_to_outside  = 0.0;
0081       double                mean_time            = 0.0;
0082       double                step_length          = 0.0;
0083       double                e_cut                = 0.0;
0084       int                   current              = -1;
0085       int                   parent               = 0;
0086       int                   combined             = 0;
0087       int                   hit_position_type    = POSITION_MIDDLE;
0088       int                   hit_flag             = 0;
0089       int                   g4ID                 = 0;
0090       EInside               last_inside          = kOutside;
0091       long long int         cell                 = 0;
0092       bool                  single_deposit_mode  = false;
0093
0094       /// Default constructor
0095       TrackerWeighted() = default;
0096
0097       /// Clear collected information and restart for new hit
0098       TrackerWeighted& clear()   {
0099         mean_pos.SetXYZ(0,0,0);
0100         distance_to_inside = 0;
0101         distance_to_outside = 0;
0102         mean_time = 0;
0103         step_length  = 0;
0104         thisPV = nullptr;
0105         post.clear();
0106         pre.clear();
0107         current  = -1;
0108         parent   = -1;
0109         combined =  0;
0110         cell     =  0;
0111         hit_flag =  0;
0112         g4ID     =  0;
0113         last_inside = kOutside;
0114         return *this;
0115       }
0116
0117       /// Start a new hit
0118       TrackerWeighted& start(const G4Step* step, const G4StepPoint* point)   {
0119         if( DEBUG == printLevel() ) {
0120           std::cout<<" DEBUG: Geant4TrackerWeightedSD::start(const G4Step* step, const G4StepPoint* point) ...."<<std::endl;
0121           Geant4StepHandler h(step);
0122           dumpStep( h, step);
0123         }
0124
0125         clear();
0126         pre.storePoint(step,point);
0127         pre.truth.deposit = 0.0;
0128         post.truth.deposit = 0.0;
0129         current = pre.truth.trackID;
0130         sensitive->mark(step->GetTrack());
0131         post.copyFrom(pre);
0132         parent = step->GetTrack()->GetParentID();
0133         g4ID = step->GetTrack()->GetTrackID();
0134
0135         Geant4StepHandler startVolume(step);
0136         thisPV = startVolume.preVolume();
0137
0138         return *this;
0139       }
0140
0141       /// Update energy and track information during hit info accumulation
0142       TrackerWeighted& update(const G4Step* step)   {
0143         if( DEBUG == printLevel() ) {
0144           std::cout<<" DEBUG: Geant4TrackerWeightedSD::update(const G4Step* step) ...."<<std::endl;
0145           Geant4StepHandler h(step);
0146           dumpStep( h, step);
0147         }
0148
0149         post.storePoint(step,step->GetPostStepPoint());
0150         Position mean    = (post.position+pre.position)*0.5;
0151         double   mean_tm = (post.truth.time+pre.truth.time)*0.5;
0152         pre.truth.deposit += post.truth.deposit;
0153         mean_pos.SetX(mean_pos.x()+mean.x()*post.truth.deposit);
0154         mean_pos.SetY(mean_pos.y()+mean.y()*post.truth.deposit);
0155         mean_pos.SetZ(mean_pos.z()+mean.z()*post.truth.deposit);
0156         mean_time += mean_tm*post.truth.deposit;
0157         step_length  += step->GetStepLength();
0158         if ( 0 == cell )   {
0159           cell = sensitive->cellID(step);
0160           if ( 0 == cell )  {
0161             cell = sensitive->volumeID(step) ;
0162             sensitive->except("+++ Invalid CELL ID for hit!");
0163           }
0164         }
0165         ++combined;
0166         return *this;
0167       }
0168
0169       /// Helper function to decide if the hit has to be extracted and saved in the collection
0170       bool mustSaveTrack(const G4Track* tr)  const   {
0171         return current > 0 && current != tr->GetTrackID();
0172       }
0173
0174       /// Extract hit information and add the created hit to the collection
0175       void extractHit(EInside ended)   {
0176         Geant4HitCollection* collection = sensitive->collection(0);
0177         extractHit(collection, ended);
0178       }
0179
0180       TrackerWeighted& calc_dist_out(const G4VSolid* solid)    {
0181         Position v(pre.momentum.unit()), &p=post.position;
0182         double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()),
0183                                            G4ThreeVector(v.X(),v.Y(),v.Z()));
0184         distance_to_outside = dist;
0185         return *this;
0186       }
0187
0188       TrackerWeighted& calc_dist_in(const G4VSolid* solid)    {
0189         Position v(pre.momentum.unit()), &p=pre.position;
0190         double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()),
0191                                            G4ThreeVector(v.X(),v.Y(),v.Z()));
0192         distance_to_inside = dist;
0193         return *this;
0194       }
0195
0196       void extractHit(Geant4HitCollection* collection, EInside ended)   {
0197         if( DEBUG == printLevel() ) {
0198           std::cout<<" DEBUG: Geant4TrackerWeightedSD::extractHit(Geant4HitCollection* collection, EInside ended) ...."<<std::endl;
0199           std::cout<<" DEBUG: =================================================="<<std::endl;
0200         }
0201
0202         double deposit  = pre.truth.deposit;
0203         if ( current != -1 )  {
0204           Position pos;
0205           Momentum mom;
0206           double   time = deposit != 0 ? mean_time / deposit : mean_time;
0207           char     dist_in[64], dist_out[64];
0208
0209           switch(hit_position_type)  {
0210           case POSITION_WEIGHTED:
0211             pos = deposit != 0 ? mean_pos / deposit : mean_pos;
0212             mom = 0.5 * (pre.momentum + post.momentum);
0213             break;
0214           case POSITION_PREPOINT:
0215             pos = pre.position;
0216             mom = pre.momentum;
0217             break;
0218           case POSITION_POSTPOINT:
0219             pos = post.position;
0220             mom = post.momentum;
0221             break;
0222           case POSITION_MIDDLE:
0223           default:
0224             pos = (post.position + pre.position) / 2.0;
0225             mom = 0.5 * (pre.momentum + post.momentum);
0226             break;
0227           }
0228
0229           if ( ended == kSurface || distance_to_outside < std::numeric_limits<float>::epsilon() )
0230             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_SURFACE;
0231           else if ( ended == kInside )
0232             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_INSIDE;
0233           else if ( ended == kOutside )
0234             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_OUTSIDE;
0235
0236           Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID,
0237                                                            pre.truth.pdgID,
0238                                                            deposit,time, step_length,
0239                                                            pos, mom);
0240           hit->flag     = hit_flag;
0241           hit->cellID   = cell;
0242           hit->g4ID     = g4ID;
0243
0244           dist_in[0] = dist_out[0] = 0;
0245           if ( !(hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) )
0246             ::snprintf(dist_in,sizeof(dist_in)," [%.2e um]",distance_to_inside/CLHEP::um);
0247           if ( !(hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE) )
0248             ::snprintf(dist_out,sizeof(dist_out)," [%.2e um]",distance_to_outside/CLHEP::um);
0249           sensitive->print("+++ G4Track:%5d CREATE hit[%03d]:%3d deps E:"
0250                            " %.2e keV Pos:%7.2f %7.2f %7.2f [mm] Start:%s%s%s%s End:%s%s%s%s",
0251                            pre.truth.trackID,int(collection->GetSize()),
0252                            combined,pre.truth.deposit/CLHEP::keV,
0253                            pos.X()/CLHEP::mm,pos.Y()/CLHEP::mm,pos.Z()/CLHEP::mm,
0254                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) ? "SURFACE" : ""),
0255                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_OUTSIDE) ? "OUTSIDE" : ""),
0256                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_INSIDE)  ? "INSIDE " : ""),
0257                            dist_in,
0258                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE)   ? "SURFACE" : ""),
0259                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_OUTSIDE)   ? "OUTSIDE" : ""),
0260                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_INSIDE)    ? "INSIDE " : ""),
0261                            dist_out);
0262           collection->add(hit);
0263         }
0264         clear();
0265       }
0266
0267       /// Method for generating hit(s) using the information of G4Step object.
0268       G4bool process(const G4Step* step, G4TouchableHistory* )   {
0269         Geant4StepHandler h(step);
0270         if( DEBUG == printLevel() ) {
0271           std::cout<<" DEBUG: Geant4TrackerWeightedSD::process(const G4Step* step, G4TouchableHistory* ) ...."<<std::endl;
0272           dumpStep( h, step);
0273         }
0274
0275         // std::cout << " process called - pre pos: " << h.prePos() << " post pos " << h.postPos()
0276         //           << " edep: " << h.deposit() << std::endl ;
0277
0278         G4VSolid*     preSolid    = h.solid(h.pre);
0279         G4VSolid*     postSolid   = h.solid(h.post);
0280         G4ThreeVector local_pre   = h.globalToLocalG4(h.prePosG4());
0281         G4ThreeVector local_post  = h.globalToLocalG4(h.postPosG4());
0282         EInside       pre_inside  = preSolid->Inside(local_pre);
0283         EInside       post_inside = postSolid->Inside(local_post);
0284
0285         const void* postPV = h.postVolume();
0286         const void* prePV  = h.preVolume();
0287         const void* postSD = h.postSD();
0288         const void* preSD  = h.preSD();
0289         G4VSolid* solid = (preSD == thisSD) ? preSolid : postSolid;
0290         // Track went into new Volume, extracted the hit in prePV, then start a new hit in thisPV.
0291         if ( current == h.trkID() && thisPV != 0 && prePV != thisPV )  {
0292           if( DEBUG == printLevel() ) {
0293             std::cout<<" DEBUG: Geant4TrackerWeightedSD: if ( current == h.trkID() && thisPV != 0 && prePV != thisPV ),"
0294                      <<" Track went into new Volume, extracted the hit in prePV, then start a new hit in thisPV."
0295                      << std::endl;
0296           }
0297           extractHit(post_inside);
0298           start(step, h.pre);
0299         }
0300         // 1) Track killed inside SD: trace incomplete. This deposition must be added as well.
0301         else if ( current == h.trkID() && !h.trkAlive() )  {
0302           hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
0303           update(step).calc_dist_out(solid).extractHit(post_inside);
0304           return true;
0305         }
0306         // 2) Track leaving SD volume. Sensitive detector changed. Store hit.
0307         else if ( current == h.trkID() && postSD != thisSD )  {
0308           update(step).calc_dist_out(solid).extractHit(kOutside);
0309           return true;
0310         }
0311         // 3) Track leaving SD volume. Store hit.
0312         else if ( current == h.trkID() && postSD == thisSD && post_inside == kSurface )  {
0313           update(step).calc_dist_out(solid).extractHit(kSurface);
0314           return true;
0315         }
0316         // 4) Track leaving SD volume. Store hit.
0317         else if ( current == h.trkID() && postSD == thisSD && post_inside == kOutside )  {
0318           update(step).calc_dist_out(solid).extractHit(post_inside);
0319           return true;
0320         }
0321         // 5) Normal update: either intermediate deposition or track is going to be killed.
0322         else if ( current == h.trkID() && postSD == thisSD && post_inside == kInside )  {
0323           last_inside = post_inside;
0324           update(step).calc_dist_out(solid);
0325           return true;
0326         }
0327         // 6) Track from secondary created in SD volume. Store hit from previous.
0328         // --> New hit will be created, to whom also this deposition belongs
0329         else if ( current != h.trkID() && current >= 0 )  {
0330           extractHit(last_inside);
0331         }
0332
0333         // If the hit got extracted, a new one must be set up
0334         if ( current < 0 )  {
0335           EInside  inside  = pre_inside;
0336           // Track entering SD volume
0337           if ( preSD != thisSD )   {
0338             start(step, h.post);
0339             inside = post_inside;
0340             sensitive->print("++++++++++ Using POST step volume to start hit -- dubious ?");
0341           }
0342           else  {
0343             start(step, h.pre);
0344           }
0345           calc_dist_in(solid);
0346           if ( inside == kSurface )
0347             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_SURFACE;
0348           else if ( inside == kInside )
0349             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_INSIDE;
0350           else if ( inside == kOutside )
0351             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_OUTSIDE;
0352
0353           // New (secondary) track created by some process starting inside the volume
0354           if ( inside == kInside )  {
0355             hit_flag |= Geant4Tracker::Hit::HIT_SECONDARY_TRACK;
0356           }
0357         }
0358
0359         // Update Data
0360         last_inside = post_inside;
0361         update(step);
0362         calc_dist_out(solid);
0363
0364         // Track killed inside SD: trace incomplete. This deposition must be added as well.
0365         if ( !h.trkAlive() )  {
0366           hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
0367           extractHit(post_inside);
0368         }
0369         // Avoid dangling hits if the track leaves the sensitive volume
0370         else if ( post_inside == kSurface )  {
0371           extractHit(post_inside);
0372         }
0373         // Avoid dangling hits if the track leaves the sensitive volume
0374         else if ( thisSD == preSD && (preSD != postSD || prePV != postPV) )  {
0375           extractHit(post_inside);
0376         }
0377         // This should simply not happen!
0378         else if ( thisSD == postSD && (preSD != postSD || prePV != postPV) )  {
0379           sensitive->error("+++++ WRONG!!! Extract. How did we get here?");
0380           extractHit(post_inside);
0381         }
0382         // In single hit mode we MUST write the hit after update
0383         else if ( single_deposit_mode )  {
0384           extractHit(post_inside);
0385         }
0386
0387         return true;
0388       }
0389
0390       /// Post-event action callback
0391       void endEvent()   {
0392         // We need to add the possibly last added hit to the collection here.
0393         // otherwise the last hit would be assigned to the next event and the
0394         // MC truth would be screwed.
0395         //
0396         if ( current > 0 )   {
0397           Geant4HitCollection* coll = sensitive->collection(0);
0398           sensitive->print("++++++++++ Found dangling hit: Is the hit extraction logic correct?");
0399           extractHit(coll,last_inside);
0400         }
0401       }
0402       /// Pre event action callback
0403       void startEvent()   {
0404         thisSD = dynamic_cast<G4VSensitiveDetector*>(&sensitive->detector());
0405       }
0406
0407       ///dumpStep
0408       void dumpStep(const Geant4StepHandler& h, const G4Step* s)  {
0409         std::stringstream str;
0410         str << " ----- step in detector " << h.sdName( s->GetPreStepPoint() )
0411             << " prePos  " << h.prePos()
0412             << " postPos " << h.postPos()
0413             << " preStatus  " << h.preStepStatus()
0414             << " postStatus  " << h.postStepStatus()
0415             << " preVolume " << h.volName( s->GetPreStepPoint() )
0416             << " postVolume " << h.volName( s->GetPostStepPoint() )
0417             << std::endl
0418             << "     momentum : "  << std::scientific
0419             <<  s->GetPreStepPoint()->GetMomentum()[0] << ", "
0420             <<  s->GetPreStepPoint()->GetMomentum()[1]<< ", "
0421             <<  s->GetPreStepPoint()->GetMomentum()[2]
0422             << " / "
0423             << s->GetPostStepPoint()->GetMomentum()[0] << ", "
0424             <<  s->GetPostStepPoint()->GetMomentum()[1] << ", "
0425             <<  s->GetPostStepPoint()->GetMomentum()[2]
0426             << ", PDG: " << s->GetTrack()->GetDefinition()->GetPDGEncoding();
0427         std::cout << str.str() << std::endl;
0428       }
0429
0430       /// GFLash processing callback
0431       G4bool process(const Geant4FastSimSpot* , G4TouchableHistory* ) {
0432         sensitive->except("GFlash/FastSim action is not implemented for SD: %s", sensitive->c_name());
0433         return false;
0434       }
0435     };
0436
0437     /// Initialization overload for specialization
0438     template <> void Geant4SensitiveAction<TrackerWeighted>::initialize() {
0439       declareProperty("HitPositionCombination", m_userData.hit_position_type);
0440       declareProperty("CollectSingleDeposits", m_userData.single_deposit_mode);
0441       m_userData.e_cut = m_sensitive.energyCutoff();
0442       m_userData.sensitive = this;
0443     }
0444
0445     /// G4VSensitiveDetector interface: Method invoked at the beginning of each event.
0446     template <> void Geant4SensitiveAction<TrackerWeighted>::begin(G4HCofThisEvent* /* hce */)   {
0447       m_userData.startEvent();
0448     }
0449
0450     /// G4VSensitiveDetector interface: Method invoked at the end of each event.
0451     template <> void Geant4SensitiveAction<TrackerWeighted>::end(G4HCofThisEvent* /* hce */)   {
0452       m_userData.endEvent();
0453     }
0454
0455     /// Define collections created by this sensitivie action object
0456     template <> void Geant4SensitiveAction<TrackerWeighted>::defineCollections() {
0457       m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>();
0458     }
0459
0460     /// Method for generating hit(s) using the information of G4Step object.
0461     template <> void Geant4SensitiveAction<TrackerWeighted>::clear(G4HCofThisEvent* /* hce */) {
0462       m_userData.clear();
0463     }
0464
0465     /// Method for generating hit(s) using the information of G4Step object.
0466     template <> G4bool
0467     Geant4SensitiveAction<TrackerWeighted>::process(const G4Step* step, G4TouchableHistory* history) {
0468       return m_userData.process(step, history);
0469     }
0470
0471     /// Method for generating hit(s) using the information of the Geant4FastSimSpot object.
0472     template <> bool
0473     Geant4SensitiveAction<TrackerWeighted>::processFastSim(const Geant4FastSimSpot* spot, G4TouchableHistory* history) {
0474       return m_userData.process(spot, history);
0475     }
0476     typedef Geant4SensitiveAction<TrackerWeighted>  Geant4TrackerWeightedAction;
0477   }
0478 }
0479
0480 using namespace dd4hep::sim;
0481
0482 #include <DDG4/Factories.h>
0483 DECLARE_GEANT4SENSITIVE(Geant4TrackerWeightedAction)