diff --git a/Calorimeter/CaloG4Sim/CaloG4Sim/CalibrationDefaultProcessing.h b/Calorimeter/CaloG4Sim/CaloG4Sim/CalibrationDefaultProcessing.h new file mode 100755 index 0000000000000000000000000000000000000000..d8eff3049605b32ef01214311e85843a975de839 --- /dev/null +++ b/Calorimeter/CaloG4Sim/CaloG4Sim/CalibrationDefaultProcessing.h @@ -0,0 +1,60 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// CaloG4::CalibrationDefaultProcessing +// 04-Mar-2004 William Seligman + +// 21-Sep-2004 WGS: Moved to CaloG4Sim. + +// The calibration studies rely on every volume in the simulation +// being made into a sensitive detector. There is a practical +// problem: What if we're still in the middle of developing code, and +// not every volume has been made sensitive yet? What if we've +// overlooked a volume? + +// This class provides a "default behavior" for all energy deposits +// that are not made in a volume that's been made sensitive for +// calibration studies. + +// Since we need to examine the particles at every G4Step, this class +// makes use of a UserAction interface class. + +#ifndef CaloG4_CalibrationDefaultProcessing_h +#define CaloG4_CalibrationDefaultProcessing_h + +// The Athena and stand-alone G4 setups have different UserAction +// classes (though they accomplish the same thing. + +#include "FadsActions/UserAction.h" + +// Forward declarations +class G4Step; +class G4VSensitiveDetector; + +namespace CaloG4 { + + class CalibrationDefaultProcessing : public FADS::UserAction { + + public: + + CalibrationDefaultProcessing(); + virtual ~CalibrationDefaultProcessing(); + + virtual void SteppingAction(const G4Step*); + + // Make the default sensitive detector available to other routines. + static G4VSensitiveDetector* GetDefaultSD() { return m_defaultSD; } + static void SetDefaultSD( G4VSensitiveDetector* ); + + private: + + // The default sensitive detector to be applied to all G4Steps + // in volumes without a CalibrationSensitiveDetector. + static G4VSensitiveDetector* m_defaultSD; + + }; + +} // namespace CaloG4 + +#endif // CaloG4_CalibrationDefaultProcessing_h diff --git a/Calorimeter/CaloG4Sim/CaloG4Sim/EscapedEnergyRegistry.h b/Calorimeter/CaloG4Sim/CaloG4Sim/EscapedEnergyRegistry.h new file mode 100755 index 0000000000000000000000000000000000000000..92023b8e53c7988960a3a0c3b1b20098440b00ed --- /dev/null +++ b/Calorimeter/CaloG4Sim/CaloG4Sim/EscapedEnergyRegistry.h @@ -0,0 +1,64 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// EscapedEnergyRegistry +// 15-Jul-2004 William Seligman + +#ifndef CaloG4_EscapedEnergyRegistry_H +#define CaloG4_EscapedEnergyRegistry_H + +// This class keeps track of which types of volumes use which +// VEscapedEnergyProcessing objects. + +// The types volumes are categorized by volume name. I anticipate that +// there will be only two entries in this registry, one for "LAr::" +// and one for "Tile::", but there may be others in the future. + +// In other words, LAr volumes will use one VEscapedEnergyProcessing +// object, and Tile volumes will another VEscapedEnergyProcessing +// object, and this class keeps track of which is which. + +// Since there's only one registry, this class uses the singleton +// pattern. + +#include "CaloG4Sim/VEscapedEnergyProcessing.h" +#include "globals.hh" + +#include <map> + +namespace CaloG4 { + + class EscapedEnergyRegistry + { + public: + + static EscapedEnergyRegistry* GetInstance(); + + ~EscapedEnergyRegistry(); + + // Add a VEscapedEnergyProcessing object to the registry. The + // name include "Adopt" because we assume responsibility for + // deleting the pointer. + void AddAndAdoptProcessing( const G4String& name, VEscapedEnergyProcessing* process ); + + // Given a volume name, + VEscapedEnergyProcessing* GetProcessing( const G4String& volumeName ) const; + + protected: + + // Constructor protected according to singleton pattern. + EscapedEnergyRegistry(); + + private: + + typedef std::map< const G4String, VEscapedEnergyProcessing* > m_processingMap_t; + typedef m_processingMap_t::iterator m_processingMap_ptr_t; + typedef m_processingMap_t::const_iterator m_processingMap_const_ptr_t; + m_processingMap_t m_processingMap; + + }; + +} // namespace CaloG4 + +#endif // CaloG4_EscapedEnergyRegistry_H diff --git a/Calorimeter/CaloG4Sim/CaloG4Sim/SimulationEnergies.h b/Calorimeter/CaloG4Sim/CaloG4Sim/SimulationEnergies.h new file mode 100755 index 0000000000000000000000000000000000000000..269c2018197fce84ca5b902d3914c51960ed391f --- /dev/null +++ b/Calorimeter/CaloG4Sim/CaloG4Sim/SimulationEnergies.h @@ -0,0 +1,133 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// SimulationEnergies.h + +// This class implements the calculations requires to categorize the +// energies deposited during the simulation. This "Monte Carlo truth" +// information is used in calibration studies and other applications. + +// 12-Aug-2009 M. Leltchouk: this code update improves performance. +// Many thanks to Zachary Marshall for introducing and coding or +// suggesting most of these changes and to Gennady Pospelov for +// supporting this code (including careful testing) and the whole +// "machinary of calibration hits" which he greatly improved and extended. + +// 07-Jan-2004 Mikhail Leltchouk with code additions and great support +// from William Glenn Seligman (WGS). +// Columbia University in the City of New York, +// Nevis Labs, 136 South Broadway, Irvington, NY 10533, USA + +// ---------------------------- + +// 29-Jun-2004 WGS: Take advantage of the new features of Geant 4.6.2: +// this class no longer has to be a singleton, and it need no longer +// inherit from G4VSteppingVerbose. + +#ifndef CaloG4_SimulationEnergies_H +#define CaloG4_SimulationEnergies_H + +#include "G4TrackStatus.hh" +#include "G4ParticleDefinition.hh" +#include "G4DynamicParticle.hh" +#include "G4ThreeVector.hh" +#include "globals.hh" + +#include <map> +#include <vector> + +// Forward declarations. +class G4Step; + +namespace CaloG4 { + + class SimulationEnergies + { + public: + + SimulationEnergies(); + virtual ~SimulationEnergies(); + + // The simple method to call from a calibration calculator: + // Examine the G4Step and return the energies required for a + // calibration hit. + + void Energies( const G4Step* , std::vector<G4double> & ); + + // Accessing detailed information: + + // It's useful to define some types here. The following enum uses + // a C++ trick: the first item in the list is set to zero. + // Therefore, "kNumberOfEnergyCategories" will be equal to the + // number of entries in the enum. + + enum eEnergyCategory + { kEm = 0, + kNonEm, + kInvisible0, + /* 12-Aug-2009 we compared three different techniques + kInvisible1, // to determine the energy lost + kInvisible2, // due to "invisible" processes in the detector + and chose one of them for production runs + */ + kEscaped, + kNumberOfEnergyCategories }; + + // This defines the results returned by the energy classification; + // these detailed results are mostly used for studies. Time is in + // [ns], energies are in [MeV]. + + typedef struct { + std::map<eEnergyCategory, G4double> energy; + } ClassifyResult_t; + + /* 12-Aug-2009 We keep only info about energy in ClassifyResult_t + for production runs (see above). + typedef struct { + G4double time; + std::map<eEnergyCategory, G4double> energy; + G4TrackStatus trackStatus; + const G4ParticleDefinition* particle; + const G4DynamicParticle* dynamicParticle; + G4String processName; + } ClassifyResult_t; + */ + + // This method performs the actual function of this class: It + // classifies the components of the energy deposited by the + // current G4Step. Allow a flag to control whether the escaped + // energy is routed to some other volume that the one in which the + // G4Step occurs. + + ClassifyResult_t Classify( const G4Step*, const G4bool processEscaped = true ); + + // Has this routine been called for this G4Step? + static G4bool StepWasProcessed() { return m_calledForStep; } + static void SetStepProcessed() { m_calledForStep = true; } + static void ResetStepProcessed() { m_calledForStep = false; } + + private: + + // Some private methods for internal calculations: + + G4double m_measurableEnergy(const G4ParticleDefinition *particleDef, + G4int PDGEncoding, + G4double totalEnergy, + G4double kineticEnergy); + + G4double m_measurableEnergyV2(const G4ParticleDefinition *particleDef, + G4int PDGEncoding, + G4double totalEnergy, + G4double kineticEnergy); + + // Escaped energy requires special processing. + G4bool m_ProcessEscapedEnergy( G4ThreeVector point, G4double energy ); + + // Used to keep track of processing state. + static G4bool m_calledForStep; + }; + +} // namespace CaloG4 + +#endif // CaloG4_SimulationEnergies_H diff --git a/Calorimeter/CaloG4Sim/CaloG4Sim/VEscapedEnergyProcessing.h b/Calorimeter/CaloG4Sim/CaloG4Sim/VEscapedEnergyProcessing.h new file mode 100755 index 0000000000000000000000000000000000000000..58d251291bc2cb63e09ea2375b1e77b6810543f4 --- /dev/null +++ b/Calorimeter/CaloG4Sim/CaloG4Sim/VEscapedEnergyProcessing.h @@ -0,0 +1,51 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// VEscapedEnergyProcessing +// 13-Jul-2004 William Seligman + +#ifndef CaloG4_VEscapedEnergyProcessing_H +#define CaloG4_VEscapedEnergyProcessing_H + +// The SimulationEnergies class provides a common procedure for +// categorizing the energy deposited in a given G4Step. However, +// different detectors have different scheme for handling one of the +// categories: escaped energy. + +// The issue is that, if a particle's energy is lost by escaping the +// simulation, you don't want to record that energy in the volume it +// escapes; you want to record that energy in the volume in which the +// particle was created. Neutrinos are a good example of this. + +// In effect, the SimulationEnergies class has to issue an "interrupt" +// to some other volume than the current G4Step, telling that other +// volume to accumulate the energy of the escaped particle. The Tile +// Simulation and the LAr simulation handle this interrupt +// differently, since they organize sensitive detectors in a different +// way. + +// This interface "hides" the different escaped-energy processing +// required by the different detectors in the simulation. + +#include "G4TouchableHandle.hh" +#include "G4ThreeVector.hh" +#include "globals.hh" + +namespace CaloG4 { + + class VEscapedEnergyProcessing + { + public: + virtual ~VEscapedEnergyProcessing() {;} + + // Method: The G4TouchableHandle for the volume in which "point" is + // located; the value of "point" itself in case additional + // processing is necessary, and the amount of escaped energy. + + virtual G4bool Process( G4TouchableHandle& handle, G4ThreeVector& point, G4double energy ) = 0; + }; + +} // namespace CaloG4 + +#endif // CaloG4_VEscapedEnergyProcessing_H diff --git a/Calorimeter/CaloG4Sim/GNUmakefile b/Calorimeter/CaloG4Sim/GNUmakefile new file mode 100755 index 0000000000000000000000000000000000000000..55076c5bb08eccb12d967c1e12ec78ac2b7755d1 --- /dev/null +++ b/Calorimeter/CaloG4Sim/GNUmakefile @@ -0,0 +1,27 @@ +# -------------------------------------------------------------- +# GNUmakefile for CaloG4Sim code +# -------------------------------------------------------------- + +name := CaloG4Sim +G4TARGET := $(name) +G4EXLIB := true + +# Note that we're only compiling a library; this GNUmakefile cannot +# create a binary on its own. + +.PHONY: all +all: lib + +# Standard G4 architecture. +include $(G4INSTALL)/config/architecture.gmk + +CPPFLAGS += $(LARG4CFLAGS) + +# Standard G4 compilation process. +include $(G4INSTALL)/config/binmake.gmk + +# Include a "realclean" to get rid of excess baggage. +.PHONY: realclean +realclean: + @echo "Also deleting old Common code backups..." + @rm -f *~ *.bak CaloG4Sim/*~ CaloG4Sim/*.bak src/*~ src/*.bak diff --git a/Calorimeter/CaloG4Sim/cmt/requirements b/Calorimeter/CaloG4Sim/cmt/requirements new file mode 100755 index 0000000000000000000000000000000000000000..f8598b4653ce149006d5a0a254b0c57915579e97 --- /dev/null +++ b/Calorimeter/CaloG4Sim/cmt/requirements @@ -0,0 +1,23 @@ +package CaloG4Sim + +author William Seligman <seligman@nevis.columbia.edu> + +# This package contains the common Geant4 Simulation code used by +# LAr and Tile. + +use AtlasPolicy AtlasPolicy-* + +# Application Action - defined in FADS. +use FadsActions FadsActions-* Simulation/G4Sim/FADS + +use Geant4 Geant4-* External + +# Build the library (and export the headers) +library CaloG4Sim *.cc +apply_pattern installed_library + +#======================================================= +private + +# Suppress warnings of unused variables. +macro_append CaloG4Sim_cppflags " -Wno-unused" diff --git a/Calorimeter/CaloG4Sim/doc/MainPage.h b/Calorimeter/CaloG4Sim/doc/MainPage.h new file mode 100644 index 0000000000000000000000000000000000000000..2c8d893d8203716b4cb88fc0aae732957e70a307 --- /dev/null +++ b/Calorimeter/CaloG4Sim/doc/MainPage.h @@ -0,0 +1,29 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/** + + +@mainpage + +This package is used by Calorimeter simulation with Calibration Hits. + +The key classes of this package are: + +1. CalibrationDefaultProcessing. This class provides a "default behavior" for all energy deposits that are not made in a volume that's been made sensitive for calibration studies. + +2. EscapedEnergyRegistry. This class keeps track of which types of volumes use which VEscapedEnergyProcessing objects. + +3. SimulationEnergies. This class implements the calculations requires to categorize the energies deposited during the simulation. It's intended for use in calibration studies. + +4. VEscapedEnergyProcessing. This interface "hides" the different escaped-energy processing required by the different detectors in the simulation. + +-------------------------------- + REQUIREMENTS +-------------------------------- + +@include requirements +@htmlinclude used_packages.html + +*/ diff --git a/Calorimeter/CaloG4Sim/src/CalibrationDefaultProcessing.cc b/Calorimeter/CaloG4Sim/src/CalibrationDefaultProcessing.cc new file mode 100755 index 0000000000000000000000000000000000000000..5160d8621633fe46c9f80cad27cc80e97d36f160 --- /dev/null +++ b/Calorimeter/CaloG4Sim/src/CalibrationDefaultProcessing.cc @@ -0,0 +1,101 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// CalibrationDefaultProcessing +// 04-Mar-2004 William Seligman + +// The calibration studies rely on every volume in the simulation +// being made into a sensitive detector. There is a practical +// problem: What if we're still in the middle of developing code, and +// not every volume has been made sensitive yet? What if we've +// overlooked a volume? + +// This class provides a "default behavior" for all energy deposits +// that are not made in a volume that's been made sensitive for +// calibration studies. + +// Since we need to examine the particles at every G4Step, this class +// makes use of a UserAction interface class. + +#include "CaloG4Sim/CalibrationDefaultProcessing.h" +#include "CaloG4Sim/SimulationEnergies.h" + +#include "G4VSensitiveDetector.hh" +#include "G4Step.hh" + +#include "FadsActions/FadsSteppingAction.h" + +using namespace FADS; + +namespace CaloG4 { + + G4VSensitiveDetector* CalibrationDefaultProcessing::m_defaultSD = 0; + + // Register this class with the UserAction class managers. + CalibrationDefaultProcessing::CalibrationDefaultProcessing() + : UserAction("CaloG4Sim:::CalibrationDefaultProcessing") + { + FadsSteppingAction::GetSteppingAction()->RegisterAction(this); + } + + + CalibrationDefaultProcessing::~CalibrationDefaultProcessing() {;} + + + void CalibrationDefaultProcessing::SteppingAction(const G4Step* a_step) + { + // Do we have a sensitive detector? + if ( m_defaultSD != 0 ) + { + // We only want to perform the default processing if no other + // calibration processing has occurred for this step. + + if ( ! CaloG4::SimulationEnergies::StepWasProcessed() ) + { + // We haven't performed any calibration processing for this + // step (probably there is no sensitive detector for the + // volume). Use the default sensitive detector to process + // this step. Note that we have to "cast away" const-ness for + // the G4Step*, due to how G4VSensitiveDetector::Hit() is + // defined. + m_defaultSD->Hit( const_cast<G4Step*>(a_step) ); + } + + // In any case, clear the flag for the next step. + CaloG4::SimulationEnergies::ResetStepProcessed(); + } + else + { + static G4bool warningPrinted = false; + if ( ! warningPrinted ) + { + warningPrinted = true; + G4cout << "CaloG4::CalibrationDefaultProcessing::SteppingAction - " + << G4endl + << " A default calibration sensitive detector was not defined." + << G4endl + << " Not all calibration energies will be recorded." + << G4endl; + } + } + } + + void CalibrationDefaultProcessing::SetDefaultSD( G4VSensitiveDetector* a_sd ) + { + if ( m_defaultSD != 0 ) + { + G4cout << "CaloG4::CalibrationDefaultProcessing::SetDefaultSD - " + << G4endl + << " The default calibration sensitive detector '" + << m_defaultSD->GetName() + << "'" << G4endl + << " is being replace by sensitive detector '" + << a_sd->GetName() + << "'" << G4endl; + } + + m_defaultSD = a_sd; + } + +} // namespace CaloG4 diff --git a/Calorimeter/CaloG4Sim/src/EscapedEnergyRegistry.cc b/Calorimeter/CaloG4Sim/src/EscapedEnergyRegistry.cc new file mode 100755 index 0000000000000000000000000000000000000000..12150d63804bfac3a04736aa401ea721f1d987af --- /dev/null +++ b/Calorimeter/CaloG4Sim/src/EscapedEnergyRegistry.cc @@ -0,0 +1,103 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// EscapedEnergyRegistry +// 15-Jul-2004 William Seligman + +#include "CaloG4Sim/EscapedEnergyRegistry.h" +#include "CaloG4Sim/VEscapedEnergyProcessing.h" + +#include "G4ios.hh" +#include "globals.hh" + +#include <map> + +namespace CaloG4 { + + // Standard implementation of a singleton pattern. + + EscapedEnergyRegistry* EscapedEnergyRegistry::GetInstance() + { + static EscapedEnergyRegistry instance; + return &instance; + } + + EscapedEnergyRegistry::EscapedEnergyRegistry() + {;} + + EscapedEnergyRegistry::~EscapedEnergyRegistry() + { + // Delete all the pointers we've adopted. + m_processingMap_ptr_t i; + for ( i = m_processingMap.begin(); i != m_processingMap.end(); i++ ) + { + delete (*i).second; + } + } + + void EscapedEnergyRegistry::AddAndAdoptProcessing( const G4String& name, VEscapedEnergyProcessing* process ) + { + // Don't bother adding a null pointer. + if ( process == 0 ) return; + + // Check that we're not adding any duplicates. + m_processingMap_ptr_t i; + for ( i = m_processingMap.begin(); i != m_processingMap.end(); i++ ) + { + if ( name == (*i).first ) + { + G4cout << "CaloG4Sim::EscapedEnergyRegistry::AddAndAdoptProcessing -" + << G4endl; + G4cout << " Trying to add a second VEscapedEnergyProcessing with the name '" + << name + << "'" << G4endl; + G4cout << " Entry is rejected!" << G4endl; + return; + } + if ( process == (*i).second ) + { + G4cout << "CaloG4Sim::EscapedEnergyRegistry::AddAndAdoptProcessing -" + << G4endl; + G4cout << " The key '" + << name + << "' has the same VEscapedEnergyProcessing object as the key '" + << (*i).first + << "'" << G4endl; + G4cout << " Entry is rejected!" << G4endl; + return; + } + } + + // There are no duplicates, so add the entry. + m_processingMap[ name ] = process; + } + + VEscapedEnergyProcessing* EscapedEnergyRegistry::GetProcessing( const G4String& volumeName ) const + { + // Search through the map. If we find an entry whose text string + // is a substring of the volume name (e.g., "LAr::" is a substring + // of "LAr::BarrelCryostat::MotherVolume"), then return that + // VEscapedEnergyProcessing object. + + // Reminder: + // m_processingMap = consists of pair< G4String, VEscapedEnergyProcessing* > + // i = iterator (pointer) to a m_processingMap entry. + // (*i) = pair< G4String, VEscapedEnergyProcessing* > + // (*i).first = a G4String + // (*i).second = a VEscapedEnergyProcessing* + + m_processingMap_const_ptr_t i; + for ( i = m_processingMap.begin(); i != m_processingMap.end(); i++ ) + { + if ( volumeName.contains( (*i).first ) ) + return (*i).second; + } + + // If we get here, then there was no entry in the map that + // matched any portion of the volume name. + + return 0; + } + +} // namespace CaloG4 diff --git a/Calorimeter/CaloG4Sim/src/SimulationEnergies.cc b/Calorimeter/CaloG4Sim/src/SimulationEnergies.cc new file mode 100755 index 0000000000000000000000000000000000000000..9efc867d8bfeb2e3929cdf1ffadb5e9bedbb476a --- /dev/null +++ b/Calorimeter/CaloG4Sim/src/SimulationEnergies.cc @@ -0,0 +1,626 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// SimulationEnergies.cc + +// This class implements the calculations requires to categorize the +// energies deposited during the simulation. This "Monte Carlo truth" +// information is used in calibration studies and other applications. + +// 12-Aug-2009 M. Leltchouk: this code update improves performance. +// Many thanks to Zachary Marshall for introducing and coding or +// suggesting most of these changes and to Gennady Pospelov and +// Vakhtang Tsulaia for supporting this code and the whole +// "machinary of calibration hits" which he greatly improved and extended. + +// 07-Jan-2004 Mikhail Leltchouk with code additions and great support +// from William Glenn Seligman (WGS). +// Columbia University in the City of New York, +// Nevis Labs, 136 South Broadway, Irvington, NY 10533, USA + +// ---------------------------- + +// 15-Jul-2004 WGS: Now use the VEscapedEnergyProcessing interface and +// EscapedEnergyRegistry from package Calorimeter/CaloG4Sim. + +// 28-Nov-2005 M. Leltchouk: protection against segfaults caused by +// geometry problems when touchableHandle->GetVolume()==0 + +// 20-Apr-2006 M. Leltchouk: internal particle mass table is used now +// in SimulationEnergies::m_measurableEnergy because in recently released +// G4 8.0 particle masses may not be available when SimulationEnergies +// constructor is looking for them, see Andrea Dell'Acqua's comment below. + +// 24-Apr-2007 M. Leltchouk: Particle Data Group (PDG) encoding for nuclei +// has been introduced in SimulationEnergies.cc and +// tagged CaloG4Sim-00-02-26 to be consistent with update in Geant4 8.2. + +// Particle Data Group (PDG) encoding for nuclei ( described in +// http://pdg.lbl.gov/2006/mcdata/mc_particle_id_contents.html ) +// was introduced in Geant4 8.2 and in CaloG4Sim-00-02-26 which +// should be used for ATLAS simulations starting from release 13.0.0. + +// CaloG4Sim-00-02-25 should be used for earlier releases where +// PDGEncoding == 0 for nuclei. + +// 18-Feb-2008 M. Leltchouk: +// Since G4AtlasApps-00-02-45 the neutrinos are killed in the 'ATLAS' +// geometries (not in Combined Test Beam geometries) by +// Simulation/G4Atlas/G4AtlasAlg/AthenaStackingAction::ClassifyNewTrack() + +// To handle both cases: with or without neutrino cut, i.e. with +// from G4Svc.G4SvcConf import G4Svc +// G4Svc.KillAllNeutrinos=True (or G4Svc.KillAllNeutrinos=False) +// the SimulationEnergies.cc code has been changed: + +// 1) Now neutrino energies are accumulated in result.energy[kEscaped] +// as soon as neutrinos appear in the list of secondaries. +// 2) These escaped energies are recorded to the cell where the escaping +// track originates (i.e. where neutrinos have been born as a result of +// some particle decay) without call of +// m_ProcessEscapedEnergy( thisTrackVertex, escapedEnergy ). +// 3) If a neutrino is tracked (was not killed) then the special early +// return from SimulationEnergies::Classify is used to avoid the second +// counting of the same neutrino energy when this neutrino escapes from +// World volume. + +#undef DEBUG_ENERGIES + +#include "CaloG4Sim/SimulationEnergies.h" +#include "CaloG4Sim/VEscapedEnergyProcessing.h" +#include "CaloG4Sim/EscapedEnergyRegistry.h" + +#include "G4EventManager.hh" +#include "G4SteppingManager.hh" +#include "G4TrackVector.hh" +#include "G4StepPoint.hh" +#include "G4Step.hh" +#include "G4Track.hh" +#include "G4TrackStatus.hh" +#include "G4ParticleDefinition.hh" +#include "G4ParticleTable.hh" +#include "G4ThreeVector.hh" +#include "G4VProcess.hh" +// #include "G4EmProcessSubType.hh" + +// Particle definitions +#include "G4Electron.hh" +#include "G4Positron.hh" +#include "G4Neutron.hh" +#include "G4Proton.hh" +#include "G4AntiNeutron.hh" +#include "G4AntiProton.hh" +#include "G4NeutrinoE.hh" +#include "G4AntiNeutrinoE.hh" +#include "G4NeutrinoMu.hh" +#include "G4AntiNeutrinoMu.hh" +#include "G4NeutrinoTau.hh" +#include "G4AntiNeutrinoTau.hh" +#include "G4Lambda.hh" +#include "G4SigmaPlus.hh" +#include "G4XiMinus.hh" +#include "G4SigmaZero.hh" +#include "G4XiZero.hh" +#include "G4OmegaMinus.hh" +#include "G4SigmaMinus.hh" +#include "G4AntiLambda.hh" +#include "G4AntiSigmaPlus.hh" +#include "G4AntiSigmaZero.hh" +#include "G4AntiXiZero.hh" +#include "G4AntiXiMinus.hh" +#include "G4AntiOmegaMinus.hh" +#include "G4AntiSigmaMinus.hh" + +// For special escaped energy processing. +#include "G4Navigator.hh" +#include "G4TransportationManager.hh" +#include "G4TouchableHandle.hh" +#include "G4TouchableHistory.hh" + +#include "G4ios.hh" +#include <vector> +#include <string> + + +namespace CaloG4 { + + /* 20-Apr-2006 M. Leltchouk + G4double SimulationEnergies::electronMass = 0; + G4double SimulationEnergies::protonMass = 0; + G4double SimulationEnergies::neutronMass = 0; + G4double SimulationEnergies::deuteronMass = 0; + G4double SimulationEnergies::tritonMass = 0; + G4double SimulationEnergies::alphaMass = 0; + G4double SimulationEnergies::helium3Mass = 0; + */ + + G4bool SimulationEnergies::m_calledForStep = false; + + + SimulationEnergies::SimulationEnergies() + { + // Initialize some static variables. + + static G4bool initialized = false; + if ( ! initialized ) + { + initialized = true; + + /* + // Constructor: initialize some useful constants. + electronMass = G4ParticleTable::GetParticleTable()->FindParticle("e-") ->GetPDGMass(); + protonMass = G4ParticleTable::GetParticleTable()->FindParticle("proton") ->GetPDGMass(); + neutronMass = G4ParticleTable::GetParticleTable()->FindParticle("neutron") ->GetPDGMass(); + deuteronMass = G4ParticleTable::GetParticleTable()->FindParticle("deuteron")->GetPDGMass(); + tritonMass = G4ParticleTable::GetParticleTable()->FindParticle("triton") ->GetPDGMass(); + alphaMass = G4ParticleTable::GetParticleTable()->FindParticle("alpha") ->GetPDGMass(); + helium3Mass = G4ParticleTable::GetParticleTable()->FindParticle("He3") ->GetPDGMass(); + +#ifdef DEBUG_ENERGIES + G4cout << "SimulationEnergies initialization: " << G4endl; + G4cout << ">>>> electronMass="<<electronMass << G4endl; + G4cout << ">>>> protonMass="<<protonMass << G4endl; + G4cout << ">>>> neutronMass="<<neutronMass << G4endl; + G4cout << ">>>> deuteronMass="<<deuteronMass << G4endl; + G4cout << ">>>> tritonMass="<<tritonMass << G4endl; + G4cout << ">>>> alphaMass="<<alphaMass << G4endl; + G4cout << ">>>> helium3Mass="<<helium3Mass << G4endl; +#endif + */ + } + } + + + SimulationEnergies::~SimulationEnergies() + {;} + + + // The "simple" call, intended for calibration calculators: + void SimulationEnergies::Energies( const G4Step* a_step , std::vector<G4double>& energies ) + { + // Call the detailed classification. Process any escaped energy. + ClassifyResult_t category = Classify( a_step, true ); + + // Extract the values we need from the result. Note that it's at + // this point we decide which of the available-energy calculations + // is to be used in a calibration hit. + if(energies.size()!=0) energies.clear(); + + energies.push_back( category.energy[SimulationEnergies::kEm] ); + energies.push_back( category.energy[SimulationEnergies::kNonEm] ); + energies.push_back( category.energy[SimulationEnergies::kInvisible0] ); + energies.push_back( category.energy[SimulationEnergies::kEscaped] ); + + // Note that we've been called for the current G4Step. + SetStepProcessed(); + } + + + // Detailed calculation. + + // -- OUTPUT (energy in MeV, time in nanosecond): + // 1) result.energy[kEm] - visible energy for electromagnetic scale + // 2) result.energy[kNonEm] - visible energy for hadronic scale + // 3) result.energy[kInvisible0] - invisible energy (Version 0) | one of these + // 3a) result.energy[kInvisible1] - invisible energy (Version 1) | 3 versions should + // 3b) result.energy[kInvisible2] - invisible energy (Version 2) | be choosen for a hit + // 4) result.energy[kEscaped] - escaped energy + // 5) result.time - "midstep" time delay compare to the light + // travel from the point (0,0,0) in World + // coordinates to the "midstep" point. + // + // For some studies, it's useful to have this method not process the + // escaped energy. If a_processEscaped==false, then the escaped + // energy will be not be routed to some other volume's hits. + + SimulationEnergies::ClassifyResult_t SimulationEnergies::Classify( const G4Step* step, + const G4bool a_processEscaped ) + { + // Initialize our result to zero. + + ClassifyResult_t result; + for (unsigned int i = 0; i != kNumberOfEnergyCategories; i++ ) + result.energy[ (eEnergyCategory) i ] = 0; + + G4Track* pTrack = step->GetTrack(); + G4ParticleDefinition* particle = pTrack->GetDefinition(); + + // If it is a step of a neutrino tracking: + if (particle == G4NeutrinoE::Definition() || particle == G4AntiNeutrinoE::Definition() || // nu_e, anti_nu_e + particle == G4NeutrinoMu::Definition() || particle == G4AntiNeutrinoMu::Definition() || // nu_mu, anti_nu_mu + particle == G4NeutrinoTau::Definition() || particle == G4AntiNeutrinoTau::Definition()) // nu_tau, anti_nu_tau + { + return result; + } + + G4TrackStatus status = pTrack->GetTrackStatus(); + G4double dEStepVisible = step->GetTotalEnergyDeposit(); + G4int processSubTypeValue=0; + if ( step->GetPostStepPoint()->GetProcessDefinedStep() != 0 ) + processSubTypeValue = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType(); //from G4VProcess.hh + + G4double EkinPreStep = step->GetPreStepPoint()->GetKineticEnergy(); + const G4DynamicParticle* dynParticle = pTrack->GetDynamicParticle(); + G4double EkinPostStep = pTrack->GetKineticEnergy(); + + // Start of calculation of invisible energy for current step. + if ( status == fStopAndKill ){ // StopAndKill stepping particle at PostStep + G4double incomingEkin = EkinPreStep - dEStepVisible; + G4double incomingEtot = dynParticle->GetMass() + incomingEkin; + + result.energy[kInvisible0] = + m_measurableEnergy(particle, + particle->GetPDGEncoding(), + incomingEtot, + incomingEkin); + } + else if (status == fAlive || status == fStopButAlive){// Alive stepping particle at PostStep + result.energy[kInvisible0] = EkinPreStep - EkinPostStep - dEStepVisible; + } + + G4SteppingManager* steppingManager = G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager(); + + // Copy internal variables from the G4SteppingManager. + G4TrackVector* fSecondary = steppingManager->GetfSecondary(); + G4int fN2ndariesAtRestDoIt = steppingManager->GetfN2ndariesAtRestDoIt(); + G4int fN2ndariesAlongStepDoIt = steppingManager->GetfN2ndariesAlongStepDoIt(); + G4int fN2ndariesPostStepDoIt = steppingManager->GetfN2ndariesPostStepDoIt(); + + G4int tN2ndariesTot = fN2ndariesAtRestDoIt + + fN2ndariesAlongStepDoIt + + fN2ndariesPostStepDoIt; + + // loop through secondary particles which were added at current step + // to the list of all secondaries of current track: + G4int loopStart = (*fSecondary).size() - tN2ndariesTot; + size_t loopEnd = (*fSecondary).size(); + if (loopStart < 0) { + loopEnd = loopEnd - loopStart; + loopStart = 0; + } + + G4ParticleDefinition *secondaryID; + G4double totalEofSecondary=0, kinEofSecondary=0, measurEofSecondary=0; + for(size_t lp1=loopStart; lp1<loopEnd; lp1++){ + + //----- extract information about each new secondary particle: + secondaryID = (*fSecondary)[lp1]->GetDefinition(); + totalEofSecondary = (*fSecondary)[lp1]->GetTotalEnergy(); + + //----- use this information: + if (secondaryID == G4NeutrinoE::Definition() || secondaryID == G4AntiNeutrinoE::Definition() || // nu_e, anti_nu_e + secondaryID == G4NeutrinoMu::Definition() || secondaryID == G4AntiNeutrinoMu::Definition() || // nu_mu, anti_nu_mu + secondaryID == G4NeutrinoTau::Definition() || secondaryID == G4AntiNeutrinoTau::Definition()) // nu_tau, anti_nu_tau + { + result.energy[kInvisible0] -= totalEofSecondary; + result.energy[kEscaped] += totalEofSecondary; + } + else { + //----- extract further information about each new secondary particle: + kinEofSecondary = (*fSecondary)[lp1]->GetKineticEnergy(); + measurEofSecondary = m_measurableEnergy(secondaryID, + secondaryID->GetPDGEncoding(), + totalEofSecondary, + kinEofSecondary); + result.energy[kInvisible0] -= measurEofSecondary; + } + } + +#ifdef DEBUG_ENERGIES + // For debugging: + G4bool allOK = true; +#endif + + if ( pTrack->GetNextVolume() == 0 ) + { + // this particle escaped World volume at this step + if ( status != fStopAndKill ) + { // checking for abnormal case + G4cerr << "SimulationEnergies::Classify particle " + << particle->GetParticleName() + << " escaped World volume with status="<<status<< G4endl; + } + + G4double escapedEnergy = + m_measurableEnergy(particle, + particle->GetPDGEncoding(), + dynParticle->GetTotalEnergy(), + EkinPostStep); + + result.energy[kInvisible0] -= escapedEnergy; + + if ( a_processEscaped ) +#ifdef DEBUG_ENERGIES + allOK = +#endif + m_ProcessEscapedEnergy( pTrack->GetVertexPosition(), escapedEnergy ); + else + result.energy[kEscaped] = escapedEnergy; + } + + // END of calculation of Invisible and Escaped energy for current step + + // Subdivision of visible energy for current step + if (particle == G4Electron::Definition() || + particle == G4Positron::Definition() || + processSubTypeValue == 12) + { + result.energy[kEm] = dEStepVisible; + } + else { + result.energy[kNonEm] = dEStepVisible; + } + +#ifdef DEBUG_ENERGIES + if ( ! allOK ) + { + std::cout << "SimulationEnergies::Classify - additional info: " + << std::endl + << "particle name=" << particle->GetParticleName() + << " volume=" << step->GetPreStepPoint()->GetPhysicalVolume()->GetName() + << " status=" << status + // 12-Aug-2009 << " process=" << processDefinedStepName + << " process SubType=" << processSubTypeValue + << std::endl; + } +#endif + + return result; + } + + + + G4double SimulationEnergies::m_measurableEnergyV2(const G4ParticleDefinition *particleDef, + G4int PDGEncoding, + G4double totalEnergy, + G4double kineticEnergy) + + // 15-Dec-2003 Mikhail Leltchouk: inspired by FORTRAN Function PrMeasE + // used by Michael Kuhlen and Peter Loch with Geant3 since 1991. + + { + G4double measurableEnergy = totalEnergy; + + if (particleDef == G4Electron::Definition() || particleDef == G4Proton::Definition() || + particleDef == G4Neutron::Definition() || PDGEncoding == 1000010020 || + PDGEncoding == 1000010030 || PDGEncoding == 1000020040 || + PDGEncoding == 1000020030 ){ + measurableEnergy = kineticEnergy; + } + else if (particleDef == G4Positron::Definition() ){ + measurableEnergy = 2.* totalEnergy - kineticEnergy; + } + else if (PDGEncoding > 1000010019 ){ //for nuclei + measurableEnergy = kineticEnergy; + } + else if (particleDef == G4Lambda::Definition() || particleDef == G4SigmaPlus::Definition() || + particleDef == G4SigmaZero::Definition() || particleDef == G4SigmaMinus::Definition() || + particleDef == G4XiMinus::Definition() || particleDef == G4XiZero::Definition() || + particleDef == G4OmegaMinus::Definition() ){ + measurableEnergy = kineticEnergy; + } + + else if (particleDef == G4AntiNeutron::Definition() ){ + measurableEnergy = 2.* totalEnergy - kineticEnergy; + } + else if (particleDef == G4AntiProton::Definition() || + particleDef == G4AntiLambda::Definition() || particleDef == G4AntiSigmaPlus::Definition() || + particleDef == G4AntiSigmaZero::Definition() || particleDef == G4AntiSigmaMinus::Definition() || // Need AntiSigmacPlus and Zero as well? + particleDef == G4AntiXiZero::Definition() || particleDef == G4AntiXiMinus::Definition() || // Need AntiXicMinus and Zero as well? + particleDef == G4AntiOmegaMinus::Definition() ){ + measurableEnergy = 2.* totalEnergy - kineticEnergy; + } + + if (measurableEnergy < -0.0001){ + G4cerr << "Calibration Hit: NEGATIVE measurableEnergyV2="<<measurableEnergy + << " < -0.0001 MeV for "<<particleDef->GetParticleName()<<G4endl; + measurableEnergy = 0.; + } + + //for tests: + //measurableEnergy = kineticEnergy; + + return measurableEnergy; + } + + + G4double SimulationEnergies::m_measurableEnergy(const G4ParticleDefinition* particleDef, + G4int PDGEncoding, + G4double totalEnergy, + G4double kineticEnergy) + + // 5-Dec-2003 Mikhail Leltchouk: extended version of FORTRAN Function PrMeasE + // used by Michael Kuhlen and Peter Loch with Geant3 since 1991. + + // Purpose: returns energy measurable in a calorimeter. + + // Rest masses of stable particles do not give rise to a signal in + // the calorimeter. The measurable energy is obtained by subtracting + // the rest mass from the total energy for these particles and those + // particles which may decay into them. For antiparticles the rest + // mass has to be added due to annihilation processes. + + // Input: particleDef - Geant4 particle definition + // PDGEncoding - Particle Data Group (PDG) particle number + // totalEnergy - particle total energy + // kineticEnergy - particle kinetic energy + + // Output: m_measurableEnergy - energy measurable in a calorimeter + + { + + // Constructor: initialize some useful constants. + // static const G4double electronMass = G4Electron::Electron()->GetPDGMass(); + //static const G4double protonMass = G4Proton::Proton()->GetPDGMass(); + //static const G4double neutronMass = G4Neutron::Neutron()->GetPDGMass(); + + const G4double electronMass = G4Electron::Definition()->GetPDGMass(); + const G4double protonMass = G4Proton::Definition()->GetPDGMass(); + const G4double neutronMass = G4Neutron::Definition()->GetPDGMass(); + +#ifdef DEBUG_ENERGIES + G4cout << "SimulationEnergies initialization: " << G4endl; + G4cout << ">>>> electronMass="<<electronMass << G4endl; + G4cout << ">>>> protonMass="<<protonMass << G4endl; + G4cout << ">>>> neutronMass="<<neutronMass << G4endl; +#endif + + G4double correctionForMass = 0.; + if (particleDef == G4Electron::Definition()){ + correctionForMass = - electronMass; + } + else if (particleDef == G4Positron::Definition()){ + correctionForMass = + electronMass; + } + else if (particleDef == G4Proton::Definition()){ + correctionForMass = - protonMass; + } + else if (particleDef == G4Neutron::Definition()){ + correctionForMass = - neutronMass; + } + else if (PDGEncoding > 1000010019 ){ //for nuclei + return kineticEnergy; + } + else if (particleDef == G4Lambda::Definition() || particleDef == G4SigmaPlus::Definition() || + particleDef == G4SigmaZero::Definition() || particleDef == G4XiMinus::Definition() || + particleDef == G4XiZero::Definition() || particleDef == G4OmegaMinus::Definition() ){ + correctionForMass = - protonMass; + } + else if (particleDef == G4SigmaMinus::Definition() ){ // Need Sigma Minus? + correctionForMass = - neutronMass; + } + + else if (particleDef == G4AntiNeutron::Definition() ){ + correctionForMass = + neutronMass; + } + else if (particleDef == G4AntiProton::Definition() || particleDef == G4AntiLambda::Definition() || + particleDef == G4AntiSigmaZero::Definition() || particleDef == G4AntiSigmaPlus::Definition() || // Need AntiSigmacPlus and Zero as well? + particleDef == G4AntiXiZero::Definition() || particleDef == G4AntiXiMinus::Definition() || // Need AntiXicMinus and Zero as well? + particleDef == G4AntiOmegaMinus::Definition() ){ + correctionForMass = + protonMass; + } + else if (particleDef == G4AntiSigmaMinus::Definition() ){ + correctionForMass = + neutronMass; + } + + G4double measurableEnergy = totalEnergy + correctionForMass; + + if (measurableEnergy < -0.0001){ + G4cerr << "Calibration Hit: NEGATIVE measurableEnergy=" + <<measurableEnergy<<" < -0.0001 MeV for "<<particleDef->GetParticleName()<<G4endl; + measurableEnergy = 0.; + } + return measurableEnergy; + } + + + + G4bool SimulationEnergies::m_ProcessEscapedEnergy( G4ThreeVector a_point, G4double a_energy ) + { + // Escaped energy requires special processing. The energy was not + // deposited in the current G4Step, but at the beginning of the + // particle's track. For example, if a neutrino escapes the + // detector, the energy should be recorded at the point where the + // neutrino was created, not at the point where it escaped the + // detector. + + // Different detectors have different procedures for handling + // escaped energy. We use the EscapedEnergyRegistry class to + // route the request for special processing to the appropriate + // procedure. + + // The first time this routine is called, we have to set up some + // Geant4 navigation pointers. + + static G4Navigator* navigator = 0; + static G4TouchableHandle touchableHandle; + + // Locate the G4TouchableHandle associated with the volume + // containing "a_point". + + if ( navigator == 0 ) + { + // Get the navigator object used by the G4TransportationManager. + navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + // Initialize the touchableHandle object. + touchableHandle = new G4TouchableHistory(); + navigator-> + LocateGlobalPointAndUpdateTouchable(a_point, + touchableHandle(), + false); + } + else + { + navigator-> + LocateGlobalPointAndUpdateTouchable(a_point, + touchableHandle(), + false); + } + + // Choose which VEscapedEnergyProcessing routine we'll use based + // on the volume name. + + CaloG4::EscapedEnergyRegistry* registry = CaloG4::EscapedEnergyRegistry::GetInstance(); + CaloG4::VEscapedEnergyProcessing* eep; + + if (touchableHandle->GetHistoryDepth()) { + // normal case: volume name exists + G4String volumeName = touchableHandle->GetVolume()->GetName(); + eep = registry->GetProcessing( volumeName ); + + if ( eep != 0 ) + // Call the appropriate routine. + return eep->Process( touchableHandle, a_point, a_energy ); + + // If we reach here, we're in a volume that has no escaped energy + // procedure (probably neither Tile nor LAr). + + // I don't know about Tile, but in LAr there is a default + // procedure for non-sensitive volumes. Let's use that (for now). + + eep = registry->GetProcessing( "LAr::" ); + if ( eep != 0 ) + return eep->Process( touchableHandle, a_point, a_energy ); + + // If we get here, the registry was never initialized for LAr. + static G4bool errorDisplayed = false; + if ( ! errorDisplayed ) + { + errorDisplayed = true; + G4cout << "SimulationEnergies::m_ProcessEscapedEnergy - " << G4endl + << " WARNING! CaloG4::EscapedEnergyRegistry was never initialized for 'LArG4::'" << G4endl + << " and LArG4Sim is the package with the code that handles CalibrationHits" << G4endl + << " in non-sensitive volumes. Not all energies deposited in this simulation" << G4endl + << " will be recorded." << G4endl; + } + + return false; + } + else { + // If we reach here, we're in an area with geometry problem + // There is a default procedure for non-sensitive volumes in LAr. + // Let's use that (for now). + + eep = registry->GetProcessing( "LAr::" ); + if ( eep != 0 ) + return eep->Process( touchableHandle, a_point, a_energy ); + + // If we get here, the registry was never initialized for LAr. + static G4bool errorDisplayed1 = false; + if ( ! errorDisplayed1 ) + { + errorDisplayed1 = true; + G4cout << "SimulationEnergies::m_ProcessEscapedEnergy - " << G4endl + <<" WARNING! touchableHandle->GetVolume()==0 geometry problem ? and also" << G4endl + << " WARNING! CaloG4::EscapedEnergyRegistry was never initialized for 'LArG4::'" << G4endl + << " and LArG4Sim is the package with the code that handles CalibrationHits" << G4endl + << " in non-sensitive volumes. Not all energies deposited in this simulation" << G4endl + << " will be recorded." << G4endl; + } + return false; + } + } + +} // namespace LArG4 +