Implement histogramming of angles of the line connecting a reference point and a point of impact

default_input_parameters.py

@@ -76,6 +76,10 @@
             'save_mp_state_time_file': -1,
             'flag_detailed_MP_info': 0,
             'flag_hist_impact_seg': 1,
+            'flag_hist_impact_angle' : False,
+            'hist_impact_angle_nbins': 360.,
+            'impact_x_ref': 0.,
+            'impact_y_ref': 0.,
             'flag_En_hist_seg': False,
             'flag_verbose_file': False,
             'flag_verbose_stdout': False,

impact_management_class.py

@@ -63,7 +63,9 @@ def __init__(
             En_hist_max, Nbin_lifetime_hist=None,
             lifetime_hist_max=None, flag_lifetime_hist=False,
             flag_seg=False, flag_En_hist_seg=False,
-            cos_angle_width=0.05, flag_cos_angle_hist=True):
+            cos_angle_width=0.05, flag_cos_angle_hist=True,
+            flag_hist_impact_angle=True, hist_impact_angle_nbins=360.,
+            impact_x_ref=0., impact_y_ref=0.):
 
         print('Start impact man. init.')
 
@@ -100,6 +102,15 @@ def __init__(
         else:
             print('Not saving cosine of angle of incident electrons.')
 
+        # recording of angle impacts (in degrees) with respect to a reference point 
+        self.flag_hist_impact_angle = flag_hist_impact_angle
+        if self.flag_hist_impact_angle:
+            self.impact_x_ref = impact_x_ref
+            self.impact_y_ref = impact_y_ref
+            self.hist_impact_angle = np.zeros(hist_impact_angle_nbins, float)
+            self.hist_impact_angle_bin_width = np.rad2deg(2 * np.pi) / hist_impact_angle_nbins
+            print(f'Saving angle of impact with respect to reference (x,y)=({self.impact_x_ref:.4f}, {self.impact_y_ref:.4f}).')
+
         self.xg_hist = xg_hist
         self.Nxg_hist = Nxg_hist
         self.bias_x_hist = bias_x_hist
@@ -167,6 +178,9 @@ def reset_energ_impact_seg(self):
     def reset_cos_angle_hist(self):
         self.cos_angle_hist *= 0
 
+    def reset_hist_impact_angle(self):
+        self.hist_impact_angle *= 0
+
     def reset_lifetime_hist_line(self):
         self.lifetime_hist_line *= 0.
 
@@ -266,10 +280,18 @@ def backtrack_and_second_emiss(self, old_pos, MP_e, tt_curr=None):
                 histf.compute_hist(x_impact, nel_impact * E_impact_eV,
                                    bias_x_hist, Dx_hist, self.energ_eV_impact_hist)
 
-                # angle histogram
+                # angle histograms
                 if self.flag_cos_angle_hist:
                     histf.compute_hist(
                         costheta_impact, nel_impact, 0., self.cos_angle_width, self.cos_angle_hist)
+                if self.flag_hist_impact_angle:
+                    x_impact_wrt_ref = x_impact - self.impact_x_ref
+                    y_impact_wrt_ref = y_impact - self.impact_y_ref
+                    impact_angle_rad = np.arctan2(y_impact_wrt_ref, x_impact_wrt_ref)
+                    impact_angle = np.rad2deg(np.mod(impact_angle_rad, 2*np.pi))
+                    histf.compute_hist(impact_angle, nel_impact, 0., 
+                                       self.hist_impact_angle_bin_width,
+                                       self.hist_impact_angle)
 
                 if flag_seg:
                     segi.update_seg_impact(

init.py

@@ -415,7 +415,9 @@ def read_input_files_and_init_components(pyecl_input_folder='./', skip_beam=Fals
             thiscloud.Dx_hist, thiscloud.scrub_en_th, cc.Nbin_En_hist, cc.En_hist_max,
             cc.Nbin_lifetime_hist, cc.lifetime_hist_max, cc.flag_lifetime_hist,
             flag_seg=flag_seg, flag_En_hist_seg=thiscloud.flag_En_hist_seg,
-            cos_angle_width=cc.cos_angle_width)
+            cos_angle_width=cc.cos_angle_width, flag_hist_impact_angle=cc.flag_hist_impact_angle,
+            impact_x_ref=cc.impact_x_ref, impact_y_ref=cc.impact_y_ref,
+            hist_impact_angle_nbins=cc.hist_impact_angle_nbins)
 
         # Init gas ionization and photoemission
         if thiscloud.gas_ion_flag == 1:

pyecloud_saver.py

@@ -196,7 +196,9 @@ def start_observing(self, Dt_ref, MP_e, beamtim, impact_man,
 
         # Init energy and cos angle histogram saving
         self._energy_cos_and_lifetime_angle_hist_init(Dt_En_hist, flag_cos_angle_hist, cos_angle_width, flag_lifetime_hist, Dt_lifetime_hist)
-
+
+        # Init impact angle histogram
+        self._impact_angle_hist_init(impact_man)
 
         #Space charge electrostatic energy
         self.t_sc_video = []
@@ -314,6 +316,9 @@ def witness(self, MP_e, beamtim, spacech_ele, impact_man,
             self._checkpoint_save(beamtim, spacech_ele, t_sc_ON, flag_presence_sec_beams,
                                   sec_beams_list, self.flag_multiple_clouds, cloud_list)
 
+            # check for impact angle hist save
+            self._impact_angle_hist_save(impact_man)
+
             self._copy_main_outp_save(beamtim)
 
         if beamtim.flag_new_bunch_pass:
@@ -462,7 +467,9 @@ def build_outp_dict(self, buildup_sim):
                     'dec_fact_out': self.dec_fact_out,
                     'chamber_area': self.area,
                     'cos_angle_hist': self.cos_angle_hist,
-                    'xg_hist_cos_angle': self.xg_hist_cos_angle
+                    'xg_hist_cos_angle': self.xg_hist_cos_angle,
+                    'hist_impact_angle': self.hist_impact_angle,
+                    'xg_hist_impact_angle': self.xg_hist_impact_angle
         }
 
         if self.flag_lifetime_hist:
@@ -1208,3 +1215,20 @@ def _energy_cos_and_lifetime_angle_hist_save(self, beamtim, impact_man, MP_e):
                 self.t_lifetime_hist.append(beamtim.tt_curr)
                 impact_man.reset_lifetime_hist_line()
                 self.t_last_lifetime_hist = beamtim.tt_curr
+
+    def _impact_angle_hist_init(self, impact_man):
+        # Angle histogram
+        self.flag_hist_impact_angle = impact_man.flag_hist_impact_angle 
+        if self.flag_hist_impact_angle:
+            N_angles = len(impact_man.hist_impact_angle)
+            dx = impact_man.hist_impact_angle_bin_width
+            self.hist_impact_angle = []
+            self.xg_hist_impact_angle = np.linspace(dx/2., np.rad2deg(2 * np.pi) - dx/2., N_angles)
+        else:
+            self.hist_impact_angle = -1
+            self.xg_hist_impact_angle = -1
+
+    def _impact_angle_hist_save(self, impact_man):
+            if self.flag_hist_impact_angle:
+                self.hist_impact_angle.append(impact_man.hist_impact_angle.copy())
+                impact_man.reset_hist_impact_angle()