OK with esynth

XRStools/XRS_swissknife.py

@@ -372,6 +372,9 @@ roiaddress   must be given in the form  roiaddress : "myfile.hdf5:/path/to/hdf5/
 but : was not found
 """)
     filename, groupname = dataadress[:pos], dataadress[pos+1:]
+    if( len(groupname) and groupname[0:1] !="/" ):
+        groupname = "/"+groupname
+    print(" RITORNO ", groupname, (groupname[0:1] !="/" ), len(groupname) )
     return filename, groupname 
 
 
@@ -1174,11 +1177,11 @@ def loadscan_2Dimages_galaxies_foilscan(mydata):
          signalfile : "filename"  # Target file for the signals
 
          isolateSpot : 0   # if different from zero selects on each image  ( when sumto1d=0 ) ROI  the spot region and sets to zero outside a radius =  isolateSpot
- 
+         scan_list : None # if given the expdata input  will be used as a template to load the scans in the scan_list
     #
     """
     
-    allowed_keys = ["roiaddress",'isolateSpot',"signalfile","expdata"   ]
+    allowed_keys = ["roiaddress",'isolateSpot',"signalfile","expdata", "scan_list"   ]
     check_allowed_keys(mydata, allowed_keys)
     
     roiaddress=None
@@ -1221,17 +1224,21 @@ def loadscan_2Dimages_galaxies_foilscan(mydata):
             
         rois_dict_group[roikey] = h5py.SoftLink(  "/"+ roigroup )
 
-
-    
     if mydata["expdata"]  not in [None, "None"]:
-        filename, dataname = split_hdf5_address( mydata["expdata"]       )
-        data =  np.array(h5py.File(filename,"r")[dataname][:])
+        if "scan_list" not in mydata or  mydata["scan_list"] in [None, "None"]:
+            filename, dataname = split_hdf5_address( mydata["expdata"]       )
+            data =  np.array(h5py.File(filename,"r")[dataname][:])
+        else:
+            data = 0
+            for scan_number in mydata["scan_list"]:
+                filename, dataname = split_hdf5_address( mydata["expdata"]%scan_number       )
+                data =  data +  np.array(h5py.File(filename,"r")[dataname][:])
+            
     else:
         # per scan puramente energetici
         data = None
         isolateSpot = 0
 
-        
     for roikey, (origin, box) in roiob.red_rois.items():
 
         roigroup=roikey[3:]+"/"  # with "ROI" removed, nly numerical part

XRStools/reponse_percussionelle.py

@@ -727,6 +727,10 @@ def get_spots_list(  filename  , groupname , filter_rois =1 ):
     xscales = {}
     enescan = None
 
+    print( " ===========================\n"*10  )
+    print( filename   )
+    print( groupname  ) 
+    
     tmp_list = list(filterRoiList(list(h5[groupname].keys())))
     tmp_list = sorted(tmp_list, key = int )
     

nonregressions/volumes/test1/01_create_rois_galaxies.yaml

 create_rois_galaxies :
-   expdata     :  "/data/scisofttmp/mirone/raffaela/kapton_00001_01.nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"  ##-- this pointsto, inside a hdf5 file, a stack of images 
-   filter_path  :  "/data/scisofttmp/mirone/raffaela/mymask.h5:/mymask"    
+   expdata     :  "/data/raffaela/kapton_00001_01.nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"  ##-- this pointsto, inside a hdf5 file, a stack of images 
+   filter_path  :  "/data/raffaela/mymask.h5:/mymask"    
    roiaddress : "results/myrois.h5:/ROIS"            # the target destination for rois
\ No newline at end of file

nonregressions/volumes/test1/02_loadscan_2Dimages_galaxies_foilscan.yaml

 loadscan_2Dimages_galaxies_foilscan :
-     roiaddress : "/data/scisofttmp/mirone/raffaela/myrois.h5:/ROIS"   # the same given in create_rois
-     expdata    :  "/data/scisofttmp/mirone/raffaela/kapton_00001_01.nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"  
+     roiaddress : "/data/raffaela/myrois.h5:/ROIS"   # the same given in create_rois
+     expdata    :  "/data/raffaela/kapton_00001_01.nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"  
      signalfile : "results/response.h5"  # Target file for the signals
      isolateSpot : 7

nonregressions/volumes/test3/batch_extraction_esynth1.py

@@ -24,6 +24,7 @@ def main():
     filter_path =  dataprefix + "mymask.h5:/mymask"
     filename_template = "GS9_%05d_01"
     data_path_template =  datadir + filename_template + ".nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"
+    reference_data_path_template =  dataprefix + "kapton_%05d_01" + ".nxs:/root_spyc_config1d_RIXS_00001/scan_data/data_07"
 
     
     monitor_path_template = None
@@ -43,14 +44,14 @@ def main():
     Zdim = 16
     Edim = 19
 
-    reference_scan_list =  list(range(1,20))
+    reference_scan_list =  [1]
     
     custom_components_file = None
     # custom_components_file = "components.h5"
-    roi_scan_num = list(range(1,20))
+    roi_scan_num = list(range(1,3))
 
     reference_clip = None
-    # reference_clip = [ 90, 180 ]
+    reference_clip = [ 0, 200 ]
 
     ## in the reference scan for each position there is a spot with a maximum. We set zero the background which is further than
     ## such radius from the maximum
@@ -66,16 +67,16 @@ def main():
     selected_rois = list(range(0,4) )
     
     steps_to_do = {
-        "do_step_make_roi":                      True,
-        "do_step_sample_extraction":             True,
+        "do_step_make_roi":                      False,
+        "do_step_sample_extraction":             False,
         
         "do_step_extract_reference_scan":        False,
-        "do_step_fit_reference_response":               False,
+        "do_step_fit_reference_response":        False,
         "do_step_resynthetise_reference":        False,
         
         "do_step_scalar":               False, 
         "do_step_interpolation_coefficients":    False,
-        "do_step_finalise_for_fit":              False
+        "do_step_finalise_for_fit":              True
     }
 
     scalar_products_target_file  = "results/scalar_prods.h5"
@@ -99,7 +100,8 @@ def main():
                       roi_scan_num          = roi_scan_num          ,
                       roi_target_path       = roi_target_path       ,
                       data_path_template    = data_path_template    ,
-
+                      reference_data_path_template = reference_data_path_template,
+                      
                       steps_to_do = steps_to_do,
                       
                       scan_interval         = scan_interval         ,                      
@@ -128,11 +130,127 @@ def main():
                     
                       energy_custom_grid    = energy_custom_grid    ,
                       custom_components_file = custom_components_file,
-                      interpolation_infos_file =  interpolation_infos_file
+                      interpolation_infos_file =  interpolation_infos_file,
+                      energy_exp_grid = energy_exp_grid
     ) 
 
+def synthetise_response(scan_address=None, target_address=None,             response_fit_options = None
+):
+    input_string = """
+    superR_fit_responses :
+       foil_scan_address : "{scan_address}"
+       nref : 7                 # the number of subdivision per pixel dimension used to 
+                                # represent the optical response function at higher resolution
+       niter_optical  :  {niter_optical}    # the number of iterations used in the optimisation of the optical
+                                # response
+       beta_optical  :  {beta_optical}     # The L1 norm factor in the regularisation 
+                                #  term for the optical functions
+       pixel_dim : 1            # The pixel response function is represented with a 
+                                #  pixel_dim**2 array
+       niter_pixel : 10        # The number of iterations in the pixel response optimisation
+                                # phase. A negative number stands for ISTA, positive for FISTA
+       beta_pixel  :  0.0    # L1 factor for the pixel response regularisation
+
+       ## The used trajectories are always written whith the calculated response 
+       ## They can be reloaded and used as initialization(and freezed with do_refine_trajectory : 0 )
+       ## Uncomment the following line if you want to reload a set of trajectories
+       ## without this options trajectories are initialised from the spots drifts
+       ##
+       #   reload_trajectories_file : "response.h5"
+
+       filter_rois : 0
+
+
+       ######
+       ## The method first find an estimation of the foil scan trajectory on each roi
+       ## then, based on this, obtain a fit of the optical response function
+       ## assuming a flat pixel response. Finally the pixel response is optimised
+       ##
+       ## There is a final phase where a global optimisation
+       ## is done in niter_global steps.
+       ##
+       ## Each step is composed of optical response fit, followed by a pixel response fit.
+       ## If do_refine_trajectory is different from zero, the trajectory is reoptimised at each step
+       ## 
+       niter_global  :  {niter_global}
+
+       ## if do_refine_trajectory=1 the start and end point of the trajectory are free
+       ##  if =2 then the start and end point are forced to a trajectory which is obtained
+       ##  from a reference scan : the foil scan may be short, then one can use the scan of
+       ##   an object to get another one : key *trajectory_reference_scan_address*
+       ##
+
+       do_refine_trajectory : 1
+
+       ## optional: only if do_refine_trajectory = 2
+
+       trajectory_reference_scansequence_address : "demo_newrois.h5:/ROI_FOIL/images/scans/"
+       trajectory_threshold   : 0.1
+
+       ## if the pixel response function is forced to be symmetrical 
+
+       simmetrizza : 1
+
+       ## where the found responses are written
+
+       target_file : {target_address}
+       # target_file : "fitted_responses.h5"
+
+    """ 
+    s=input_string.format(  scan_address=scan_address ,
+                           target_address=target_address,
+                           niter_optical = response_fit_options[ "niter_optical"],
+                           beta_optical=response_fit_options["beta_optical"],
+                           niter_global=response_fit_options["niter_global"]
+    )
+    process_input( s , exploit_slurm_mpi = 1, stop_omp = True) 
+
+
+    
+def    resynthetise_scan(
+        old_scan_address= None,
+        response_file  =  None ,
+        target_address =  None,
+        original_roi_path = None,
+        resynth_z_square = None
+    ):
+
+
+    input_string = """
+     superR_recreate_rois :
+     ### we have calculated the responses in responsefilename
+         ### and we want to enlarge the scan  by a margin of 3 times
+         ### the original scan on the right and on the left 
+         ###  ( so for a total of a 7 expansion factor )
+
+         responsefilename :  {response_file}
+         nex : 0
+
+         ## the old scan covered by the old rois
+         old_scan_address : {old_scan_address}
+
+         ## where new rois and bnew scan are written
+         target_filename : {target_address}
+         filter_rois      : 0
+
+         original_roi_path : {original_roi_path}
+
+         resynth_z_square : {resynth_z_square}
+
+"""
+    s=input_string.format(  response_file = response_file ,
+                            target_address = target_address,
+                            old_scan_address=old_scan_address,
+                            original_roi_path = original_roi_path +"/rois_definition",
+                            resynth_z_square = resynth_z_square)
+    
+    process_input( s , exploit_slurm_mpi = 0, stop_omp = True) 
+
+
 
 
+
+    
 def process_input(s, go=0, exploit_slurm_mpi = 0, stop_omp = False):
     open("input_tmp_%d.par"%go, "w").write(s)
     background_activator = ""
@@ -254,6 +372,8 @@ def InterpInfo_Esynt_components(peaks_shifts , energy_exp_grid = None,    custom
     
 def  InterpInfo_Esynt( peaks_shifts , energy_exp_grid = None,    custom_ene_list =  None):
 
+    print(energy_exp_grid)
+    print(peaks_shifts)
     info_dict = {"energy_exp_grid":list(energy_exp_grid), "de_list": list(peaks_shifts)}
 
     N_custom = len(custom_ene_list)
@@ -365,18 +485,35 @@ def  InterpInfo_Esynt( peaks_shifts , energy_exp_grid = None,    custom_ene_list
 
         return info_dict    
 
-    
-def get_reference( roi_path = None, reference_target_file = None ):
-    inputstring = """
 
-    loadscan_2Dimages_galaxies_foilscan :
-        roiaddress : {roi_path}           
-        expdata    :  None  
-        signalfile : {reference_target_file}  # Target file for the signals
-    """ .format( roi_path = roi_path,  reference_target_file = reference_target_file ) 
+def get_reference(   roi_path =  None,
+                     data_path_template = None,
+                     monitor_path_template = None ,
+                     reference_scan_list = None,
+                     extracted_reference_target_file = None,
+                     isolate_spot_by = None
+):
+    signal_path = extracted_reference_target_file + ":/"
 
-    process_input( inputstring , exploit_slurm_mpi = 0) 
+    input_string = """
+    loadscan_2Dimages_galaxies_foilscan :
+       roiaddress :  {roi_path} 
+       expdata    :  {data_path_template} 
+       signalfile :  "{extracted_reference_target_file}"
+       isolateSpot :  {isolate_spot_by}
+       scan_list : {reference_scan_list}  
+    """
+    s=input_string.format(
+        data_path_template = data_path_template,
+        reference_scan_list = reference_scan_list,
+        
 
+        roi_path = roi_path,
+        isolate_spot_by = isolate_spot_by,
+        signal_path = signal_path,
+        extracted_reference_target_file = extracted_reference_target_file
+    )
+    process_input( s , exploit_slurm_mpi = 0) 
 
 
 def get_scalars( iE = None,
@@ -388,7 +525,7 @@ def get_scalars( iE = None,
     inputstring = """
     superR_scal_deltaXimages_Esynt :
          sample_address : {signals_file}:/E{iE}
-         delta_address : {reference_file}:/Scan0
+         delta_address : {reference_file}:/rois_and_reference/Scan0
          load_factors_from : 
          nbin : 1
          target_address : {target_file}:/{iE}/scal_prods
@@ -408,9 +545,7 @@ def get_volume_Esynt( scalarprods_file = None, interpolation_file = None):
     inputstring = """    
      superR_getVolume_Esynt :
         scalprods_address : {scalarprods_file}:/
-        target_address : {scalarprods_file}:/data_for_volumes
         dict_interp : {interpolation_file}
-        debin : [1, 1]
         output_prefix : DATASFORCC/test0_
     """.format( scalarprods_file = scalarprods_file ,
                 interpolation_file = interpolation_file
@@ -456,7 +591,8 @@ def tools_sequencer(  peaks_shifts = None,
                       roi_scan_num = None,
                       roi_target_path          = None,
                       data_path_template = None,
-
+                      reference_data_path_template = None, 
+                      
                       steps_to_do = None,
                       
                       scan_interval = None,
@@ -467,7 +603,6 @@ def tools_sequencer(  peaks_shifts = None,
                       monitor_path_template = None,                      
                       signals_target_file      = None,
 
-
                       reference_scan_list = None,
                       reference_clip = None,
                       extracted_reference_target_file    = None  ,
@@ -482,28 +617,24 @@ def tools_sequencer(  peaks_shifts = None,
                       
                       scalar_products_target_file  = None,
 
-                    
                       energy_custom_grid    = None    ,
                       custom_components_file = None,
                       interpolation_infos_file =  None,
 
-                      energy_exp_grid  = None,
+                      energy_exp_grid  = None
 
-                      reference_target_file    = None,
-                      scalarprods_target_file  = None,
-                      custom_components = None
 ) :
 
 
     
     if(steps_to_do["do_step_make_roi"]):   # ROI selection and reference scan
-        select_rois(  data_path_template = data_path_template, filter_path = filter_path, roi_target_path = roi_target_path, scans_to_use = roi_scan_num   )
+        select_rois(  data_path_template = reference_data_path_template, filter_path = filter_path, roi_target_path = roi_target_path, scans_to_use = roi_scan_num   )
     roi_path = roi_target_path
 
 
     if("do_step_make_reference" in steps_to_do and steps_to_do["do_step_make_reference"]):  
-        get_reference( roi_path = roi_path ,    reference_target_file = reference_target_file   )
-    reference_file = reference_target_file
+        get_reference( roi_path = roi_path ,    reference_target_file = resynthetised_reference_and_roi_target_file    )
+    reference_file = resynthetised_reference_and_roi_target_file
 
 
 
@@ -522,23 +653,66 @@ def tools_sequencer(  peaks_shifts = None,
     signals_file = signals_target_file    
 
 
+    if(steps_to_do["do_step_extract_reference_scan"]):  # of course we need the REFERENCE SCAN
+        get_reference(
+                       roi_path = roi_path,
+                       data_path_template = reference_data_path_template,
+                       monitor_path_template = monitor_path_template  ,
+                       extracted_reference_target_file = extracted_reference_target_file,
+                       isolate_spot_by = isolate_spot_by,
+                       reference_scan_list = reference_scan_list
+        )        
+        if reference_clip is not None:
+            clip1, clip2= reference_clip
+            print(extracted_reference_target_file ) 
+            ftarget = h5py.File( extracted_reference_target_file  ,"r+")
+            target_group = ftarget["Scan0" ]
+            for k in target_group.keys():
+                if k != "motorDict":
+                    for dsn in ["matrix"]:
+                        mat = target_group[k][dsn][()]
+                        del target_group[k][dsn]
+                        target_group[k][dsn] = mat[clip1:clip2]
+            ftarget.close()
+
+    if(steps_to_do["do_step_fit_reference_response"]):  
+        synthetise_response(
+            scan_address=  extracted_reference_target_file +":Scan0" ,
+            target_address =  +":/FIT",
+            response_fit_options = response_fit_options
+        )
+
+
+        
+    if(steps_to_do["do_step_resynthetise_reference"]):  
+        resynthetise_scan(
+            old_scan_address=  extracted_reference_target_file +":/Scan0" ,
+            response_file  = response_target_file +":/FIT",
+            target_address =  resynthetised_reference_and_roi_target_file +  ":/rois_and_reference",
+            original_roi_path = roi_path,
+            resynth_z_square = resynth_z_square
+        )
+    
 
 
     if(steps_to_do["do_step_scalar"]):
-        os.system("rm %s"%scalarprods_target_file)
+        os.system("rm %s"%scalar_products_target_file)
+
+
+
         for iE in range(Edim) :
             get_scalars(  iE = iE,
                           signals_file = signals_file,
-                          reference_file = reference_file,
-                          target_file =  scalarprods_target_file
+                          reference_file = resynthetised_reference_and_roi_target_file,
+                          target_file =  scalar_products_target_file
             )
 
-    scalarprods_file = scalarprods_target_file
+    scalarprods_file = scalar_products_target_file
 
 
     interpolation_infos_file = "interpolation_infos.json"
     if(steps_to_do["do_step_interpolation_coefficients"]):    # INTERPOLATION  ESYNTH
-        if use_custom_components is None:        
+        if custom_components_file is None:        
             info_dict = InterpInfo_Esynt(  peaks_shifts ,
                                            energy_exp_grid = energy_exp_grid,
                                            custom_ene_list =  energy_custom_grid
@@ -547,7 +721,7 @@ def tools_sequencer(  peaks_shifts = None,
             info_dict = InterpInfo_Esynt_components(  peaks_shifts,
                                                       energy_exp_grid = energy_exp_grid,
                                                       custom_ene_list =  energy_custom_grid,
-                                                      custom_components = use_custom_components
+                                                      custom_components = custom_components_file
             )
         json.dump(info_dict,open( interpolation_infos_file,"w"), indent=4)