Plot dE/dx from rat

plot_stopping_power.py

#!/usr/bin/env python                                                                                              
# Plot the stopping power versus energy, and deposited energy versus distance travelled                               
# Author P G Jones - 13/09/2012 <[email protected]> 
import ROOT
import rat
import time
import sys
import math

def plot_file(file_name):
    """ Plot the stopping power."""
    depositation_hist = ROOT.TH1D("EnergyDepositation", "Energy Depositation", 600, 0.0, 60.0)
    stopping_plot = ROOT.TGraph()
    plot_point = 0
    # Now loop over the event data                                                                               
    for ds, run in rat.dsreader(file_name):
        mc = ds.GetMC()
        deposit_total = 0
        print mc.GetMCTrackCount(), "track(s)",
        for itrack in range(0, mc.GetMCTrackCount()):
            mc_track = mc.GetMCTrack(itrack)
            distance = 0.0 # Distance track has travelled                                                             
            deposit_total += mc_track.GetMCTrackStep(0).GetDepositedEnergy()
            if not mc_track.GetParticleName() == "e-":
		continue
            for istep in range(1, mc_track.GetMCTrackStepCount()):
		mc_track_step = mc_track.GetMCTrackStep(istep)
		distance += mc_track_step.GetLength() / 10.0
		dEdx = mc_track_step.GetDepositedEnergy() / ( mc_track_step.GetLength() / 10.0 ) # Convert to cm
		stopping_plot.SetPoint(plot_point, mc_track_step.GetKE(), dEdx)
		depositation_hist.Fill(distance, mc_track_step.GetDepositedEnergy())
		plot_point += 1
		deposit_total += mc_track_step.GetDepositedEnergy()
	print deposit_total, "MeV deposited."
    stopping_plot.GetXaxis().SetTitle("Energy Kinetic [MeV]")
    stopping_plot.GetYaxis().SetTitle("dE/dx [MeV/cm]")
    depositation_hist.GetXaxis().SetTitle("Distance [cm]")
    depositation_hist.GetYaxis().SetTitle("Total Energy Deposited [MeV]")
    return (stopping_plot, depositation_hist, dedx_hist)