diff --git a/src/modules/moduleCollisions.f95 b/src/modules/moduleCollisions.f95
index 7367b29..5af2644 100644
--- a/src/modules/moduleCollisions.f95
+++ b/src/modules/moduleCollisions.f95
@@ -135,6 +135,7 @@ MODULE moduleCollisions
       vRel = species(1)%obj%m*NORM2(part_i%v - part_j%v)**2
       sigmaVrel = self%crossSec%get(vRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
+      PRINT *, sigmaVrelMaxNew/sigmaVrelMax, RAND()
       IF (sigmaVrelMaxNew/sigmaVrelMax > RAND()) THEN
         !Applies the collision
         v_i = NORM2(part_i%v)
diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95
index e48e3c6..df637f2 100644
--- a/src/modules/moduleMesh.f95
+++ b/src/modules/moduleMesh.f95
@@ -95,6 +95,8 @@ MODULE moduleMesh
     INTEGER(KIND=OMP_LOCK_KIND):: lock
     !Number of collisions per volume
     INTEGER:: nColl = 0
+    !Total weight of particles inside cell
+    REAL(8):: totalWeight = 0.D0
     CONTAINS
       PROCEDURE(initVol_interface),     DEFERRED, PASS:: init
       PROCEDURE(scatter_interface),     DEFERRED, PASS:: scatter
diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index ad9c29e..9291f8f 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -454,6 +454,7 @@ MODULE moduleMeshCyl
         !Assign particle to listPart_in
         CALL OMP_SET_LOCK(self%lock)
         CALL self%listPart_in%add(part)
+        self%totalWeight = self%totalWeight + part%weight
         CALL OMP_UNSET_LOCK(self%lock)
 
       ELSE
@@ -525,43 +526,48 @@ MODULE moduleMeshCyl
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
+      use moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVolCyl), INTENT(inout):: self
-      INTEGER:: Npart !Number of particles inside the cell
-      REAL(8):: Fn !Specific weight
-      REAL(8):: Pmax !Maximum probability of collision
+      INTEGER:: nPart !Number of particles inside the cell
+      REAL(8):: pMax !Maximum probability of collision
       INTEGER:: rnd !random index
       TYPE(particle), POINTER:: part_i, part_j
       INTEGER:: n !collision
       INTEGER:: ij, k
       REAL(8):: sigmaVrelMaxNew
 
-      Fn = species(1)%obj%weight!TODO: Check how to do this for multiple species
-      Npart = self%listPart_in%amount
-      Pmax = Fn*self%sigmaVrelMax*tau/self%volume
-      self%nColl = INT(REAL(Npart*(Npart-1))*Pmax*0.5D0)
+      self%nColl = 0
+      nPart = self%listPart_in%amount
+      IF (nPart > 1) THEN
+        pMax = self%totalWeight*self%sigmaVrelMax*tau/self%volume
+        self%nColl = INT(REAL(nPart)*pMax*0.5D0)
 
-      DO n = 1, self%NColl
-        !Select random numbers
-        rnd = 1 + FLOOR(Npart*RAND())
-        part_i => self%listPart_in%get(rnd)
-        rnd = 1 + FLOOR(Npart*RAND())
-        part_j => self%listPart_in%get(rnd)
-        ij = interactionIndex(part_i%pt, part_j%pt)
-        sigmaVrelMaxNew = 0.D0
-        DO k = 1, interactionMatrix(ij)%amount
-          CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
+        DO n = 1, self%nColl
+          !Select random numbers
+          rnd = 1 + FLOOR(nPart*RAND())
+          part_i => self%listPart_in%get(rnd)
+          rnd = 1 + FLOOR(nPart*RAND())
+          part_j => self%listPart_in%get(rnd)
+          ij = interactionIndex(part_i%pt, part_j%pt)
+          sigmaVrelMaxNew = 0.D0
+          DO k = 1, interactionMatrix(ij)%amount
+            CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
 
+          END DO
+
+          !Update maximum cross section*v_rel per each collision
+          IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
+            self%sigmaVrelMax = sigmaVrelMaxNew
+
+          END IF
+            
         END DO
 
-        !Update maximum cross section*v_rel per each collision
-        IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
-          self%sigmaVrelMax = sigmaVrelMaxNew
+      END IF
 
-        END IF
-          
-      END DO
+      self%totalWeight = 0.D0
 
       !Reset output in nodes
       CALL self%resetOutput()