##############################################################
# Check chromaticity of DAB from calib spectra (1204-window) #
##############################################################

set user     "AST"
set toppath  "/sps/baoradio/AmasNancay/${user}"
set source   "Abell1205"
set srclower "abell1205"
set date     "20110415"
set cycle    "8"
set mode     "Off"
set path "${toppath}/${source}/${date}${srclower}/${mode}/calibcycle${cycle}"

#
# Power vs. run (calib monitoring)
#-----------------------------------
#shell ls ${path}/*.fits 
shell ls ${path}/*.fits | wc -l > pp
vecfrascii mm pp
vec2var mm nfiles

set mean0 " "
set mean1 " "
for i 0:$nfiles
  openfits ${path}/medfiltmtx${i}.fits
  del g0 g1 s0 s1 n0 n1
  objaoper /home/medfiltmtx${i} row 0 g0
  objaoper /home/medfiltmtx${i} row 1 g1
  c++exec sa_size_t n0=g0.NElts()-1; TVector<r_4>s0=g0(Range(1,n0)); KeepObj(s0); \
          sa_size_t n1=g1.NElts()-1; TVector<r_4>s1=g1(Range(1,n1)); KeepObj(s1);
  m0 = ${s0.sum}/${s0.size}
  m1 = ${s1.sum}/${s1.size}
  set mean0 "${mean0} ${m0}"
  set mean1 "${mean1} ${m1}"
end
echo ${mean0}
echo ${mean1}
line2vec vmean0 ${mean0}
line2vec vmean1 ${mean1}

maxplot = ${vmean1.max}+5
newwin 1 1
graphicatt "xylimits=0,116,0,3"
plot2d vmean1 n val 1 "red cpts nsta notit"
plot2d vmean0 n val 1 "same blue cpts nsta notit"
settitle "Mean calib $source ${date} ${mode} Ch 0 (blue) Ch 1 (red)"
setaxelabels "File #" "I (a.u)"

#
# Identify baseline, DAB_1, DAB_2 from monitoring plot
#------------------------------------------------------
## Baseline
del ifirst count nspec sum0 sum1
ifirst = 1
count = 0
for i 0:34
#  echo "^^^^^^^^" $i $count
  count = $count+1
  del g0 g1
  objaoper /home/medfiltmtx${i} row 0 g0
  objaoper /home/medfiltmtx${i} row 1 g1
  if ( $ifirst == 1 ) then
    ifirst = 0
    c++exec TVector<r_4> sum0=g0; KeepObj(sum0); \
            TVector<r_4> sum1=g1; KeepObj(sum1);
  else
    c++exec sum0+=g0; sum1+=g1;
  endif
end
for i 70:$nfiles
  count = $count+1
#  echo "^^^^^^^^" $i $count
  del g0 g1
  objaoper /home/medfiltmtx${i} row 0 g0
  objaoper /home/medfiltmtx${i} row 1 g1
  c++exec sum0+=g0; sum1+=g1;
end
echo $count
line2vec nspec $count
del base0 base1
c++exec TVector<r_4> base0=sum0; KeepObj(base0); base0/=nspec[0]; \
        TVector<r_4> base1=sum1; KeepObj(base1); base1/=nspec[0];

## DAB Level 1
idablev1 = 36
idablev2 = 51

del ifirst count nspec sum0 sum1
ifirst = 1
count = 0
for i $idablev1:$idablev2 
  count = $count+1
  echo "^^^^^^^^" $i $count
  del g0 g1
  objaoper /home/medfiltmtx${i} row 0 g0
  objaoper /home/medfiltmtx${i} row 1 g1
  if ( $ifirst == 1 ) then
    ifirst = 0
    c++exec TVector<r_4> sum0=g0; KeepObj(sum0); \
            TVector<r_4> sum1=g1; KeepObj(sum1);
  else
    c++exec sum0+=g0; sum1+=g1;
  endif
end
echo $count
line2vec nspec $count
del dab01 dab11
c++exec TVector<r_4> dab01=sum0; KeepObj(dab01); dab01/=nspec[0]; \
        TVector<r_4> dab11=sum1; KeepObj(dab11); dab11/=nspec[0];

## DAB Level 2
idablev3 = 53
idablev4 = 67

del ifirst count nspec sum0 sum1
ifirst = 1
count = 0
for i $idablev3:$idablev4
  count = $count+1
  echo "^^^^^^^^" $i $count
  del g0 g1
  objaoper /home/medfiltmtx${i} row 0 g0
  objaoper /home/medfiltmtx${i} row 1 g1
  if ( $ifirst == 1 ) then
    ifirst = 0
    c++exec TVector<r_4> sum0 = g0; KeepObj(sum0); \
            TVector<r_4> sum1 = g1; KeepObj(sum1);
  else
    c++exec sum0+=g0; sum1+=g1;
  endif
end
echo $count
line2vec nspec $count
del dab02 dab12
c++exec TVector<r_4> dab02=sum0; KeepObj(dab02); dab02/=nspec[0]; \
        TVector<r_4> dab12=sum1; KeepObj(dab12); dab12/=nspec[0];

## Plot results
newwin 1 1
plot2d base0 (n/8192)*250+1250 val 1 "blue cpts notit nstat"
plot2d base1 (n/8192)*250+1250 val 1 "same red cpts notit nstat"
settitle "Mean baseline $source $date Ch0 blue Ch1 red"

newwin 1 1
graphicatt "xylimits=1250,1500,0,3"
plot2d dab01 (n/8192)*250+1250 val 1 "blue cpts notit nstat"
plot2d dab11 (n/8192)*250+1250 val 1 "same red cpts notit nstat"
settitle "Mean First DAB $source $date Ch0 blue Ch1 red"

newwin 1 1
graphicatt "xylimits=1250,1500,0,3"
plot2d dab02 (n/8192)*250+1250 val 1 "blue cpts notit nstat"
plot2d dab12 (n/8192)*250+1250 val 1 "same red cpts notit nstat"
settitle "Mean Second DAB $source $date Ch0 blue Ch1 red"

## Save results in ppf files
del mtxbase mtxdab1 mtxdab2

c++exec TMatrix<r_4> mtxbase(2,base0.NElts()); mtxbase.Row(0)=base0; mtxbase.Row(1)=base1; KeepObj(mtxbase); \
        TMatrix<r_4> mtxdab1(2,dab01.NElts()); mtxdab1.Row(0)=dab01; mtxdab1.Row(1)=dab11; KeepObj(mtxdab1); \
        TMatrix<r_4> mtxdab2(2,dab02.NElts()); mtxdab2.Row(0)=dab02; mtxdab2.Row(1)=dab12; KeepObj(mtxdab2);

saveppf mtxbase ${toppath}/${source}/${date}${srclower}/baseline_${date}_${srclower}_${mode}.ppf
saveppf mtxdab1 ${toppath}/${source}/${date}${srclower}/dab1_${date}_${srclower}_${mode}.ppf
saveppf mtxdab2 ${toppath}/${source}/${date}${srclower}/dab2_${date}_${srclower}_${mode}.ppf