Hi yt, When making projection plots of total gas density and cold gas density, averaged over the projection length (i.e. with a uniform weighting), and using cut_region() to filter out hot gas to make the cold gas plots, it appears that some areas in the cold gas plot have a higher density than the corresponding areas in the total gas plot, which shouldn't happen. If you look at the top-middle of the two plots (link below) and compare the densities, the cold gas appears to have a higher density there. The code is below: (link to plots: https://drive.google.com/drive/folders/1UJ6fmKNFQzGr8cfjov5-GwOzfS8nt2e8?usp...) import yt ds = yt.load("/tigress/cen/CZ15/redshift0100") z = 0.62 # the box size is 120/0.7/(1+z) Mpc, where z is redshift boxSize = 120e3/0.7/(1+z) # virial radius in code units, r_v is virial radius in kpc codeWidth = r_v / boxSize center = [xCoord, yCoord, zCoord] left_edge = [xCoord-2.1*codeWidth, yCoord-2.1*codeWidth, zCoord-2.1*codeWidth] right_edge = [xCoord+2.1*codeWidth, yCoord+2.1*codeWidth, zCoord+2.1*codeWidth] # want a box with full width about 4 times the virial radius reg = ds.region(center, left_edge,right_edge) # total gas plot prj = yt.ProjectionPlot(ds, "x", "density", center=center, width=(4*r_v, 'kpc'), data_source=reg, method='integrate', weight_field='ones') u = 100 / 0.048 l = 0.07 / 0.048 prj.set_zlim('density', l,u) prj.set_unit('density','mbd') prj.set_cmap(field='density', cmap='CMRmap') prj.show() # cold gas plot cold_reg = reg.cut_region(["obj['temperature'] < 1e5"]) prj = yt.ProjectionPlot(ds, "x", "density", center=center, width=(4*r_v, 'kpc'), data_source=cold_reg, method='integrate', weight_field='ones') prj.set_zlim('density', l,u) prj.set_unit('density','mbd') prj.set_cmap(field='density', cmap='CMRmap') prj.show() Thanks, Andrew