diff --git a/src/simulation.cpp b/src/simulation.cpp
index 763add5d7..d31bcad3d 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -201,86 +201,18 @@ void process_calculate_xs_events(QueueItem * queue, int n)
 {
   // Sort queue by energy
   std::sort(queue, queue+n);
-  //std::sort(queue, queue+n, by_energy);
-
-  // Then, stable sort by material (so as to preserve energy ordering)
-  //std::stable_sort(queue, queue+n, by_material);
 
   // Save last_ members, find grid index
-  int lost_particles = 0;
-  #pragma omp parallel for reduction(+:lost_particles) schedule(dynamic,DYNAMIC_SIZE)
+  #pragma omp parallel for schedule(dynamic,DYNAMIC_SIZE)
   for (int i = 0; i < n; i++) {
 	  Particle *p = particles + queue[i].idx; 
-	  //std::cout << "particle offset = " << queue[i] << std::endl;
-    // Set the random number stream
-	// TODO: Move RNG seeds to particle storage
-  /*
-    if (p->type_ == Particle::Type::neutron) {
-		p->stream_ = STREAM_TRACKING;
-    } else {
-		p->stream_ = STREAM_PHOTON;
-    }
-
-    // Store pre-collision particle properties
-    p->wgt_last_ = p->wgt_;
-    p->E_last_ = p->E_;
-    p->u_last_ = p->u();
-    p->r_last_ = p->r();
-    
-    // Reset event variables
-    p->event_ = EVENT_KILL;
-    p->event_nuclide_ = NUCLIDE_NONE;
-    p->event_mt_ = REACTION_NONE;
-
-    // If the cell hasn't been determined based on the particle's location,
-    // initiate a search for the current cell. This generally happens at the
-    // beginning of the history and again for any secondary particles
-    if (p->coord_[p->n_coord_ - 1].cell == C_NONE) {
-      if (!find_cell(p, false)) {
-        p->mark_as_lost("Could not find the cell containing particle "
-          + std::to_string(p->id_));
-        //return;
-		lost_particles += 1;
-		continue;
-      }
-
-      // set birth cell attribute
-      if (p->cell_born_ == C_NONE) p->cell_born_ = p->coord_[p->n_coord_ - 1].cell;
-    }
-
-    // Write particle track.
-    if (p->write_track_) write_particle_track(*p);
-
-    if (settings::check_overlaps) check_cell_overlap(p);
-    */
     p->event_calculate_xs_I();
   }
 
-  if( lost_particles > 0 )
-	  exit(1);
-
   #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
   for( int i = 0; i < n; i++ )
   {
 	  Particle * p = particles + queue[i].idx;
-    /*
-	  // Calculate microscopic and macroscopic cross sections
-	  if (p->material_ != MATERIAL_VOID) {
-		  if (settings::run_CE) {
-			  if (p->material_ != p->material_last_ || p->sqrtkT_ != p->sqrtkT_last_) {
-				  // If the material is the same as the last material and the
-				  // temperature hasn't changed, we don't need to lookup cross
-				  // sections again.
-				  model::materials[p->material_]->calculate_xs(*p);
-			  }
-		  } // else MG not supported
-	  } else {
-		  p->macro_xs_.total      = 0.0;
-		  p->macro_xs_.absorption = 0.0;
-		  p->macro_xs_.fission    = 0.0;
-		  p->macro_xs_.nu_fission = 0.0;
-	  }
-    */
     p->event_calculate_xs_II();
   }
 
@@ -300,75 +232,9 @@ void process_calculate_xs_events(QueueItem * queue, int n)
 
 void process_advance_particle_events()
 {
-  //for (auto& p : advance_particle_queue) {
   #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
   for (int i = 0; i < advance_particle_queue_length; i++) {
 	  Particle * p = particles + advance_particle_queue[i].idx;
-    //p->trace_ == (p->id_ == 0);
-
-    /*
-    // Sample a distance to collision
-    double d_collision;
-    if (p->type_ == Particle::Type::electron ||
-        p->type_ == Particle::Type::positron) {
-      d_collision = 0.0;
-    } else if (p->macro_xs_.total == 0.0) {
-      d_collision = INFINITY;
-    } else {
-      d_collision = -std::log(prn(p->current_seed())) / p->macro_xs_.total;
-    }
-
-    // -------------- break here? -------------------
-
-    // Find the distance to the nearest boundary
-    p->boundary_ = distance_to_boundary(p);
-
-    // Select smaller of the two distances
-    double distance;
-	int idx;
-    if (p->boundary_.distance < d_collision) {
-		#pragma omp atomic capture
-		idx = surface_crossing_queue_length++;
-      surface_crossing_queue[idx].idx = advance_particle_queue[i].idx;
-      surface_crossing_queue[idx].E = p->E_;
-      surface_crossing_queue[idx].material = p->material_;
-      surface_crossing_queue[idx].type = p->type_;
-      distance = p->boundary_.distance;
-    } else {
-		#pragma omp atomic capture
-		idx = collision_queue_length++;
-		collision_queue[idx].idx = advance_particle_queue[i].idx;
-		collision_queue[idx].E = p->E_;
-		collision_queue[idx].material = p->material_;
-		collision_queue[idx].type = p->type_;
-      distance = d_collision;
-    }
-
-    // -------------- break here? -------------------
-
-    // Advance particle
-    for (int j = 0; j < p->n_coord_; ++j) {
-      p->coord_[j].r += distance * p->coord_[j].u;
-    }
-
-    // -------------- break here? -------------------
-
-    // Score track-length tallies
-    if (!model::active_tracklength_tallies.empty()) {
-      score_tracklength_tally(p, distance);
-    }
-
-    // Score track-length estimate of k-eff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
-        p->type_ == Particle::Type::neutron) {
-      p->tally_tracklength_ += p->wgt_ * distance * p->macro_xs_.nu_fission;
-    }
-
-    // Score flux derivative accumulators for differential tallies.
-    if (!model::active_tallies.empty()) {
-      score_track_derivative(p, distance);
-    }
-  */
     p->event_advance();
     if( p->collision_distance_ > p->boundary_.distance ) 
     {
@@ -396,48 +262,13 @@ void process_advance_particle_events()
 
 void process_surface_crossing_events()
 {
-  //for (auto& p : surface_crossing_queue) {
   #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
   for (int i = 0; i < surface_crossing_queue_length; i++) {
 	  Particle * p = particles + surface_crossing_queue[i].idx;
-    /*
-    // Set surface that particle is on and adjust coordinate levels
-    p->surface_ = p->boundary_.surface_index;
-    p->n_coord_ = p->boundary_.coord_level;
-
-    // Saving previous cell data
-    for (int j = 0; j < p->n_coord_; ++j) {
-      p->cell_last_[j] = p->coord_[j].cell;
-    }
-    p->n_coord_last_ = p->n_coord_;
-
-    if (p->boundary_.lattice_translation[0] != 0 ||
-        p->boundary_.lattice_translation[1] != 0 ||
-        p->boundary_.lattice_translation[2] != 0) {
-      // Particle crosses lattice boundary
-      cross_lattice(p, p->boundary_);
-      p->event_ = EVENT_LATTICE;
-    } else {
-      // Particle crosses surface
-      p->cross_surface();
-      p->event_ = EVENT_SURFACE;
-    }
-    // Score cell to cell partial currents
-    if (!model::active_surface_tallies.empty()) {
-      score_surface_tally(p, model::active_surface_tallies);
-    }
-
-    //if (!p->alive_ && !simulation::secondary_bank.empty()) {
-    if (!p->alive_ && !p->secondary_bank_.empty()) {
-      revive_particle_from_secondary(p);
-    }
-    */
     p->event_cross_surface();
-
+    p->event_revive_from_secondary();
     if (p->alive_)
-    {
       dispatch_xs_event(surface_crossing_queue[i].idx);
-    }
   }
 
   surface_crossing_queue_length = 0;
@@ -445,166 +276,18 @@ void process_surface_crossing_events()
 
 void process_collision_events()
 {
-  //for (auto& p : collision_queue) {
   #pragma omp parallel for schedule(dynamic,DYNAMIC_SIZE)
   for (int i = 0; i < collision_queue_length; i++) {
 	  Particle * p = particles + collision_queue[i].idx;
-    /*
-	  //std::cout << "Beginning collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-    // Score collision estimate of keff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
-        p->type_ == Particle::Type::neutron) {
-      p->tally_collision_ += p->wgt_ * p->macro_xs_.nu_fission
-        / p->macro_xs_.total;
-    }
-
-    // Score surface current tallies -- this has to be done before the collision
-    // since the direction of the particle will change and we need to use the
-    // pre-collision direction to figure out what mesh surfaces were crossed
-
-    if (!model::active_meshsurf_tallies.empty())
-      score_surface_tally(p, model::active_meshsurf_tallies);
-	    
-	//std::cout << "After surface tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Clear surface component
-    p->surface_ = 0;
-
-    if (settings::run_CE) {
-      collision(p);
-    } else {
-      collision_mg(p);
-    }
-
-	//std::cout << "After collision() of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Score collision estimator tallies -- this is done after a collision
-    // has occurred rather than before because we need information on the
-    // outgoing energy for any tallies with an outgoing energy filter
-    if (!model::active_collision_tallies.empty()) score_collision_tally(p);
-    if (!model::active_analog_tallies.empty()) {
-      if (settings::run_CE) {
-        score_analog_tally_ce(p);
-      } else {
-        score_analog_tally_mg(p);
-      }
-    }
-	//std::cout << "After analog tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Reset banked weight during collision
-    p->n_bank_ = 0;
-    p->n_bank_second_ = 0;
-    p->wgt_bank_ = 0.0;
-    for (int& v : p->n_delayed_bank_) v = 0;
-
-    // Reset fission logical
-    p->fission_ = false;
-
-    // Save coordinates for tallying purposes
-    p->r_last_current_ = p->r();
-
-    // Set last material to none since cross sections will need to be
-    // re-evaluated
-    p->material_last_ = C_NONE;
-
-    // Set all directions to base level -- right now, after a collision, only
-    // the base level directions are changed
-    for (int j = 0; j < p->n_coord_ - 1; ++j) {
-      if (p->coord_[j + 1].rotated) {
-        // If next level is rotated, apply rotation matrix
-        const auto& m {model::cells[p->coord_[j].cell]->rotation_};
-        const auto& u {p->coord_[j].u};
-        p->coord_[j + 1].u = u.rotate(m);
-      } else {
-        // Otherwise, copy this level's direction
-        p->coord_[j+1].u = p->coord_[j].u;
-      }
-    }
-
-    // Score flux derivative accumulators for differential tallies.
-    if (!model::active_tallies.empty()) score_collision_derivative(p);
-    */
     p->event_collide();
-
-    //if (!p->alive_ && !simulation::secondary_bank.empty()) {
-    /*
-    if (!p->alive_ && !p->secondary_bank_.empty()) {
-      revive_particle_from_secondary(p);
-    }
-    */
-
     p->event_revive_from_secondary();
-
     if (p->alive_)
-    {
       dispatch_xs_event(collision_queue[i].idx);
-      //std::cout << "Ended collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-      //assert(std::isfinite(p->E_) );
-    }
   }
 
   collision_queue_length = 0;
 }
 
-/*
-void check_energies(void)
-{
-	int * Q;
-	int n;
-
-
-	Q = calculate_fuel_xs_queue;
-	n = calculate_fuel_xs_queue_length;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index xs FUEL " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = calculate_nonfuel_xs_queue;
-	n = calculate_nonfuel_xs_queue_length ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index xs Non fuel " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = advance_particle_queue;
-	n = advance_particle_queue_length     ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index advance particle " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = surface_crossing_queue;
-	n = surface_crossing_queue_length     ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index surface crossing " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = collision_queue;
-	n = collision_queue_length            ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index collision " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-}
-*/
 
 double get_time()
 {
@@ -648,35 +331,17 @@ void transport()
 	while (remaining_work > 0) {
     start = get_time();
 
-
 		// Figure out work for this subiteration
 		int n_particles = MAX_PARTICLES_IN_FLIGHT;
 		if( n_particles > remaining_work)
 			n_particles = remaining_work;
 
-		//std::cout << "Initializing particle histories..." << std::endl;
-		// Initialize all histories
-		// TODO: Parallelize
+    #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
+    for (int i = 0; i < n_particles; i++) {
+      initialize_history(particles + i, source_offset + i + 1);
+    }
 
-		/*
-		for( int p = 0; p < mpi::n_procs; p++)
-		{
-			MPI_Barrier(mpi::intracomm);
-			if( p == mpi::rank )
-			{
-			*/
-        #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
-				for (int i = 0; i < n_particles; i++) {
-					initialize_history(particles + i, source_offset + i + 1);
-				}
-				/*
-			}
-		}
-		*/
-
-		//std::cout << "Enqueing particles for XS Lookups..." << std::endl;
 		// Add all particles to advance particle queue
-		// TODO: Parallelize
      #pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
 		for (int i = 0; i < n_particles; i++) {
 			dispatch_xs_event(i);
@@ -688,62 +353,32 @@ void transport()
 		int event_kernel_executions = 0;
 		while (true) {
 			event_kernel_executions++;
-      /*
-			std::cout << "Fuel XS Lookups = " << calculate_fuel_xs_queue_length << std::endl;
-			std::cout << "Non Fuel XS Lookups = " << calculate_nonfuel_xs_queue_length << std::endl;
-			std::cout << "Advance Particles = " << advance_particle_queue_length << std::endl;
-			std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl;
-			std::cout << "Collisions = " << collision_queue_length << std::endl;
-      */
-			/*
-			Particle * p = particles +  1;
-		   std::cout << "E = " << p->E_ << " and Position {" <<
-			   p->r().x << ", " <<
-			   p->r().y << ", " <<
-			   p->r().z << "}" << std::endl;
-			   */
 			int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length});
-			//check_energies();
 			if (max == 0) {
 				break;
 			} else if (max == calculate_fuel_xs_queue_length) {
-				//std::cout << "pre fuel XS check..." << std::endl;
-				//check_energies(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
-				//std::cout << "Performing Fuel XS Lookups..." << std::endl;
 				start = get_time();
 				process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
 				stop = get_time();
 				time_fuel_xs += (stop-start);
 				calculate_fuel_xs_queue_length = 0;
 			} else if (max == calculate_nonfuel_xs_queue_length) {
-				//std::cout << "pre non fuel XS check..." << std::endl;
-				//check_energies(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
-			//	std::cout << "Performing Non Fuel XS Lookups..." << std::endl;
 				start = get_time();
 				process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
 				stop = get_time();
 				time_nonfuel_xs += (stop-start);
 				calculate_nonfuel_xs_queue_length = 0;
 			} else if (max == advance_particle_queue_length) {
-				//std::cout << "pre advancing check..." << std::endl;
-				//check_energies(advance_particle_queue, advance_particle_queue_length);
-				//std::cout << "Advancing Particles..." << std::endl;
 				start = get_time();
 				process_advance_particle_events();
 				stop = get_time();
 				time_advance += (stop-start);
 			} else if (max == surface_crossing_queue_length) {
-				//std::cout << "pre surface crossing check..." << std::endl;
-				//check_energies(surface_crossing_queue, surface_crossing_queue_length);
-				//std::cout << "Surface Crossings..." << std::endl;
 				start = get_time();
 				process_surface_crossing_events();
 				stop = get_time();
 				time_surf += (stop-start);
 			} else if (max == collision_queue_length) {
-				//std::cout << "pre Colliding check..." << std::endl;
-				//check_energies(collision_queue, collision_queue_length);
-				//std::cout << "Colliding..." << std::endl;
 				start = get_time();
 				process_collision_events();
 				stop = get_time();