diff --git a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt
index b90d5b9..8a62bf1 100644
--- a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt
+++ b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt
@@ -1,14 +1,17 @@
 package inr.numass.trapping
 
 import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
-import org.apache.commons.math3.random.SynchronizedRandomGenerator
+import org.apache.commons.math3.random.JDKRandomGenerator
 import org.apache.commons.rng.UniformRandomProvider
 import org.apache.commons.rng.simple.RandomSource
 import java.io.File
 import java.io.PrintStream
 import java.nio.file.Files
 import java.nio.file.StandardOpenOption
-import java.util.stream.Stream
+import java.util.stream.LongStream
+
+private val seedGenerator = JDKRandomGenerator()
+
 
 /**
  * Created by darksnake on 04-Jun-16.
@@ -20,7 +23,10 @@ class SimulationManager() {
 
     var comment = ""
 
-    var generator: UniformRandomProvider = RandomSource.create(RandomSource.SPLIT_MIX_64)
+    /**
+     * A supplier for random generator. Each track has its own generator
+     */
+    var generatorFactory: (Long) -> UniformRandomProvider = { RandomSource.create(RandomSource.MT_64, seedGenerator.nextInt()) }
     var reportFilter: (Simulator.SimulationResult) -> Boolean = { it.state == Simulator.EndState.ACCEPTED }
 
     var initialE = 18000.0
@@ -85,8 +91,8 @@ class SimulationManager() {
             outputDirectory.mkdirs()
         }
 
-        val generator = SynchronizedRandomGenerator(RandomGeneratorBridge(generator))
-        val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField, generator)
+
+        val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField)
 
         val counter = Counter()
 
@@ -100,18 +106,19 @@ class SimulationManager() {
 
 
         val outputPath = outputDirectory.toPath().resolve("$fileName.out")
+
         PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output ->
             output.println("# " + header.replace("\n", "\n# "))//adding comment symbols
 
             System.out.printf("%nStarting simulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong())
-            output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state")
-            Stream
-                    .generate {
+            output.printf("%s\t%s\t%s\t%s\t%s\t%s\t%s%n", "id", "E", "theta", "theta_start", "colNum", "L", "state")
+            LongStream.rangeClosed(1, num.toLong()).parallel()
+                    .mapToObj { id ->
+                        val generator = RandomGeneratorBridge(generatorFactory(id))
                         val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi
                         val z = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH
-                        simulator.simulate(initialE, theta, z).also { counter.count(it) }
+                        simulator.simulate(id, generator, initialE, theta, z).also { counter.count(it) }
                     }
-                    .limit(num.toLong()).parallel()
                     .filter(reportFilter)
                     .forEach { res ->
                         printOne(output, res)
@@ -148,7 +155,7 @@ class SimulationManager() {
     }
 
     private fun printOne(out: PrintStream, res: Simulator.SimulationResult) {
-        out.printf("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString())
+        out.printf("%d\t%g\t%g\t%g\t%d\t%g\t%s%n", res.id, res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString())
         out.flush()
     }
 
diff --git a/src/main/kotlin/inr/numass/trapping/Simulator.kt b/src/main/kotlin/inr/numass/trapping/Simulator.kt
index 90f135c..6ce096b 100644
--- a/src/main/kotlin/inr/numass/trapping/Simulator.kt
+++ b/src/main/kotlin/inr/numass/trapping/Simulator.kt
@@ -23,8 +23,8 @@ class Simulator(
         val thetaPinch: Double,
         val gasDensity: Double,// m^-3
         val bSource: Double,
-        val magneticField: ((Double) -> Double)?,
-        val generator: RandomGenerator) {
+        val magneticField: ((Double) -> Double)?
+) {
 
     enum class EndState {
 
@@ -55,20 +55,12 @@ class Simulator(
         //проверяем нормировку
         assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2))
 
-        val alpha = generator.nextDouble()
+        val alpha = pos.generator.nextDouble()
 
-        val delta: Pair<Double, Double>
-        if (alpha > sigmaEl) {
-            if (alpha > sigmaEl + sigmaExc) {
-                //ionization case
-                delta = Scatter.randomion(pos.e, generator)
-            } else {
-                //excitation case
-                delta = Scatter.randomexc(pos.e, generator)
-            }
-        } else {
-            // elastic
-            delta = Scatter.randomel(pos.e, generator)
+        val delta: Pair<Double, Double> = when {
+            alpha < sigmaEl -> Scatter.randomel(pos.e, pos.generator)// elastic
+            alpha > sigmaEl + sigmaExc -> Scatter.randomion(pos.e, pos.generator)//ionization case
+            else -> Scatter.randomexc(pos.e, pos.generator)//excitation case
         }
 
         //Обновляем значени угла и энергии независимо ни от чего
@@ -85,7 +77,7 @@ class Simulator(
      * @param e
      * @return
      */
-    private fun freePath(e: Double): Double {
+    private fun freePath(generator: RandomGenerator, e: Double): Double {
         //FIXME redundant cross-section calculation
         //All cross-sections are in m^2
         return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample()
@@ -196,15 +188,15 @@ class Simulator(
      * Симулируем один пробег электрона от начального значения и до вылетания из
      * иточника или до того момента, как энергия становится меньше eLow.
      */
-    fun simulate(initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult {
+    fun simulate(id: Long, generator: RandomGenerator, initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult {
         assert(initEnergy > 0)
         assert(initTheta > 0 && initTheta < Math.PI)
         assert(abs(initZ) <= SOURCE_LENGTH / 2.0)
 
-        val pos = State(initEnergy, initTheta, initZ)
+        val pos = State(generator, initEnergy, initTheta, initZ)
 
         while (!pos.isFinished) {
-            val dl = freePath(pos.e) // path to next scattering
+            val dl = freePath(generator, pos.e) // path to next scattering
             // propagate to next scattering position
             propagate(pos, dl)
 
@@ -220,7 +212,7 @@ class Simulator(
             }
         }
 
-        return SimulationResult(pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l)
+        return SimulationResult(id, pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l)
     }
 
     fun resetDebugCounters() {
@@ -236,7 +228,7 @@ class Simulator(
         }
     }
 
-    class SimulationResult(var state: EndState, var E: Double, var theta: Double, var initTheta: Double, var collisionNumber: Int, var l: Double)
+    data class SimulationResult(val id: Long, val state: EndState, val E: Double, val theta: Double, val initTheta: Double, val collisionNumber: Int, val l: Double)
 
     /**
      * Current electron position in simulation. Not thread safe!
@@ -245,7 +237,7 @@ class Simulator(
      * @property theta Current theta recalculated to the field in the center of the source
      * @property z current z. Zero is the center of the source
      */
-    private inner class State(internal var e: Double, internal var theta: Double, internal var z: Double) {
+    private inner class State(val generator: RandomGenerator, internal var e: Double, internal var theta: Double, internal var z: Double) {
         /**
          * Current total path
          */
@@ -349,7 +341,7 @@ class Simulator(
         /**
          * Reverse electron direction
          */
-        internal fun flip() {
+        fun flip() {
             if (theta < 0 || theta > PI) {
                 throw Error()
             }
@@ -361,7 +353,7 @@ class Simulator(
          *
          * @return
          */
-        internal fun field(): Double {
+        fun field(): Double {
             return this@Simulator.field(z)
         }
 
@@ -370,7 +362,7 @@ class Simulator(
          *
          * @return
          */
-        internal fun realTheta(): Double {
+        fun realTheta(): Double {
             if (magneticField == null) {
                 return theta
             } else {
@@ -391,7 +383,7 @@ class Simulator(
          * @param dTheta
          * @return resulting angle
          */
-        internal fun addTheta(dTheta: Double): Double {
+        fun addTheta(dTheta: Double): Double {
             //Генерируем случайный фи
             val phi = generator.nextDouble() * 2.0 * Math.PI
 
@@ -406,7 +398,7 @@ class Simulator(
 
             val result = rot.applyTo(init.cartesian)
 
-            val newtheta = acos(result.z)
+            val newTheta = acos(result.z)
 
             //            //следим чтобы угол был от 0 до Pi
             //            if (newtheta < 0) {
@@ -418,10 +410,10 @@ class Simulator(
 
             //change back to virtual angles
             if (magneticField == null) {
-                theta = newtheta
+                theta = newTheta
             } else {
-                theta = asin(sin(newtheta) * sqrt(bSource / field()))
-                if (newtheta > PI / 2) {
+                theta = asin(sin(newTheta) * sqrt(bSource / field()))
+                if (newTheta > PI / 2) {
                     theta = PI - theta
                 }
             }
@@ -435,8 +427,8 @@ class Simulator(
 
     companion object {
 
-        val SOURCE_LENGTH = 3.0
-        private val DELTA_L = 0.1 //step for propagate calculation
+        const val SOURCE_LENGTH = 3.0
+        private const val DELTA_L = 0.1 //step for propagate calculation
     }
 
 
diff --git a/src/main/kotlin/inr/numass/trapping/Trapping.kt b/src/main/kotlin/inr/numass/trapping/Trapping.kt
index eb1a61f..9ed354d 100644
--- a/src/main/kotlin/inr/numass/trapping/Trapping.kt
+++ b/src/main/kotlin/inr/numass/trapping/Trapping.kt
@@ -1,7 +1,5 @@
 package inr.numass.trapping
 
-import org.apache.commons.rng.simple.RandomSource
-import java.io.File
 import java.time.Duration
 import java.time.Instant
 
@@ -21,12 +19,11 @@ fun main(args: Array<String>) {
         SimulationManager().apply {
             comment = "Out of the box cross-sections"
             fileName = "trap[$e]"
-            generator = RandomSource.create(RandomSource.MWC_256)
             setFields(0.6, 3.7, 7.2)
             gasDensity = 1e19
             initialE = e
             range = 4000.0
-        }.simulateAll(1e7)
+        }.simulateAll(1_000_000)
     }
 
     val finishTime = Instant.now()