DDG4/plugins/Geant4MaterialScanner.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/Objects.h>
0016 #include <DDG4/Defs.h>
0017 #include <DDG4/Geant4SteppingAction.h>
0018
0019 // Forward declarations
0020 class G4LogicalVolume;
0021
0022 /// Namespace for the AIDA detector description toolkit
0023 namespace dd4hep {
0024
0025   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0026   namespace sim   {
0027
0028     /// Class to perform directional material scans using Geantinos.
0029     /**
0030      *
0031      *  \author  M.Frank
0032      *  \version 1.0
0033      *  \ingroup DD4HEP_SIMULATION
0034      */
0035     class Geant4MaterialScanner : public Geant4SteppingAction  {
0036     protected:
0037       /// Structure to hold the information of one simulation step.
0038       class StepInfo {
0039       public:
0040         /// Pre-step and Post-step position
0041         Position pre, post;
0042         /// Reference to the logical volue
0043         const G4LogicalVolume* volume;
0044
0045         /// Initializing constructor
0046         StepInfo(const Position& pre, const Position& post, const G4LogicalVolume* volume);
0047         /// Copy constructor
0048         StepInfo(const StepInfo& c);
0049         /// Default destructor
0050         ~StepInfo() {}
0051         /// Assignment operator
0052         StepInfo& operator=(const StepInfo& c);
0053       };
0054       typedef std::vector<StepInfo*> Steps;
0055
0056       double m_sumX0     = 0E0;
0057       double m_sumLambda = 0E0;
0058       double m_sumPath   = 0E0;
0059       Steps  m_steps;
0060
0061     public:
0062       /// Standard constructor
0063       Geant4MaterialScanner(Geant4Context* context, const std::string& name);
0064       /// Default destructor
0065       virtual ~Geant4MaterialScanner();
0066       /// User stepping callback
0067       virtual void operator()(const G4Step* step, G4SteppingManager* mgr);
0068       /// Begin-of-tracking callback
0069       virtual void begin(const G4Track* track);
0070       /// End-of-tracking callback
0071       virtual void end(const G4Track* track);
0072       /// Registered callback on Begin-event
0073       void beginEvent(const G4Event* event);
0074     };
0075   }
0076 }
0077
0078 //====================================================================
0079 //  AIDA Detector description implementation
0080 //--------------------------------------------------------------------
0081 //
0082 //  Author     : M.Frank
0083 //
0084 //====================================================================
0085
0086 // Framework include files
0087 #include <DD4hep/InstanceCount.h>
0088 #include <DD4hep/Printout.h>
0089 #include <DDG4/Geant4TouchableHandler.h>
0090 #include <DDG4/Geant4StepHandler.h>
0091 #include <DDG4/Geant4EventAction.h>
0092 #include <DDG4/Geant4TrackingAction.h>
0093 #include <CLHEP/Units/SystemOfUnits.h>
0094 #include <G4LogicalVolume.hh>
0095 #include <G4Material.hh>
0096
0097 using namespace dd4hep::sim;
0098
0099 #include <DDG4/Factories.h>
0100 DECLARE_GEANT4ACTION(Geant4MaterialScanner)
0101
0102 /// Initializing constructor
0103 Geant4MaterialScanner::StepInfo::StepInfo(const Position& prePos, const Position& postPos, const G4LogicalVolume* vol)
0104 : pre(prePos), post(postPos), volume(vol)
0105 {
0106 }
0107
0108 /// Copy constructor
0109 Geant4MaterialScanner::StepInfo::StepInfo(const StepInfo& c)
0110   : pre(c.pre), post(c.post), volume(c.volume)
0111 {
0112 }
0113
0114 /// Assignment operator
0115 Geant4MaterialScanner::StepInfo& Geant4MaterialScanner::StepInfo::operator=(const StepInfo& c)  {
0116   pre = c.pre;
0117   post = c.post;
0118   volume = c.volume;
0119   return *this;
0120 }
0121
0122 /// Standard constructor
0123 Geant4MaterialScanner::Geant4MaterialScanner(Geant4Context* ctxt, const std::string& nam)
0124   : Geant4SteppingAction(ctxt,nam)
0125 {
0126   m_needsControl = true;
0127   eventAction().callAtBegin(this,&Geant4MaterialScanner::beginEvent);
0128   trackingAction().callAtEnd(this,&Geant4MaterialScanner::end);
0129   trackingAction().callAtBegin(this,&Geant4MaterialScanner::begin);
0130   InstanceCount::increment(this);
0131 }
0132
0133 /// Default destructor
0134 Geant4MaterialScanner::~Geant4MaterialScanner() {
0135   InstanceCount::decrement(this);
0136 }
0137
0138 /// User stepping callback
0139 void Geant4MaterialScanner::operator()(const G4Step* step, G4SteppingManager*) {
0140   Geant4StepHandler h(step);
0141 #if 0
0142   Geant4TouchableHandler pre_handler(step);
0143   std::string prePath = pre_handler.path();
0144   Geant4TouchableHandler post_handler(step);
0145   std::string postPath = post_handler.path();
0146 #endif
0147   G4LogicalVolume* logVol = h.logvol(h.pre);
0148   m_steps.emplace_back(new StepInfo(h.prePos(), h.postPos(), logVol));
0149 }
0150
0151 /// Registered callback on Begin-event
0152 void Geant4MaterialScanner::beginEvent(const G4Event* /* event */)   {
0153   for_each(m_steps.begin(),m_steps.end(),detail::DestroyObject<StepInfo*>());
0154   m_steps.clear();
0155   m_sumX0 = 0;
0156   m_sumLambda = 0;
0157   m_sumPath = 0;
0158 }
0159
0160 /// Begin-of-tracking callback
0161 void Geant4MaterialScanner::begin(const G4Track* track) {
0162   printP2("Starting tracking action for track ID=%d",track->GetTrackID());
0163   for_each(m_steps.begin(),m_steps.end(),detail::DestroyObject<StepInfo*>());
0164   m_steps.clear();
0165   m_sumX0 = 0;
0166   m_sumLambda = 0;
0167   m_sumPath = 0;
0168 }
0169
0170 /// End-of-tracking callback
0171 void Geant4MaterialScanner::end(const G4Track* track) {
0172   if ( !m_steps.empty() )  {
0173     constexpr const char* line = " +--------------------------------------------------------------------------------------------------------------------------------------------------\n";
0174     constexpr const char* fmt1 = " | %5d %-20s %3.0f %8.3f %8.4f %11.4f  %11.4f %10.3f %8.2f %11.6f %11.6f  (%7.2f,%7.2f,%7.2f)\n";
0175     constexpr const char* fmt2 = " | %5d %-20s %3.0f %8.3f %8.4f %11.6g  %11.6g %10.3f %8.2f %11.6f %11.6f  (%7.2f,%7.2f,%7.2f)\n";
0176     const Position& pre = m_steps[0]->pre;
0177     const Position& post = m_steps[m_steps.size()-1]->post;
0178
0179     ::printf("%s + Material scan between: x_0 = (%7.2f,%7.2f,%7.2f) [cm] and x_1 = (%7.2f,%7.2f,%7.2f) [cm]  TrackID:%d: \n%s",
0180              line,pre.X()/CLHEP::cm,pre.Y()/CLHEP::cm,pre.Z()/CLHEP::cm,post.X()/CLHEP::cm,post.Y()/CLHEP::cm,post.Z()/CLHEP::cm,track->GetTrackID(),line);
0181     ::printf(" |     \\   %-11s        Atomic                 Radiation   Interaction               Path   Integrated  Integrated    Material\n","Material");
0182     ::printf(" | Num. \\  %-11s   Number/Z   Mass/A  Density    Length       Length    Thickness   Length      X0        Lambda      Endpoint  \n","Name");
0183     ::printf(" | Layer \\ %-11s            [g/mole]  [g/cm3]     [cm]        [cm]          [cm]      [cm]     [cm]        [cm]     (     cm,     cm,     cm)\n","");
0184     ::printf("%s",line);
0185     int count = 1;
0186     for(Steps::const_iterator i=m_steps.begin(); i!=m_steps.end(); ++i, ++count)  {
0187       const G4LogicalVolume* logVol = (*i)->volume;
0188       G4Material* material = logVol->GetMaterial();
0189       const Position& prePos  = (*i)->pre;
0190       const Position& postPos = (*i)->post;
0191       Position direction = postPos - prePos;
0192       double length  = direction.R()/CLHEP::cm;
0193       double intLen  = material->GetNuclearInterLength()/CLHEP::cm;
0194       double radLen  = material->GetRadlen()/CLHEP::cm;
0195       double density = material->GetDensity()/(CLHEP::gram/CLHEP::cm3);
0196       double nLambda = length / intLen;
0197       double nx0     = length / radLen;
0198       double Aeff    = 0.0;
0199       double Zeff    = 0.0;
0200       const char* fmt = radLen >= 1e5 ? fmt2 : fmt1;
0201       const double* fractions = material->GetFractionVector();
0202       for(size_t j=0; j<material->GetNumberOfElements(); ++j)  {
0203         Zeff += fractions[j]*(material->GetElement(j)->GetZ());
0204         Aeff += fractions[j]*(material->GetElement(j)->GetA())/CLHEP::gram;
0205       }
0206       m_sumX0     += nx0;
0207       m_sumLambda += nLambda;
0208       m_sumPath   += length;
0209       ::printf(fmt,count,material->GetName().c_str(),
0210                Zeff, Aeff, density, radLen, intLen, length,
0211                m_sumPath,m_sumX0,m_sumLambda,
0212                postPos.X()/CLHEP::cm,postPos.Y()/CLHEP::cm,postPos.Z()/CLHEP::cm);
0213       //cout << *m << endl;
0214     }
0215     for_each(m_steps.begin(),m_steps.end(),detail::DestroyObject<StepInfo*>());
0216     m_steps.clear();
0217   }
0218 }