diff --git a/examples/make_binned_analysis.py b/examples/make_binned_analysis.py
index 01f0f91..dccba4b 100755
--- a/examples/make_binned_analysis.py
+++ b/examples/make_binned_analysis.py
@@ -28,6 +28,7 @@ from gammalib import *
 from math import *
 import os
 import glob
+import sys
 
 
 # ============================= #
@@ -82,7 +83,7 @@ def pipeline_v1():
 	sim["emin"].real(emin)
 	sim["emax"].real(emax)
 	sim.execute()
-	print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
+	sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += sim.telapse()
@@ -104,7 +105,7 @@ def pipeline_v1():
 	bin["yref"].real(dec)
 	bin["proj"].string(proj)
 	bin.execute()
-	print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
+	sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += bin.telapse()
@@ -119,15 +120,15 @@ def pipeline_v1():
 	like["caldb"].string(caldb)
 	like["irf"].string(irf)
 	like.execute()
-	print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
+	sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += like.telapse()
 	cpu_seconds  += like.celapse()
 	
 	# Show total times
-	print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
-	print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
+	sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n")
+	sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n")
 
 	# Return
 	return
@@ -181,7 +182,7 @@ def pipeline_v2():
 	sim["emin"].real(emin)
 	sim["emax"].real(emax)
 	sim.run()
-	print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
+	sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += sim.telapse()
@@ -200,7 +201,7 @@ def pipeline_v2():
 	bin["yref"].real(dec)
 	bin["proj"].string(proj)
 	bin.run()
-	print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
+	sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += bin.telapse()
@@ -209,18 +210,18 @@ def pipeline_v2():
 	# Perform maximum likelihood fitting
 	like = ctlike(bin.obs())
 	like.run()
-	print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
+	sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n")
 
 	# Update timing
 	wall_seconds += like.telapse()
 	cpu_seconds  += like.celapse()
 		
 	# Show total times
-	print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
-	print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
+	sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n")
+	sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n")
 
 	# Show model fitting results
-	#print like.obs().models()
+	#sys.stdout.write(like.obs().models()+"\n")
 
 	# Plot counts
 	plot_counts(bin.obs())
@@ -249,9 +250,9 @@ def plot_counts(observations):
 		styles = ['b-', 'g-', 'y-', 'n-']
 
 		# Dump header
-		print ""
-		print "Make plots (using matplotlib):"
-		print "=============================="
+		sys.stdout.write("\n")
+		sys.stdout.write("Make plots (using matplotlib):\n")
+		sys.stdout.write("==============================\n")
 		
 		# Create figure 1
 		plt.figure(1,figsize=(12,6))
@@ -274,7 +275,7 @@ def plot_counts(observations):
 				energy.append(ebounds.elogmean(i).TeV())
 
 			# Create spectrum
-			print "Extract data:"
+			sys.stdout.write("Extract data:\n")
 			counts = [0.0 for i in range(ebounds.size())]
 			for bin in cube:
 				index         = ebounds.index(bin.energy())
@@ -288,10 +289,10 @@ def plot_counts(observations):
 			plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
 
 			# Extract models
-			print "Extract models:"
+			sys.stdout.write("Extract models:\n")
 			sum_model = [0.0 for i in range(ebounds.size())]
 			for k, m in enumerate(observations.models()):
-				print "- "+m.name()
+				sys.stdout.write("- "+m.name()+"\n")
 				model = [0.0 for i in range(ebounds.size())]
 				for bin in cube:
 					index        = ebounds.index(bin.energy())
@@ -322,7 +323,7 @@ def plot_counts(observations):
 			cube = cast_GCTAEventCube(obs.events())
 			
 			# Create offset histogram
-			print "Extract data:"
+			sys.stdout.write("Extract data:\n")
 			nx       = 30
 			doffset2 = 0.01
 			offset2  = [(i+0.5)*doffset2 for i in range(nx)]
@@ -341,10 +342,10 @@ def plot_counts(observations):
 			plt.semilogy(offset2, counts, 'ro', label='data')
 
 			# Extract models
-			print "Extract models:"
+			sys.stdout.write("Extract models:\n")
 			sum_model = [0.0 for i in range(nx)]
 			for k, m in enumerate(observations.models()):
-				print "- "+m.name()
+				sys.stdout.write("- "+m.name()+"\n")
 				model = [0.0 for i in range(nx)]
 				for bin in cube:
 					off   = bin.dir().dist_deg(crab)
@@ -368,7 +369,7 @@ def plot_counts(observations):
 		plt.show()
 
 	except:
-		print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
+		sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n")
 
 	# Return
 	return
@@ -398,10 +399,10 @@ if __name__ == '__main__':
 	around the Crab has been performed.
 	"""
 	# Dump header
-	print "*************************************"
-	print "*    CTA binned analysis scripts    *"
-	print "*************************************"
-	print "... this script will take 1-2 minutes"
+	sys.stdout.write("*************************************\n")
+	sys.stdout.write("*    CTA binned analysis scripts    *\n")
+	sys.stdout.write("*************************************\n")
+	sys.stdout.write("... this script will take 1-2 minutes\n")
 
 	# Remove any existing result files
 	list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
@@ -410,14 +411,14 @@ if __name__ == '__main__':
 			os.remove(file)
 	
 	# Save intermediate results on disk
-	print ""
-	print "Executable analysis pipeline:"
-	print "============================="
+	sys.stdout.write("\n")
+	sys.stdout.write("Executable analysis pipeline:\n")
+	sys.stdout.write("=============================\n")
 	pipeline_v1()
 
 	# Do analysis in memory
-	print ""
-	print "In memory analysis pipeline:"
-	print "============================"
+	sys.stdout.write("\n")
+	sys.stdout.write("In memory analysis pipeline:\n")
+	sys.stdout.write("============================\n")
 	pipeline_v2()
-	print ""
+	sys.stdout.write("\n")
diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py
index ea7cc58..042384d 100755
--- a/examples/show_pull_histogram.py
+++ b/examples/show_pull_histogram.py
@@ -58,11 +58,10 @@ def read_pull(filename, parname):
 			try:
 				index = row.index(parname)
 			except:
-				print 'ERROR: Parameter "'+parname+'" not found in list:'
+				sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n')
 				for p in row:
-					print '       "'+p+'"'
+					sys.stdout.write('       "'+p+'"\n')
 				raise NameError(parname)
-			#print index
 
 		# Handle data rows
 		else:
@@ -73,7 +72,6 @@ def read_pull(filename, parname):
 	
 	# Create numpy array
 	a = np.array(values)
-	#print a
 	
 	# Return array
 	return a
@@ -89,7 +87,7 @@ if __name__ == '__main__':
 	# Print usage information
 	usage = "Usage: show_pull_histogram filename parname [bins]"
 	if len(sys.argv) < 3:
-		print usage
+		sys.stdout.write(usage+"\n")
 		sys.exit()
 
 	# Extract parameters
diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py
index 050d28d..13dcdb8 100755
--- a/examples/show_ts_distribution.py
+++ b/examples/show_ts_distribution.py
@@ -53,10 +53,9 @@ def read_ts(filename, tsname="TS"):
 			try:
 				index = row.index(tsname)
 			except:
-				print 'ERROR: Column "TS" not found in file'
-				print row
+				sys.stdout.write('ERROR: Column "TS" not found in file\n')
+				sys.stdout.write(row+"\n")
 				raise NameError("TS")
-			#print index
 
 		# Handle data rows
 		else:
@@ -227,7 +226,7 @@ if __name__ == '__main__':
 	usage = "Usage: show_ts_distribution filename" \
 	        " [-n bins] [-c column] [-t title] [-p plot]"
 	if len(sys.argv) < 2:
-		print usage
+		sys.stdout.write(usage+"\n")
 		sys.exit()
 
 	# Extract parameters
@@ -257,10 +256,10 @@ if __name__ == '__main__':
 					try:
 						par['value'] = sys.argv[i]
 					except:
-						print usage
+						sys.stdout.write(usage+"\n")
 						sys.exit()
 				else:
-					print usage
+					sys.stdout.write(usage+"\n")
 					sys.exit()
 		
 		# Next item
@@ -274,7 +273,7 @@ if __name__ == '__main__':
 	
 	# Read values from CSV file
 	values = read_ts(filename, tsname=tsname)
-	print len(values), "values read."
+	sys.stdout.write(len(values), "values read.\n")
 	
 	# Show histogram
 	if plot == "pdf":
diff --git a/scripts/cspull.py b/scripts/cspull.py
index 7d55a16..159a590 100755
--- a/scripts/cspull.py
+++ b/scripts/cspull.py
@@ -90,7 +90,7 @@ class cspull(GApplication):
 			pars = GPars(parfile)
 		except:
 			# Signal if parfile was not found
-			print "Parfile "+parfile+" not found. Create default parfile."
+			sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
 			
 			# Create default parfile
 			pars = GPars()
diff --git a/scripts/cssens.py b/scripts/cssens.py
index f5b1888..10c4b4e 100755
--- a/scripts/cssens.py
+++ b/scripts/cssens.py
@@ -101,7 +101,7 @@ class cssens(GApplication):
 			pars = GPars(parfile)
 		except:
 			# Signal if parfile was not found
-			print "Parfile "+parfile+" not found. Create default parfile."
+			sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
 			
 			# Create default parfile
 			pars = GPars()
diff --git a/scripts/cstsdist.py b/scripts/cstsdist.py
index 578e2bf..cd89c5f 100755
--- a/scripts/cstsdist.py
+++ b/scripts/cstsdist.py
@@ -91,7 +91,7 @@ class cstsdist(GApplication):
 			pars = GPars(parfile)
 		except:
 			# Signal if parfile was not found
-			print "Parfile "+parfile+" not found. Create default parfile."
+			sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
 			
 			# Create default parfile
 			pars = GPars()
diff --git a/scripts/cta_make_nodes.py b/scripts/cta_make_nodes.py
index 4fa9384..6100304 100755
--- a/scripts/cta_make_nodes.py
+++ b/scripts/cta_make_nodes.py
@@ -86,7 +86,8 @@ def create_nodes(emin, emax, enumbins):
 		write_node(file, energy, scale)
 		
 		# Debug
-		#print energy, math.pow(10.0, int(math.log10(scale))-1.0)
+		#sys.stdout.write(energy)
+        #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n")
 	
 	# Close file
 	file.close()
@@ -105,7 +106,7 @@ if __name__ == '__main__':
 	# Check command line
 	usage = "Usage: cta_make_nodes emin emax enumbins"
 	if len(sys.argv) < 3:
-		print usage
+		sys.stdout.write(usage+"\n")
 		sys.exit()
 	
 	# Extract parameters
diff --git a/scripts/cta_root2irf.py b/scripts/cta_root2irf.py
index 8600989..5afebab 100755
--- a/scripts/cta_root2irf.py
+++ b/scripts/cta_root2irf.py
@@ -149,7 +149,7 @@ def make_files(filename, prfname, bkgname):
 	
 	# Write header
 	fprf.write("log(E) Area  r68  r80  ERes. BG Rate  Diff Sens\n")
-	print "log(E) Area  Area80  r68  r80  ERes. BG Rate  Diff Sens  BGRateSqDeg  BGInt  (BGIntControl)"
+	#sys.stdout.write("log(E) Area  Area80  r68  r80  ERes. BG Rate  Diff Sens  BGRateSqDeg  BGInt  (BGIntControl)\n")
 
 	# Get histograms. Konrad's root files do not have an EffectiveArea80
 	# histogram, but his effective area is for 80% containment radius
@@ -222,8 +222,8 @@ def make_files(filename, prfname, bkgname):
 		fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
 
 		# Show performance file
-		print "%.1f  %.1f  %.1f  %.4f  %.4f %.4f %.7f  %.5e  %.5e %.5e (%.5e)" % \
-		      (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
+		#print "%.1f  %.1f  %.1f  %.4f  %.4f %.4f %.7f  %.5e  %.5e %.5e (%.5e)" % \
+		#      (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
 
 	# Write trailer
 	fprf.write("---------------------------------------------\n")
@@ -285,7 +285,7 @@ if __name__ == '__main__':
 		# Build output filenames
 		prfname = name + ".dat"
 		bkgname = "bkg_"+name+".txt"
-		print prfname, bkgname
+		#print prfname, bkgname
 
 		# Create files
 		make_files(file, prfname, bkgname)
diff --git a/test/test_python.py b/test/test_python.py
index c5a496a..19c5302 100755
--- a/test/test_python.py
+++ b/test/test_python.py
@@ -177,10 +177,10 @@ if __name__ == '__main__':
 	Test ctools Python modules.
 	"""
 	# Dump header
-	print ""
-	print "*************************************"
-	print "* ctools Python interface testing *"
-	print "************************************"
+	sys.stdout.write("\n")
+	sys.stdout.write("***********************************\n")
+	sys.stdout.write("* ctools Python interface testing *\n")
+	sys.stdout.write("***********************************\n")
 	
 	# Set PFILES environment variable
 	try:
@@ -197,11 +197,11 @@ if __name__ == '__main__':
 			os.remove(file)
 	
 	# Save intermediate results on disk
-	print "Test executable analysis: ",
+	sys.stdout.write("Test executable analysis: ")
 	pipeline_v1()
-	print " ok"
+	sys.stdout.write(" ok\n")
 
 	# Keep results in memory
-	print "Test in-memory analysis: ",
+	sys.stdout.write("Test in-memory analysis: ")
 	pipeline_v2()
-	print " ok"
+	sys.stdout.write(" ok\n")