Mass swapping: fhot using Mvir_infall to compute the BH accretion in …

…hot halo mode and small modifications

src/gas_cooling.cpp

@@ -392,6 +392,7 @@ double GasCooling::cooling_rate(Subhalo &subhalo, Galaxy &galaxy, double z, doub
 	if(subhalo.subhalo_type == Subhalo::SATELLITE && subhalo.Vvir_infall != 0 &&
 	               dark_matter_params.apply_fix_to_mass_swapping_events){
 		vvir = subhalo.Vvir_infall;
+		fhot = mhot / subhalo.Mvir_infall;
 	}
 
 	double zhot = 0;

src/merger_tree_reader.cpp

@@ -145,6 +145,7 @@ const std::vector<SubhaloPtr> SURFSReader::read_subhalos(unsigned int batch)
 	        const auto fname_transients = get_transients_filename(batch);
 		hdf5::Reader batch_transients(fname_transients);
 		transientsIndex = batch_transients.read_dataset_v<Subhalo::id_t>("nodeIndex");
+		//		transientsHostIndex = batch_transients.read_dataset_v<Subhalo::id_t>("hostIndex");
 	}
 
 	auto n_subhalos = Mvir.size();
@@ -307,6 +308,10 @@ const std::vector<HaloPtr> SURFSReader::read_halos(unsigned int batch)
 		}
 		subhalo->host_halo = halo;
 		halo->add_subhalo(std::move(subhalo));
+//		if (!transients_prefix.empty()){
+//			// create flag to indicate this halo is transient
+//		        halo->transient = std::find(std::begin(transientsHostIndex), std::end(transientsHostIndex), hostIndex[i]) != std::end(transientsHostIndex);
+//		}
 	}
 	subhalos.clear();
 
@@ -315,14 +320,14 @@ const std::vector<HaloPtr> SURFSReader::read_halos(unsigned int batch)
 	os << "This should take another ~" << memory_amount(halos.size() * (sizeof(Halo) + sizeof(HaloPtr))) << " of memory";
 	LOG(info) << os.str();
 
-	// Calculate halos' vvir and concentration
-	t = Timer();
-	omp_dynamic_for(halos, threads, 10000, [&](const HaloPtr &halo, unsigned int thread_idx) {
-		auto z = simulation_params.redshifts[halo->snapshot];
-		halo->Vvir = dark_matter_halos->halo_virial_velocity(halo->Mvir, z);
-		halo->concentration = dark_matter_halos->nfw_concentration(halo->Mvir,z);
-	});
-	LOG(info) << "Calculated Vvir and concentration for new Halos in " << t;
+//	// Calculate halos' vvir and concentration
+//	t = Timer();
+//	omp_dynamic_for(halos, threads, 10000, [&](const HaloPtr &halo, unsigned int thread_idx) {
+//		auto z = simulation_params.redshifts[halo->snapshot];
+//		halo->Vvir = dark_matter_halos->halo_virial_velocity(halo->Mvir, z);
+//		halo->concentration = dark_matter_halos->nfw_concentration(halo->Mvir,z);
+//	});
+//	LOG(info) << "Calculated Vvir and concentration for new Halos in " << t;
 
 	return halos;
 }

src/tree_builder.cpp

@@ -222,19 +222,20 @@ SubhaloPtr TreeBuilder::define_central_subhalo(HaloPtr &halo, SubhaloPtr &subhal
 		mvir = halo->Mvir;
 	}
 
-	double z= sim_params.redshifts[subhalo->snapshot];
-	double npart = mvir/sim_params.particle_mass;
+	double z = sim_params.redshifts[subhalo->snapshot];
+	double npart = subhalo->Mvir/sim_params.particle_mass;
 
 	subhalo->concentration = darkmatterhalos->nfw_concentration(mvir, z);
-
+	
 	if (subhalo->concentration < 1) {
-		throw invalid_argument("concentration is <1, cannot continue. Please check input catalogue");
+	        throw invalid_argument("concentration is <1, cannot continue. Please check input catalogue");
 	}
 	subhalo->lambda = darkmatterhalos->halo_lambda(*subhalo, mvir, z, npart);
 	subhalo->Vvir = darkmatterhalos->halo_virial_velocity(mvir, z);
 
-	halo->concentration = subhalo->concentration;
 	halo->lambda = subhalo->lambda;
+	halo->Vvir = darkmatterhalos->halo_virial_velocity(halo->Mvir, z);
+	halo->concentration = darkmatterhalos->nfw_concentration(halo->Mvir,z);
 
 	/** If virial velocity of halo (which is calculated from the total mass
 	and redshift) is smaller than the virial velocity of the central subhalo, which is