[update] Flux is coming up

analysis/etc/FC_time_align.yaml

-GAMMAFLASH_DETECTION_THRESHOLD_FRACTION: 0.333
\ No newline at end of file
+GAMMAFLASH_DETECTION_THRESHOLD_FRACTION: 0.333
+
+bidim_tmin: -200.0
+bidim_tmax: 14000.0
+bidim_tnbins: 7100
\ No newline at end of file

analysis/etc/fc_energy_spect.yaml

+min_neutron_energy: 5.0
+max_neutron_energy: 35000.0
+
+min_time: 150.0 # nanosec
+
+time_shift: 0.
\ No newline at end of file

analysis/etc/fc_flux_spect.yaml

+min_neutron_energy: 10.0
+max_neutron_energy: 27000.0
+
+number_of_bins: 26990
+
+nf_eval: 'ENDF/B-VIII.0: U-235(N,F)'
\ No newline at end of file

analysis/run_on_cc.sh

@@ -35,6 +35,14 @@ echo "# Running tasks"
 #./run_task tasks/do_raw2histo.py
 echo "#... FC time align"
 ./run_task tasks/do_fc_timealign.py
+echo "#... FC find energy cut"
+./run_task tasks/do_fc_ecut.py
+echo "#... FC get pile up rate"
+./run_task tasks/do_fc_pileup.py
+echo "#... FC cut 2D and project on time"
+./run_task tasks/do_fc_cut_and_project.py
+echo "#... FC -> neutron flux [MC]"
+./run_task tasks/do_fc_to_flux.py
 # echo "#... ge_time_align"
 #./run_task tasks/do_ge_timealign.py
 #echo "#... ge_calibrate"

analysis/run_task

@@ -2,5 +2,5 @@
 source `hostname`_env/bin/activate
 #export PYTHONPATH=/iphc-work/dnr/usr/ghenning/183W/analysis/lib:/iphc-work/dnr/usr/ghenning/183W/analysis/scripts
 
-time ./env.run -m doit --verbosity 2 --reporter executed-only -n 3 -f $@ --dir .
+time ./env.run -m doit --verbosity 2 --reporter executed-only -n 1 -f $@ --dir .
 

analysis/scripts/time_to_energy_1D.py

@@ -28,10 +28,19 @@ def time_to_energy_func(tof: float, dof: float, mass: float) -> float:
     E = MASS_OF_NEUTRON * (gamma - 1.0)
     return E
 
+def energy_to_time_func(en: float, dof:float, mass: float) -> float:
+    EoverM = en / mass
+    lambda_ = 1. / (1. + EoverM) ** 2.0
+    beta = sqrt(1. - lambda_)
+    Tof = dof/beta/SPEED_OF_LIGHT
+    return Tof
+
 def time_to_energy_1D(h: histo,
                       dof: float = 1.0,
                       mass: float = MASS_OF_NEUTRON,
-                      tmin: float = 0.0, tmax: float = 1e30) -> histo:
+                      tmin: float = 0.0, tmax: float = 1e30,
+                      emin: float = 0.0, emax: float = 1e30,
+                      tshift: float = 0.0) -> histo:
     '''time_to_energy_1D : convert histogram bin edges from time (ns) to energy (keV)
     
     params:
@@ -43,15 +52,22 @@ def time_to_energy_1D(h: histo,
     
     return an histogram with all bins recalibrated (and reordered).
     '''
-    new_h = h.copy()
+    
+    tmin_frome = energy_to_time_func(emax, dof, mass)
+    tmax_frome = energy_to_time_func(emin, dof, mass)
+    print(f"# time limits by energy constraints : {tmin_frome}:{tmax_frome}")
+    tmin, tmax = max(tmin, tmin_frome), min(tmax_frome, tmax)
+    true_emin, true_emax = max(emin, time_to_energy_func(tmax + tshift, dof, mass)), min(emax, time_to_energy_func(tmin + tshift, dof, max))
+    new_h = h.empty_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)
+        new_h[i].upedge = time_to_energy_func(h[i].lowedge + tshift, dof, mass)
+        new_h[i].lowedge = time_to_energy_func(h[i].upedge + tshift, dof, mass)
+        new_h[i].count = h[i].count
     #enf for enumerate bind
-    new_h._frombins(new_h.bins)
+    new_h._frombins( list( b for b in new_h.bins if (min(b.lowedge, b.upedge) > true_emin) and (max(b.lowedge, b.upedge) <true_emax)))
     return new_h
 
 if __name__=="__main__":
@@ -67,14 +83,24 @@ if __name__=="__main__":
                         help="Minimum time (in ns)")
     parser.add_argument('--tmax', type=float, 
                         default=1e30, nargs='?',
-                        help="Maximum time (in ns)")    
+                        help="Maximum time (in ns)")   
+    parser.add_argument('--tshift', type=float,
+                        default=0.0, nargs='?',
+                        help="time shift")
+    parser.add_argument('--emax', type=float,
+                        default=1e30, nargs='?',
+                        help="Maximum energy")
+    parser.add_argument('--emin', type=float,
+                        default=0., nargs='?',
+                        help="Minimum energy")
     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)
+                              tmin=args.tmin, tmax=args.tmax,
+                              emin=args.emin, emax=args.emax)
         print(h)
 
 # EOF

analysis/tasks/do_fc_cut_and_project.py

@@ -81,6 +81,7 @@ def task_FC_get_background_count():
                                           'src_file': source_file,
                                           'target_file': target_file,
                                                  'FC_det': FC_det}),
+                   
                 ]
             }
 
@@ -122,7 +123,7 @@ def task_FC_remove_background():
             }
 
 
-# Multiply by !+pileup rate
+# Multiply by 1+pileup rate
 def task_FC_pileup_scaleup():
     '''get the bg count per bin from data before the gamma flash'''
     def pileup_scaleup(info_file: str, 
@@ -133,7 +134,8 @@ def task_FC_pileup_scaleup():
         h_puscu = h * (1.0 + pileup) 
         open(target_file, 'w').write(str(h_puscu))
         return True
-
+    
+    for FC_det in FC_detectors:
         # source file
         source_file = pathlib.Path(f'./output/FC_time_projection/{FC_det}_clean.egated_tproject_bgsub.h1.txt')
         # info on pileup
@@ -142,17 +144,17 @@ def task_FC_pileup_scaleup():
         target_file = f'./output/FC_time_projection/{FC_det}_clean.egated_tproject_bgsub_puscu.h1.txt'
                                           
         yield {
-                'name': target_file,
-                'file_dep': [source_file,
-                             info_file],
-                'targets': [target_file],
-                'actions': [
-                    (create_folder, (str(source_file.parent),)),
-                    (pileup_scaleup, None, {'info_file': info_file,
-                                          'src_file': source_file,
-                                          'target_file': target_file}),
-                ]
-            }
+            'name': target_file,
+            'file_dep': [source_file,
+                         info_file],
+            'targets': [target_file],
+            'actions': [
+                (create_folder, (str(source_file.parent),)),
+                (pileup_scaleup, None, {'info_file': info_file,
+                                        'src_file': source_file,
+                                        'target_file': target_file}),
+                        ]
+        }
         
         
 def task_draw_1D():

analysis/tasks/do_fc_timealign.py

@@ -188,6 +188,7 @@ def task_calibrate_FC_2D():
                     
 def task_merge_calibrated_bidim():
     '''Merge the claibrated H2 into one file'''
+    yml_conf = yaml.load(open('etc/FC_time_align.yaml', 'r'), Loader=yaml.Loader)
     for FC_det in FC_detectors:
         for flavor in ('all', 'clean'):
             # define target params
@@ -197,15 +198,16 @@ def task_merge_calibrated_bidim():
             for week in validated_data_weeks:
                 for h2_file in pathlib.Path().glob(f'./output/calibrated_2D/{week}/{FC_det}_*.h2_{flavor}.txt'):
                     files_list.append(str(h2_file))
+            # end for
             command_ = " ".join(["./bin/fill_into_h2.`hostname`",
-                             "--xnbins=7000 --xmin=-200.0 --xmax=6800.",
+                             f"--xnbins={yml_conf['bidim_tnbins']} --xmin={yml_conf['bidim_tmin']} --xmax={yml_conf['bidim_tmax']}",
                              "--ynbins=10000 --ymin=0.0 --ymax=20000.0",
                              " ".join(files_list),
                              f"> {target_file}"])
             
             yield {
                 'name': f"{FC_det}, {flavor}",
-                'file_dep': files_list,
+                'file_dep': ['etc/FC_time_align.yaml'] + files_list,
                 'targets': [target_file],
                 'actions': [
                     (create_folder, (target_dir,)),

analysis/tasks/do_fc_to_flux.py

+'''
+This file contain all the necessary steps for a raw .evt to histogram conversion
+'''
+
+from doit.tools import create_folder
+import yaml
+import pathlib
+from itertools import product, chain
+from random import uniform as rand_uniform
+
+import scipy
+from scipy.interpolate import interp1d
+
+import json
+
+from pyhisto import Histogram
+
+
+the_mapping = yaml.load(open('etc/mapping.yaml'),
+                        Loader=yaml.Loader)
+detectors = the_mapping['detectors']
+FC_detectors = list(filter(lambda x:x.startswith('U'), detectors))
+Ge_detectors = list(filter(lambda x: not x.startswith('U'), detectors))
+
+validated_data_weeks = yaml.load(open('data/validated_data_list.yaml'),
+                               Loader=yaml.Loader)['data_weeks']
+
+
+
+
+
+def task_FC_time_to_en():
+    '''cut the FC 2D according to ecut and prject the time'''    
+    for FC_det in FC_detectors:
+        source_file = pathlib.Path(f'./output/FC_time_projection/{FC_det}_clean.egated_tproject_bgsub_puscu.h1.txt')
+        target_dir = "./output/flux"
+        target_file = f"{target_dir}/{FC_det}_nevents.inenergy.h1.txt"
+        # log parameters
+        this_dof = yaml.load(open('./etc/fc_dof.yaml'), Loader=yaml.Loader)['FC_DOF'][FC_det]
+        spect_yml = yaml.load(open('./etc/fc_energy_spect.yaml'), Loader=yaml.Loader)
+        emin, emax = spect_yml['min_neutron_energy'], spect_yml['max_neutron_energy']
+        tmin = spect_yml['min_time']
+        tshift = spect_yml['time_shift']
+        #      prepr call to script time_to_energy_1D
+        yield {
+                'name': FC_det,
+                'actions': [
+                            (create_folder, (target_dir,)),
+                            f"./env.run scripts/time_to_energy_1D.py --dof={this_dof} --tmin={tmin} --emin={emin} --emax={emax} --tshift={tshift} {source_file} > {target_file}",
+                            f"./env.run scripts/Draw1Dh.py {target_file}"
+                           ],
+                'targets': [target_file],
+                'file_dep': [
+                             source_file, 
+                             pathlib.Path(f'./etc/fc_energy_spect.yaml'),
+                             pathlib.Path(f'./etc/fc_dof.yaml')
+                            ],
+            }
+
+def task_FC_spect_to_perkeV():
+    '''Convert FC Nevents in energy spectraum into the same, but tih Neve/keV in Y'''
+    def spect_to_perkeV(source_file: str,
+                        target_file: str) -> bool:
+        hsrc = Histogram(fromfile=source_file)
+        htrgt = hsrc.empty_copy()
+        for i, b in enumerate(hsrc):
+            htrgt[i].count = b.count / b.width()
+        open(target_file, 'w').write(str(htrgt))
+        return True
+    for FC_det in FC_detectors:
+        source_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy.h1.txt')
+        target_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy_perkeV.h1.txt')
+        yield {
+            'name': FC_det,
+            'file_dep': [source_file],
+            'targets': [target_file],
+            'actions': [
+                (spect_to_perkeV, None, {'source_file': source_file,
+                                          'target_file': target_file}),
+                 f"./env.run scripts/Draw1Dh.py {target_file}"
+            ],
+        }
+        
+def task_interpolate_perkeV():
+    '''interpolate the spect per keV into a spectra in keV'''
+    def interp_to_keVbins(source_file: str, 
+                          target_file: str,
+                          spect_yaml: str) -> bool:
+        hsrc = Histogram(fromfile=source_file)
+        spect_yml = yaml.load(open(spect_yaml), Loader=yaml.Loader)
+        htrgt = Histogram(nbins=spect_yml['number_of_bins'],
+                          xmin=spect_yml['min_neutron_energy'],
+                          xmax=spect_yml['max_neutron_energy'])
+        # Keeping original number of hits
+        Sorigin = sum(hsrc)
+        # interpolating
+        hinterpolator = interp1d(*hsrc.xy(),
+                                 kind='quadratic',
+                                 fill_value="extrapolate")
+        
+        for i, b in enumerate(htrgt):
+            htrgt.bins[i].count = max(0, hinterpolator(b.center()))
+        # restoring sum
+        htrgt *= Sorigin / sum(htrgt)
+        # exporting the histogram
+        open(target_file, 'w').write(str(htrgt))
+        return True
+        
+    for FC_det in FC_detectors:
+        source_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy_perkeV.h1.txt')
+        target_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy_perkeV_interp.h1.txt')
+        yield {
+            'name': FC_det,
+            'file_dep': [source_file, pathlib.Path('./etc/fc_flux_spect.yaml')],
+            'targets': [target_file],
+            'actions': [
+                (interp_to_keVbins, None, {'source_file': source_file,
+                                           'target_file': target_file,
+                                           'spect_yaml': './etc/fc_flux_spect.yaml'}),
+                 f"./env.run scripts/Draw1Dh.py {target_file}"
+            ],
+        }
+        
+def task_divide_per_sigmanf():
+    '''divide the per keV histogram bins count per the value of the cross section at the venter of the bin'''
+    def divide_per_sigmanf(source_file: str, 
+                           target_file: str,
+                           spect_yaml: str,
+                           nf_db: str) -> bool:
+        hsrc = Histogram(fromfile=source_file)
+        eval_to_use = yaml.load(open(spect_yaml), Loader=yaml.Loader)['nf_eval']
+        htrgt = hsrc.empty_copy()
+        # getting the evaluation
+        eval_db = json.load(open(nf_db,'r'))
+        the_eval_pts = next(filter(lambda x: x.get('fName')==eval_to_use,
+                                   eval_db['funcs']))['pts']
+        evalinterpolator = interp1d(list([pts['x']*1000. for pts in the_eval_pts]),
+                                    list([pts['y']*1.0 for pts in the_eval_pts]),
+                                    kind='linear',
+                                    fill_value="extrapolate")
+        
+        for i, b in enumerate(hsrc):
+            htrgt.bins[i].count = b.count / evalinterpolator(b.center())
+        open(target_file, 'w').write(str(htrgt))
+        return True
+    
+    for FC_det in FC_detectors:
+        source_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy_perkeV_interp.h1.txt')
+        target_file = pathlib.Path(f'./output/flux/{FC_det}_nevents.inenergy_perkeV_interp_divpersigma.h1.txt')
+        yield {
+            'name': FC_det,
+            'file_dep': [source_file, pathlib.Path('./etc/fc_flux_spect.yaml'), pathlib.Path('./etc/nf_crosssections/db.json')],
+            'targets': [target_file],
+            'actions': [
+                (divide_per_sigmanf, None, {'source_file': source_file,
+                                            'target_file': target_file,
+                                            'spect_yaml': './etc/fc_flux_spect.yaml',
+                                            'nf_db': './etc/nf_crosssections/db.json'}),
+                 f"./env.run scripts/Draw1Dh.py --yscale=log --xscale=log {target_file}"
+            ],
+        }
+        
+# LAST TASK : multiply by all the characteristics of the target : mass, target density, ...
\ No newline at end of file