Patch to avoid sgroup returning bad results when the number of requested groups is too large.

lamp_mac/TOF/sgroup.pro

@@ -140,36 +140,44 @@ COMMON c_lamp_access, inst
 COMMON printing, iprint, outstring
 
 take_datp, datp
-par =datp.p
-x_in=datp.x
-y_in=datp.y
-e_in=datp.e
+par  = datp.p
+x_in = datp.x
+y_in = datp.y
+e_in = datp.e
 
-iprint=0
+iprint = 0
 
 ; -------------------------------------------
 ; Donnees experimentales
 
 ;Print,'SGROUP ## S.ROLS 2002 ##'
 
-nchannels=n_elements(x_in) & nspectra=N_ELEMENTS(y_in)
+nchannels = n_elements(x_in) 
+nspectra  = n_elements(y_in)
 
 if (nspectra eq 1) then begin
-    str_mess='Data already grouped... we assume ng=1 and no /sint'
-    y_buf=y_in
-    w_buf=w_in
-    e_buf=e_in
+    str_mess = 'Data already grouped... we assume ng=1 and no /sint'
+    y_buf = y_in
+    w_buf = w_in
+    e_buf = e_in
     goto, nogrp
 endif
 
+if ngroups gt nspectra/2 then begin
+    ngroups = nspectra/2
+    print, "Current algorithm limits the number of groups to be <= number of spectra/2."
+    print, "Grouping to ngroups = ", ngroups
+endif
+
 ; -------------------------------------------
 ; tests sur la dimensions des erreurs
 
-sein=n_elements(e_in) & swin=n_elements(w_in)
+sein = n_elements(e_in) 
+swin = n_elements(w_in)
 
 if (sein ne swin) then begin
-    e_in=SQRT(w_in)
-    print,"SGROUP:warning - errors were not calculated ... put them to sqrt(w)"
+    e_in = SQRT(w_in)
+    print, "SGROUP: warning - errors were not calculated ... put them to sqrt(w)"
 endif
 
 
@@ -178,19 +186,19 @@ IF (NOT KEYWORD_SET(scale)) THEN scale=1.
 ; ------------------------------------------
 ; Redimensionnement si THETA_MIN et THETA_MAX indiqu_s
 
-IF n_elements(theta_min) eq 1 THEN tmin=theta_min ELSE tmin=y_in[0]
-IF n_elements(theta_max) eq 1 THEN tmax=theta_max ELSE tmax=y_in[nspectra-1]
+IF n_elements(theta_min) eq 1 THEN tmin = theta_min ELSE tmin = y_in[0]
+IF n_elements(theta_max) eq 1 THEN tmax = theta_max ELSE tmax = y_in[nspectra-1]
 
-ind_in=Where((y_in le tmax) and (y_in ge tmin),ncount)
-y_ini =y_in[ind_in]   & nspectra=ncount
-w_ini =w_in[*,ind_in] & e_ini=e_in[*,ind_in]
-tmin  =y_ini[0]       & tmax=y_ini[ncount-1]
+ind_in = where((y_in le tmax) and (y_in ge tmin), ncount)
+y_ini = y_in[ind_in]   & nspectra = ncount
+w_ini = w_in[*,ind_in] & e_ini = e_in[*,ind_in]
+tmin  = y_ini[0]       & tmax = y_ini[ncount-1]
 
-if KEYWORD_SET(sint) THEN sintt=SIN(y_ini*!pi/180.) ELSE sintt=y_ini*0.+1.
+if KEYWORD_SET(sint) THEN sintt = SIN(y_ini*!pi/180.) ELSE sintt = y_ini*0.+1.
 
-for iy=0,ncount-1 do begin
-    w_ini[*,iy]=w_ini[*,iy]*sintt[iy]
-    e_ini[*,iy]=e_ini[*,iy]*sintt[iy]
+for iy = 0, ncount-1 do begin
+    w_ini[*,iy] = w_ini[*,iy]*sintt[iy]
+    e_ini[*,iy] = e_ini[*,iy]*sintt[iy]
 endfor
 
 ; -------------------------------------------
@@ -200,66 +208,68 @@ IF n_elements(ngroups) eq 1 THEN BEGIN
 
    itest=FIX(ngroups) - ngroups
 
-   if(ngroups le 0) then begin
-       PRINT,'SGROUP: Error - Number of groups MUST be > 0'
-       w_out=w_buf
+   if (ngroups le 0) then begin
+       PRINT, 'SGROUP: Error - Number of groups MUST be > 0'
+       w_out = w_buf
    endif
-   if(itest ne 0) then begin
-       PRINT,'SGROUP: Error - Number of groups MUST be integer'
-       w_out=w_buf
+   if (itest ne 0) then begin
+       PRINT, 'SGROUP: Error - Number of groups MUST be integer'
+       w_out = w_buf
    endif
    if(nspectra lt ngroups) then begin
-       PRINT,'SGROUP: Error - Number of groups MUST be <= initial number of spectra'
-       w_out=w_buf
+       PRINT, 'SGROUP: Error - Number of groups MUST be <= initial number of spectra'
+       w_out = w_buf
    endif
 ENDIF ELSE ngroups=1
 
-ratio=FIX(nspectra/ngroups) & rest=nspectra-ratio*ngroups
+ratio = FIX(nspectra/ngroups) 
+rest = nspectra - ratio*ngroups
 
-if iprint ne 0 then print,inst,nspectra,ngroups,ratio,rest
+if iprint ne 0 then print, inst, nspectra, ngroups, ratio, rest
 
 imin=0 & imax=0
 
-w_buf=FLTARR(nchannels,ngroups)
-e_buf=FLTARR(nchannels,ngroups)
-y_buf=FLTARR(ngroups)
+w_buf = FLTARR(nchannels,ngroups)
+e_buf = FLTARR(nchannels,ngroups)
+y_buf = FLTARR(ngroups)
 
-for i=0,ngroups-1 do begin
-  ajou=0
-  if rest gt 0 then ajou=ratio else ajou=ratio-1
-  if imin ne imax then imin=imax+1
-  imax=imin+ajou
+for i = 0, ngroups-1 do begin
+  ajou = 0
+  if rest gt 0 then ajou = ratio else ajou = ratio-1
+  if imin ne imax then imin = imax+1
+  imax = imin+ajou
   if imax gt imin then begin
-    w_buf[*,i]=total(w_ini[*,imin:imax],2)/(ajou+1)
-    e_buf[*,i]=SQRT(total((e_ini[*,imin:imax])^2,2))/(ajou+1)
-    y_buf[i]=total(y_ini[imin:imax])/(ajou+1)
-    rest=rest-1
+    w_buf[*,i] = total(w_ini[*,imin:imax],2)/(ajou+1)
+    e_buf[*,i] = SQRT(total((e_ini[*,imin:imax])^2,2))/(ajou+1)
+    y_buf[i] = total(y_ini[imin:imax])/(ajou+1)
+    rest = rest-1
   endif else begin
-    w_buf[*,i]=w_ini[*,imin]
-    e_buf[*,i]=e_ini[*,imin]
-    y_buf[i]=y_ini[imin]
-    rest=0
+    w_buf[*,i] = w_ini[*,imin]
+    e_buf[*,i] = e_ini[*,imin]
+    y_buf[i] = y_ini[imin]
+    rest = 0
   endelse
 endfor
 
-nogrp:e_out=e_buf*scale
-w_out=w_buf*scale
-if NOT(KEYWORD_SET(file_out)) THEN file_out='input_mupho.txt'
-if (n_elements(nsp) ne 1) THEN nsp=1
-if KEYWORD_SET(mupho)     THEN BEGIN
-HS_DOS_WRITE, w_out, DATP=datp,NSPECIES=nsp,T2MIN=tmin,T2MAX=tmax,SANGLES=sangles   $
-                                 , CHANNELS=Channel,CNT_EFF=icnt_eff,DW=dw,AMASS=amass,XS=xs,NPH=n_ph,UNTER=unt,$
-                                 FILE_OUT=file_out,CONC2=conc2,EMAX=emax
-ENDIF
+nogrp: e_out = e_buf*scale
+w_out = w_buf*scale
+if not(KEYWORD_SET(file_out)) then file_out = 'input_mupho.txt'
+if (n_elements(nsp) ne 1) then nsp = 1
+if KEYWORD_SET(mupho) then begin
+  HS_DOS_WRITE, w_out, DATP=datp, NSPECIES=nsp, T2MIN=tmin, T2MAX=tmax, SANGLES=sangles $
+                     , CHANNELS=Channel, CNT_EFF=icnt_eff, DW=dw, AMASS=amass, XS=xs $
+                     , NPH=n_ph, UNTER=unt, FILE_OUT=file_out, CONC2=conc2, EMAX=emax
+endif
 
-str_mess='GROUPING Successfull !: '+'NG = '+strcompress(string(ngroups),/remove_all)+' Theta_min = ' $
-         +strcompress(string(tmin),/remove_all)+' Theta_max = '+strcompress(string(tmax),/remove_all)
-if KEYWORD_SET(sint)  then str_mess=str_mess+' With sin(theta) weighting'
-if KEYWORD_SET(mupho) then str_mess=str_mess+' muphocor input file is : input_mupho.txt'
+str_mess = 'GROUPING Successfull !: ' + 'NG = ' + strcompress(string(ngroups),/remove_all) $
+         + ' Theta_min = ' + strcompress(string(tmin),/remove_all) + ' Theta_max = ' $
+         + strcompress(string(tmax),/remove_all)
+if KEYWORD_SET(sint)  then str_mess = str_mess+' With sin(theta) weighting'
+if KEYWORD_SET(mupho) then str_mess = str_mess+' muphocor input file is : input_mupho.txt'
 
-fin:print,str_mess
-x_out=x_in
-y_out=y_buf
+fin: print,str_mess
+x_out = x_in
+y_out = y_buf
 
 mod_datp, datp, "e", e_out
 mod_datp, datp, "y", y_out