diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index fcefb8c..5314cd1 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -39,7 +39,7 @@ jobs:
         uses: actions/cache@v4
         with:
           path: ${{ env.G4DATA_DIR }}
-          key: g4data-G4EMLOW-8.2-G4ENSDFSTATE-2.3
+          key: g4data-G4EMLOW-8.6.1-G4ENSDFSTATE-3.0
 
       - name: Download Geant4 EMLOW & ENSDFSTATE datasets
         if: steps.cache-g4data.outputs.cache-hit != 'true'
@@ -47,19 +47,19 @@ jobs:
           set -euo pipefail
           mkdir -p "$G4DATA_DIR"
           cd "$G4DATA_DIR"
-          curl -sSLO https://cern.ch/geant4-data/datasets/G4EMLOW.8.2.tar.gz
-          tar -xzf G4EMLOW.8.2.tar.gz
-          rm G4EMLOW.8.2.tar.gz
-          curl -sSLO https://cern.ch/geant4-data/datasets/G4ENSDFSTATE.2.3.tar.gz
-          tar -xzf G4ENSDFSTATE.2.3.tar.gz
-          rm G4ENSDFSTATE.2.3.tar.gz
+          curl -sSLO https://cern.ch/geant4-data/datasets/G4EMLOW.8.6.1.tar.gz
+          tar -xzf G4EMLOW.8.6.1.tar.gz
+          rm G4EMLOW.8.6.1.tar.gz
+          curl -sSLO https://cern.ch/geant4-data/datasets/G4ENSDFSTATE.3.0.tar.gz
+          tar -xzf G4ENSDFSTATE.3.0.tar.gz
+          rm G4ENSDFSTATE.3.0.tar.gz
 
       - name: Pull OpenTOPAS Docker image
         run: docker pull $TOPAS_DOCKER_IMAGE
 
       - name: Download TOPAS docker launcher
         run: |
-          curl -sSfL https://raw.githubusercontent.com/OpenTOPAS/OpenTOPAS/fix-bugs/docker/topas-docker -o docker-run
+          curl -sSfL https://raw.githubusercontent.com/OpenTOPAS/OpenTOPAS/main/docker/topas-docker -o docker-run
           chmod +x docker-run
 
       - name: Install nrtest tooling
diff --git a/README.md b/README.md
index 6c21a9a..eb0360b 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# TOPAS-nBio
+# TOPAS-nBio (Version 4.1.0)
 This is the TOPAS-nBio extension repository, a Monte Carlo simulation framework for (sub-) cellular radiobiology.
 
 TOPAS-nBio is described here: https://topas-nbio.readthedocs.io/. 
diff --git a/geometry/dna/TsLinearDNA.cc b/geometry/dna/TsLinearDNA.cc
index eb7eaeb..792b7d3 100755
--- a/geometry/dna/TsLinearDNA.cc
+++ b/geometry/dna/TsLinearDNA.cc
@@ -12,7 +12,6 @@
 // A strand of DNA.
 
 #include "TsLinearDNA.hh"
-
 #include "TsParameterManager.hh"
 
 #include "G4VPhysicalVolume.hh"
diff --git a/processes/TsEmDNAPhysics.cc b/processes/TsEmDNAPhysics.cc
index 35158db..5f654b1 100644
--- a/processes/TsEmDNAPhysics.cc
+++ b/processes/TsEmDNAPhysics.cc
@@ -241,7 +241,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBElasticModel* modelDNAPTBElastic = new G4DNAPTBElasticModel("THF/TMP/PY", particle);
                         
                         G4DNAModelInterface* e_elasticInteraction = new G4DNAModelInterface("e-_elastic_interaction");
-                        e_elasticInteraction->RegisterModel(e_modelDNARutherfordElastic, particle);
+                        e_elasticInteraction->RegisterModel(e_modelDNARutherfordElastic);
                         e_elasticInteraction->RegisterModel(modelDNAPTBElastic);
                         e_elasticInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -255,7 +255,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBElasticModel* modelDNAPTBElastic = new G4DNAPTBElasticModel("THF/TMP/PY", particle);
                         
                         G4DNAModelInterface* e_elasticInteraction = new G4DNAModelInterface("e-_elastic_interaction");
-                        e_elasticInteraction->RegisterModel(e_modelDNACPA100Model, particle);
+                        e_elasticInteraction->RegisterModel(e_modelDNACPA100Model);
                         e_elasticInteraction->RegisterModel(modelDNAPTBElastic);
                         e_elasticInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -270,7 +270,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBElasticModel* modelDNAPTBElastic = new G4DNAPTBElasticModel("THF/TMP/PY", particle);
                         
                         G4DNAModelInterface* e_elasticInteraction = new G4DNAModelInterface("e-_elastic_interaction");
-                        e_elasticInteraction->RegisterModel(e_modelDNAUeharaScreenedRutherfordElastic, particle);
+                        e_elasticInteraction->RegisterModel(e_modelDNAUeharaScreenedRutherfordElastic);
                         e_elasticInteraction->RegisterModel(modelDNAPTBElastic);
                         e_elasticInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -347,7 +347,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBExcitationModel* modelDNAPTBExcitation = new G4DNAPTBExcitationModel("THF/TMP/PY",particle);
                         G4DNAModelInterface* e_excitationInteraction = new G4DNAModelInterface("e-_excitation_interaction");
                         
-                        e_excitationInteraction->RegisterModel(e_modelDNAEmfietzoglouExcitation,particle);
+                        e_excitationInteraction->RegisterModel(e_modelDNAEmfietzoglouExcitation);
                         e_excitationInteraction->RegisterModel(modelDNAPTBExcitation);
                         e_excitationInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -360,7 +360,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBExcitationModel* modelDNAPTBExcitation = new G4DNAPTBExcitationModel("THF/TMP/PY",particle);
                         G4DNAModelInterface* e_excitationInteraction = new G4DNAModelInterface("e-_excitation_interaction");
                         
-                        e_excitationInteraction->RegisterModel(e_modelDNABornExcitation,particle);
+                        e_excitationInteraction->RegisterModel(e_modelDNABornExcitation);
                         e_excitationInteraction->RegisterModel(modelDNAPTBExcitation);
                         e_excitationInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -454,7 +454,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBIonisationModel* modelDNAPTBIonisation = new G4DNAPTBIonisationModel("THF/TMP/PY",particle);
                         
                         G4DNAModelInterface* e_ionisationInteraction = new G4DNAModelInterface("e-_ionisation_interaction");
-                        e_ionisationInteraction->RegisterModel(e_modelDNAEmfietzoglouIonisation,particle);
+                        e_ionisationInteraction->RegisterModel(e_modelDNAEmfietzoglouIonisation);
                         e_ionisationInteraction->RegisterModel(modelDNAPTBIonisation);
                         e_ionisationInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -482,7 +482,7 @@ void TsEmDNAPhysics::ConstructProcess()
                         G4DNAPTBIonisationModel* modelDNAPTBIonisation = new G4DNAPTBIonisationModel("THF/TMP/PY",particle);
                         
                         G4DNAModelInterface* e_ionisationInteraction = new G4DNAModelInterface("e-_ionisation_interaction");
-                        e_ionisationInteraction->RegisterModel(e_modelDNABornIonisation,particle);
+                        e_ionisationInteraction->RegisterModel(e_modelDNABornIonisation);
                         e_ionisationInteraction->RegisterModel(modelDNAPTBIonisation);
                         e_ionisationInteraction->RegisterModel(new G4DNAVacuumModel());
                         
@@ -904,3 +904,4 @@ void TsEmDNAPhysics::ConstructProcess()
     G4LossTableManager::Instance()->SetAtomDeexcitation(de);
 }
 
+
diff --git a/processes/TsEmDNAPhysics.hh b/processes/TsEmDNAPhysics.hh
index 9c2624d..e2338bd 100644
--- a/processes/TsEmDNAPhysics.hh
+++ b/processes/TsEmDNAPhysics.hh
@@ -36,4 +36,4 @@ private:
     G4int  verbose;
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/processes/TsIRT.cc b/processes/TsIRT.cc
index ff46f50..8a5c232 100755
--- a/processes/TsIRT.cc
+++ b/processes/TsIRT.cc
@@ -26,7 +26,6 @@
 #include <set>
 
 TsIRT::TsIRT(TsParameterManager* pM, G4String parmName): TsVIRTProcedure(pM,parmName){
-	fReactionConf = new TsIRTConfiguration(parmName, pM);
 
 	fSpeciesIndex = 0;
 	fMoleculesName = fReactionConf->GetMoleculeNames();
@@ -123,7 +122,7 @@ void TsIRT::AddMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4Thre
 TsIRTConfiguration::TsMolecule TsIRT::ConstructMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	pdg = fMoleculesIDs[name];
 
 	TsIRTConfiguration::TsMolecule aMol;
@@ -371,7 +370,8 @@ void TsIRT::ConductReactions() {
 				
 			} else {
 				fReactionConf->Diffuse(fChemicalSpecies[iM],irt-fChemicalSpecies[iM].time);
-				for ( int ip = 0; ip < 3; ip++ )
+                //positions = fReactionConf->GetBackgroundPositionOfProducts(fChemicalSpecies[iM], indexOfReaction);
+                for ( int ip = 0; ip < 3; ip++ )
 					positions.push_back(fChemicalSpecies[iM].position);
 				
 				binReaction = fReactionConf->GetReaction(indexOfReaction);
diff --git a/processes/TsIRT.hh b/processes/TsIRT.hh
index 442edcb..b27e711 100755
--- a/processes/TsIRT.hh
+++ b/processes/TsIRT.hh
@@ -49,7 +49,6 @@ protected:
 	std::unordered_map<G4int,std::unordered_map<G4int,std::unordered_map<G4int,std::unordered_map<G4int,G4bool>>>> fSpaceBinned;
 	std::map<G4int, G4String> fMoleculesName;
 
-	TsIRTConfiguration* fReactionConf;
 	TsParameterManager* fPm;
 	G4String fName;
 
@@ -59,4 +58,3 @@ protected:
 	
 };
 #endif
-
diff --git a/processes/TsIRTConfiguration.cc b/processes/TsIRTConfiguration.cc
index cbd31e5..3a7ff26 100755
--- a/processes/TsIRTConfiguration.cc
+++ b/processes/TsIRTConfiguration.cc
@@ -47,6 +47,16 @@ fKick(false), fAllTotallyDiffusionControlled(false)
 	fExistingMolecules["trioxide"]         = "O3^-1";
 	fExistingMolecules["ozone"]            = "O3^0";
 	fExistingMolecules["none"]             = "None";
+
+	// Normalize Geant4 molecule naming to the canonical TOPAS-nBio names.
+	fGeant4NameOverrides["\u00b0OH^0"]   = "OH^0";
+	fGeant4NameOverrides["H_O2\u00b0^0"] = "HO2^0";
+	fGeant4NameOverrides["O_2^0"] = "O2^0";
+	fGeant4NameOverrides["O_2^-1"] = "O2^-1";
+	fGeant4NameOverrides["O_3^0"] = "O3^0";
+	fGeant4NameOverrides["O_3^-1"] = "O3^-1";
+	fGeant4NameOverrides["\u00b0O^0"] = "O^0";
+	fGeant4NameOverrides["O\u00b0^-1"] = "O^-1";
 	
 	AddMolecule("H^0",     7.0e9*nm*nm/s,   0, 0.19*nm);
 	AddMolecule("OH^0",    2.2e9*nm*nm/s,   0, 0.22*nm);
@@ -65,6 +75,8 @@ fKick(false), fAllTotallyDiffusionControlled(false)
 	AddMolecule("O3^0",    2.0e9*nm*nm/s,  0, 0.20*nm);
 	AddMolecule("None",    0.0e9*nm*nm/s,  0, 0.00*nm);
 	
+	RegisterGeant4Aliases();
+	
 	// Re-set diffusion coefficient, radius and/or charge.
 	std::vector<G4String>* moleculeNames = new std::vector<G4String>;
 	fPm->GetParameterNamesStartingWith("Mo/", moleculeNames);
@@ -104,7 +116,7 @@ fKick(false), fAllTotallyDiffusionControlled(false)
 				G4int molID =  fMoleculesID[fExistingMolecules[molName]];
 				fMoleculesDefinition[molID].charge = fPm->GetUnitlessParameter(fullName);
 				G4cout << molName << ": Re-set charge to "
-				<< fMoleculesDefinition[molID].radius/nm << " nm" << std::endl;
+				<< fMoleculesDefinition[molID].charge << std::endl;
 				moleculeExists = true;
 			}
 		}
@@ -323,6 +335,9 @@ fKick(false), fAllTotallyDiffusionControlled(false)
 		InsertBackgroundReaction(reactorA, reactorB,vProduct,reactionRate,concentration,false);
 	}
 	
+    ResolveReactionParameters();
+    PrintReactionsInformation();
+    
 	// Re-scale chemistry parameters based on temperature
 	parName = "Ch/" + chemistryList + "/Temperature";
 	fScaleForTemperature = false;
@@ -336,33 +351,39 @@ fKick(false), fAllTotallyDiffusionControlled(false)
 			fKick = fPm->GetBooleanParameter(parName);
 		
 		AdjustReactionAndDiffusionRateForTemperature();
+        G4cout << "======================================== performed Temperature Scaling ===========================================" << G4endl;
+        ResolveReactionParameters();
+        PrintReactionsInformation();
 	}
-	
-	if ( fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance") ) {
-		fpHSolventConcentration = 0.0;
-		fpHValue                = 7.1;
-		
-		fpHSolvent = fPm->GetStringParameter("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance");
-		G4StrUtil::to_lower(fpHSolvent);
-		if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration") &&
-			fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH")) {
-			G4String message = "Cannot be defined when parameter: Ch/" + chemistryList + "/ModelAcidPropertiesWithConcentration is used.";
-			Quit("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH",message);
-		}
-		
-		if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration") ) {
-			fpHSolventConcentration = fPm->GetDoubleParameter("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration","molar concentration");
-			AdjustReactionRateForPH("Concentration");
-		}
-		
-		if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH") ) {
-			fpHValue = fPm->GetUnitlessParameter("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH");
-			AdjustReactionRateForPH("PH");
-		}
-	}
-
-	ResolveReactionParameters();
-	PrintReactionsInformation();
+    
+    
+    
+    if ( fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance") ) {
+        fpHSolventConcentration = 0.0;
+        fpHValue                = 7.1;
+        
+        fpHSolvent = fPm->GetStringParameter("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance");
+        G4StrUtil::to_lower(fpHSolvent);
+        if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration") &&
+            fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH")) {
+            G4String message = "Cannot be defined when parameter: Ch/" + chemistryList + "/ModelAcidPropertiesWithConcentration is used.";
+            Quit("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH",message);
+        }
+        
+        if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration") ) {
+            fpHSolventConcentration = fPm->GetDoubleParameter("Ch/"+chemistryList+"/ModelAcidPropertiesWithConcentration","molar concentration");
+            AdjustReactionRateForPH("Concentration");
+        }
+        
+        if (fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH") ) {
+            fpHValue = fPm->GetUnitlessParameter("Ch/"+chemistryList+"/ModelAcidPropertiesWithpH");
+            AdjustReactionRateForPH("PH");
+        }
+        
+        G4cout << "======================================== performed ph Scaling ===========================================" << G4endl;
+        ResolveReactionParameters();
+        PrintReactionsInformation();
+    }
 }
 
 
@@ -384,7 +405,7 @@ void TsIRTConfiguration::AddMolecule(G4String name, G4double diffusionCoefficien
 		fMoleculesDefinition[moleculeID] = aMolecule;
 		fMoleculesID[name] = moleculeID;
 		fMoleculesName[moleculeID] = name;
-		fLastMoleculeID = moleculeID;
+		fLastMoleculeID = std::max(fLastMoleculeID, moleculeID);
 	}else{
 		fLastMoleculeID++;
 		fMoleculesDefinition[fLastMoleculeID] = aMolecule;
@@ -404,6 +425,33 @@ void TsIRTConfiguration::AddMolecule(G4String name, G4double diffusionCoefficien
 }
 
 
+G4String TsIRTConfiguration::NormalizeMoleculeName(const G4String& name) {
+	const std::string key = name;
+	auto it = fGeant4NameOverrides.find(key);
+	if (it != fGeant4NameOverrides.end()) {
+		if (fLoggedOverrides.find(name) == fLoggedOverrides.end()) {
+			//G4cout << "TsIRTConfiguration: remapping Geant4 molecule name '" << name
+			//	   << "' to '" << it->second << "'" << G4endl;
+			fLoggedOverrides.insert(name);
+		}
+		return it->second.c_str();
+	}
+	return name;
+}
+
+void TsIRTConfiguration::RegisterGeant4Aliases() {
+	for (const auto& overridePair : fGeant4NameOverrides) {
+		const G4String alias = overridePair.first.c_str();
+		const G4String canonical = overridePair.second.c_str();
+		auto idIt = fMoleculesID.find(canonical);
+		if (idIt != fMoleculesID.end()) {
+			fMoleculesID[alias] = idIt->second;
+			fExistingMolecules[alias] = canonical;
+		}
+	}
+}
+
+
 G4bool TsIRTConfiguration::MoleculeExists(G4String name) {
 	if ( fMoleculesID.find(name) == fMoleculesID.end() ) return false;
 	return true;
@@ -487,6 +535,32 @@ void TsIRTConfiguration::ResolveReactionParameters() {
 
 			kdiff = 4 * pi * sumDiffCoeff * effectiveReactionRadius * Avogadro;
 
+            /*
+            if ( kobs >= kdiff ) {
+                G4cout << "TOPAS-nBio is exiting due to an error in Chemistry setup for reaction type " << reactionType << ":" <<G4endl;
+                G4cout << "  " << GetMoleculeNameFromMoleculeID(molA) << " + "
+                << GetMoleculeNameFromMoleculeID(molB) << G4endl;
+                G4cout << "This reaction has kobs >= kdiff after pH scaling (kobs=" << kobs
+                << ", kdiff=" << kdiff << ")." << G4endl;
+                G4cout << "This causes (kdif-kobs) negative or zero. When defining kact, it set the reaction probability to unrealistic values. Adjust rate or type." << G4endl;
+                fPm->AbortSession(1);
+            }
+        
+            G4double ratio = kobs / kdiff;
+            if (ratio > 0.5) {
+                G4cout << "IRT: type " << reactionType << " "
+                << GetMoleculeNameFromMoleculeID(molA) << " + "
+                << GetMoleculeNameFromMoleculeID(molB)
+                << " kobs=" << kobs << " kdiff=" << kdiff
+                << " ratio=" << ratio << G4endl;
+            }
+            // optionally treat ratio >= 0.99 as fatal
+            if (ratio >= 0.99) {
+                G4cout << GetMoleculeNameFromMoleculeID(molA) << " " << GetMoleculeNameFromMoleculeID(molB) << G4endl;
+                fPm->AbortSession(1);
+            }
+            */
+            
 			kact = kdiff * kobs / (kdiff - kobs);
 
             G4double sigmaEffRs = reactionRadius + Rs;
@@ -645,9 +719,39 @@ std::vector<G4ThreeVector> TsIRTConfiguration::ResampleReactantsPosition(TsMolec
 	G4double dtA = std::fabs(time-molA.time);
 	G4double dtB = std::fabs(time-molB.time);
 	G4double dt  = dtA > dtB ? dtA : dtB;
-	
-	if ( D1 == 0 ) molB.position = r1;
-	else if ( D2 == 0 )	molA.position = r2;
+    std::vector<G4ThreeVector> result;
+    G4int nOfProducts = fReactions[index].products.size();
+    G4ThreeVector position;
+    
+    if ( D1 == 0 ) {
+        molB.position = r1;
+        molB.time = time;
+        molA.time = time;
+        result.push_back(r1);
+        result.push_back(r1);
+        if (nOfProducts <= 2) {
+            return result;
+        } else {
+            result.push_back(r1);
+            result.push_back(r1);
+            return result;
+        }
+        
+    } else if ( D2 == 0 ) {
+        molA.position = r2;
+        molA.time = time;
+        molB.time = time;
+        result.push_back(r2);
+        result.push_back(r2);
+        if (nOfProducts <= 2) {
+            return result;
+        } else {
+            result.push_back(r2);
+            result.push_back(r2);
+            return result;
+        }
+    }
+    
 	if (dt == 0) dt = 1*ps;
 	
 	G4ThreeVector S1 = r1-r2;
@@ -678,16 +782,13 @@ std::vector<G4ThreeVector> TsIRTConfiguration::ResampleReactantsPosition(TsMolec
 	molB.position = R2;
 	molB.time = time;
 	
-	std::vector<G4ThreeVector> result;
-	G4int nOfProducts = fReactions[index].products.size();
-	
 	result.push_back(R1);
 	result.push_back(R2);
 
 	if(nOfProducts <= 2) return result;
 
 	// weighted-position
-	G4ThreeVector position = R1*std::sqrt(D2)/(std::sqrt(D1) + std::sqrt(D2)) + R2*std::sqrt(D1)/(std::sqrt(D1) + std::sqrt(D2));
+    position = R1*std::sqrt(D2)/(std::sqrt(D1) + std::sqrt(D2)) + R2*std::sqrt(D1)/(std::sqrt(D1) + std::sqrt(D2));
 
 	result.push_back(position);
 	result.push_back(R1 + 0.5 * (R1 + R2));
@@ -702,11 +803,14 @@ std::vector<G4ThreeVector> TsIRTConfiguration::GetBackgroundPositionOfProducts(T
 	
 	G4double D1   = fMoleculesDefinition[molA.id].diffusionCoefficient;
 	G4double kobs = fReactions[index].kobs;
+    if (kobs < 0 )
+        kobs = 1e10 * 1.0e-3*m*m*m/(mole*s); // dummy kobs for now
+        
 	G4ThreeVector CurrentPos = molA.position;
 	std::vector<G4int> Products = fReactions[index].products;
 
 	for (size_t i = 0; i < Products.size(); i++) {
-		G4double r  = kobs / (4*CLHEP::pi*(2*D1));
+        G4double r  = kobs / (4*CLHEP::pi*(2*D1)*CLHEP::Avogadro);
 		G4double theta = 2 * CLHEP::pi * G4UniformRand();
 		G4double phi   = acos(1 - 2 * G4UniformRand());
 		G4double x = r * sin(phi) * cos(theta);
@@ -784,7 +888,7 @@ G4int  TsIRTConfiguration::ContactFirstOrderAndBackgroundReactions(TsMolecule mo
 		G4double R3 = R*R*R;
 		G4double Cs = fReactions[u].concentration;
 		if (Cs > 50) {continue;}          // No Contact Reactions with water molecules
-		Cs *= CLHEP::Avogadro;  // nm3 to M multiply by 10^-24 Nav = 6.02214076×10^23x10^-24 = 0.602214076
+		Cs *= CLHEP::Avogadro;  // nm3 to M multiply by 10^-24 Nav = 6.02214076×10^23x10^-24 = 0.602214076. 
 		G4double prob1 = std::exp(-1.33333333*CLHEP::pi*R3*Cs);
 
 		if ( G4UniformRand() < 1. - prob1 ) return (int)u;
@@ -939,7 +1043,7 @@ std::vector<G4double> TsIRTConfiguration::GetH2SO4ComponentsConcentrationP(G4dou
 	
 	for (size_t i = 0; i < Result.size(); i++) {
 		if ((Result[i] > 0 and Result[i] < 3.0*_C) and i % 2 != 1) {
-			if (Result[i + 1] == 0) {
+			if ( i+1< Result.size() && Result[i + 1] == 0) {
 				_H_pos = Result[i];
 			}
 		}
@@ -1206,122 +1310,200 @@ void TsIRTConfiguration::AdjustReactionAndDiffusionRateForTemperature() {
 }
 
 void TsIRTConfiguration::AdjustReactionRateForPH(G4String pHOrConcentration) {
-	std::vector<G4double> AcidComponents;
-	G4double Ionic   = 0.0;
-	G4double HCon    = 0.0;
-	G4double HCon25  = 1.00E-7;
-	G4double OHCon   = 0.0;
-	G4double OHCon25 = 1.00E-7;
-	G4double HSO4Con = 0.0;
-
-	if (pHOrConcentration == "PH" && fpHSolvent == "h2so4") {
-		AcidComponents = GetH2SO4ComponentsConcentrationPH(fpHValue);
-		Ionic          = GetIonicStrength(AcidComponents);
-		HCon           = AcidComponents[0];
-		HSO4Con        = AcidComponents[1];
-		OHCon          = AcidComponents[3];
-		G4cout << "-- Adjust for PH of H2SO4 " << G4endl;
-
-	}
-
-	else if (pHOrConcentration == "Concentration" && fpHSolvent == "h2so4") {
-		AcidComponents = GetH2SO4ComponentsConcentrationP(fpHSolventConcentration/fPm->GetUnitValue("M"));
-		Ionic          = GetIonicStrength(AcidComponents);
-		HCon           = AcidComponents[0];
-		HSO4Con        = AcidComponents[1];
-		OHCon          = AcidComponents[3];
-		G4cout << "-- Adjust for Concentration of H2SO4 " << G4endl;
-	}
-
-	else if (fpHSolvent == "generic") {
-		HCon  = pow(10,-fpHValue);
-		OHCon = 1E-14 / HCon;
-		AcidComponents = {HCon, 0.0, 0.0, OHCon, 0.0, 0.0};
-		Ionic          = GetIonicStrength(AcidComponents);
-		G4cout << "-- Adjust for a generic substance " << G4endl;
-
-	} else {
-		G4cout << "-- Is doing nothing " << pHOrConcentration << " " << fpHSolvent << G4endl;
-	}
-
-	G4cout << G4endl;
-	G4cout << " ###-------- pH Scaling Starts ---------###" << G4endl;
-
-	for(size_t i = 0; i < fReactions.size(); i++) {
-		G4int chargeA   = fMoleculesDefinition[fReactions[i].reactorA].charge;
-		G4int chargeB   = fMoleculesDefinition[fReactions[i].reactorB].charge;
-		G4String ReactA = fMoleculesName[fReactions[i].reactorA];
-		G4String ReactB = fMoleculesName[fReactions[i].reactorB];
-		std::vector<G4String> products;
-		G4double k_Before = 0;
-
-		for (size_t vsize = 0; vsize < fReactions[i].products.size(); vsize++){
-			products.push_back(fMoleculesName[fReactions[i].products[vsize]]);
-		}
-
-		fReactions[i].kobs /= fPm->GetUnitValue("/M/s");
-		fReactions[i].scavengingCapacity /= 1/s;
-
-		if ((chargeA != 0 && chargeB != 0)) {
-			if (fReactions[i].reactionType != 6) {
-				k_Before = fReactions[i].kobs;
-				fReactions[i].kobs = IonicRate(Ionic, fReactions[i]);
-			}
-
-			else if (fReactions[i].reactionType == 6 && ReactB == "H3O^1") {
-				k_Before = fReactions[i].scavengingCapacity;
-				fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * HCon;
-			}
-
-			else if (fReactions[i].reactionType == 6 && ReactB == "OH^-1") {
-				k_Before = fReactions[i].scavengingCapacity;
-				fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * OHCon;
-			}
-
-			else {
-				G4cout << "========= "<< fReactions[i].reactionType << G4endl;
-				k_Before = fReactions[i].kobs;
-//				fReactions[i].kobs = fReactions[i].kobs * 1E-7;
-				fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i]);
-			}
-		}
-
-		else if ((fReactions[i].reactionType == 6) && (ReactB == "H3O^1")) {
-			k_Before = fReactions[i].scavengingCapacity;
-			fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * HCon / HCon25;
-		}
-
-		else if ((fReactions[i].reactionType == 6) && (ReactB == "OH^-1")) {
-			k_Before = fReactions[i].scavengingCapacity;
-			fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * OHCon / OHCon25;
-		}
-
-		else if ((fReactions[i].reactionType == 6) && (ReactB == "HSO4^-1")) {
-			k_Before = fReactions[i].scavengingCapacity;
-			fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * HSO4Con;
-		}
-
-		if (k_Before > 0) {
-			G4cout << " Reaction Type: " << fReactions[i].reactionType << " | Reaction: " << ReactA << " + " << ReactB;
-			for (int prod = 0; prod < (int)products.size(); prod++) {
-				if (prod == 0)
-					G4cout << " -> ";
-				G4cout << products[prod];
-				if (prod < (int)products.size() - 1) {
-					G4cout << " + ";
-				}
-			}
-			if ( fReactions[i].reactionType == 6 )
-				G4cout << G4endl << " ---- scav: " << k_Before << " ---> "  << fReactions[i].scavengingCapacity << G4endl << G4endl;
-			else
-				G4cout << G4endl << " ---- kobs: " << k_Before << " ---> "  << fReactions[i].kobs << G4endl << G4endl;
-		}
-		fReactions[i].kobs *= fPm->GetUnitValue("/M/s");
-		fReactions[i].scavengingCapacity *= 1/s;
-	}
-
-	G4cout << " ###-------- pH Scaling Ends ---------###" << G4endl;
-	G4cout << G4endl;
+    std::vector<G4double> AcidComponents;
+    G4double Ionic   = 0.0;
+    G4double HCon    = 0.0;
+    G4double HCon25  = 1.00E-7;
+    G4double OHCon   = 0.0;
+    G4double OHCon25 = 1.00E-7;
+    G4double HSO4Con = 0.0;
+    
+    if (pHOrConcentration == "PH" && fpHSolvent == "h2so4") {
+        AcidComponents = GetH2SO4ComponentsConcentrationPH(fpHValue);
+        Ionic          = GetIonicStrength(AcidComponents);
+        HCon           = AcidComponents[0];
+        HSO4Con        = AcidComponents[1];
+        OHCon          = AcidComponents[3];
+        G4cout << "-- Adjust for PH of H2SO4 " << G4endl;
+        
+    }
+    
+    else if (pHOrConcentration == "Concentration" && fpHSolvent == "h2so4") {
+        AcidComponents = GetH2SO4ComponentsConcentrationP(fpHSolventConcentration/fPm->GetUnitValue("M"));
+        Ionic          = GetIonicStrength(AcidComponents);
+        HCon           = AcidComponents[0];
+        HSO4Con        = AcidComponents[1];
+        OHCon          = AcidComponents[3];
+        G4cout << "-- Adjust for Concentration of H2SO4 " << G4endl;
+    }
+    
+    else if (fpHSolvent == "generic") {
+        HCon  = pow(10,-fpHValue);
+        OHCon = 1E-14 / HCon;
+        AcidComponents = {HCon, 0.0, 0.0, OHCon, 0.0, 0.0};
+        Ionic          = GetIonicStrength(AcidComponents);
+        G4cout << "-- Adjust for a generic substance " << G4endl;
+        
+    } else {
+        G4cout << "-- Is doing nothing " << pHOrConcentration << " " << fpHSolvent << G4endl;
+    }
+    
+    G4cout << G4endl;
+    G4cout << " ###-------- pH Scaling Starts ---------###" << G4endl;
+    
+    for(size_t i = 0; i < fReactions.size(); i++) {
+        G4int chargeA   = fMoleculesDefinition[fReactions[i].reactorA].charge;
+        G4int chargeB   = fMoleculesDefinition[fReactions[i].reactorB].charge;
+        
+        if ( chargeA == 0 && chargeB == 0 ) // No scaling in neutral  molecules
+            continue;
+        /*
+        if ( fReactions[i].reactionType != 6 ) { // scaling only for charged molecules, no scavegers
+            if (chargeA * chargeB == 0 ) {
+                continue;
+            }
+        }
+        */
+        G4String ReactA = fMoleculesName[fReactions[i].reactorA];
+        G4String ReactB = fMoleculesName[fReactions[i].reactorB];
+        std::vector<G4String> products;
+        G4double k_Before = 0;
+        
+        for (size_t vsize = 0; vsize < fReactions[i].products.size(); vsize++){
+            products.push_back(fMoleculesName[fReactions[i].products[vsize]]);
+        }
+        
+        fReactions[i].kobs /= fPm->GetUnitValue("/M/s");
+        fReactions[i].scavengingCapacity /= 1/s;
+        
+        /*
+        if ( fReactions[i].reactionType != 6 ) {
+            k_Before = fReactions[i].kobs;
+            fReactions[i].kobs = IonicRate(Ionic, fReactions[i]);
+        } else { // fReactions[i].reactionType == 6
+            if  (ReactB == "H3O^1" ) {
+                k_Before = fReactions[i].scavengingCapacity;
+                fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * HCon;
+            } else if (ReactB == "OH^-1" ) {
+                k_Before = fReactions[i].scavengingCapacity;
+                fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * OHCon;
+            } else if (ReactB == "HSO4^-1" ) {
+                k_Before = fReactions[i].scavengingCapacity;
+                fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * HCon / HCon25;
+            } else {
+                G4cout << "========= "<< fReactions[i].reactionType << G4endl;
+                k_Before = fReactions[i].kobs;
+                //                fReactions[i].kobs = fReactions[i].kobs * 1E-7;
+                fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i]);
+            }
+        }*/
+        /*
+         if ((chargeA != 0 && chargeB != 0)) {
+         if (fReactions[i].reactionType != 6) {
+         k_Before = fReactions[i].kobs;
+         fReactions[i].kobs = IonicRate(Ionic, fReactions[i]);
+         }
+         
+         else if (fReactions[i].reactionType == 6 && ReactB == "H3O^1") {
+         k_Before = fReactions[i].scavengingCapacity;
+         fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * HCon;
+         }
+         
+         else if (fReactions[i].reactionType == 6 && ReactB == "OH^-1") {
+         k_Before = fReactions[i].scavengingCapacity;
+         fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i])/ 1e-7 * OHCon;
+         }
+         
+         else {
+         G4cout << "========= "<< fReactions[i].reactionType << G4endl;
+         k_Before = fReactions[i].kobs;
+         //                fReactions[i].kobs = fReactions[i].kobs * 1E-7;
+         fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i]);
+         }
+         }
+         
+         else if ((fReactions[i].reactionType == 6) && (ReactB == "H3O^1")) {
+         k_Before = fReactions[i].scavengingCapacity;
+         fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * HCon / HCon25;
+         }
+         
+         else if ((fReactions[i].reactionType == 6) && (ReactB == "OH^-1")) {
+         k_Before = fReactions[i].scavengingCapacity;
+         fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * OHCon / OHCon25;
+         }
+         
+         else if ((fReactions[i].reactionType == 6) && (ReactB == "HSO4^-1")) {
+         k_Before = fReactions[i].scavengingCapacity;
+         fReactions[i].scavengingCapacity = fReactions[i].scavengingCapacity * HSO4Con;
+         }
+         */
+        const bool chargedPair  = (chargeA != 0 && chargeB != 0);
+        const bool isScavenging = (fReactions[i].reactionType == 6);
+        auto setScavenging = [&](G4double newValue) {
+            k_Before = fReactions[i].scavengingCapacity;
+            fReactions[i].scavengingCapacity = newValue;
+        };
+
+        if (chargedPair) {
+            if (!isScavenging) {
+                k_Before = fReactions[i].kobs;
+                fReactions[i].kobs = IonicRate(Ionic, fReactions[i]);
+            } else if (ReactB == "H3O^1") {
+                setScavenging(IonicRate(Ionic, fReactions[i]) / 1e-7 * HCon);
+            } else if (ReactB == "OH^-1") {
+                setScavenging(IonicRate(Ionic, fReactions[i]) / 1e-7 * OHCon);
+            } else {
+                G4cout << "========= "<< fReactions[i].reactionType << G4endl;
+                k_Before = fReactions[i].kobs;
+                //                fReactions[i].kobs = fReactions[i].kobs * 1E-7;
+                fReactions[i].scavengingCapacity = IonicRate(Ionic, fReactions[i]);
+            }
+        } else if (isScavenging) {
+            if (ReactB == "H3O^1") {
+                setScavenging(fReactions[i].scavengingCapacity * HCon / HCon25);
+            } else if (ReactB == "OH^-1") {
+                setScavenging(fReactions[i].scavengingCapacity * OHCon / OHCon25);
+            } else if (ReactB == "HSO4^-1") {
+                setScavenging(fReactions[i].scavengingCapacity * HSO4Con);
+            }
+        }
+        
+        if (k_Before > 0) {
+            G4cout << " Reaction Type: " << fReactions[i].reactionType << " | Reaction: " << ReactA << " + " << ReactB;
+            for (int prod = 0; prod < (int)products.size(); prod++) {
+                if (prod == 0)
+                    G4cout << " -> ";
+                G4cout << products[prod];
+                if (prod < (int)products.size() - 1) {
+                    G4cout << " + ";
+                }
+            }
+            if ( fReactions[i].reactionType == 6 )
+                G4cout << G4endl << " ---- scav: " << k_Before << " ---> "  << fReactions[i].scavengingCapacity << G4endl << G4endl;
+            else
+                G4cout << G4endl << " ---- kobs: " << k_Before << " ---> "  << fReactions[i].kobs << G4endl << G4endl;
+        }
+        // add this guard before units are restored
+        if (fReactions[i].reactionType == 2 || fReactions[i].reactionType == 4) {
+            G4double kobsUnitless  = fReactions[i].kobs; // still divided by /M/s
+            G4double kdiffUnitless = fReactions[i].kdif / fPm->GetUnitValue("/M/s");
+            if (kobsUnitless >= 0.99 * kdiffUnitless) {
+                G4cout << "TOPAS-nBio is exiting due to an error in Chemistry setup for reaction type " << fReactions[i].reactionType << ":" <<G4endl;
+                G4cout << "  " << GetMoleculeNameFromMoleculeID(fReactions[i].reactorA) << " + "
+                << GetMoleculeNameFromMoleculeID(fReactions[i].reactorB) << G4endl;
+                G4cout << "This reaction has kobs >= kdiff after pH scaling (kobs=" << kobsUnitless
+                << ", kdiff=" << kdiffUnitless << ")." << G4endl;
+                G4cout << "This causes (kdif-kobs) negative or zero. When defining kact, it set the reaction probability to unrealistic values. Adjust rate or type." << G4endl;
+                //fPm->AbortSession(1);
+            }
+        }
+        
+        fReactions[i].kobs *= fPm->GetUnitValue("/M/s");
+        fReactions[i].scavengingCapacity *= 1/s;
+    }
+    
+    G4cout << " ###-------- pH Scaling Ends ---------###" << G4endl;
+    G4cout << G4endl;
 }
 
 
@@ -1391,6 +1573,9 @@ G4double TsIRTConfiguration::GetOnsagerRadius(G4int molA, G4int molB) {
 
 
 void TsIRTConfiguration::PrintMoleculesInformation() {
+    if (!G4Threading::IsMasterThread())
+        return;
+    
 	for ( auto& molecules : fMoleculesID ) {
 		G4String name = molecules.first;
 		G4int id = molecules.second;
@@ -1410,6 +1595,9 @@ void TsIRTConfiguration::PrintMoleculesInformation() {
 
 void TsIRTConfiguration::PrintReactionsInformation() {
 
+    if (!G4Threading::IsMasterThread())
+        return;
+    
 	G4IosFlagsSaver iosfs(G4cout);
 	std::map<size_t, std::vector<TsMolecularReaction> > temporal;
 	for ( size_t i = 0; i < fReactions.size(); i++) {
@@ -1431,8 +1619,12 @@ void TsIRTConfiguration::PrintReactionsInformation() {
 
 	G4cout << G4endl;
 	G4int n = 0;
-	for ( size_t i = 1; i <= temporal.size(); i++ ) {
-		for ( auto& reactions : temporal[i] ) {
+	
+    for (const auto& kv : temporal) { //for ( size_t i = 1; i <= temporal.size(); i++ ) {
+        size_t i = kv.first;
+        const auto& reactionsInType = kv.second;
+        
+        for ( auto& reactions : reactionsInType) { //temporal[i] ) {
 			G4int type = reactions.reactionType;
 			G4int molA = reactions.reactorA;
 			G4int molB = reactions.reactorB;
diff --git a/processes/TsIRTConfiguration.hh b/processes/TsIRTConfiguration.hh
index b78d504..9fd54f4 100755
--- a/processes/TsIRTConfiguration.hh
+++ b/processes/TsIRTConfiguration.hh
@@ -8,6 +8,7 @@
 #include <vector>
 #include <map>
 #include <unordered_map>
+#include <set>
 
 class TsParameterManager;
 class TsIRTUtils;
@@ -24,7 +25,7 @@ public:
 	void InsertBackgroundReaction(G4String, G4String, std::vector<G4int>,
 								  G4double, G4double, G4bool);
 	void SetTimeLimits(G4double, G4double);
-
+	G4String NormalizeMoleculeName(const G4String& name);
 	
 
 	void QuitIfMoleculeNotFound(G4String mol);
@@ -95,6 +96,8 @@ private:
 	std::map<G4String, G4int> fMoleculesID;
 	std::map<G4int, G4String> fMoleculesName;
 	std::map<G4String, G4String> fExistingMolecules;
+	std::unordered_map<std::string, std::string> fGeant4NameOverrides;
+	std::set<G4String> fLoggedOverrides;
 	
 	std::map<G4int, TsMolecularReaction > fReactions;
 	std::map<G4int, std::vector<std::pair<G4int,G4int>>> fMoleculeCanReactWith;
@@ -115,6 +118,7 @@ private:
 	
 	G4bool fAllTotallyDiffusionControlled;
 	
+	void RegisterGeant4Aliases();
 public:
 	
 	G4double GetIndependentReactionTime(TsMolecule molA, TsMolecule molB, G4int indexOfReaction);
@@ -165,4 +169,3 @@ public:
 
 };
 #endif
-
diff --git a/processes/TsIRTContinuous.cc b/processes/TsIRTContinuous.cc
index 5918b50..1de81c0 100755
--- a/processes/TsIRTContinuous.cc
+++ b/processes/TsIRTContinuous.cc
@@ -26,8 +26,6 @@
 #include <set>
 
 TsIRTContinuous::TsIRTContinuous(TsParameterManager* pM, G4String parmName): TsVIRTProcedure(pM,parmName){
-	fReactionConf = new TsIRTConfiguration(parmName, pM);
-
 	fSpeciesIndex = 0;
 	fMoleculesName = fReactionConf->GetMoleculeNames();
 
@@ -148,7 +146,7 @@ void TsIRTContinuous::AddMolecule(G4Track* aTrack, G4double time, G4int molecule
 TsIRTConfiguration::TsMolecule TsIRTContinuous::ConstructMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	pdg = fMoleculesIDs[name];
 
 	TsIRTConfiguration::TsMolecule aMol;
diff --git a/processes/TsIRTContinuous.hh b/processes/TsIRTContinuous.hh
index 2586bf9..4ff1c3c 100755
--- a/processes/TsIRTContinuous.hh
+++ b/processes/TsIRTContinuous.hh
@@ -54,7 +54,6 @@ protected:
 	std::map<G4int, G4String> fMoleculesName;
 	std::vector<G4double> fPulseTime;
 
-	TsIRTConfiguration* fReactionConf;
 	TsParameterManager* fPm;
 	G4String fName;
 
@@ -72,4 +71,3 @@ protected:
 
 };
 #endif
-
diff --git a/processes/TsVIRTProcedure.cc b/processes/TsVIRTProcedure.cc
index 59cddf0..a333886 100644
--- a/processes/TsVIRTProcedure.cc
+++ b/processes/TsVIRTProcedure.cc
@@ -200,7 +200,7 @@ void TsVIRTProcedure::AddMolecule(TsIRTConfiguration::TsMolecule aMol) {
 void TsVIRTProcedure::AddMolecule(G4Step* aStep, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aStep->GetPreStepPoint()->GetPosition();
-	const G4String& name = GetMolecule(aStep->GetTrack())->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aStep->GetTrack())->GetName());
 	
 	pdg = fMoleculesIDs[name];
 	
@@ -228,7 +228,7 @@ void TsVIRTProcedure::AddMolecule(G4Step* aStep, G4double time, G4int moleculeID
 void TsVIRTProcedure::AddMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset, G4bool isDNA) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	
 	pdg = fMoleculesIDs[name];
 	      
@@ -255,7 +255,7 @@ void TsVIRTProcedure::AddMolecule(G4Track* aTrack, G4double time, G4int molecule
 void TsVIRTProcedure::AddMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	
 	pdg = fMoleculesIDs[name];
 	
@@ -282,7 +282,7 @@ void TsVIRTProcedure::AddMolecule(G4Track* aTrack, G4double time, G4int molecule
 void TsVIRTProcedure::AddMolecule(const G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	
 	pdg = fMoleculesIDs[name];
 	
@@ -309,7 +309,7 @@ void TsVIRTProcedure::AddMolecule(const G4Track* aTrack, G4double time, G4int mo
 TsIRTConfiguration::TsMolecule TsVIRTProcedure::ConstructMolecule(G4Track* aTrack, G4double time, G4int moleculeID, G4ThreeVector offset) {
 	G4int pdg = -1;
 	G4ThreeVector position = aTrack->GetPosition();
-	const G4String& name = GetMolecule(aTrack)->GetName();
+	G4String name = fReactionConf->NormalizeMoleculeName(GetMolecule(aTrack)->GetName());
 	pdg = fMoleculesIDs[name];
 
 	TsIRTConfiguration::TsMolecule aMol;
diff --git a/scorers/TsScoreDNADamageWithIRT.cc b/scorers/TsScoreDNADamageWithIRT.cc
index e2a7d6c..31ec806 100755
--- a/scorers/TsScoreDNADamageWithIRT.cc
+++ b/scorers/TsScoreDNADamageWithIRT.cc
@@ -213,8 +213,8 @@ void TsScoreDNADamageWithIRT::UserHookForEndOfEvent() {
 		G4int totalSSB = totalDirectSSB + totalIndirectSSB;
 		G4cout << " --- IRT ends " << GetEventID() << G4endl;
 
-		G4int numberOfBasePairs = G4int(totalSpecies/4);
-		G4double perDA = numberOfBasePairs * 660/0.34;
+		G4int numberOfBasePairs = G4int(totalSpecies/2);
+		G4double perDA = numberOfBasePairs * 660;
 		G4double per100eV_to_umolPerJoule = 0.1036;
 
 		G4cout << " -- total SSB: " << totalSSB << G4endl;
@@ -222,8 +222,8 @@ void TsScoreDNADamageWithIRT::UserHookForEndOfEvent() {
 		G4cout << "=========== Total indirect SSB      : " << totalIndirectSSB << G4endl;
 		G4cout << "=========== Total dose in envelope  : " << totalDose/gray<< " Gy " <<G4endl;
 		G4cout << "=========== Total edep in envelope  : " << fEnergyDepositPerEvent/MeV << " MeV " << G4endl;
-		G4cout << "=========== SSB in plasmid       /Gy: " <<  totalSSB/totalDose /gray       << " /Gy" << G4endl;
-		G4cout << "=========== SSB in plasmid    /Gy/Da: " <<  totalSSB/totalDose/gray/perDA << " /Gy/Da" << G4endl;
+		G4cout << "=========== SSB in plasmid       /Gy: " <<  totalSSB/(totalDose/gray)       << " /Gy" << G4endl;
+		G4cout << "=========== SSB in plasmid    /Gy/Da: " <<  totalSSB/(totalDose/gray)/perDA << " /Gy/Da" << G4endl;
 		G4cout << "=========== SSB in plasmid    /100eV: " << 100*totalSSB/(fEnergyDepositPerEvent/eV) << " /100 eV" << G4endl;
 		G4cout << "=========== SSB in plasmid    umol/J: " << per100eV_to_umolPerJoule*100*totalSSB/(fEnergyDepositPerEvent/eV) << " umol/J" << G4endl;
 		G4cout << "=========== Total events to reach a dose of " << fPrescribedDose/gray << " (Gy) = " << GetEventID() << G4endl;