Fix checkreq warning. (CaloG4Sim-00-02-33)

Calorimeter/CaloG4Sim/CaloG4Sim/CalibrationDefaultProcessing.h

+/*
+  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

Calorimeter/CaloG4Sim/CaloG4Sim/EscapedEnergyRegistry.h

+/*
+  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

Calorimeter/CaloG4Sim/CaloG4Sim/SimulationEnergies.h

+/*
+  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

Calorimeter/CaloG4Sim/CaloG4Sim/VEscapedEnergyProcessing.h

+/*
+  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

Calorimeter/CaloG4Sim/GNUmakefile

+# --------------------------------------------------------------
+# 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

Calorimeter/CaloG4Sim/cmt/requirements

+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"

Calorimeter/CaloG4Sim/doc/MainPage.h

+/*
+  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 
+
+*/

Calorimeter/CaloG4Sim/src/CalibrationDefaultProcessing.cc

+/*
+  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

Calorimeter/CaloG4Sim/src/EscapedEnergyRegistry.cc

+/*
+  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

Calorimeter/CaloG4Sim/src/SimulationEnergies.cc

+/*
+  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
+