Merge pull request #55 from jbrage/feature/43-reorganize-hadrons

Feature/43 reorganize electrons

electrons/Boag_theory.py

@@ -1,46 +1,47 @@
 import numpy as np
+
 # import matplotlib.pyplot as plt
 
 
 alpha = 1.6e-6  # recombination coefficient [cm^3 / s]
-k_1 = 1.36      # positive charge carriers [cm^2 / (V s)]
-k_2 = 2.1       # negative charge carriers [cm^2 / (V s)]
+k_1 = 1.36  # positive charge carriers [cm^2 / (V s)]
+k_2 = 2.1  # negative charge carriers [cm^2 / (V s)]
 e_charge = 1.60217662e-19  # [C]
 
 
-'''
+"""
 References:
     [1] Boag and Currant (1980) "Current collection and ionic recombination in small cylindrical ionization chambers
         exposed to pulsed radiation"
     [2] Andreo, Burns, Nahum et al, Fundamentals of Ionizing Radiation
     [3] Boutillon (1998) "Volume recombination parameter in ionization chambers"
     [4] Liszka et al (2018) "Ion recombination and polarity correction factors for a plane–parallel ionization chamber
         in a proton scanning beam"
-'''
+"""
 
 
 def Boag_pulsed(Qdensity_C_cm3, d_cm, V):
-    '''
+    """
     The Boag theory (ref [1]) for recombination in a parallel-plate ionization chamber in pulsed beams
 
     The theory assumes:
         - uniform charge carrier distribution
         - uniform electric field
         - enables the inclusion of a free electron component (not included here)
-    '''
+    """
 
     mu = alpha / (e_charge * (k_1 + k_2))
     # ICRU report 34 suggests mu = 3.02e10 V m /C, see ref [2] page 518
 
     r = Qdensity_C_cm3
 
-    u = mu * r * d_cm*d_cm / V
-    f = 1./u * np.log(1 + u)
+    u = mu * r * d_cm * d_cm / V
+    f = 1.0 / u * np.log(1 + u)
     return f
 
 
 def Boag_Continuous(Qdensity_C_cm3_s, d_cm, V):
-    '''
+    """
     The theory (ref [3, 4]) for recombination in a parallel-plate ionization chamber in continuous beams
 
     The theory assumes:
@@ -51,11 +52,11 @@ def Boag_Continuous(Qdensity_C_cm3_s, d_cm, V):
         - f
         - f minus 1 std
         - f plus 1 std
-    '''
+    """
 
     def f_c(mu):
-        ksi_squared = mu * (d_cm**4 / (V*V)) * Qdensity_C_cm3_s
-        return 1./(1 + ksi_squared)
+        ksi_squared = mu * (d_cm**4 / (V * V)) * Qdensity_C_cm3_s
+        return 1.0 / (1 + ksi_squared)
 
     mu_c = alpha / (6 * e_charge * k_1 * k_2)  # [V^2 s / (cm C)]
     # mu_c = 6.73e11 # [V^2 s / (cm C)]

electrons/cython/run_simulation.py

@@ -0,0 +1,69 @@
+from electrons.cython.pulsed_e_beam import pulsed_beam_PDEsolver
+from electrons.cython.continuous_e_beam import continuous_beam_PDEsolver
+
+import argparse
+
+SOLVER_MAP = {"continous": continuous_beam_PDEsolver, "pulsed": pulsed_beam_PDEsolver}
+
+
+def run_simulation(
+    solver_name="continous",
+    voltage_V=300,
+    electrode_gap=0.1,
+    electron_density_per_cm3=1e9,
+    verbose=True,
+):
+    parameters = {
+        "voltage_V": voltage_V,
+        "d_cm": electrode_gap,
+        "elec_per_cm3": electron_density_per_cm3,
+        "show_plot": False,
+        "print_parameters": False,
+    }
+
+    if solver_name not in SOLVER_MAP.keys():
+        print(f'Invalid solver type "{solver_name}", defaulting to Continous solver.')
+        solver_name = "continous"
+
+    if verbose:
+        print(f"Running the simulation using the {solver_name} solver.")
+        print(f"Voltage: {voltage_V} [V]")
+        print(f"Electrode gap: {electrode_gap} [cm]")
+        print(f"Electron density per cm3: {electron_density_per_cm3}")
+
+    solver = SOLVER_MAP[solver_name]
+
+    # return the collection efficiency
+    result = solver(parameters)
+
+    if solver_name == "continous":
+        f = result[1]["f"]
+    else:
+        f = result
+
+    return f
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "solver_name",
+        type=str,
+        default="continous",
+        help="The type of the solver to use",
+    )
+
+    parser.add_argument(
+        "--verbose",
+        "-v",
+        type=bool,
+        default=True,
+    )
+
+    args = parser.parse_args()
+
+    result = run_simulation(**vars(args))
+
+    print("calculated f is", result)

electrons/example_continuous_electron_beam.py

@@ -7,92 +7,118 @@
 import itertools
 
 
-def IonTracks_continuous(voltage_V=300, electrode_gap_cm=0.1, elec_per_cm3=1e9, show_plot=False, print_parameters=False):
-
+def IonTracks_continuous(
+    voltage_V=300,
+    electrode_gap_cm=0.1,
+    elec_per_cm3=1e9,
+    show_plot=False,
+    print_parameters=False,
+):
     # convert to a dict
-    parameters = {"voltage_V": voltage_V,
-    	          "d_cm": electrode_gap_cm, 
-    	          "elec_per_cm3": elec_per_cm3,
-    	          "show_plot": show_plot,  
-                  "print_parameters" : print_parameters,
-                  }      
+    parameters = {
+        "voltage_V": voltage_V,
+        "d_cm": electrode_gap_cm,
+        "elec_per_cm3": elec_per_cm3,
+        "show_plot": show_plot,
+        "print_parameters": print_parameters,
+    }
 
-    result_df, time_variation_df = continuous_beam_PDEsolver(parameters)    
+    result_df, time_variation_df = continuous_beam_PDEsolver(parameters)
 
     return result_df, time_variation_df
 
 
-
-
 if __name__ == "__main__":
-
-    """"
+    """
     # calculate an example for default parameters
     result_df, time_variation_df = IonTracks_continuous()
     print(result_df)
-    
-    # plot the result to show how the charge collection converges after the time it takes for 
+
+    # plot the result to show how the charge collection converges after the time it takes for
     # a partcile to move between the two electrodes. This is the recombination in the continuous situaiton
     fig, ax = plt.subplots()
     ax.set_title("Charge collection versus time")
     sns.lineplot(ax=ax, data=time_variation_df, x="time_us", y="f", label="Charge collection")
     ax.set_xlabel("time (us)")
-    ax.set_ylabel("Collection efficiency")    
+    ax.set_ylabel("Collection efficiency")
     ax.axvline(x=result_df["convergence_time_s"].values * 1e6, label="Drift time b/w gap", c="r", ls=":")
-    ax.axhline(y=result_df["f"].values, label="Converged efficiency", c="r", ls="--")    
+    ax.axhline(y=result_df["f"].values, label="Converged efficiency", c="r", ls="--")
     ax.ticklabel_format(style="plain")
     ax.legend()
     fig.tight_layout()
-    fig.savefig("fig_charge_collection_versus_time.pdf", bbox_inches="tight")
-"""   
-
+    fig.savefig("fig_charge_collection_versus_time.pdf", bbox_inches="tight")"""
+
     # compare the IonTracks results to the Boag theory for a continuous beam
-    
+
     # set parameters for the IonTracks calculation
     d_cm = 0.1
     charge_density_C_cm3_s = np.asarray([0.02, 0.08]) * 1e-6
     data_dict = dict(
         electrode_gap_cm=[d_cm],
-        elec_per_cm3_s = charge_density_C_cm3_s/ e_charge,
-        voltage_V = np.asarray([60, 120, 200]),
+        elec_per_cm3_s=charge_density_C_cm3_s / e_charge,
+        voltage_V=np.asarray([60, 120, 200]),
     )
-        
+
     # create a data frame with all the variables
-    data_df = pd.DataFrame.from_records(data=itertools.product(*data_dict.values()), columns=data_dict.keys())
+    data_df = pd.DataFrame.from_records(
+        data=itertools.product(*data_dict.values()), columns=data_dict.keys()
+    )
 
     # Calculate the recombination using IonTracks
     IonTracks_df = pd.DataFrame()
     for idx, data in data_df.iterrows():
-        result_df, _ = IonTracks_continuous(voltage_V=data.voltage_V, electrode_gap_cm=data.electrode_gap_cm, elec_per_cm3=data.elec_per_cm3_s)
-        IonTracks_df = pd.concat([IonTracks_df, result_df], ignore_index=True)    
+        result_df, _ = IonTracks_continuous(
+            voltage_V=data.voltage_V,
+            electrode_gap_cm=data.electrode_gap_cm,
+            elec_per_cm3=data.elec_per_cm3_s,
+        )
+        IonTracks_df = pd.concat([IonTracks_df, result_df], ignore_index=True)
     print(IonTracks_df)
 
     # rename for plot
-    IonTracks_df["electron_density"] = IonTracks_df["elec_per_cm3"].map("IonTracks: {:0.2E} e$^{{-}}$/cm$^3$".format)
-
-
+    IonTracks_df["electron_density"] = IonTracks_df["elec_per_cm3"].map(
+        "IonTracks: {:0.2E} e$^{{-}}$/cm$^3$".format
+    )
+
     # generate the data for the Boag theory
-    data_dict["charge_density_C_cm3_s"] = charge_density_C_cm3_s # replace by the charge density, not electron
-    data_dict["voltage_V"] = 10 **np.linspace(np.log10(IonTracks_df.voltage_V.min()), np.log10(IonTracks_df.voltage_V.max()), 100)
-    Boag_df = pd.DataFrame.from_records(data=itertools.product(*data_dict.values()), columns=data_dict.keys())
+    data_dict[
+        "charge_density_C_cm3_s"
+    ] = charge_density_C_cm3_s  # replace by the charge density, not electron
+    data_dict["voltage_V"] = 10 ** np.linspace(
+        np.log10(IonTracks_df.voltage_V.min()),
+        np.log10(IonTracks_df.voltage_V.max()),
+        100,
+    )
+    Boag_df = pd.DataFrame.from_records(
+        data=itertools.product(*data_dict.values()), columns=data_dict.keys()
+    )
 
     # calculate the Boag collection efficiency for each point
     for idx, data in Boag_df.iterrows():
-        f, _, _ = Boag_continuous(Qdensity_C_cm3_s=data.charge_density_C_cm3_s, d_cm=data.electrode_gap_cm, V=data.voltage_V)
+        f, _, _ = Boag_continuous(
+            Qdensity_C_cm3_s=data.charge_density_C_cm3_s,
+            d_cm=data.electrode_gap_cm,
+            V=data.voltage_V,
+        )
         Boag_df.loc[idx, "f"] = f
     Boag_df["ks"] = 1 / Boag_df["f"]
-    Boag_df["charge_density_C_cm3_s"] = Boag_df["charge_density_C_cm3_s"].map("Boag: {:0.2E} C/cm$^3$".format)
-
+    Boag_df["charge_density_C_cm3_s"] = Boag_df["charge_density_C_cm3_s"].map(
+        "Boag: {:0.2E} C/cm$^3$".format
+    )
+
     # plot the results
     fig, ax = plt.subplots()
-
-    sns.lineplot(ax=ax, data=Boag_df, x="voltage_V", y="ks", hue="charge_density_C_cm3_s")
-    sns.scatterplot(ax=ax, data=IonTracks_df, x="voltage_V", y="ks", hue="electron_density")    
+
+    sns.lineplot(
+        ax=ax, data=Boag_df, x="voltage_V", y="ks", hue="charge_density_C_cm3_s"
+    )
+    sns.scatterplot(
+        ax=ax, data=IonTracks_df, x="voltage_V", y="ks", hue="electron_density"
+    )
 
     ax.set_title("Continuous beam: $d = {:0.2g}$ cm air gap".format(d_cm))
     ax.set_xlabel("Voltage [V]")
     ax.set_ylabel("Recombination correction factor ($k_s$)")
     ax.set_xscale("log")
-    ax.legend(loc='best', frameon=False)
+    ax.legend(loc="best", frameon=False)
     fig.savefig("fig_example_continuous.pdf")
-