adapting the energy scan retrieval to the new fit scheme

XRStools/XRS_swissknife.py

@@ -1195,7 +1195,7 @@ def loadscan_2Dimages_galaxies(mydata):
 
     Ydim = mydata['Ydim']
     Zdim = mydata['Zdim']
-    Edim = mydata['Zdim']
+    Edim = mydata['Edim']
 
     hf = h5py.File(mydata["signalfile"],"w")
     for iE in range(Edim):
@@ -1323,15 +1323,26 @@ def loadscan_2Dimages_galaxies_foilscan(mydata):
             
         rois_dict_group[roikey] = h5py.SoftLink(  "/"+ roigroup )
 
-    filename, dataname = split_hdf5_address( mydata["expdata"]       )
-
-    data =  np.array(h5py.File(filename,"r")[dataname][:])
 
     
+    if mydata["expdata"]  not in [None, "None"]:
+        filename, dataname = split_hdf5_address( mydata["expdata"]       )
+        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
-        sliced = data[:, origin[0]:origin[0]+box.shape[0], origin[1]:origin[1]+box.shape[1]] * box
+
+        if data is not None:
+            sliced = data[:, origin[0]:origin[0]+box.shape[0], origin[1]:origin[1]+box.shape[1]] * box
+        else:
+            sliced    = np.ones([1, box.shape[0], box.shape[1]],"f")
+            sliced[:] = box
 
         if isolateSpot:
             imageLines = np.sum(sliced,axis=1)
@@ -1887,12 +1898,17 @@ def create_rois_galaxies(mydata):
 
     repertorio = mydata["expdata"]
 
-    filename, dataname = split_hdf5_address(repertorio)
-
-
-    
-    image4roi =  (np.array(h5py.File(filename,"r")[dataname])).sum(axis=0)
-    
+    if "scans" not in mydata:
+        filename, dataname = split_hdf5_address(repertorio)
+        image4roi =  (np.array(h5py.File(filename,"r")[dataname])).sum(axis=0)
+    else:
+        scans = mydata["scans"]
+        image4roi = 0.0
+        for sc in scans:
+            filename, dataname = split_hdf5_address(repertorio%sc)
+            image4roi +=  (np.array(h5py.File(filename,"r")[dataname])).sum(axis=0)
+        
+        
 
     if image4roi.shape == (512,768) :
         layout="2X3-12"

nonregressions/volumes/test2/batch_extraction_esynth1.py

+#  newfactors e' gia pronto , dopo verifica fuzionalita di roi_sel in ximages
+
+
+import numpy as np
+import h5py
+import glob
+import json
+
+
+BATCH_PARALLELISM = 1
+
+import os
+
+def process_input(s, go=0, exploit_slurm_mpi = 0, stop_omp = False):
+    open("input_tmp_%d.par"%go, "w").write(s)
+    background_activator = ""
+    if (go % BATCH_PARALLELISM ):
+        background_activator = "&"
+
+    prefix=""
+    if stop_omp:
+        prefix = prefix +"export OMP_NUM_THREADS=1 ;"
+        
+    if (  exploit_slurm_mpi==0  ):
+        os.system(prefix +"mpirun -n 1 XRS_swissknife  input_tmp_%d.par  %s"%(go, background_activator))
+    elif (  exploit_slurm_mpi>0  ):
+        os.system(prefix + "mpirun XRS_swissknife  input_tmp_%d.par  %s"%(go, background_activator) )
+    else:
+        os.system(prefix + "mpirun -n %d XRS_swissknife  input_tmp_%d.par  %s"%(abs( exploit_slurm_mpi  ), go, background_activator) )
+ 
+
+def select_rois(  data_path_template=None, filter_path=None, roi_target_path=None, scans_to_use=None   ):
+    
+    inputstring = """
+
+    create_rois_galaxies :
+       expdata     :  {data_path_template}
+       filter_path  :  {filter_path}
+       roiaddress :       {roi_target_path}      # the target destination for rois
+       scans  : {scans_to_use}
+
+    """ .format(data_path_template = data_path_template,
+                filter_path = filter_path,
+                roi_target_path =  roi_target_path,
+                scans_to_use = scans_to_use
+    ) 
+    process_input( inputstring , exploit_slurm_mpi = 0 )
+
+def extract_sample_givenrois(
+        roi_path = None,
+        data_path_template = None,
+        monitor_path_template = None,
+        scan_interval = None,
+        Ydim = None,
+        Zdim = None,
+        Edim = None,
+        signal_target_file = None
+         ):
+
+    inputstring = """
+        loadscan_2Dimages_galaxies :
+
+             roiaddress : {roi_path}
+
+             expdata  :  {data_path_template}
+
+             monitor_address : {monitor_path_template}
+
+             scan_interval    :  {scan_interval}
+
+             Ydim : {Ydim}
+             Zdim : {Zdim}
+             Edim : {Edim}
+
+             signalfile : {signal_target_file}
+
+    """.format( roi_path = roi_path,
+                data_path_template = data_path_template,
+                monitor_path_template = monitor_path_template,
+                scan_interval = scan_interval,
+                Ydim = Ydim,
+                Zdim = Zdim,
+                Edim = Edim,
+                signal_target_file = signal_target_file) 
+    
+
+    process_input( inputstring, exploit_slurm_mpi = 0)
+
+
+def  InterpInfo_Esynt( cenom , energy_exp_grid = None,    custom_ene_list =  None):
+
+    info_dict = {}
+
+    N_custom = len(custom_ene_list)
+    N_data   = len( energy_exp_grid ) 
+    
+    for i_interval in range(len(custom_ene_list)):
+        info_dict[str(i_interval)]      =   {}            
+        info_dict[str(i_interval)]["E"] = custom_ene_list[ i_interval   ]
+        info_dict[str(i_interval)]["coefficients"]={}
+
+        for i_n in range(len(energy_exp_grid)):
+
+            info_dict[str(i_interval)]["coefficients"  ][  str(i_n)  ]={}
+            for roi_num, de in enumerate(     cenom   ):
+                info_dict[str(i_interval)]["coefficients"  ][ str(i_n) ][ str(roi_num) ] = 0
+
+
+    for i_interval in range( N_custom -1):
+        
+        cE1 =  custom_ene_list[ i_interval   ]
+        cE2 =  custom_ene_list[ i_interval+1 ]
+        
+        for i_ene, t_ene  in enumerate( energy_exp_grid)   :
+            
+            for roi_num, de in enumerate(     cenom   ):
+                if  t_ene+de < cE1 or t_ene+de > cE2:
+                    continue
+                alpha = (cE2-(t_ene+de) )/(cE2-cE1)
+
+                info_dict[str(i_interval)]["coefficients"  ][  str(i_ene)  ][  str(roi_num)  ]   = alpha
+                info_dict[str(i_interval+1)]["coefficients"][  str(i_ene)  ][  str(roi_num)  ] = 1-alpha
+                
+    print( info_dict)
+    return info_dict    
+
+
+
+    def __init__(self, cenom, interp_file, source,  custom_ene_list = None):
+        
+        volum_list = list(interp_file[source].keys())
+        scan_num_list = np.array([ int( t.split("_") [1]) for t in volum_list])
+        
+        ene_list      = np.array([    interp_file[source][vn]["scans"]["Scan%03d"%sn ]["motorDict"]["energy"][()] for vn,sn in zip(volum_list, scan_num_list   )   ])
+
+        print ( " ecco la scannumlist " , scan_num_list)
+        print (" ecco ene_list", ene_list)
+        
+        
+        self.volum_list    =  volum_list
+        self.scan_num_list =  scan_num_list
+        self.ene_list      =  ene_list
+
+        order = np.argsort(  self.ene_list    )
+        
+        self.ene_list  = self.ene_list [order]
+
+        if custom_ene_list is None:
+            self.custom_ene_list      = self.ene_list
+        else:
+            self.custom_ene_list   = custom_ene_list
+
+        self.scan_num_list  = self.scan_num_list [order]
+        self.volum_list  = [ self.volum_list [ii]  for ii in order  ] 
+        
+        self.interp_file=interp_file
+        self.source= source
+        self.cenom=cenom
+        
+    def interpola_Esynt(self,  roi_sel=roi_sel ):
+        print ( " ECCO I DATI ")
+        print (  self.ene_list  ) 
+        print (  self.cenom   )
+
+        info_dict = {}
+
+
+        for i_intervallo in range(len(self.custom_ene_list)):
+            info_dict[str(i_intervallo)]      =   {}            
+            info_dict[str(i_intervallo)]["E"] = self.custom_ene_list[ i_intervallo   ]
+            info_dict[str(i_intervallo)]["coefficients"]={}
+            for t_vn, t_sn, t_ene in list(zip(self.volum_list,  self.scan_num_list, self.ene_list    )):
+                info_dict[str(i_intervallo)]["coefficients"  ][  t_vn  ]={}
+
+            
+        for i_intervallo in range(len(self.custom_ene_list)-1):
+            cE1 =  self.custom_ene_list[ i_intervallo   ]
+            cE2 =  self.custom_ene_list[ i_intervallo+1 ]
+            for t_vn, t_sn, t_ene in list(zip(self.volum_list,  self.scan_num_list, self.ene_list    ))[0:]:
+                for roi_num, de in enumerate(     self.cenom   ):
+                    if roi_num not in roi_sel:
+                        continue
+                    if  t_ene+de < cE1 or t_ene+de > cE2:
+                        continue
+                    alpha = (cE2-(t_ene+de) )/(cE2-cE1)
+
+                    info_dict[str(i_intervallo)]["coefficients"  ][  str(t_vn)  ][  str(roi_num)  ]   = alpha
+                    info_dict[str(i_intervallo+1)]["coefficients"][  str(t_vn)  ][  str(roi_num)  ] = 1-alpha
+
+        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 ) 
+
+    process_input( inputstring , exploit_slurm_mpi = 1) 
+
+
+
+
+
+def get_scalars( iE = None,
+                 signals_file      = None,
+                 reference_file    = None,
+                 target_file    = None 
+):
+
+    inputstring = """
+    superR_scal_deltaXimages_Esynt :
+         sample_address : {signals_file}:/E{iE}
+         delta_address : {reference_file}:/Scan0
+         load_factors_from :
+         nbin : 1
+         target_address : {target_file}:/scals_{iE}/scal_prods
+    """ % format( iE = iE,
+                  signals_file      = signals_file      ,
+                  reference_file    = reference_file    ,
+                  target_file    = target_file
+    )
+    process_input( inputstring, exploit_slurm_mpi = 0)    
+
+
+    
+def get_volume_Esynt( scalarprods_file = None, interpolation_file = None):
+    
+    os.system("mkdir DATASFORCC")
+
+    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
+    ) 
+    
+    s=inputstring.format(   nick=nick , volumes_file=volumes_file)
+    print ( " INPUT ", s)
+    process_input(s,  exploit_slurm_mpi = 1)
+
+    
+def myOrder(tok):
+    if("volume" not  in tok):
+        tokens = tok.split("_")
+        print( tokens)
+        return int(tokens[1])*10000+ int(tokens[2])
+    else:
+        return 0
+    
+def reshuffle(   volumefile  = "volumes.h5",   nick = None    ):
+
+    h5file_root = h5py.File( volumefile ,"r+" )
+    h5file = h5file_root[nick]
+    scankeys = list( h5file.keys())
+    scankeys.sort(key=myOrder)
+    print( scankeys) 
+    
+    volumes = []
+    for k in scankeys:
+        if k[:1]!="_":
+            continue
+        print( k)
+        if "volume" in h5file[k]:
+            volumes.append( h5file[k]["volume"]  )
+    # volume = np.concatenate(volumes,axis=0)
+    volume = np.array(volumes)
+    if "concatenated_volume" in h5file:
+        del h5file["concatenated_volume"]
+    h5file["concatenated_volume"]=volume
+    h5file_root.close()
+    
+## THE FOLLOWING PART IS THE RELEVANT ONE
+
+cenom = np.array(
+                 [  6471.983002,
+                    6470.612314,
+                    6471.318152,
+                    6471.217828,
+                 ])
+
+Enominal = 6470.0
+
+
+cenom -= Enominal
+
+datadir =  "/data/scisofttmp/mirone/Loic1/data/"
+
+filter_path =  datadir + "../mask.h5:/filter"
+filename_template = "HC-D_2Dmap_%04d"
+data_path_template =  datadir + filename_template + ".nxs:/root.spyc.config1d_RIXS_0024/scan_data/data_03"
+monitor_path_template = datadir + filename_template +"_monitor.hd5:/monitor"
+
+
+
+energy_custom_grid = np.array([6746.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6790.5])
+energy_exp_grid  = h5py.File( datadir +(filename_template%1)+".nxs" ,"r")["/root.spyc.config1d_RIXS_0024/scan_data/actuator_1_1"][()]
+
+
+scan_interval   =  [1,476]   # from 1 to 475 included 
+scans_to_use_for_roi = list(range(1,20))
+
+first_scan_n = 1
+
+
+Ydim = 25
+Zdim = 19
+Edim = 9
+
+
+
+
+roi_target_path  = "results/myrois.h5:/ROIS"
+reference_target_file = "results/response.h5"
+signal_target_file = "results/extracted.h5" 
+scalarprods_target_file  = "results/scalar_prods.h5"
+
+
+if(0):   # 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=scans_to_use_for_roi   )
+roi_path = roi_target_path
+
+
+if(0):  
+    get_reference( roi_path = roi_path ,    reference_target_file = reference_target_file   )
+reference_file = reference_target_file
+
+
+    
+if(1): # SAMPLE extraction
+    extract_sample_givenrois(
+        roi_path              = roi_path               ,
+        data_path_template    = data_path_template     ,
+        monitor_path_template = monitor_path_template  ,
+        
+        scan_interval      = scan_interval      ,
+        Ydim               = Ydim               ,
+        Zdim               = Zdim               ,
+        Edim               = Edim               ,
+        signal_target_file = signal_target_file 
+    )
+signal_file = signal_target_file    
+    
+if(0):    # INTERPOLATION  ESYNTH
+    
+    info_dict = InterpInfo_Esynt(  cenom , energy_exp_grid = energy_exp_grid,    custom_ene_list =  energy_custom_grid      )
+    json.dump(info_dict,open("interpolation.json","w"), indent=4)
+        
+interpolation_file = "interpolation.json"
+if(1):
+    for iE in range(Edim) :
+        get_scalars(  iE = iE,
+                      signals_file = signals_file,
+                      reference_file = reference_file,
+                      target_file =  scalarprods_target_file  
+        )
+        
+scalarprods_file = scalarprods_target_file
+
+
+        
+# ### ESYNTH
+if(0):
+    get_volume_Esynt( scalarprods_file = scalarprods_file,
+                      interpolation_file = interpolation_file):
+
+
+