Changeset 553

Timestamp:

Oct 6, 2011, 7:32:54 AM (13 years ago)

Author:

campagne

Message:

include image (jec)

Location:

BAORadio/AmasNancay/trunk

Files:

• mergeRawOnOff.cc (modified) (6 diffs)
• plotDiffOnOff.pic (modified) (1 diff)

BAORadio/AmasNancay/trunk/mergeRawOnOff.cc

@@ -40 +40 @@
 const sa_size_t NUMBER_OF_CHANNELS = 2;
 const sa_size_t NUMBER_OF_FREQ     = 8192;
-const sa_size_t TOTAL_NUM_CYCLES   = 10000;
+const sa_size_t TOTAL_NUM_CYCLES   = 500;
 const r_4    LOWER_FREQUENCY = 1250.0; //MHz
 const r_4    TOTAL_BANDWIDTH = 250.0; //MHz
@@ -56 +56 @@
 //Utility functions
 
+sa_size_t freqTochan(r_4 f){
+  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
+}
+//---------
   class median_of_empty_list_exception:public std::exception{
     virtual const char* what() const throw() {
@@ -309 +313 @@
     listOfSpectra.sort(stringCompare);
     //Loop of spectra matrix
-    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
+    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ;  ++iSpec){
       if (para.debuglev_>90){
         cout << " spactra <" << *iSpec << ">" << endl;
@@ -316 +320 @@
       fin.GetObject(aSpec,*iSpec);
 
-      if ((sa_size_t)nSpectra < TOTAL_NUM_CYCLES ) {
-        tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
-      } else {
-        stringstream tmp;
-        tmp<<nSpectra;
-        throw "FATAL: SIZE ERROR: too many cycles:"+tmp.str();
-      }
+      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
+
       nSpectra++;
     }//eo loop on spectra in a file
   }//eo loop on files
   
-  //Resize
 
  
@@ -347 +345 @@
   reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
 
+
+  sa_size_t f1320=freqTochan(1320.);
+  sa_size_t f1345=freqTochan(1345.);
+  sa_size_t f1355=freqTochan(1355.);
+  sa_size_t f1380=freqTochan(1380.);
+  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
+  if (para.debuglev_>9){
+    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
+  }
+  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    double meanMed1;
+    double sigmaMed1;
+    TVector<r_4> band1;
+    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
+    MeanSigma(band1,meanMed1,sigmaMed1);
+    double meanMed2;
+    double sigmaMed2;
+    TVector<r_4> band2;
+    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
+    MeanSigma(band2,meanMed2,sigmaMed2);
+    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
+  }  
+  meanMed.Print(cout);
+  
+  //Compute the sigma in the range 1320MHz-1380MHz
+  if (para.debuglev_>9){
+    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
+  }
+  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
+  sa_size_t redf1320=(1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
+  sa_size_t redf1380=(1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
+
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    double meanSigma;
+    double sigmaSigma;
+    TVector<r_4> band;
+    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
+    MeanSigma(band,meanSigma,sigmaSigma);
+    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
+    sigmaMed(iCh) = meanSigma;
+  }
+  sigmaMed.Print(cout);
+  
+  
+  if (para.debuglev_>9){
+    cout << "Compute medianOfSpectraNorm" << endl;
+  }
+  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
+    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
+  }
+
+  
+
   {//Save the result
     stringstream tmp;
@@ -356 +410 @@
     string tag = "median";
     fos << PPFNameTag(tag) << medianOfSpectra;
+
+    tag = "medianNorm";
+    fos << PPFNameTag(tag) << medianOfSpectraNorm;
+    
+
     tag = "meanmedred";
     fos << PPFNameTag(tag) << meanRedMtx;

BAORadio/AmasNancay/trunk/plotDiffOnOff.pic

@@ -90 +90 @@
 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
+
+
+##Normalized General Median
+del s0 s1
+objaoper medianNorm row 0 s0
+objaoper medianNorm row 1 s1
+graphicatt "xylimits=1250,1500,-2,2"
+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)" "Norm Median[(ON-OFF)[r,c]/Gain[r]]"  $defatt
+