diff --git a/CMakeLists.txt b/CMakeLists.txt
index ea11ce3..55f1da0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,9 +1,8 @@
 #----------------------------------------------------------------------------
 # Setup the project
-cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
-project(OpenCMSG4ver1)
+cmake_minimum_required(VERSION 3.16...3.27)
+project(OpenCMSG4)
 
-list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS})
 
 #----------------------------------------------------------------------------
 # Find Geant4 package, activating all available UI and Vis drivers by default
@@ -37,31 +36,11 @@ if(ROOT_FOUND)
 endif()
 include(${ROOT_USE_FILE})
 
-#----------------------------------------------------------------------------
-# Find Pythia6 (optional package)
-#
-find_package(Pythia6 QUIET)
-if(PYTHIA6_FOUND)
-  message(STATUS "G4 Examples: Pythia6 found. --> HepMCEx01 example with Pythia6 enabled.") 
-  add_definitions(-DG4LIB_USE_PYTHIA) 
-else()
-  set(PYTHIA6_LIBRARIES "")  
-endif()
-
-#----------------------------------------------------------------------------
-# GCC COMPILER FLAGS
-#
-SET(CMAKE_CXX_STANDARD 14)
-SET(CMAKE_CXX_STANDARD_REQUIRED ON)
-
-SET(GCC_COMPILE_FLAGS "-Wshadow")
-SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${GCC_COMPILE_FLAGS}")
-
 #----------------------------------------------------------------------------
 # Locate sources and headers for this project
 # NB: headers are included so they will show up in IDEs
 #
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include 
+include_directories(${PROJECT_SOURCE_DIR}/include 
                     ${Geant4_INCLUDE_DIR}
 					${HEPMC_INCLUDE_DIR})
 file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc)
@@ -72,8 +51,9 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
 #
 add_executable(OpenCMSG4 OpenCMSG4.cc ${sources} ${headers})
 target_link_libraries(OpenCMSG4 ${Geant4_LIBRARIES}
-						${HEPMC_LIBRARIES} ${HEPMC_FIO_LIBRARIES}
-						${PYTHIA6_LIBRARIES} ${ROOT_LIBRARIES} gfortran)	
+						${HEPMC_LIBRARIES} 
+						${ROOT_LIBRARIES})	
+
 
 #----------------------------------------------------------------------------
 # Copy all scripts to the build directory, i.e. the directory in which we
@@ -110,7 +90,7 @@ endforeach()
 # Add program to the project targets
 # (this avoids the need of typing the program name after make)
 #
-add_custom_target(OpenCMSG4ver1 DEPENDS OpenCMSG4)
+# add_custom_target(OpenCMSG4ver1 DEPENDS OpenCMSG4)
 
 #----------------------------------------------------------------------------
 # Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
diff --git a/GNUmakefile b/GNUmakefile
index 79a9740..5dedcd4 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -8,13 +8,32 @@ G4TARGET := $(name)
 G4EXLIB := true
 
 ifndef G4INSTALL
-  G4INSTALL = ../../..
+  G4INSTALL = ../../../../..
 endif
 
+
 .PHONY: all
-all: lib bin
 
-include $(G4INSTALL)/config/binmake.gmk
+
+ifdef HEPMC_DIR
+all : lib bin
+
+  include $(G4INSTALL)/config/binmake.gmk
+
+  INCFLAGS  += -I$(HEPMC_DIR)/include
+
+  LDLIBS1  += -L$(HEPMC_DIR)/lib -lHepMC
+
+FCFLAGS += -c
+
+
+else
+
+all:
+	@echo 'ERROR - HEPMC_DIR not defined in the environment !'
+	@echo '        Requires HepMC release 1.27.'
+endif
+
 
 visclean:
 	rm -f g4*.prim g4*.eps g4*.wrl
diff --git a/exampleRun_hepmc.mac b/exampleRun_hepmc.mac
index 0f818da..ef823d1 100644
--- a/exampleRun_hepmc.mac
+++ b/exampleRun_hepmc.mac
@@ -1,7 +1,9 @@
 #/vis/ogl/set/displayListLimit 10000000000
 #/vis/viewer/set/viewpointThetaPhi 65 315 deg
 
-/detector/sensitiveMaterialOnly false
+/run/numberOfThreads 1
+
+#/detector/sensitiveMaterialOnly false
 /detector/trackerMode 111111
 /detector/ecalMode 111
 
@@ -10,6 +12,6 @@
 /generator hepmc
 /hepmcAscii/verbose 0
 /hepmcAscii/open example_MyPythia.dat
-/analysis/setFileName OpenCMSG4_hepmc.root
+/root/setFileName OpenCMSG4_particle.root
 
 /run/beamOn 10
diff --git a/exampleRun_particle.mac b/exampleRun_particle.mac
index 44a34c0..3e87dfe 100644
--- a/exampleRun_particle.mac
+++ b/exampleRun_particle.mac
@@ -1,6 +1,8 @@
 #/vis/ogl/set/displayListLimit 10000000000
 #/vis/viewer/set/viewpointThetaPhi 65 315 deg
 
+/run/numberOfThreads 1
+
 /detector/trackerSensitiveMaterialOnly false
 /detector/trackerMode 111111
 /detector/ecalMode 111
@@ -20,4 +22,4 @@
 /generator/particle/sigmaEta 0
 /root/setFileName OpenCMSG4_particle.root
 
-/run/beamOn 100
+/run/beamOn 10
diff --git a/include/Analysis.hh b/include/Analysis.hh
index 2f40a79..3ec94e2 100644
--- a/include/Analysis.hh
+++ b/include/Analysis.hh
@@ -1,7 +1,8 @@
 #ifndef Analysis_h
 #define Analysis_h 1
 
-#include "g4root.hh"
+#include "G4AnalysisManager.hh"
+// #include "g4root.hh"
 //#include "g4xml.hh"
 //#include "g4csv.hh"
 
diff --git a/include/CustomParticleFactory.hh b/include/CustomParticleFactory.hh
index ed6651e..fb39d6c 100755
--- a/include/CustomParticleFactory.hh
+++ b/include/CustomParticleFactory.hh
@@ -12,8 +12,8 @@
 class CustomParticleFactory 
 {
  private:
-  static bool loaded;
-  static std::set<G4ParticleDefinition *> m_particles;
+  static thread_local bool loaded;
+  static thread_local std::set<G4ParticleDefinition *> m_particles;
 
  public:
   static void loadCustomParticles();
diff --git a/include/FullModelHadronicProcess.hh b/include/FullModelHadronicProcess.hh
index 5db3ed9..26c9938 100755
--- a/include/FullModelHadronicProcess.hh
+++ b/include/FullModelHadronicProcess.hh
@@ -41,7 +41,7 @@ private:
 
   G4double GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*);
   
-  void CalculateMomenta( G4FastVector<G4ReactionProduct,256> &vec,
+  void CalculateMomenta( std::vector<G4ReactionProduct*> &vec,
 			 G4int &vecLen,
 			 const G4HadProjectile *originalIncident,
 			 const G4DynamicParticle *originalTarget,
@@ -58,7 +58,7 @@ private:
 				    const G4ReactionProduct &targetParticle,
 				    G4ReactionProduct &leadParticle );
     
-  void Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen);
+  void Rotate(std::vector<G4ReactionProduct*> &vec, G4int &vecLen);
 
   const G4DynamicParticle* FindRhadron(G4ParticleChange*);
 
diff --git a/include/HepMCAsciiReader.hh b/include/HepMCAsciiReader.hh
index 6fbf31d..279b906 100644
--- a/include/HepMCAsciiReader.hh
+++ b/include/HepMCAsciiReader.hh
@@ -1,51 +1,44 @@
-#ifndef _HEPMC_ASCII_READER_H
-#define _HEPMC_ASCII_READER_H
+#ifndef HEPMC_ASCII_READER_H
+#define HEPMC_ASCII_READER_H
 
-#include "G4VPrimaryGenerator.hh"
 #include "HepMC/IO_GenEvent.h"
+#include "HepMCG4Interface.hh"
 
 class HepMCAsciiReaderMessenger;
 
-class HepMCAsciiReader : public G4VPrimaryGenerator
+class HepMCAsciiReader : public HepMCG4Interface
 {
-	public:
-		HepMCAsciiReader();
-		~HepMCAsciiReader();
+  protected:
+    G4String filename;
+    HepMC::IO_GenEvent* asciiInput;
 
-	protected:
-		G4String filename;
-		G4int verbose;
-		HepMCAsciiReaderMessenger* messenger;
+    G4int verbose;
+    HepMCAsciiReaderMessenger* messenger;
 
-		HepMC::IO_GenEvent* asciiInput;
-		HepMC::GenEvent* hepmcevt;
-		virtual G4bool CheckVertexInsideWorld(const G4ThreeVector& pos) const;
-		virtual HepMC::GenEvent* GenerateHepMCEvent();
-		void HepMC2G4(const HepMC::GenEvent* hepmcevt, G4Event* event);
+    virtual HepMC::GenEvent* GenerateHepMCEvent();
 
-	public:
-		virtual void GeneratePrimaryVertex(G4Event * event);
+  public:
+    HepMCAsciiReader();
+    ~HepMCAsciiReader();
 
-	// set/get methods
-		void SetFileName(G4String name);
-		G4String GetFileName() const;
+    // set/get methods
+    void SetFileName(G4String name);
+    G4String GetFileName() const;
 
-		void SetVerboseLevel(G4int i);
-		G4int GetVerboseLevel() const; 
-
-		// methods...
-		void Initialize();
+    void SetVerboseLevel(G4int i);
+    G4int GetVerboseLevel() const;
 
+    // methods...
+    void Initialize();
 };
 
-
 // ====================================================================
 // inline functions
 // ====================================================================
 
 inline void HepMCAsciiReader::SetFileName(G4String name)
 {
-  filename= name;
+  filename = name;
 }
 
 inline G4String HepMCAsciiReader::GetFileName() const
@@ -55,7 +48,7 @@ inline G4String HepMCAsciiReader::GetFileName() const
 
 inline void HepMCAsciiReader::SetVerboseLevel(G4int i)
 {
-  verbose= i;
+  verbose = i;
 }
 
 inline G4int HepMCAsciiReader::GetVerboseLevel() const
diff --git a/include/HepMCAsciiReaderMessenger.hh b/include/HepMCAsciiReaderMessenger.hh
index 7063077..47e276d 100644
--- a/include/HepMCAsciiReaderMessenger.hh
+++ b/include/HepMCAsciiReaderMessenger.hh
@@ -1,5 +1,5 @@
-#ifndef HEPMC_G4_ASCII_READER_MESSENGER_H
-#define HEPMC_G4_ASCII_READER_MESSENGER_H
+#ifndef HEPMC_ASCII_READER_MESSENGER_H
+#define HEPMC_ASCII_READER_MESSENGER_H
 
 #include "G4UImessenger.hh"
 
@@ -9,21 +9,21 @@ class G4UIcmdWithoutParameter;
 class G4UIcmdWithAString;
 class G4UIcmdWithAnInteger;
 
-class HepMCAsciiReaderMessenger : public G4UImessenger {
-public:
-  HepMCAsciiReaderMessenger(HepMCAsciiReader* agen);
-  ~HepMCAsciiReaderMessenger();
+class HepMCAsciiReaderMessenger : public G4UImessenger
+{
+  public:
+    HepMCAsciiReaderMessenger(HepMCAsciiReader* agen);
+    ~HepMCAsciiReaderMessenger();
 
-  void SetNewValue(G4UIcommand* command, G4String newValues);
-  G4String GetCurrentValue(G4UIcommand* command);
+    void SetNewValue(G4UIcommand* command, G4String newValues);
+    G4String GetCurrentValue(G4UIcommand* command);
 
-private:
-  HepMCAsciiReader* gen;
-
-  G4UIdirectory* dir;
-  G4UIcmdWithAnInteger* verbose;
-  G4UIcmdWithAString* open;
+  private:
+    HepMCAsciiReader* gen;
 
+    G4UIdirectory* dir;
+    G4UIcmdWithAnInteger* verbose;
+    G4UIcmdWithAString* open;
 };
 
 #endif
diff --git a/include/HepMCG4Interface.hh b/include/HepMCG4Interface.hh
new file mode 100644
index 0000000..9473b36
--- /dev/null
+++ b/include/HepMCG4Interface.hh
@@ -0,0 +1,50 @@
+#ifndef HEPMC_G4_INTERFACE_H
+#define HEPMC_G4_INTERFACE_H
+
+#include "HepMC/GenEvent.h"
+
+#include "G4VPrimaryGenerator.hh"
+
+/// A base class for primary generation via HepMC object.
+/// This class is derived from G4VPrimaryGenerator.
+
+class HepMCG4Interface : public G4VPrimaryGenerator
+{
+  protected:
+    // Note that the life of HepMC event object must be handled by users.
+    // In the default implementation, a current HepMC event will be
+    // deleted at GeneratePrimaryVertex() in the next event.
+    HepMC::GenEvent* hepmcEvent;  // (care for single event case only)
+
+    // We  have to take care for the position of primaries because
+    // primary vertices outside the world voulme give rise to G4Execption.
+    // If the default implementation is not adequate, an alternative
+    // can be implemented in your own class.
+    virtual G4bool CheckVertexInsideWorld(const G4ThreeVector& pos) const;
+
+    // service method for conversion from HepMC::GenEvent to G4Event
+    void HepMC2G4(const HepMC::GenEvent* hepmcevt, G4Event* g4event);
+
+    // Implement this method in his/her own concrete class.
+    // An empty event will be created in default.
+    virtual HepMC::GenEvent* GenerateHepMCEvent();
+
+  public:
+    HepMCG4Interface();
+    virtual ~HepMCG4Interface();
+
+    // set/get methods
+    HepMC::GenEvent* GetHepMCGenEvent() const;
+
+    // The default behavior is that a single HepMC event generated by
+    // GenerateHepMCEvent() will be converted to G4Event through HepMC2G4().
+    virtual void GeneratePrimaryVertex(G4Event* anEvent);
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+inline HepMC::GenEvent* HepMCG4Interface::GetHepMCGenEvent() const
+{
+  return hepmcEvent;
+}
+
+#endif
diff --git a/include/MyG4ReactionDynamics.hh b/include/MyG4ReactionDynamics.hh
index 22f1e9c..0443623 100755
--- a/include/MyG4ReactionDynamics.hh
+++ b/include/MyG4ReactionDynamics.hh
@@ -37,7 +37,7 @@
 #include "G4DynamicParticle.hh"
 #include "G4ReactionProduct.hh"
 #include "G4Nucleus.hh"
-#include "G4FastVector.hh"
+// #include "G4FastVector.hh"
 #include "G4HadProjectile.hh"
 
 enum{ MYGHADLISTSIZE=256};
@@ -57,7 +57,7 @@ enum{ MYGHADLISTSIZE=256};
     { return 0.0; }
     
     G4bool GenerateXandPt(                   // derived from GENXPT
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effect included
      const G4HadProjectile *originalIncident,
@@ -70,7 +70,7 @@ enum{ MYGHADLISTSIZE=256};
      G4ReactionProduct &leadingStrangeParticle );
     
     void SuppressChargedPions(
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      const G4ReactionProduct &modifiedOriginal,
      G4ReactionProduct &currentParticle,
@@ -80,7 +80,7 @@ enum{ MYGHADLISTSIZE=256};
      G4bool &targetHasChanged );
       
     G4bool TwoCluster(                       // derived from TWOCLU
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effect included
      const G4HadProjectile *originalIncident,
@@ -93,7 +93,7 @@ enum{ MYGHADLISTSIZE=256};
      G4ReactionProduct &leadingStrangeParticle );
     
     void TwoBody(                         // derived from TWOB
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      G4ReactionProduct &modifiedOriginal,
      const G4DynamicParticle *originalTarget,
@@ -107,11 +107,11 @@ enum{ MYGHADLISTSIZE=256};
     G4double GenerateNBodyEvent(                // derived from PHASP
      const G4double totalEnergy,
      const G4bool constantCrossSection,
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen );
     
     void ProduceStrangeParticlePairs(
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      const G4ReactionProduct &modifiedOriginal,
      const G4DynamicParticle *originalTarget, 
@@ -121,7 +121,7 @@ enum{ MYGHADLISTSIZE=256};
      G4bool &targetHasChanged );
     
     void NuclearReaction(                     // derived from NUCREC
-     G4FastVector<G4ReactionProduct,4> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen,
      const G4HadProjectile *originalIncident,
      const G4Nucleus &aNucleus,
@@ -138,14 +138,14 @@ enum{ MYGHADLISTSIZE=256};
      const G4Nucleus &targetNucleus,
      G4ReactionProduct &currentParticle,
      G4ReactionProduct &targetParticle,
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen );
     
     void Defs1(
      const G4ReactionProduct &modifiedOriginal,
      G4ReactionProduct &currentParticle,
      G4ReactionProduct &targetParticle,
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen );
     
     void AddBlackTrackParticles(
@@ -159,14 +159,14 @@ enum{ MYGHADLISTSIZE=256};
      const G4ReactionProduct &modifiedOriginal,
      G4double spall,
      const G4Nucleus &aNucleus,
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen );
     
     void MomentumCheck(
      const G4ReactionProduct &modifiedOriginal,
      G4ReactionProduct &currentParticle,
      G4ReactionProduct &targetParticle,
-     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+     std::vector<G4ReactionProduct*> &vec,
      G4int &vecLen );
     
     G4double normal();
diff --git a/src/B4CustomPhysics.cc b/src/B4CustomPhysics.cc
index e6b3512..34811df 100755
--- a/src/B4CustomPhysics.cc
+++ b/src/B4CustomPhysics.cc
@@ -69,11 +69,12 @@ void B4CustomPhysics::ConstructExotics()
 #include "FullModelHadronicProcess.hh"
 //#include "ToyModelHadronicProcess.hh"
 
-#include "G4PionMinusInelasticProcess.hh"
+#include "G4HadronPhysicsFTFP_BERT.hh"
+// #include "G4PionMinusInelasticProcess.hh"
 //#include "G4LEPionMinusInelastic.hh"
 //#include "G4HEPionMinusInelastic.hh"
 
-#include "G4PionPlusInelasticProcess.hh"
+// #include "G4PionPlusInelasticProcess.hh"
 //#include "G4LEPionPlusInelastic.hh"
 //#include "G4HEPionPlusInelastic.hh"
 #include "G4DecayTable.hh"
diff --git a/src/B4PhysicsList.cc b/src/B4PhysicsList.cc
index f9ab7ab..0604b81 100755
--- a/src/B4PhysicsList.cc
+++ b/src/B4PhysicsList.cc
@@ -85,7 +85,7 @@ B4PhysicsList::B4PhysicsList() : G4VModularPhysicsList() {
   RegisterPhysics( new G4NeutronTrackingCut());
   
   //Custom Physics
-  RegisterPhysics( new B4CustomPhysics());
+//   RegisterPhysics( new B4CustomPhysics());
 
 }
 
diff --git a/src/CustomParticleFactory.cc b/src/CustomParticleFactory.cc
index 7c401c6..bbdd0ea 100755
--- a/src/CustomParticleFactory.cc
+++ b/src/CustomParticleFactory.cc
@@ -13,8 +13,8 @@
 #include "G4PhaseSpaceDecayChannel.hh"
 #include "G4SystemOfUnits.hh"
 
-bool CustomParticleFactory::loaded = false;
-std::set<G4ParticleDefinition *> CustomParticleFactory::m_particles;
+thread_local bool CustomParticleFactory::loaded = false;
+thread_local std::set<G4ParticleDefinition *> CustomParticleFactory::m_particles;
 
 bool CustomParticleFactory::isCustomParticle(G4ParticleDefinition *particle)
 {
diff --git a/src/FullModelHadronicProcess.cc b/src/FullModelHadronicProcess.cc
index 5d8c3a1..4ec91b4 100755
--- a/src/FullModelHadronicProcess.cc
+++ b/src/FullModelHadronicProcess.cc
@@ -4,7 +4,7 @@
 #include "G4ParticleTable.hh"
 #include "Decay3Body.hh"
 #include "MyG4ReactionDynamics.hh"
-#include "G4HadReentrentException.hh"
+// #include "G4HadReentrentException.hh"
 #include "CustomPDGParser.hh"
 #include "CustomParticle.hh"
 #include "G4GenericMessenger.hh"
@@ -260,9 +260,9 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
   if (!TargetSurvives) targetHasChanged = true; //Ditto here
   G4bool quasiElastic = false;
   if (rp.size()==2) quasiElastic = true; //Oh well...
-  G4FastVector<G4ReactionProduct,256> vec;  // vec will contain the secondary particles
+  std::vector<G4ReactionProduct*> vec;  // vec will contain the secondary particles
   G4int vecLen = 0;
-  vec.Initialize( 0 );
+  vec.clear();
 
   
   
@@ -275,7 +275,8 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
 	G4ReactionProduct* pa = new G4ReactionProduct;
 	pa->SetDefinition( theParticleDefinitions[i] );
 	(G4UniformRand() < 0.5) ? pa->SetSide( -1 ) : pa->SetSide( 1 );
-	vec.SetElement( vecLen++, pa );
+	vec.push_back( pa );
+  vecLen++;
       }
   }
 
@@ -460,7 +461,7 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
 
 
 void FullModelHadronicProcess::CalculateMomenta(
-					       G4FastVector<G4ReactionProduct,256> &vec,
+					       std::vector<G4ReactionProduct*> &vec,
 					       G4int &vecLen,
 					       const G4HadProjectile *originalIncident,   // the original incident particle
 					       const G4DynamicParticle *originalTarget,
@@ -599,10 +600,10 @@ void FullModelHadronicProcess::CalculateMomenta(
 							   targetHasChanged, leadFlag,
 							   leadingStrangeParticle );
 	}
-      catch(G4HadReentrentException aC)
+      catch(...)
 	{
-	  aC.Report(G4cout);
-	  throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
+	  G4Exception("FullModelHadronicProcess", "Hadronic001", FatalException,
+      "Failing to calculate momenta");
 	}
     }
   if( finishedTwoClu )
@@ -666,7 +667,7 @@ G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(
   return lead;
 }
 
-void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
+void FullModelHadronicProcess::Rotate(std::vector<G4ReactionProduct*> &vec, G4int &vecLen)
 {
   G4double rotation = 2.*pi*G4UniformRand();
   cache = rotation;
diff --git a/src/HCalConstruction.cc b/src/HCalConstruction.cc
index 11d2e8f..45accff 100644
--- a/src/HCalConstruction.cc
+++ b/src/HCalConstruction.cc
@@ -5,8 +5,8 @@
 #include "G4PVPlacement.hh"
 
 // HCal Barrel Variables
-constexpr G4int nofHcalBarEta = 2*15;  // 2*17 - default value
-constexpr G4int nofHcalBarPhi = 54;    // 72 - default value
+constexpr G4int nofHcalBarEta = 2*17;  // 2*17 - default value
+constexpr G4int nofHcalBarPhi = 72;    // 72 - default value
 constexpr G4int nofHcalBarR = 18;       // 18 - default value
 constexpr G4int nofHcalBarCells = nofHcalBarEta*nofHcalBarPhi*nofHcalBarR;
 constexpr G4double diffRAbsBar = 55*mm;
diff --git a/src/HepMCAsciiReader.cc b/src/HepMCAsciiReader.cc
index 8273736..6f1afba 100644
--- a/src/HepMCAsciiReader.cc
+++ b/src/HepMCAsciiReader.cc
@@ -1,113 +1,40 @@
 #include "HepMCAsciiReader.hh"
-#include "HepMCAsciiReaderMessenger.hh"
 
-#include "G4RunManager.hh"
-#include "G4LorentzVector.hh"
-#include "G4Event.hh"
-#include "G4PrimaryParticle.hh"
-#include "G4PrimaryVertex.hh"
-#include "G4TransportationManager.hh"
-#include "G4PhysicalConstants.hh"
-#include "G4SystemOfUnits.hh"
+#include "HepMCAsciiReaderMessenger.hh"
 
-#include <iostream>
 #include <fstream>
+#include <iostream>
 
-HepMCAsciiReader::HepMCAsciiReader():
-  filename("xxx.dat"), verbose(0), hepmcevt(0)
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+HepMCAsciiReader::HepMCAsciiReader() : filename("xxx.dat"), verbose(0)
 {
-	asciiInput = new HepMC::IO_GenEvent("example_MyPythia.dat", std::ios::in);
-	messenger= new HepMCAsciiReaderMessenger(this);
+  asciiInput = new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);
+
+  messenger = new HepMCAsciiReaderMessenger(this);
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 HepMCAsciiReader::~HepMCAsciiReader()
 {
-	delete asciiInput;
-	delete hepmcevt;
-	delete messenger;
+  delete asciiInput;
+  delete messenger;
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void HepMCAsciiReader::Initialize()
 {
-	delete asciiInput;
+  delete asciiInput;
 
-	asciiInput= new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);
+  asciiInput = new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 HepMC::GenEvent* HepMCAsciiReader::GenerateHepMCEvent()
 {
-	auto hepmcEvt = asciiInput->read_next_event();
-	if(!hepmcEvt) return 0;
-//	else hepmcEvt->print();
-
-	if(verbose>0) hepmcEvt-> print();
-
-	return hepmcEvt;
-}
-
-G4bool HepMCAsciiReader::CheckVertexInsideWorld(const G4ThreeVector& pos) const
-{
-	auto navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
-	auto world = navigator->GetWorldVolume();
-	auto solid = world->GetLogicalVolume()->GetSolid();
-
-	EInside qinside = solid->Inside(pos);
-	if(qinside != kInside) return false;
-	else return true;
-}
-
-void HepMCAsciiReader::HepMC2G4(const HepMC::GenEvent* hepmcEvt, G4Event* event)
-{
-	for(HepMC::GenEvent::vertex_const_iterator vitr = hepmcEvt->vertices_begin();
-				vitr != hepmcEvt->vertices_end(); ++vitr)
-	{
-		G4bool qvtx = false;
-		for(HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
-				pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
-		{
-			if(!(*pitr)->end_vertex() && (*pitr)->status()==1)
-			{
-				qvtx = true;
-				break;
-			}
-		}
-		if(!qvtx) continue;
-		
-		HepMC::FourVector pos = (*vitr)->position();
-		G4LorentzVector xvtx(pos.x(), pos.y(), pos.z(), pos.t());
-		if (! CheckVertexInsideWorld(xvtx.vect()*mm)) continue;
+  HepMC::GenEvent* evt = asciiInput->read_next_event();
+  if (!evt) return 0;  // no more event
 
-		G4PrimaryVertex* g4vtx = 
-				new G4PrimaryVertex(xvtx.x()*mm, xvtx.y()*mm, xvtx.z()*mm, xvtx.t()*mm/c_light);
-		
-		for (HepMC::GenVertex::particle_iterator vpitr= (*vitr)->particles_begin(HepMC::children);
-						vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr)
-		{
-			if( (*vpitr)->status() != 1 ) continue;
+  if (verbose > 0) evt->print();
 
-			G4int pdgcode= (*vpitr)-> pdg_id();
-			pos= (*vpitr)-> momentum();
-			G4LorentzVector p(pos.px(), pos.py(), pos.pz(), pos.e());
-			G4PrimaryParticle* g4prim = 
-						new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
-//std::cout<<g4prim->GetParticleDefinition()->GetParticleName()<<std::endl;
-			g4vtx-> SetPrimary(g4prim);
-		}
-		event->AddPrimaryVertex(g4vtx);
-	}
+  return evt;
 }
-
-void HepMCAsciiReader::GeneratePrimaryVertex(G4Event * event)
-{
-	delete hepmcevt;
-G4cout<<"Check coming here"<<G4endl;
-	hepmcevt = GenerateHepMCEvent();
-	if(! hepmcevt) 
-	{
-		G4cout << "HepMCInterface: no generated particles. run terminated..."<< G4endl;
-		G4RunManager::GetRunManager()-> AbortRun();
-		return;
-	}
-	HepMC2G4(hepmcevt, event);
-}
-
diff --git a/src/HepMCAsciiReaderMessenger.cc b/src/HepMCAsciiReaderMessenger.cc
index 13785ee..116004a 100644
--- a/src/HepMCAsciiReaderMessenger.cc
+++ b/src/HepMCAsciiReaderMessenger.cc
@@ -1,28 +1,26 @@
-#include "G4UIdirectory.hh"
-#include "G4UIcmdWithoutParameter.hh"
-#include "G4UIcmdWithAString.hh"
-#include "G4UIcmdWithAnInteger.hh"
 #include "HepMCAsciiReaderMessenger.hh"
+
 #include "HepMCAsciiReader.hh"
 
+#include "G4UIcmdWithAString.hh"
+#include "G4UIcmdWithAnInteger.hh"
+#include "G4UIcmdWithoutParameter.hh"
+#include "G4UIdirectory.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-HepMCAsciiReaderMessenger::HepMCAsciiReaderMessenger
-                             (HepMCAsciiReader* agen)
-  : gen(agen)
+HepMCAsciiReaderMessenger::HepMCAsciiReaderMessenger(HepMCAsciiReader* agen) : gen(agen)
 {
-  dir= new G4UIdirectory("/hepmcAscii/");
-  dir-> SetGuidance("Reading HepMC event from an Ascii file");
+  dir = new G4UIdirectory("/hepmcAscii/");
+  dir->SetGuidance("Reading HepMC event from an Ascii file");
 
-  verbose=
-    new G4UIcmdWithAnInteger("/hepmcAscii/verbose", this);
-  verbose-> SetGuidance("Set verbose level");
-  verbose-> SetParameterName("verboseLevel", false, false);
-  verbose-> SetRange("verboseLevel>=0 && verboseLevel<=1");
+  verbose = new G4UIcmdWithAnInteger("/hepmcAscii/verbose", this);
+  verbose->SetGuidance("Set verbose level");
+  verbose->SetParameterName("verboseLevel", false, false);
+  verbose->SetRange("verboseLevel>=0 && verboseLevel<=1");
 
-  open= new G4UIcmdWithAString("/hepmcAscii/open", this);
-  open-> SetGuidance("(re)open data file (HepMC Ascii format)");
-  open-> SetParameterName("input ascii file", true, true);
+  open = new G4UIcmdWithAString("/hepmcAscii/open", this);
+  open->SetGuidance("(re)open data file (HepMC Ascii format)");
+  open->SetParameterName("input ascii file", true, true);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -35,30 +33,29 @@ HepMCAsciiReaderMessenger::~HepMCAsciiReaderMessenger()
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void HepMCAsciiReaderMessenger::SetNewValue(G4UIcommand* command,
-                                              G4String newValues)
+void HepMCAsciiReaderMessenger::SetNewValue(G4UIcommand* command, G4String newValues)
 {
-  if (command==verbose) {
-    int level= verbose-> GetNewIntValue(newValues);
-    gen-> SetVerboseLevel(level);
-  } else if (command==open) {
-    gen-> SetFileName(newValues);
-    G4cout << "HepMC Ascii inputfile: "
-           << gen-> GetFileName() << G4endl;
-    gen-> Initialize();
+  if (command == verbose) {
+    int level = verbose->GetNewIntValue(newValues);
+    gen->SetVerboseLevel(level);
+  }
+  else if (command == open) {
+    gen->SetFileName(newValues);
+    G4cout << "HepMC Ascii inputfile: " << gen->GetFileName() << G4endl;
+    gen->Initialize();
   }
 }
 
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4String HepMCAsciiReaderMessenger::GetCurrentValue(G4UIcommand* command)
 {
   G4String cv;
 
   if (command == verbose) {
-    cv= verbose-> ConvertToString(gen-> GetVerboseLevel());
-  } else  if (command == open) {
-    cv= gen-> GetFileName();
+    cv = verbose->ConvertToString(gen->GetVerboseLevel());
+  }
+  else if (command == open) {
+    cv = gen->GetFileName();
   }
   return cv;
 }
diff --git a/src/HepMCG4Interface.cc b/src/HepMCG4Interface.cc
new file mode 100644
index 0000000..e8e0365
--- /dev/null
+++ b/src/HepMCG4Interface.cc
@@ -0,0 +1,103 @@
+#include "HepMCG4Interface.hh"
+
+#include "G4Event.hh"
+#include "G4LorentzVector.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4PrimaryParticle.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4RunManager.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4TransportationManager.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+HepMCG4Interface::HepMCG4Interface() : hepmcEvent(0) {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+HepMCG4Interface::~HepMCG4Interface()
+{
+  delete hepmcEvent;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool HepMCG4Interface::CheckVertexInsideWorld(const G4ThreeVector& pos) const
+{
+  G4Navigator* navigator =
+    G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
+
+  G4VPhysicalVolume* world = navigator->GetWorldVolume();
+  G4VSolid* solid = world->GetLogicalVolume()->GetSolid();
+  EInside qinside = solid->Inside(pos);
+
+  if (qinside != kInside)
+    return false;
+  else
+    return true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void HepMCG4Interface::HepMC2G4(const HepMC::GenEvent* hepmcevt, G4Event* g4event)
+{
+  for (HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
+       vitr != hepmcevt->vertices_end(); ++vitr)
+  {  // loop for vertex ...
+
+    // real vertex?
+    G4bool qvtx = false;
+    for (HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
+         pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
+    {
+      if (!(*pitr)->end_vertex() && (*pitr)->status() == 1) {
+        qvtx = true;
+        break;
+      }
+    }
+    if (!qvtx) continue;
+
+    // check world boundary
+    HepMC::FourVector pos = (*vitr)->position();
+    G4LorentzVector xvtx(pos.x(), pos.y(), pos.z(), pos.t());
+    if (!CheckVertexInsideWorld(xvtx.vect() * mm)) continue;
+
+    // create G4PrimaryVertex and associated G4PrimaryParticles
+    G4PrimaryVertex* g4vtx =
+      new G4PrimaryVertex(xvtx.x() * mm, xvtx.y() * mm, xvtx.z() * mm, xvtx.t() * mm / c_light);
+
+    for (HepMC::GenVertex::particle_iterator vpitr = (*vitr)->particles_begin(HepMC::children);
+         vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr)
+    {
+      if ((*vpitr)->status() != 1) continue;
+
+      G4int pdgcode = (*vpitr)->pdg_id();
+      pos = (*vpitr)->momentum();
+      G4LorentzVector p(pos.px(), pos.py(), pos.pz(), pos.e());
+      G4PrimaryParticle* g4prim =
+        new G4PrimaryParticle(pdgcode, p.x() * GeV, p.y() * GeV, p.z() * GeV);
+
+      g4vtx->SetPrimary(g4prim);
+    }
+    g4event->AddPrimaryVertex(g4vtx);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+HepMC::GenEvent* HepMCG4Interface::GenerateHepMCEvent()
+{
+  HepMC::GenEvent* aevent = new HepMC::GenEvent();
+  return aevent;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void HepMCG4Interface::GeneratePrimaryVertex(G4Event* anEvent)
+{
+  // delete previous event object
+  delete hepmcEvent;
+
+  // generate next event
+  hepmcEvent = GenerateHepMCEvent();
+  if (!hepmcEvent) {
+    G4cout << "HepMCG4Interface: no generated particles. run terminated..." << G4endl;
+    G4RunManager::GetRunManager()->AbortRun();
+    return;
+  }
+  HepMC2G4(hepmcEvent, anEvent);
+}
diff --git a/src/MyG4ReactionDynamics.cc b/src/MyG4ReactionDynamics.cc
index 8ef97d0..40e6fa9 100755
--- a/src/MyG4ReactionDynamics.cc
+++ b/src/MyG4ReactionDynamics.cc
@@ -53,7 +53,7 @@
 #include "G4AntiNeutron.hh"
 #include "Randomize.hh"
 #include <iostream>
-#include "G4HadReentrentException.hh"
+// #include "G4HadReentrentException.hh"
 #include <signal.h>
 #include "G4ParticleTable.hh"
 
@@ -102,7 +102,7 @@
 */
  
  G4bool MyG4ReactionDynamics::GenerateXandPt(
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen,
    G4ReactionProduct &modifiedOriginal,   // Fermi motion & evap. effects included
    const G4HadProjectile *originalIncident,   // the original incident particle
@@ -314,7 +314,8 @@
           pVec->SetSide( -1 );                // backside particle, but not a nucleon
         }
         pVec->SetNewlyAdded( true );                // true is the same as IPA(i)<0
-        vec.SetElement( vecLen++, pVec );    
+        vec.push_back( pVec );    
+        vecLen++;
         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
         backwardEnergy -= pVec->GetMass()/GeV;
         ++backwardCount;
@@ -467,7 +468,7 @@
             {
               for(G4int ctr=0; ctr<vecLen; ctr++) delete vec[ctr];
    	      vecLen = 0;
-              throw G4HadReentrentException(__FILE__, __LINE__,
+              G4Exception("MyG4ReactionDynamics", "Hadronic", FatalException,
 	      "MyG4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
             }
             vec[i]->SetSide( -3 );
@@ -1013,13 +1014,18 @@
       else
         pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
       
-      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;  // tempV contains the backward nucleons
-      tempV.Initialize( backwardNucleonCount );
+      std::vector<G4ReactionProduct*> tempV;  // tempV contains the backward nucleons
       G4int tempLen = 0;
-      if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
+      if( targetParticle.GetSide() == -3 ) {
+        tempV.push_back( &targetParticle );
+        tempLen++;
+      }
       for( i=0; i<vecLen; ++i )
       {
-        if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
+        if( vec[i]->GetSide() == -3 ) {
+          tempV.push_back( vec[i] );
+          tempLen++;
+        }
       }
       if( tempLen != backwardNucleonCount )
       {
@@ -1233,10 +1239,12 @@
       tempR[0] = currentParticle;
       tempR[1] = targetParticle;
       for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
-      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
-      tempV.Initialize( vecLen+2 );
+      std::vector<G4ReactionProduct*> tempV;
       G4int tempLen = 0;
-      for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
+      for( i=0; i<vecLen+2; ++i ) {
+        tempV.push_back( &tempR[i] );
+        tempLen++;
+      }
       constantCrossSection = true;
 
       wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
@@ -1382,7 +1390,7 @@
   }
  
  void MyG4ReactionDynamics::SuppressChargedPions(
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen,
    const G4ReactionProduct &modifiedOriginal,
    G4ReactionProduct &currentParticle,
@@ -1457,7 +1465,7 @@
   }
  
  G4bool MyG4ReactionDynamics::TwoCluster(
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen,
    G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effects included
    const G4HadProjectile *originalIncident, // the original incident particle
@@ -1512,7 +1520,7 @@
       {
         for(G4int ctr=0; ctr<vecLen; ctr++) delete vec[ctr];
 	vecLen = 0;
-        throw G4HadReentrentException(__FILE__, __LINE__,
+        G4Exception("MyG4ReactionDynamics", "Hadronic", FatalException,
 	"MyG4ReactionDynamics::TwoCluster: Negative number of particles");
       }
       delete vec[vecLen-1];
@@ -1613,7 +1621,8 @@
         pVec->SetSide( -2 );    // backside particle
         extraMass += pVec->GetMass()/GeV;
         pVec->SetNewlyAdded( true );
-        vec.SetElement( vecLen++, pVec );
+        vec.push_back( pVec );
+        vecLen++;
         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
       }
     }
@@ -1891,20 +1900,25 @@
       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
     if( forwardCount > 1 )     // tempV will contain the forward particles
     {
-      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
-      tempV.Initialize( forwardCount );
+      std::vector<G4ReactionProduct*> tempV;
       G4bool constantCrossSection = true;
       G4int tempLen = 0;
-      if( currentParticle.GetSide() == 1 )
-        tempV.SetElement( tempLen++, &currentParticle );
-      if( targetParticle.GetSide() == 1 )
-        tempV.SetElement( tempLen++, &targetParticle );
+      if( currentParticle.GetSide() == 1 ) {
+        tempV.push_back( &currentParticle );
+        tempLen++;
+      }
+      if( targetParticle.GetSide() == 1 ) {
+        tempV.push_back( &targetParticle );
+        tempLen++;
+      }
       for( i=0; i<vecLen; ++i )
       {
         if( vec[i]->GetSide() == 1 )
         {
-          if( tempLen < 18 )
-            tempV.SetElement( tempLen++, vec[i] );
+          if( tempLen < 18 ) {
+            tempV.push_back( vec[i] );
+            tempLen++;
+          }
           else
           {
             vec[i]->SetSide( -1 );
@@ -1929,20 +1943,25 @@
       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
     if( backwardCount > 1 )   //  tempV will contain the backward particles,
     {                         //  but not those created from the intranuclear cascade
-      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
-      tempV.Initialize( backwardCount );
+      std::vector<G4ReactionProduct*> tempV;
       G4bool constantCrossSection = true;
       G4int tempLen = 0;
-      if( currentParticle.GetSide() == -1 )
-        tempV.SetElement( tempLen++, &currentParticle );
-      if( targetParticle.GetSide() == -1 )
-        tempV.SetElement( tempLen++, &targetParticle );
+      if( currentParticle.GetSide() == -1 ) {
+        tempV.push_back( &currentParticle );
+        tempLen++;
+      }
+      if( targetParticle.GetSide() == -1 ) {
+        tempV.push_back( &targetParticle );
+        tempLen++;
+      }
       for( i=0; i<vecLen; ++i )
       {
         if( vec[i]->GetSide() == -1 )
         {
-          if( tempLen < 18 )
-            tempV.SetElement( tempLen++, vec[i] );
+          if( tempLen < 18 ) {
+            tempV.push_back( vec[i] );
+            tempLen++;
+          }
           else
           {
             vec[i]->SetSide( -2 );
@@ -2133,11 +2152,13 @@
       tempR[1] = targetParticle;
       for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
 
-      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
-      tempV.Initialize( vecLen+2 );
+      std::vector<G4ReactionProduct*> tempV;
       G4bool constantCrossSection = true;
       G4int tempLen = 0;
-      for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
+      for( i=0; i<vecLen+2; ++i ){
+        tempV.push_back( &tempR[i] );
+        tempLen++;
+      }
 
       if( tempLen >= 2 )
       {
@@ -2287,7 +2308,7 @@
   }
  
  void MyG4ReactionDynamics::TwoBody(
-  G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+  std::vector<G4ReactionProduct*> &vec,
   G4int &vecLen,
   G4ReactionProduct &modifiedOriginal,
   const G4DynamicParticle */*originalTarget*/,
@@ -2611,7 +2632,7 @@
  G4double MyG4ReactionDynamics::GenerateNBodyEvent(
   const G4double totalEnergy,                // MeV
   const G4bool constantCrossSection,
-  G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+  std::vector<G4ReactionProduct*> &vec,
   G4int &vecLen )
   {
 //      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
@@ -2890,7 +2911,7 @@
    const G4ReactionProduct &modifiedOriginal,
    G4ReactionProduct &currentParticle,
    G4ReactionProduct &targetParticle,
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen )
   {
     const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
@@ -2951,7 +2972,7 @@
   const G4Nucleus &targetNucleus,
   G4ReactionProduct &currentParticle,
   G4ReactionProduct &targetParticle,
-  G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+  std::vector<G4ReactionProduct*> &vec,
   G4int &vecLen )
   {
     // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
@@ -3223,7 +3244,7 @@
    const G4ReactionProduct &modifiedOriginal,
    G4double spall,
    const G4Nucleus &targetNucleus,
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen )
   {
     // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
@@ -3279,7 +3300,12 @@
           p1->SetDefinition( aProton );
         else
           p1->SetDefinition( aNeutron );
-        vec.SetElement( vecLen, p1 );
+        if(vecLen < vec.size()) {
+          vec[vecLen] = p1;          
+        }
+        else{
+          vec.push_back( p1 );
+        }
         ++spall;
         G4double cost = G4UniformRand() * 2.0 - 1.0;
         G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
@@ -3351,7 +3377,12 @@
 	  delete p2;
 	  break;
 	}
-        vec.SetElement( vecLen, p2 );
+        if(vecLen < vec.size()) {
+          vec[vecLen] = p2;          
+        }
+        else{
+          vec.push_back( p2 );
+        }
         vec[vecLen]->SetNewlyAdded( true );
         vec[vecLen]->SetKineticEnergy( kinetic*GeV );
 	kinCreated+=kinetic;
@@ -3369,7 +3400,7 @@
    const G4ReactionProduct &modifiedOriginal,
    G4ReactionProduct &currentParticle,
    G4ReactionProduct &targetParticle,
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen )
   {
     const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
@@ -3406,7 +3437,7 @@
   }
 
  void MyG4ReactionDynamics::ProduceStrangeParticlePairs(
-   G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen,
    const G4ReactionProduct &modifiedOriginal,
    const G4DynamicParticle *originalTarget, 
@@ -3544,7 +3575,8 @@
 	  vec[0]->SetMayBeKilled(false);
 	  p1->SetMayBeKilled(false);
         }
-        vec.SetElement( vecLen++, p1 );
+        vec.push_back( p1 );
+        vecLen++;
       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
       }
       else
@@ -3615,7 +3647,8 @@
            break;
         }
         (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
-        vec.SetElement( vecLen++, p1 );
+        vec.push_back( p1 );
+        vecLen++;
       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
       }
       else                                        // replace
@@ -3797,7 +3830,7 @@
 
  void
   MyG4ReactionDynamics::NuclearReaction(
-   G4FastVector<G4ReactionProduct,4> &vec,
+   std::vector<G4ReactionProduct*> &vec,
    G4int &vecLen,
    const G4HadProjectile *originalIncident,
    const G4Nucleus &targetNucleus,
@@ -3953,12 +3986,15 @@
     //
     // use phase space routine in centre of mass system
     //
-    G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
-    tempV.Initialize( nt );
+    std::vector<G4ReactionProduct*> tempV;
     G4int tempLen = 0;
-    tempV.SetElement( tempLen++, v[0] );
-    tempV.SetElement( tempLen++, v[1] );
-    if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
+    tempV.push_back( v[0] );
+    tempV.push_back( v[1] );
+    tempLen += 2;
+    if( nt == 3 ) {
+      tempV.push_back( v[2] );
+      tempLen++;
+    }
     G4bool constantCrossSection = true;
     GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
     v[0]->Lorentz( *v[0], pseudo2 );
@@ -4072,16 +4108,19 @@
     vecLen = 0;
     if( particleIsDefined )
     {
-      vec.SetElement( vecLen++, v[0] );
+      vec.push_back( v[0] );
+      vecLen++;
     }
     else
     {
       delete v[0];
     }
-    vec.SetElement( vecLen++, v[1] );
+    vec.push_back( v[1] );
+    vecLen++;
     if( nt == 3 )
     {
-      vec.SetElement( vecLen++, v[2] );
+      vec.push_back( v[2] );
+      vecLen++;
     }
     else
     {
diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc
index 4981542..ca6a051 100644
--- a/src/PrimaryGeneratorAction.cc
+++ b/src/PrimaryGeneratorAction.cc
@@ -7,7 +7,7 @@
 #include "G4ParticleGun.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PrimaryGeneratorAction::PrimaryGeneratorAction()
+PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction()
 {
   // default generator is particle gun.
   currentGenerator = particleGun= new ParticleGun();
diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc
index 9b74be7..f1dcbdf 100644
--- a/src/SteppingAction.cc
+++ b/src/SteppingAction.cc
@@ -24,8 +24,8 @@ constexpr long rhbit = 0b1000000000000000000;
 constexpr long phbit = 0b1000000000000000000000;
 constexpr long etbit = 0b10000000000000000000000;
 constexpr double rhbitres = 10*um;
-constexpr double phbitres = 4*pow(10,-6);
-constexpr double etbitres = 2*pow(10,-6);
+constexpr double phbitres = 4e-6;
+constexpr double etbitres = 2e-6;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
diff --git a/src/TrackerConstruction.cc b/src/TrackerConstruction.cc
index 9a69498..4b2b293 100644
--- a/src/TrackerConstruction.cc
+++ b/src/TrackerConstruction.cc
@@ -6,7 +6,7 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-constexpr G4double trPhiAng = 2;
+constexpr G4double trPhiAng = 360;
 constexpr G4int PIB_num = 3;
 constexpr G4double PIB_rMin[] = {44,73,102};
 constexpr G4double PIB_thick = 380;