How to use the statistics.varf2d_MPI function in statistics

total_iter += 1
		avg = EMAN2_cppwrap.EMData(nx, ny, 1, False)
		for im in data:  EMAN2_cppwrap.Util.add_img(avg, im)

		utilities.reduce_EMData_to_root(avg, myid, main_node)

		if myid == main_node:
			if CTF:
				tavg = EMAN2_cppwrap.Util.divn_filter(avg, ctf_2_sum)
			else:	 tavg = EMAN2_cppwrap.Util.mult_scalar(avg, 1.0/float(nima))
		else:
			tavg = EMAN2_cppwrap.EMData(nx, ny, 1, False)                               

		if Fourvar:
			utilities.bcast_EMData_to_all(tavg, myid, main_node)
			vav, rvar = statistics.varf2d_MPI(myid, data, tavg, mask, "a", CTF)

		if myid == main_node:
			if Fourvar:
				tavg    = fundamentals.fft(EMAN2_cppwrap.Util.divn_img(fundamentals.fft(tavg), vav))
				vav_r	= EMAN2_cppwrap.Util.pack_complex_to_real(vav)

			# normalize and mask tavg in real space
			tavg = fundamentals.fft(tavg)
			stat = EMAN2_cppwrap.Util.infomask( tavg, mask, False ) 
			tavg -= stat[0]
			EMAN2_cppwrap.Util.mul_img(tavg, mask)
			# For testing purposes: shift tavg to some random place and see if the centering is still correct
			#tavg = rot_shift3D(tavg,sx=3,sy=-4)
			tavg = fundamentals.fft(tavg)

		if Fourvar:  del vav

total_iter += 1
		avg = EMData(nx, ny, 1, False)
		for im in data:  Util.add_img(avg, im)

		reduce_EMData_to_root(avg, myid, main_node)

		if myid == main_node:
			if CTF:
				tavg = Util.divn_filter(avg, ctf_2_sum)
			else:	 tavg = Util.mult_scalar(avg, 1.0/float(nima))
		else:
			tavg = EMData(nx, ny, 1, False)                               

		if Fourvar:
			bcast_EMData_to_all(tavg, myid, main_node)
			vav, rvar = varf2d_MPI(myid, data, tavg, mask, "a", CTF)

		if myid == main_node:
			if Fourvar:
				tavg    = fft(Util.divn_img(fft(tavg), vav))
				vav_r	= Util.pack_complex_to_real(vav)

			# normalize and mask tavg in real space
			tavg = fft(tavg)
			stat = Util.infomask( tavg, mask, False ) 
			tavg -= stat[0]
			Util.mul_img(tavg, mask)
			# For testing purposes: shift tavg to some random place and see if the centering is still correct
			#tavg = rot_shift3D(tavg,sx=3,sy=-4)
			tavg = fft(tavg)

		if Fourvar:  del vav

total_iter += 1
		avg = EMAN2_cppwrap.EMData(nx, ny, 1, False)
		for im in range(ldata):
			EMAN2_cppwrap.Util.add_img(avg, fundamentals.fshift(data[im], shift_x[im]))

		utilities.reduce_EMData_to_root(avg, myid, main_node)

		if myid == main_node:
			if CTF:  tavg = EMAN2_cppwrap.Util.divn_filter(avg, ctf_2_sum)
			else:    tavg = EMAN2_cppwrap.Util.mult_scalar(avg, 1.0/float(nima))
		else:
			tavg = utilities.model_blank(nx,ny)

		if Fourvar:
			utilities.bcast_EMData_to_all(tavg, myid, main_node)
			vav, rvar = statistics.varf2d_MPI(myid, data, tavg, mask, "a", CTF)

		if myid == main_node:
			if Fourvar:
				tavg    = fundamentals.fft(EMAN2_cppwrap.Util.divn_img(fundamentals.fft(tavg), vav))
				vav_r	= EMAN2_cppwrap.Util.pack_complex_to_real(vav)
			# normalize and mask tavg in real space
			tavg = fundamentals.fft(tavg)
			stat = EMAN2_cppwrap.Util.infomask( tavg, mask, False )
			tavg -= stat[0]
			EMAN2_cppwrap.Util.mul_img(tavg, mask)
			tavg.write_image("tavg.hdf",Iter)
			# For testing purposes: shift tavg to some random place and see if the centering is still correct
			#tavg = rot_shift3D(tavg,sx=3,sy=-4)

		if Fourvar:  del vav
		utilities.bcast_EMData_to_all(tavg, myid, main_node)