Changeset 552

Timestamp:

Oct 5, 2011, 3:13:36 PM (13 years ago)

Author:

campagne

Message:

improvment (jec)

Location:

BAORadio/AmasNancay/trunk

Files:

• mergeRawOnOff.cc (modified) (5 diffs)
• plotDiffOnOff.pic (modified) (4 diffs)

BAORadio/AmasNancay/trunk/mergeRawOnOff.cc

@@ -55 +55 @@
 //--------------------------------------------------------------
 //Utility functions
+
+  class median_of_empty_list_exception:public std::exception{
+    virtual const char* what() const throw() {
+    return "Attempt to take the median of an empty list of numbers.  "
+      "The median of an empty list is undefined.";
+  }
+  };
+  template<class RandAccessIter>
+    double median(RandAccessIter begin, RandAccessIter end) 
+    throw(median_of_empty_list_exception){
+    if(begin == end){ throw median_of_empty_list_exception(); }
+    std::size_t size = end - begin;
+    std::size_t middleIdx = size/2;
+    RandAccessIter target = begin + middleIdx;
+    std::nth_element(begin, target, end);
+    
+    if(size % 2 != 0){ //Odd number of elements
+      return *target;
+    }else{            //Even number of elements
+      double a = *target;
+      RandAccessIter targetNeighbor= target-1;
+      std::nth_element(begin, targetNeighbor, end);
+      return (a+*targetNeighbor)/2.0;
+    }
+  }
+
 //--------------------------------------------------------------
 char *regexp (const char *string, const char *patrn, int *begin, int *end) {    
@@ -291 +317 @@
 
       if ((sa_size_t)nSpectra < TOTAL_NUM_CYCLES ) {
-        //STORE THE spectra
-        //      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)).Show();
-        //      aSpec.Show();
         tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
-//      aSpec.Show();
-//      cout << "deb 1" << endl;
-//      tmp = aSpec;
-//      cout << "deb 2" << endl;
       } else {
         stringstream tmp;
@@ -310 +329 @@
   //Resize
 
-  nSpectra--;//For further selection in Range
-  if (para.debuglev_>90){
-    cout << "Cycles < 0," <<  nSpectra << ">" << endl;
-  }
  
-
   TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
   //Compute the median for each freq. and Channel
   for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
     for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
-      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra)).CompactAllDimensions());
+      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
       vector<r_4> tmp;
       tmp0.FillTo(tmp);
-      nth_element(tmp.begin(),tmp.begin()+tmp.size()/2,tmp.end());
-      medianOfSpectra(iCh,freq) = *(tmp.begin()+tmp.size()/2);
+      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
     }
   }
@@ -338 +351 @@
     tmp << nSpectra;
     string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
-    cout << "Save median based on " <<  (nSpectra+1) << " cycles " << endl;
+    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
     POutPersist fos(fileName);
 
@@ -347 +360 @@
     tag = "sigmamedred";
     fos << PPFNameTag(tag) << sigmaRedMtx;
+    tag = "cycleVsfreq";
+    
+    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
+    fos << PPFNameTag(tag) << tarr;
   }
 }

BAORadio/AmasNancay/trunk/plotDiffOnOff.pic

@@ -6 +6 @@
 set ncycles $2
 set nFreqBin $3
+set action $4
+
+set defatt "font=helvetica,bold,20"
+
+
+
+if ( ${action} == "mean" ) then
+
 openppf meanDiffOnOffRaw_${source}-${ncycles}Cycles.ppf
 
@@ -13 +21 @@
 objaoper mean row 0 s0
 objaoper mean row 1 s1
-graphicatt "xylimits=1250,1500,-0.01,0.01"
+graphicatt "${defatt} xylimits=1250,1500,-0.01,0.01"
 newwin 1 1
 plot2d s0 (n/8192)*250+1250 val n>0 "cpts blue nsta notit"
 plot2d s1 (n/8192)*250+1250 val n>0 "same cpts red nsta notit"
-settitle "<(ON-OFF)[r,c]/Gain[r]>_{r,c} $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt
+setaxelabels "Freq. (MHz)" "<(ON-OFF)[r,c]/Gain[r]>_{r,c}" $defatt
+
 
 
@@ -27 +37 @@
 plot2d rs0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit"
 plot2d rs1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit"
-settitle "<(ON-OFF)[r,c]/Gain[r]>_{r,c} Rebin(${nFreqBin}) $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt
+setaxelabels "Freq. (MHz)" "<(ON-OFF)[r,c]/Gain[r]>_{r,c} Rebin(${nFreqBin})" $defatt
+
 
 ##Reduced Mean
@@ -36 +48 @@
 plot2d rss0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit"
 plot2d rss1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit"
-settitle "Sigma[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin}) $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 
+setaxelabels "Freq. (MHz)" "Sigma[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt
+
+else
+
+openppf medianDiffOnOffRaw_${source}-${ncycles}Cycles.ppf
 
 
+##General Median
+del s0 s1
+objaoper median row 0 s0
+objaoper median row 1 s1
+graphicatt "xylimits=1250,1500,-0.01,0.01"
+newwin 1 1
+plot2d s0 (n/8192)*250+1250 val n>0 "cpts blue nsta notit"
+plot2d s1 (n/8192)*250+1250 val n>0 "same cpts red nsta notit"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt
+setaxelabels "Freq. (MHz)" "Median[(ON-OFF)[r,c]/Gain[r]]"  $defatt
+
+
+##Reduced Mean of the median
+newwin 1 1
+del rs0 rs1
+objaoper meanmedred row 0 rs0
+objaoper meanmedred row 1 rs1
+graphicatt "xylimits=1250,1500,-0.01,0.01"
+plot2d rs0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit"
+plot2d rs1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt
+setaxelabels "Freq. (MHz)" "Mean Median[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt
+
+
+##Reduced Sigma of the median
+newwin 1 1
+del rss0 rss1
+objaoper sigmamedred row 0 rss0
+objaoper sigmamedred row 1 rss1
+graphicatt "xylimits=1250,1500,0,0.01"
+plot2d rss0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit"
+plot2d rss1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit"
+settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt
+setaxelabels "Freq. (MHz)" "Sigma Median[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt
+ 
+
+endif
+