[update] all the way to ganmma yield

analysis/etc/ge_count_smoothing.yaml

+#define the width of the moving window used for smoothing of the ngsub_noneg hsitpogram
+
+side: 60 

analysis/etc/ge_region_of_interest.yaml

+
+min_gamma_energy: 20.0
+max_gamma_energy: 1199.0
+
+min_time: 400.0
+max_time: 6799.0
\ No newline at end of file

analysis/etc/ge_time_to_energy.yaml

+# parameters for time to energy converion
+
+distance_of_flight: 28.7998
+tof_limits: [400., 10000.]

analysis/etc/transitions_to_look_at.yaml

+# List all transitions to look at
+inel259: 
+   energy: 259.48
+   name: L13L05
+   signal: [258.0, 261.0]
+   background: [[257, 258], [261, 262]]
+   parasits: []
+
+
+    
+inel209:
+    energy: 209.8
+    name: L06LL02
+    signal: [209.01, 210.6001]
+    background: [[206.001, 209.001], [210.7, 212.001]]
+    parasits: []
+
+inel84:
+    energy: 84.7
+    name: L05L03
+    signal: [84.0200001, 85.890000001]
+    background: [[83.3000001, 83.70000001], [86.0000001, 86.900000000009]]
+    parasits: []
+    
+n2n229:
+    energy: 229.3
+    name: 182L02L01
+    signal: [228.5, 230.0]
+    background: [[227.5, 228.2], [230.5, 231.0]]
+    parasits: []
+    
+n2n100:
+    energy: 100.1
+    name: 182L01L00
+    signal: [99.7, 101.0]
+    background: [[101.1001, 101.701], [99.001, 99.55001]]
+    parasits: []
+    
+n2n351:
+    energy: 350.99
+    name: 182L03L02
+    signal: [350.0001, 352.0001]
+    background: [[348.000, 349.50001], [352.3001, 252.90001]]
+    parasits: []
+    
+n2n1121:
+    energy: 1121.29
+    name: 182L06L01
+    signal: [119.50001, 1123.750001]
+    background: [[1117.4001, 1118.700012], [1124.3001, 1126.60005]]
+    parasits: []
+    
+
+what122:
+    energy: 122.6
+    name: what122
+    signal: [121.70001, 123.673001]
+    background: [[120.00001, 121.396001], [124.800001, 126.000001]]
+    parasits: []
\ No newline at end of file

analysis/run_on_cc.sh

@@ -37,9 +37,12 @@ echo "# Running tasks"
 #./run_task tasks/do_ge_timealign.py
 #echo "#... ge_calibrate"
 #./run_task tasks/do_ge_calibrate.py
-echo "#... ge_calibrated_2D"
-./run_task tasks/do_ge_calibrate_2D.py
-
+#echo "#... ge_calibrated_2D"
+#./run_task tasks/do_ge_calibrate_2D.py
+#echo "#... cut the bidim to the interesting region"
+#./run_task tasks/do_ge_inel.py
+echo "#... getting the transitions yields"
+./run_task tasks/do_extract_transitions_gnum.py
 
 echo "# Removing job venv"
 rm -rf `hostname`_env

analysis/scripts/min1D.py

+#!/usr/env python
+# -*- coding: utf-8 -*-
+# -*- format: python -*-
+# -*- author: G. Henning -*-
+# -*- created: Tue Mar 13 14:09:01 CET 2018 -*-
+# -*- copyright: GH/IPHC 2018 -*-
+# -*- file: add1D -*-
+# -*- purpose: -*-
+
+'''
+provide one function to impose a minimum value to all bins in an 1d histogram
+'''
+
+import sys,os
+
+import argparse
+
+from pyhisto.histogram import Histogram as histo
+
+def min1D(h: histo,
+                min_val: float = 0.0) -> histo:
+    '''Return an histogram with all bin counts equal to min(min_vval, bin_count)
+    
+    params:
+    - h (pyhisto.Histogram), required
+    - min_val (float), default = 0.0
+    '''
+    new_h = h.copy()
+    for b in new_h:
+        b.count = max(min_val, b.count)
+    new_h._frombins(new_h.bins)
+    return new_h
+
+if __name__=="__main__":
+    parser = argparse.ArgumentParser(description='calibrate a 1D hisotgrams')
+    parser.add_argument('--min_val', type=float, nargs='?',
+                        default=0.0,
+                        help="minimum value")
+    parser.add_argument('file', type=str,
+                        nargs='?',
+                        help="histogram.1d.txt file to process")
+    args = parser.parse_args()
+    if os.path.exists(args.file) and os.path.isfile(args.file):
+        h = min1D(histo(fromfile=args.file),
+                        args.min_val)
+        print(h)
+
+# EOF

analysis/scripts/smooth_1D.py

+#!/usr/env python
+# -*- coding: utf-8 -*-
+# -*- format: python -*-
+# -*- author: G. Henning -*-
+# -*- created: Tue Mar 13 14:09:01 CET 2018 -*-
+# -*- copyright: GH/IPHC 2018 -*-
+# -*- file: add1D -*-
+# -*- purpose: -*-
+
+'''
+provide one function to calibrate a 1D histogram
+'''
+
+import sys,os
+
+import argparse
+
+from pyhisto.histogram import Histogram as histo
+
+
+
+def smooth_1D(h: histo,
+              side: int = 1) -> histo:
+    '''smoothy_1D : smooth the bin content using a moving window average
+    
+    params:
+    - h (pyhisto.Histogram), required
+    - side (float), default = 1: number of bin to average over on each side.
+                                 The total size of the moving window will be 2 * side + 1
+    '''
+    new_h = h.copy()
+    len_h = len(h)
+    for i, b in enumerate(h):
+        first = max(0, i - side)
+        last = min(len_h, i + side + 1)
+        width = float(last-first)
+        new_h[i].count = 1. / width * sum(h[first:last])
+        
+    #enf for enumerate bind
+    new_h._frombins(new_h.bins)
+    new_h *= sum(h)/sum(new_h)
+    return new_h
+
+if __name__=="__main__":
+    parser = argparse.ArgumentParser(description='calibrate a 1D hisotgrams')
+    parser.add_argument('--side', type=int, nargs='?',
+                        default=1,
+                        help="Side (half width) of the average window")   
+    parser.add_argument('file', type=str,
+                        help="histogram.1d.txt file to read and calibrate")
+    args = parser.parse_args()
+    if os.path.exists(args.file) and os.path.isfile(args.file):
+        h = smooth_1D(histo(fromfile=args.file),
+                              side=args.side)
+        print(h)
+
+# EOF

analysis/scripts/sub2D.py

+#!/usr/env python
+# -*- coding: utf-8 -*-
+# -*- format: python -*-
+# -*- author: G. Henning -*-
+# -*- created: Tue Mar 13 14:09:01 CET 2018 -*-
+# -*- copyright: GH/IPHC 2018 -*-
+# -*- file: add1D -*-
+# -*- purpose: -*-
+
+'''
+Module docstring
+'''
+import sys,os
+
+import argparse
+
+from pyhisto.histogram2D import Histogram2D as histo2D
+
+def sub2D(h2_in: histo2D,
+          xmin: float = -1.0, xmax: float = -1.0,
+          ymin: float = -1.0, ymax: float = -1.0,
+         ) -> histo2D:
+
+    # Retrun a new histogram 
+    return h2_in.subhisto(xmin, xmax, ymin, ymax)
+
+if __name__=="__main__":
+    parser = argparse.ArgumentParser(description='calibrate a 2D hisotgrams')
+    parser.add_argument('--xmin', type=float, nargs='?',
+                        default=-1.,
+                        help="X min")
+    parser.add_argument('--xmax', type=float, 
+                        default=-1., nargs='?',
+                        help="X max")
+    parser.add_argument('--ymin', type=float, nargs='?',
+                        default=-1.,
+                        help="Y min")
+    parser.add_argument('--ymax', type=float, 
+                        default=-1., nargs='?',
+                        help="Y max")
+    parser.add_argument('file', type=str,
+                        help="histogram.2d.txt file to read ")
+    args = parser.parse_args()
+    if os.path.exists(args.file) and os.path.isfile(args.file):
+        print(sub2D(histo2D(fromstring=open(args.file, 'r').read()),
+                          args.xmin, args.xmax,
+                          args.ymin, args.ymax))
+    sys.exit()
+
+# EOF

analysis/scripts/time_to_energy_1D.py

+#!/usr/env python
+# -*- coding: utf-8 -*-
+# -*- format: python -*-
+# -*- author: G. Henning -*-
+# -*- created: Tue Mar 13 14:09:01 CET 2018 -*-
+# -*- copyright: GH/IPHC 2018 -*-
+# -*- file: add1D -*-
+# -*- purpose: -*-
+
+'''
+provide one function to calibrate a 1D histogram
+'''
+
+import sys,os
+
+import argparse
+
+from math import sqrt
+
+from pyhisto.histogram import Histogram as histo
+
+SPEED_OF_LIGHT = 0.299792458 # meters / nanoseconds 
+MASS_OF_NEUTRON = 939565.4133 # keV/c^2
+
+def time_to_energy_func(tof: float, dof: float, mass: float) -> float:
+    beta = (dof / tof) / SPEED_OF_LIGHT
+    gamma = 1.0 / sqrt(1.0 - beta ** 2.)
+    E = MASS_OF_NEUTRON * (gamma - 1.0)
+    return E
+
+def time_to_energy_1D(h: histo,
+                      dof: float = 1.0,
+                      mass: float = MASS_OF_NEUTRON,
+                      tmin: float = 0.0, tmax: float = 1e30) -> histo:
+    '''time_to_energy_1D : convert histogram bin edges from time (ns) to energy (keV)
+    
+    params:
+    - h (pyhisto.Histogram), required
+    - dof (float), default = 1.0: Distance of fligth in meters
+    - mass (float), default = MASS_OF_NEUTRON: mass of the particle in keV/c^2 (default = neutron) 
+    - tmin (float), default = 0.0: minimum time to process, will ignore others
+    - tmax (float), default = 1e30: maximum time to process
+    
+    return an histogram with all bins recalibrated (and reordered).
+    '''
+    new_h = h.copy()
+    for i, b in enumerate(h):
+        if (min(b.lowedge, b.upedge) < tmin
+            or max(b.lowedge, b.upedge)> tmax):
+            continue 
+        new_h[i].upedge = time_to_energy_func(b.lowedge, dof, mass)
+        new_h[i].lowedge = time_to_energy_func(b.upedge, dof, mass)
+    #enf for enumerate bind
+    new_h._frombins(new_h.bins)
+    return new_h
+
+if __name__=="__main__":
+    parser = argparse.ArgumentParser(description='calibrate a 1D hisotgrams')
+    parser.add_argument('--dof', type=float, nargs='?',
+                        default=1.0,
+                        help="Distance of flight in meters")
+    parser.add_argument('--mass', type=float, 
+                        default=MASS_OF_NEUTRON, nargs='?',
+                        help="Particle mass in keV/c^2 (default=neutron)")
+    parser.add_argument('--tmin', type=float, 
+                        default=0., nargs='?',
+                        help="Minimum time (in ns)")
+    parser.add_argument('--tmax', type=float, 
+                        default=1e30, nargs='?',
+                        help="Maximum time (in ns)")    
+    parser.add_argument('file', type=str,
+                        help="histogram.1d.txt file to read and calibrate")
+    args = parser.parse_args()
+    if os.path.exists(args.file) and os.path.isfile(args.file):
+        h = time_to_energy_1D(histo(fromfile=args.file),
+                              dof=args.dof, mass=args.mass,
+                              tmin=args.tmin, tmax=args.tmax)
+        print(h)
+
+# EOF

analysis/tasks/do_extract_transitions_gnum.py

+'''
+This file uses energy time and energy calibration coefficients and performs the calibration on bidims
+'''
+
+
+from doit.tools import create_folder
+from doit.action import CmdAction
+import yaml
+import pathlib
+from itertools import product, chain
+
+from pyhisto import Histogram
+
+#
+# Loading the mapping
+#
+
+the_mapping = yaml.load(open('etc/mapping.yaml'),
+                        Loader=yaml.Loader)
+detectors = the_mapping['detectors']
+Ge_detectors = list(filter(lambda x: not x.startswith('U'), detectors))
+
+the_transitions = yaml.load(open('etc/transitions_to_look_at.yaml'),
+                        Loader=yaml.Loader)
+
+
+def task_project_transitions_signal():
+    '''project transitions bidim according to signal'''
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt"
+                dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.signal.h1.txt"
+            
+                yield {
+                    'name': dest_file,
+                    'file_dep': [source_file],
+                    'actions': [
+                        (create_folder, (f"./output/transitions/{this_transition}/",)),
+                        " ".join(["./env.run scripts/project2D.py",
+                                  "--axis=x",
+                                  f"--min={the_transitions[this_transition]['signal'][0]}",
+                                  f"--max={the_transitions[this_transition]['signal'][1]}",
+                                  source_file,
+                                  f"> {dest_file}"])
+                    ],
+                    'targets' : [dest_file],
+                }
+                
+def task_project_transitions_bg():
+    '''project transitions bidim according to background'''
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt"
+                
+                for bg_num, bg_region in enumerate(the_transitions[this_transition]['background']):
+                    dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bg_{bg_num}.h1.txt"           
+                    yield {
+                        'name': dest_file,
+                        'file_dep': [source_file],
+                        'actions': [
+                            (create_folder, (f"./output/transitions/{this_transition}/",)),
+                            " ".join(["./env.run scripts/project2D.py",
+                                      "--axis=x",
+                                      f"--min={bg_region[0]}",
+                                      f"--max={bg_region[1]}",
+                                      source_file,
+                                      f"> {dest_file}"])
+                    ],
+                    'targets' : [dest_file],
+                }
+
+def task_compute_factors():
+    '''compute the factors to apply to the backgrounds to do the bg substraction'''
+    def compute_factor(src_file, signal, BGs, dest_file):
+        source_h = Histogram(fromfile=src_file)
+        nbin_signal = len(source_h.slice(signal[0],signal[1]))
+        nbin_bg = sum([len(source_h.slice(bg_region[0],bg_region[1])) for bg_region in BGs])
+        open(dest_file, 'w').write(f"bg: {nbin_signal/nbin_bg} \n")
+        return True
+    
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt_py.txt"
+                dest_file = f"./output/transitions/{this_transition}/factor_{ge_det}_{flavor}.yaml"
+                yield {
+                    'name': dest_file,
+                    'file_dep': [source_file],
+                    'targets': [dest_file],
+                    'actions': [
+                        (compute_factor, [], {'src_file': source_file, 
+                                              'signal': the_transitions[this_transition]['signal'],
+                                              'BGs': the_transitions[this_transition]['background'],
+                                              'dest_file': dest_file}),
+                    ]
+                }
+                
+def task_bg_remove_h1():
+    '''remove the background from h1'''
+    def substract_bg(signal, bg_files, factor_file, dest_file):
+        h_signal = Histogram(fromfile=signal)
+        bg_factor = yaml.load(open(factor_file),
+                              Loader=yaml.Loader)['bg']
+        for bgf in bg_files:
+            h_bg = Histogram(fromfile=bgf)
+            h_bg *= bg_factor
+            h_signal -= h_bg
+        open(dest_file, 'w').write(str(h_signal))
+        return True
+        
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.signal.h1.txt"
+                
+                n_bg = len(the_transitions[this_transition]['background'])
+                bg_files = list([f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bg_{bg_num}.h1.txt" for bg_num in range(n_bg)])
+                
+                factor_file = f"./output/transitions/{this_transition}/factor_{ge_det}_{flavor}.yaml"
+                
+                dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub.h1.txt"
+                
+                yield {
+                    'name': source_file,
+                    'file_dep': [source_file, factor_file] + bg_files,
+                    'targets': [dest_file],
+                    'actions': [
+                        (substract_bg, [], {'signal': source_file, 
+                                            'bg_files': bg_files,
+                                            'factor_file': factor_file,
+                                            'dest_file': dest_file}),
+                    ]
+                }
+                
+def task_nonegative_counts():
+    '''set negative bins to 0'''
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub.h1.txt"
+                dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub_noneg.h1.txt"
+                yield {
+                    'name': source_file,
+                    'file_dep': [source_file],
+                    'targets': [dest_file],
+                    'actions': [
+                        f"./env.run scripts/min1D.py {source_file} > {dest_file}",
+                    ]
+                }
+
+                
+def task_smooth_gcount():
+    '''apply smoothing to data'''
+    def _smooth_gcount_cmd(input_h1_file, smooth_param_yaml, target_h):
+        '''return command string'''
+        smooth_params = yaml.load(open(smooth_param_yaml), Loader=yaml.Loader)
+        smooth_side = smooth_params.get('side', 1)
+        return " ".join([
+            "./env.run scripts/smooth_1D.py",
+            f"--side={smooth_side}", 
+            f"{input_h1_file}",
+            f"> {target_h}"
+        ])
+    
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub_noneg.h1.txt"
+                dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub_noneg_smooth.h1.txt"
+                
+                yield {
+                    'name': source_file,
+                    'file_dep': [source_file, "etc/ge_count_smoothing.yaml"],
+                    'targets': [dest_file],
+                    'actions': [
+                        CmdAction((_smooth_gcount_cmd, 
+                                  [], {
+                                      'input_h1_file': source_file,
+                                      'smooth_param_yaml': "etc/ge_time_to_energy.yaml",
+                                      'target_h': dest_file
+                                  }))
+                    ]
+                }
+        
+def task_gcount_time_to_energy():
+    '''convert an histogram in time into an histogram in energy'''
+    def _gcount_time_to_energy_cmd(source_h, t2E_param_file, target_h):
+        '''returns the comannd string'''
+        t2E_params = yaml.load(open(t2E_param_file), Loader=yaml.Loader)
+        distance_of_flight = t2E_params['distance_of_flight']
+        tmin, tmax = t2E_params['tof_limits']
+        return " ".join(["./env.run scripts/time_to_energy_1D.py",
+                         f"--dof={distance_of_flight}", 
+                         f"--tmin={tmin} --tmax={tmax}",
+                         f"{source_h}",
+                         f"> {target_h}"])
+    for this_transition in the_transitions.keys():
+        for ge_det in Ge_detectors:
+            for flavor in ('all', 'clean'):
+                source_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.bgsub_noneg_smooth.h1.txt"
+                dest_file = f"./output/transitions/{this_transition}/{ge_det}_{flavor}.gamcount_perenergy.h1.txt"
+                
+                yield {
+                    'name': source_file,
+                    'file_dep': [source_file, "etc/ge_time_to_energy.yaml"],
+                    'targets': [dest_file],
+                    'actions': [
+                        CmdAction((_gcount_time_to_energy_cmd, 
+                                  [], {
+                                      'source_h': source_file,
+                                      't2E_param_file': "etc/ge_time_to_energy.yaml",
+                                      'target_h': dest_file
+                                  }))
+                    ]
+                }
+                
+                
+def task_draw_1D():
+    '''draw all the h1'''
+    for one_histo in pathlib.Path().glob('output/transitions/*/*h1.txt'):
+        yield {
+            'name': str(one_histo),
+            'file_dep': [one_histo],
+            'actions':[
+                f"./env.run scripts/Draw1Dh.py {str(one_histo)}"
+            ],
+            'targets': [f"{str(one_histo)}.png"]
+        }
+
+                    
\ No newline at end of file

analysis/tasks/do_ge_inel.py

+'''
+This file uses energy time and energy calibration coefficients and performs the calibration on bidims
+'''
+
+
+from doit.tools import create_folder
+import yaml
+import pathlib
+from itertools import product, chain
+
+
+#
+# Loading the mapping
+#
+
+the_mapping = yaml.load(open('etc/mapping.yaml'),
+                        Loader=yaml.Loader)
+detectors = the_mapping['detectors']
+Ge_detectors = list(filter(lambda x: not x.startswith('U'), detectors))
+
+
+
+def task_interestingbidim():
+    '''get only the interesting part of the bidim'''
+    for ge_det in Ge_detectors:
+        for flavor in ('all', 'clean'):
+            source_file = f"./output/calibrated_2D/{ge_det}_merged.h2_{flavor}.txt"
+            dest_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt"
+            region_of_interest = yaml.load(open("etc/ge_region_of_interest.yaml"), Loader=yaml.Loader)
+            yield {
+                'name': dest_file,
+                'file_dep': [source_file, "etc/ge_region_of_interest.yaml"],
+                'actions': [
+                    (create_folder, (f"./output/interesting_2D/",)),
+                    " ".join(["./env.run scripts/sub2D.py",
+                              f"--ymin={region_of_interest['min_gamma_energy']}",
+                              f"--ymax={region_of_interest['max_gamma_energy']}",
+                              f"--xmin={region_of_interest['min_time']}",
+                              f"--xmax={region_of_interest['max_time']}",
+                              source_file,
+                              f"> {dest_file}"])
+                ],
+                'targets' : [dest_file],
+            }
+            
+def task_project_intersting():
+    for ge_det in Ge_detectors:
+        for flavor in ('all', 'clean'):
+            source_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt"
+            yield {
+                'name': source_file,
+                'file_dep': [source_file,],
+                'actions': [
+                    f'./env.run scripts/project2D.py --axis=x {source_file} > {source_file}_px.txt',
+                    f'./env.run scripts/project2D.py --axis=y {source_file} > {source_file}_py.txt',
+                ],
+                'targets' : [f'{source_file}_px.txt', f'{source_file}_py.txt'],
+            }
+
+def task_draw_interesting_bidim():
+    '''prooject and draw the h2'''
+    for ge_det in Ge_detectors:
+        for flavor in ('all', 'clean'):
+            source_file = f"./output/interesting_2D/{ge_det}_interest.h2_{flavor}.txt"
+            yield {
+                'name': source_file,
+                'file_dep': [source_file, f'{source_file}_px.txt', f'{source_file}_py.txt'],
+                'actions': [
+                    f'./env.run scripts/Draw1Dh.py {source_file}_py.txt',
+                    f'./env.run scripts/Draw1Dh.py --yscale=log {source_file}_px.txt',
+                    f'./env.run scripts/Draw2Dh.py {source_file}',
+                ],
+                'targets': [
+                    f'{source_file}_px.txt.png', 
+                    f'{source_file}_py.txt.png',
+                    f'{source_file}.txt.png'
+                ]
+            }
+