Modifications by T. Saerbeck to polarized treatment, adding a calculation and...

Modifications by T. Saerbeck to polarized treatment, adding a calculation and output of spin asymmetry

lamp_mac/FIGARO/cosmos_analysis.pro

+;TS 2019-01-11 Added a print output for upper and lower end of the wavelength of each time channel
 ;TS 2018-06-04 Removed the restriction that om_fwhm is disregarded for polarized runs since the resolution calculation now should be based only on the UP polarization. 
 ;PG 2017-08-07 Added wavelength resolution smearing of wide beam at mid-chopper position (incoherent case), updated gravity correction, mainly openoffset correction, added chopper radius into c_params.cr 
 ;TS 2016-10-25 Deleted om_fwhm^2 since this messes up in comparison between different runs when intensities change and 
@@ -2272,7 +2273,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
            cosmos_logmessage, ' Unable to fit direct peak. This may be due to poor statistics.'
            if total(foregddir) le 0 then begin
               cosmos_logmessage, ' Defaulting to (foreground) peak width of 13 (Auto)'
-              peakdir = cosmos_anal_peaksearchrefined(total(cntdiri, 2), xrangedir, [6,6])
+              peakdir = cosmos_anal_peaksearchrefined(total(cntdiri, 2), xrangedir, [6,6]) ;TS 2018-10-01
               peakdirrange = indgen(13) + peakdir - 6
            endif else peakdir = cosmos_anal_peaksearchrefined(total(cntdiri, 2), xrangedir, foregddir)
         endif
@@ -2293,7 +2294,7 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
      if (peakref lt 0L) then begin
         cosmos_logmessage, ' Unable to fit reflect peak. This may be due to poor statistics.'
         if total(foregdref) le 0 then begin
-           peakref = cosmos_anal_peaksearchrefined(total(cntrefi, 2), xrangeref, [6,6])
+           peakref = cosmos_anal_peaksearchrefined(total(cntrefi, 2), xrangeref, [6,6]) ;TS 2018-10-01
            peakrefrange = indgen(13) + peakref - 6
         endif else peakref = cosmos_anal_peaksearchrefined(total(cntrefi, 2), xrangeref, foregdref)
      endif
@@ -2636,7 +2637,13 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
            if (dir_ref_diff / ref_total_tofd) gt 0.01 then cosmos_logmessage, ' Run no. ' + cosmos_str_make(runno) + ': Different TOF distances from direct and reflect runs.', /WARNING
         endif
         lambda = 1e10 * (c_params.planckperkg * ((findgen(tsize) + 0.5) * parref.channelwidth + parref.delay) / ref_total_tofd)
-    
+        ;TS 2019-01-11 Lower and upper bound of wavelength for each time channel.
+        ;lambdaLowerbound = 1e10 * (c_params.planckperkg * ((findgen(tsize)) * parref.channelwidth + parref.delay) / ref_total_tofd)
+        ;lambdaUpperbound = 1e10 * (c_params.planckperkg * ((findgen(tsize)+1) * parref.channelwidth + parref.delay) / ref_total_tofd)
+        ;print,'array', lambda
+        ;print,'lambdaLowerbound', lambdaLowerbound
+        ;print,'lambdaUpperbound', lambdaUpperbound
+   
         cosmos_logmessage, ' '
         cosmos_logmessage, ' Constants used in calculating lambda :'
         cosmos_logmessage, '  poff         = ' + cosmos_str_make(c_params.poff)
@@ -2955,10 +2962,15 @@ pro cosmos_anal_run, result, groupedresult, direct, reflect, instr, theta, runno
   	    tempa=tempratio*abs(parref.s2w-parref.s3w)+parref.s2w
   	    tempb=tempratio*(parref.s2w+parref.s3w)+parref.s2w
   	    tempwidthfwhm=0.49*(tempb^3-tempa^3)/(tempb^2-tempa^2)
+        ;print, 'new' , tempratio , tempa , tempb 
 
-  	    width_res = tempwidthfwhm/1000.0 * parref.period / (2 * !pi * c_params.cr)* c_params.planckperkg / $;calculate beam width lambda resolution
-           ( lambda[lambda_range] * 1.e-10) / parref.tofd
-         err_res = sqrt(err_res^2+width_res^2)
+        ;calculate beam width lambda resolution
+  	    width_res = tempwidthfwhm/1000.0 * parref.period / $
+  	                (2 * !pi * c_params.cr)* c_params.planckperkg / $
+                    ( lambda[lambda_range] * 1.e-10) / parref.tofd
+        err_res = sqrt(err_res^2+width_res^2)
+        ;print, width_res
+  	    ;err_res = sqrt(err_res^2)
   
   	    ;s3_fwhm=(0.68*parref.s3w)/c_params.interslit  ; fwhm angles of second slit  PG on 16/02/15: corrected for interslit distance
 		    s3_fwhm=(0.68*parref.s3w)/(parref.slit3*1000.0+parref.sdetd*1000.0)  ; fwhm angles of second slit  PG on 28/04/15: corrected for interslit distance
@@ -5507,7 +5519,7 @@ pro cosmos_run_background
   
        endfor ;end for islice = 0,slices-1
      
-     endif else begin ;merging and writing of polarized data
+     endif else begin ;polopt>0 merging and writing of polarized data
 
        ;Check if result1 is a kinetic run. If so check if result2 and result3 contain the same number of time slices.
        res1size = size(result1)
@@ -5532,16 +5544,16 @@ pro cosmos_run_background
          kinetic = 0
          slices = 1
        endelse
-
-       polstring = ['00', '01', '10', '11']
-
+       ;TS 2019-07 Added a field in polstring for calculating and outputting the Spin Asymmetry
+       polstring = ['00', '01', '10', '11', 'SA']
        for ipol = 0, 3 do begin
-
+          
          ;TS 2016-10-21 COSMOS was writing too many files if only two channels were measured. For THREE it is OK since I am calculating the remaining in the correction procedure.
          if polopt le 2 then begin
            if ipol ne 0 then ipol=3
            if ipol ne 0 then cosmos_logmessage, 'Skipped writing of Pol[01] and POL[10]', /WARNING
          endif
+         
          ;TS 2016-10-21 END: COSMOS was writing too many files if only two channels were measured
 
          for islice = 0,slices-1 do begin
@@ -5799,6 +5811,7 @@ pro cosmos_run_background
              if c_settings.outputfmt[1] eq 1 then cosmos_file_write, outfile, result, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /dat
              ;MFT
              if c_settings.outputfmt[2] eq 1 then cosmos_file_write, outfile, result, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /mft
+             
              ;Lambda format --> one different file per angle
              if c_settings.outputfmt[4] eq 1 then begin
                if c_settings.group eq 1 then begin
@@ -5867,6 +5880,17 @@ pro cosmos_run_background
                  if ipol eq 3 then $
                    cosmos_lamp_write, outlamp[0]+3, resultlog, file_basename(outfile)
                endelse
+               if polopt le 2 then begin
+                 if ipol eq 0 then $
+                   ref00 = result
+                 if ipol eq 3 then $
+                   ref11 = result
+               endif else begin
+                 if ipol eq 0 then $
+                   ref00 = result
+                 if ipol eq 3 then $
+                   ref11 = result
+               endelse
                ;TS 2016-10-21 END:: Changed where the pol files get written
              endif else begin
                cosmos_logmessage, 'Log dataset for Lamp empty. --> POL:: ' + cosmos_str_make(ipol) , /WARNING ;TS 2016-10-21 Identify the pol channels that are empty since the message pops up if only two or three channels measured
@@ -5875,13 +5899,48 @@ pro cosmos_run_background
            endif ;end if n_elements(result) gt 1 --> output files
 
          endfor ;end for islice = 0,slices-1
+         
        endfor ;end for ipol = 0,3
-
-     endelse
+       
+       ;TS 2019-07 Added output for Spin Asymmetry for Polarized Data
+       print, 'calculate and write the Spin Asymmetry (SA)'
+       Flipratio = ref00 ;fltarr(5, n_elements(result[0,*]) ) ;[Q, dQ, I, dI, lambda]
+       Flipratio[2,*] = (ref00[2,*]-ref11[2,*]) / (ref00[2,*]+ref11[2,*])
+       print, 'Calculate SA Error'
+       Flipratio[3,*] = 2*sqrt((ref11[2,*]*ref11[2,*]*ref00[3,*]*ref00[3,*]+ref00[2,*]*ref00[2,*]*ref11[3,*]*ref11[3,*])/(ref00[2,*]+ref11[2,*])^4);
+       if polopt ge 1 then begin
+          ;FRlog = cosmos_anal_log(Flipratio)
+          if polopt le 2 then begin
+             cosmos_lamp_write, outlamp[0]+2, Flipratio, file_basename(outfile)
+             outfile = file_dirname(outfile0)+path_sep()+'pol'+polstring[4]+'_'+file_basename(outfile0)
+             ;print, outfile
+             ;Write bundled data to file(s) and logged bundled data to lamp
+             ;AFT
+             if c_settings.outputfmt[0] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /aft
+             ;DAT
+             if c_settings.outputfmt[1] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /dat
+             ;MFT
+             if c_settings.outputfmt[2] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /mft
+          endif else begin
+             cosmos_lamp_write, outlamp[0]+4, Flipratio, file_basename(outfile)
+             outfile = file_dirname(outfile0)+path_sep()+'pol'+polstring[4]+'_'+file_basename(outfile0)
+             ;print, outfile
+             ;Write bundled data to file(s) and logged bundled data to lamp
+             ;AFT
+             if c_settings.outputfmt[0] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /aft
+             ;DAT
+             if c_settings.outputfmt[1] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /dat
+             ;MFT
+             if c_settings.outputfmt[2] eq 1 then cosmos_file_write, outfile, Flipratio, direct1, direct2, direct3, rrr1, rrr2, rrr3, instr1, instr2, instr3, theta1, theta2, theta3, water_file, /mft
+          endelse    
+       endif 
+       ;TS 2019-07 END Added output for Spin Asymmetry for Polarized Data
+     endelse ;end polopt>0 section
      
      t_elapsed = systime(1) - t_ini
      cosmos_logmessage, 'Row treated in ' + cosmos_str_make(t_elapsed) + ' seconds'
-     ;print, c_params.row, t_elapsed, systime(1) - 1427188740.
+     ;TS 2019-07 Added a reminder in which order the polarizations are output
+     if polopt ge 1 then cosmos_logmessage, 'Polarized data output:  [00, 01, 10, 11, FR]'
     
   endif ;end if (c_params.row lt c_params.rowfinish)