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  = 0.5 * (pre.momentum + post.momentum);
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             break;
0213           case POSITION_PREPOINT:
0214             pos = pre.position;
0215             break;
0216           case POSITION_POSTPOINT:
0217             pos = post.position;
0218             break;
0219           case POSITION_MIDDLE:
0220           default:
0221             pos = (post.position + pre.position) / 2.0;
0222             break;
0223           }
0224
0225           if ( ended == kSurface || distance_to_outside < std::numeric_limits<float>::epsilon() )
0226             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_SURFACE;
0227           else if ( ended == kInside )
0228             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_INSIDE;
0229           else if ( ended == kOutside )
0230             hit_flag |= Geant4Tracker::Hit::HIT_ENDED_OUTSIDE;
0231
0232           Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID,
0233                                                            pre.truth.pdgID,
0234                                                            deposit,time, step_length,
0235                                pos, mom);
0236           hit->flag     = hit_flag;
0237           hit->cellID   = cell;
0238           hit->g4ID     = g4ID;
0239
0240           dist_in[0] = dist_out[0] = 0;
0241           if ( !(hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) )
0242             ::snprintf(dist_in,sizeof(dist_in)," [%.2e um]",distance_to_inside/CLHEP::um);
0243           if ( !(hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE) )
0244             ::snprintf(dist_out,sizeof(dist_out)," [%.2e um]",distance_to_outside/CLHEP::um);
0245           sensitive->print("+++ G4Track:%5d CREATE hit[%03d]:%3d deps E:"
0246                            " %.2e keV Pos:%7.2f %7.2f %7.2f [mm] Start:%s%s%s%s End:%s%s%s%s",
0247                            pre.truth.trackID,int(collection->GetSize()),
0248                            combined,pre.truth.deposit/CLHEP::keV,
0249                            pos.X()/CLHEP::mm,pos.Y()/CLHEP::mm,pos.Z()/CLHEP::mm,
0250                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) ? "SURFACE" : ""),
0251                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_OUTSIDE) ? "OUTSIDE" : ""),
0252                            ((hit_flag&Geant4Tracker::Hit::HIT_STARTED_INSIDE)  ? "INSIDE " : ""),
0253                            dist_in,
0254                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE)   ? "SURFACE" : ""),
0255                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_OUTSIDE)   ? "OUTSIDE" : ""),
0256                            ((hit_flag&Geant4Tracker::Hit::HIT_ENDED_INSIDE)    ? "INSIDE " : ""),
0257                            dist_out);
0258           collection->add(hit);
0259         }
0260         clear();
0261       }
0262
0263       /// Method for generating hit(s) using the information of G4Step object.
0264       G4bool process(const G4Step* step, G4TouchableHistory* )   {
0265         Geant4StepHandler h(step);
0266         if( DEBUG == printLevel() ) {
0267           std::cout<<" DEBUG: Geant4TrackerWeightedSD::process(const G4Step* step, G4TouchableHistory* ) ...."<<std::endl;
0268           dumpStep( h, step);
0269         }
0270
0271         // std::cout << " process called - pre pos: " << h.prePos() << " post pos " << h.postPos()
0272         //    << " edep: " << h.deposit() << std::endl ;
0273
0274         G4VSolid*     preSolid    = h.solid(h.pre);
0275         G4VSolid*     postSolid   = h.solid(h.post);
0276         G4ThreeVector local_pre   = h.globalToLocalG4(h.prePosG4());
0277         G4ThreeVector local_post  = h.globalToLocalG4(h.postPosG4());
0278         EInside       pre_inside  = preSolid->Inside(local_pre);
0279         EInside       post_inside = postSolid->Inside(local_post);
0280
0281         const void* postPV = h.postVolume();
0282         const void* prePV  = h.preVolume();
0283         const void* postSD = h.postSD();
0284         const void* preSD  = h.preSD();
0285         G4VSolid* solid = (preSD == thisSD) ? preSolid : postSolid;
0286         // Track went into new Volume, extracted the hit in prePV, then start a new hit in thisPV.
0287         if ( current == h.trkID() && thisPV != 0 && prePV != thisPV )  {
0288           if( DEBUG == printLevel() ) {
0289             std::cout<<" DEBUG: Geant4TrackerWeightedSD: if ( current == h.trkID() && thisPV != 0 && prePV != thisPV ),"
0290                      <<" Track went into new Volume, extracted the hit in prePV, then start a new hit in thisPV."
0291                      << std::endl;
0292           }
0293           extractHit(post_inside);
0294           start(step, h.pre);
0295         }
0296         // 1) Track killed inside SD: trace incomplete. This deposition must be added as well.
0297         else if ( current == h.trkID() && !h.trkAlive() )  {
0298           hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
0299           update(step).calc_dist_out(solid).extractHit(post_inside);
0300           return true;
0301         }
0302         // 2) Track leaving SD volume. Sensitive detector changed. Store hit.
0303         else if ( current == h.trkID() && postSD != thisSD )  {
0304           update(step).calc_dist_out(solid).extractHit(kOutside);
0305           return true;
0306         }
0307         // 3) Track leaving SD volume. Store hit.
0308         else if ( current == h.trkID() && postSD == thisSD && post_inside == kSurface )  {
0309           update(step).calc_dist_out(solid).extractHit(kSurface);
0310           return true;
0311         }
0312         // 4) Track leaving SD volume. Store hit.
0313         else if ( current == h.trkID() && postSD == thisSD && post_inside == kOutside )  {
0314           update(step).calc_dist_out(solid).extractHit(post_inside);
0315           return true;
0316         }
0317         // 5) Normal update: either intermediate deposition or track is going to be killed.
0318         else if ( current == h.trkID() && postSD == thisSD && post_inside == kInside )  {
0319           last_inside = post_inside;
0320           update(step).calc_dist_out(solid);
0321           return true;
0322         }
0323         // 6) Track from secondary created in SD volume. Store hit from previous.
0324         // --> New hit will be created, to whom also this deposition belongs
0325         else if ( current != h.trkID() && current >= 0 )  {
0326           extractHit(last_inside);
0327         }
0328
0329         // If the hit got extracted, a new one must be set up
0330         if ( current < 0 )  {
0331           EInside  inside  = pre_inside;
0332           // Track entering SD volume
0333           if ( preSD != thisSD )   {
0334             start(step, h.post);
0335             inside = post_inside;
0336             sensitive->print("++++++++++ Using POST step volume to start hit -- dubious ?");
0337           }
0338           else  {
0339             start(step, h.pre);
0340           }
0341           calc_dist_in(solid);
0342           if ( inside == kSurface )
0343             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_SURFACE;
0344           else if ( inside == kInside )
0345             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_INSIDE;
0346           else if ( inside == kOutside )
0347             hit_flag |= Geant4Tracker::Hit::HIT_STARTED_OUTSIDE;
0348
0349           // New (secondary) track created by some process starting inside the volume
0350           if ( inside == kInside )  {
0351             hit_flag |= Geant4Tracker::Hit::HIT_SECONDARY_TRACK;
0352           }
0353         }
0354
0355         // Update Data
0356         last_inside = post_inside;
0357         update(step);
0358         calc_dist_out(solid);
0359
0360         // Track killed inside SD: trace incomplete. This deposition must be added as well.
0361         if ( !h.trkAlive() )  {
0362           hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK;
0363           extractHit(post_inside);
0364         }
0365         // Avoid dangling hits if the track leaves the sensitive volume
0366         else if ( post_inside == kSurface )  {
0367           extractHit(post_inside);
0368         }
0369         // Avoid dangling hits if the track leaves the sensitive volume
0370         else if ( thisSD == preSD && (preSD != postSD || prePV != postPV) )  {
0371           extractHit(post_inside);
0372         }
0373         // This should simply not happen!
0374         else if ( thisSD == postSD && (preSD != postSD || prePV != postPV) )  {
0375           sensitive->error("+++++ WRONG!!! Extract. How did we get here?");
0376           extractHit(post_inside);
0377         }
0378         // In single hit mode we MUST write the hit after update
0379         else if ( single_deposit_mode )  {
0380           extractHit(post_inside);
0381         }
0382
0383         return true;
0384       }
0385
0386       /// Post-event action callback
0387       void endEvent()   {
0388         // We need to add the possibly last added hit to the collection here.
0389         // otherwise the last hit would be assigned to the next event and the
0390         // MC truth would be screwed.
0391         //
0392         if ( current > 0 )   {
0393           Geant4HitCollection* coll = sensitive->collection(0);
0394           sensitive->print("++++++++++ Found dangling hit: Is the hit extraction logic correct?");
0395           extractHit(coll,last_inside);
0396         }
0397       }
0398       /// Pre event action callback
0399       void startEvent()   {
0400         thisSD = dynamic_cast<G4VSensitiveDetector*>(&sensitive->detector());
0401       }
0402
0403       ///dumpStep
0404       void dumpStep(const Geant4StepHandler& h, const G4Step* s)  {
0405         std::stringstream str;
0406         str << " ----- step in detector " << h.sdName( s->GetPreStepPoint() )
0407             << " prePos  " << h.prePos()
0408             << " postPos " << h.postPos()
0409             << " preStatus  " << h.preStepStatus()
0410             << " postStatus  " << h.postStepStatus()
0411             << " preVolume " << h.volName( s->GetPreStepPoint() )
0412             << " postVolume " << h.volName( s->GetPostStepPoint() )
0413             << std::endl
0414             << "     momentum : "  << std::scientific
0415             <<  s->GetPreStepPoint()->GetMomentum()[0] << ", "
0416             <<  s->GetPreStepPoint()->GetMomentum()[1]<< ", "
0417             <<  s->GetPreStepPoint()->GetMomentum()[2]
0418             << " / "
0419             << s->GetPostStepPoint()->GetMomentum()[0] << ", "
0420             <<  s->GetPostStepPoint()->GetMomentum()[1] << ", "
0421             <<  s->GetPostStepPoint()->GetMomentum()[2]
0422             << ", PDG: " << s->GetTrack()->GetDefinition()->GetPDGEncoding();
0423         std::cout << str.str() << std::endl;
0424       }
0425
0426       /// GFLash processing callback
0427       G4bool process(const Geant4FastSimSpot* , G4TouchableHistory* ) {
0428         sensitive->except("GFlash/FastSim action is not implemented for SD: %s", sensitive->c_name());
0429         return false;
0430       }
0431     };
0432
0433     /// Initialization overload for specialization
0434     template <> void Geant4SensitiveAction<TrackerWeighted>::initialize() {
0435       declareProperty("HitPositionCombination", m_userData.hit_position_type);
0436       declareProperty("CollectSingleDeposits", m_userData.single_deposit_mode);
0437       m_userData.e_cut = m_sensitive.energyCutoff();
0438       m_userData.sensitive = this;
0439     }
0440
0441     /// G4VSensitiveDetector interface: Method invoked at the begining of each event.
0442     template <> void Geant4SensitiveAction<TrackerWeighted>::begin(G4HCofThisEvent* /* hce */)   {
0443       m_userData.startEvent();
0444     }
0445
0446     /// G4VSensitiveDetector interface: Method invoked at the end of each event.
0447     template <> void Geant4SensitiveAction<TrackerWeighted>::end(G4HCofThisEvent* /* hce */)   {
0448       m_userData.endEvent();
0449     }
0450
0451     /// Define collections created by this sensitivie action object
0452     template <> void Geant4SensitiveAction<TrackerWeighted>::defineCollections() {
0453       m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>();
0454     }
0455
0456     /// Method for generating hit(s) using the information of G4Step object.
0457     template <> void Geant4SensitiveAction<TrackerWeighted>::clear(G4HCofThisEvent* /* hce */) {
0458       m_userData.clear();
0459     }
0460
0461     /// Method for generating hit(s) using the information of G4Step object.
0462     template <> G4bool
0463     Geant4SensitiveAction<TrackerWeighted>::process(const G4Step* step, G4TouchableHistory* history) {
0464       return m_userData.process(step, history);
0465     }
0466
0467     /// Method for generating hit(s) using the information of the Geant4FastSimSpot object.
0468     template <> bool
0469     Geant4SensitiveAction<TrackerWeighted>::processFastSim(const Geant4FastSimSpot* spot, G4TouchableHistory* history) {
0470       return m_userData.process(spot, history);
0471     }
0472     typedef Geant4SensitiveAction<TrackerWeighted>  Geant4TrackerWeightedAction;
0473   }
0474 }
0475
0476 using namespace dd4hep::sim;
0477
0478 #include <DDG4/Factories.h>
0479 DECLARE_GEANT4SENSITIVE(Geant4TrackerWeightedAction)