diff --git a/examples/13.IAXO/Neutrons.rml b/examples/13.IAXO/Neutrons.rml
index 6d83fbb..b2179b6 100644
--- a/examples/13.IAXO/Neutrons.rml
+++ b/examples/13.IAXO/Neutrons.rml
@@ -16,6 +16,7 @@ Author: Luis Obis (lobis@unizar.es)
         -->
 
         <parameter name="seed" value="17022"/>
+        <parameter name="storeHadronicTargetInfo" value="true"/>
         <parameter name="nRequestedEntries" value="1"/>
 
         <generator type="cosmic">
diff --git a/include/Application.h b/include/Application.h
index 0408982..a8b7a44 100644
--- a/include/Application.h
+++ b/include/Application.h
@@ -26,6 +26,8 @@ struct Options {
     int nRequestedEntries = 0;
     int timeLimitSeconds = 0;
 
+    int runNumber = -1;
+
     // reference to original argc and argv necessary to pass to G4UIExecutive
     int argc;
     char** argv;
diff --git a/include/PrimaryGeneratorAction.h b/include/PrimaryGeneratorAction.h
index 4888275..a9069cf 100644
--- a/include/PrimaryGeneratorAction.h
+++ b/include/PrimaryGeneratorAction.h
@@ -37,6 +37,7 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     SimulationManager* fSimulationManager;
 
     static std::mutex fDistributionFormulaMutex;
+    static std::mutex fPrimaryGenerationMutex;
 
     std::vector<TRestGeant4Particle> fTempParticles;
 
diff --git a/src/Application.cxx b/src/Application.cxx
index 658b553..6f40343 100644
--- a/src/Application.cxx
+++ b/src/Application.cxx
@@ -244,6 +244,18 @@ Options ProcessCommandLineOptions(int argc, char* const argv[]) {
                 cerr << "--time option requires one argument." << endl;
                 exit(1);
             }
+        } else if (arg == "--runNumber") {
+            if (i + 1 < argc) {  // Make sure we aren't at the end of argv!
+                options.runNumber =
+                    stoi(argv[++i]);  // Increment 'i' so we don't get the argument as the next argv[i].
+                if (options.runNumber < 0) {
+                    cout << "--runNumber option error: runNumber cannot be negative" << endl;
+                    exit(1);
+                }
+            } else {
+                cerr << "--runNumber option requires one argument." << endl;
+                exit(1);
+            }
         } else {
             const string argument = argv[i];
             if (argument[0] == '-') {
@@ -347,6 +359,10 @@ void Application::Run(const CommandLineOptions::Options& options) {
 
     run->LoadConfigFromFile(inputRmlClean);
 
+    if (options.runNumber >= 0) {
+        run->SetRunNumber(options.runNumber);
+    }
+
     if (!options.outputFile.empty()) {
         run->SetOutputFileName(options.outputFile);
     }
@@ -362,7 +378,6 @@ void Application::Run(const CommandLineOptions::Options& options) {
 
     run->AddMetadata(fSimulationManager.GetRestMetadata());
     run->AddMetadata(fSimulationManager.GetRestPhysicsLists());
-
     run->PrintMetadata();
 
     run->FormOutputFile();
diff --git a/src/DataModel.cxx b/src/DataModel.cxx
index da441ee..2e18db9 100644
--- a/src/DataModel.cxx
+++ b/src/DataModel.cxx
@@ -211,6 +211,29 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) {
     fKineticEnergy.emplace_back(track->GetKineticEnergy() / CLHEP::keV);
     fMomentumDirection.emplace_back(momentum.x(), momentum.y(), momentum.z());
 
+    string isotopeName;
+    int atomicNumber = 0;
+    int atomicMassNumber = 0;
+
+    if (metadata->GetStoreHadronicTargetInfo() && track->GetCurrentStepNumber() != 0 &&
+        process->GetProcessType() == G4ProcessType::fHadronic) {
+        auto hadronicProcess = dynamic_cast<const G4HadronicProcess*>(process);
+        G4Nucleus nucleus = *(hadronicProcess->GetTargetNucleus());
+
+        auto isotope = nucleus.GetIsotope();
+        if (isotope) {
+            isotopeName = isotope->GetName();
+            atomicNumber = isotope->GetZ();
+            atomicMassNumber = isotope->GetN();
+        }
+    }
+
+    if (metadata->GetStoreHadronicTargetInfo()) {
+        fHadronicTargetIsotopeName.emplace_back(isotopeName);
+        fHadronicTargetIsotopeZ.emplace_back(atomicNumber);
+        fHadronicTargetIsotopeA.emplace_back(atomicMassNumber);
+    }
+
     SimulationManager::GetOutputManager()->AddEnergyToVolumeForParticleForProcess(energy, volumeName,
                                                                                   particleName, processName);
 }
diff --git a/src/PrimaryGeneratorAction.cxx b/src/PrimaryGeneratorAction.cxx
index fdec6fe..f53e732 100644
--- a/src/PrimaryGeneratorAction.cxx
+++ b/src/PrimaryGeneratorAction.cxx
@@ -3,12 +3,11 @@
 
 #include <TF1.h>
 #include <TF2.h>
-#include <TRestGeant4Event.h>
 #include <TRestGeant4Metadata.h>
+#include <TRestGeant4ParticleSourceCosmics.h>
 #include <TRestGeant4PrimaryGeneratorInfo.h>
 
 #include <G4Event.hh>
-#include <G4Geantino.hh>
 #include <G4IonTable.hh>
 #include <G4ParticleDefinition.hh>
 #include <G4ParticleTable.hh>
@@ -39,7 +38,8 @@ G4ThreeVector ComputeCosmicPosition(const G4ThreeVector& direction, double radiu
     return position;
 }
 
-std::mutex PrimaryGeneratorAction::fDistributionFormulaMutex = std::mutex();
+std::mutex PrimaryGeneratorAction::fDistributionFormulaMutex;
+std::mutex PrimaryGeneratorAction::fPrimaryGenerationMutex;
 TF1* PrimaryGeneratorAction::fEnergyDistributionFunction = nullptr;
 TF1* PrimaryGeneratorAction::fAngularDistributionFunction = nullptr;
 TF2* PrimaryGeneratorAction::fEnergyAndAngularDistributionFunction = nullptr;
@@ -231,6 +231,7 @@ void PrimaryGeneratorAction::SetGeneratorSpatialDensity(TString str) {
 }
 
 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) {
+    lock_guard<mutex> lock(fPrimaryGenerationMutex);
     auto simulationManager = fSimulationManager;
     TRestGeant4Metadata* restG4Metadata = simulationManager->GetRestMetadata();
 
@@ -241,6 +242,51 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) {
     // GeneratePrimaries is called first, and we want to store event origin and position inside
     // we should have already written the information from previous event to disk (in endOfEventAction)
 
+    if (restG4Metadata->GetNumberOfSources() == 1 &&
+        string(restG4Metadata->GetParticleSource(0)->GetName()) == "TRestGeant4ParticleSourceCosmics") {
+        auto source = dynamic_cast<TRestGeant4ParticleSourceCosmics*>(restG4Metadata->GetParticleSource(0));
+
+        if (string(restG4Metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorType().Data()) !=
+            "Cosmic") {
+            cerr << "PrimaryGeneratorAction - ERROR: cosmic generator only supports cosmic type: "
+                 << restG4Metadata->GetGeant4PrimaryGeneratorInfo().GetSpatialGeneratorType().Data() << endl;
+            exit(1);
+        }
+
+        source->Update();
+
+        const auto particle = source->GetParticles().at(0);
+
+        fParticleGun.SetParticleDefinition(
+            G4ParticleTable::GetParticleTable()->FindParticle(particle.GetParticleName().Data()));
+
+        if (fCosmicCircumscribedSphereRadius == 0.) {
+            // radius in mm
+            fCosmicCircumscribedSphereRadius = fSimulationManager->GetRestMetadata()
+                                                   ->GetGeant4PrimaryGeneratorInfo()
+                                                   .GetSpatialGeneratorCosmicRadius();
+        }
+
+        const auto position = ComputeCosmicPosition(fParticleGun.GetParticleMomentumDirection(),
+                                                    fCosmicCircumscribedSphereRadius);
+
+        fParticleGun.SetParticlePosition(position);
+        fParticleGun.SetParticleEnergy(particle.GetEnergy() * keV);
+        fParticleGun.SetParticleMomentumDirection({particle.GetMomentumDirection().X(),
+                                                   particle.GetMomentumDirection().Y(),
+                                                   particle.GetMomentumDirection().Z()});
+        fParticleGun.GeneratePrimaryVertex(event);
+
+        /*
+        cout << "PrimaryGeneratorAction - INFO: cosmic particle generated: "
+             << fParticleGun.GetParticleDefinition()->GetParticleName()
+             << " energy: " << fParticleGun.GetParticleEnergy() / keV << " keV"
+             << " position: " << fParticleGun.GetParticlePosition() / mm << " mm"
+             << " direction: " << fParticleGun.GetParticleMomentumDirection() << endl;
+        */
+        return;
+    }
+
     for (int i = 0; i < restG4Metadata->GetNumberOfSources(); i++) {
         restG4Metadata->GetParticleSource(i)->Update();
     }