Changeset 544

Timestamp:

Oct 5, 2011, 9:51:55 AM (13 years ago)

Author:

campagne

Message:

introduce median calculous (jec)

File:

• BAORadio/AmasNancay/mergeRawOnOff.cc (modified) (10 diffs)

BAORadio/AmasNancay/mergeRawOnOff.cc

@@ -39 +39 @@
 //-----------------------------------------------
 const sa_size_t NUMBER_OF_CHANNELS = 2;
-const sa_size_t  NUMBER_OF_FREQ = 8192;
+const sa_size_t NUMBER_OF_FREQ     = 8192;
+const sa_size_t TOTAL_NUM_CYCLES   = 10000;
 const r_4    LOWER_FREQUENCY = 1250.0; //MHz
 const r_4    TOTAL_BANDWIDTH = 250.0; //MHz
 //-----------------------------------------------
+//Input parameters
+struct Param {
+  int debuglev_;      //debug
+  string inPath_;     //root directory of the input files
+  string outPath_;    //output files are located here
+  string sourceName_; //source name & subdirectory of the input files
+  string ppfFile_;    //generic name of the input files
+  int nSliceInFreq_;  //used by reduceSpectra() fnc
+} para;
+//--------------------------------------------------------------
+//Utility functions
+//--------------------------------------------------------------
 char *regexp (const char *string, const char *patrn, int *begin, int *end) {    
         int i, w=0, len;                 
@@ -62 +75 @@
         return word;
 }
-//------------------------------------------------ 
-
-
-struct Param {
-  int debuglev_;
-  string inPath_;
-  string outPath_;
-  string sourceName_;
-  string ppfFile_;
-} para;
-//-----
+//------- 
 sa_size_t round_sa(r_4 r) {
   return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
@@ -145 +148 @@
     
     if (b != 0) return false;
-    if (e != aStr.size()) return false;
+    if (e != (int)aStr.size()) return false;
     return true;
 
@@ -153 +156 @@
 };
 //-------------------------------------------------------
+//Rebin in frequence + compute mean and sigma
+void reduceSpectra(const TMatrix<r_4>& specMtxInPut, 
+                    TMatrix<r_4>& meanMtx, 
+                    TMatrix<r_4>& sigmaMtx) {
+  sa_size_t nSliceFreq = para.nSliceInFreq_;
+  sa_size_t deltaFreq =  NUMBER_OF_FREQ/nSliceFreq;
+  for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
+    sa_size_t freqLow= iSlice*deltaFreq;
+    sa_size_t freqHigh= freqLow + deltaFreq -1;
+    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+      TVector<r_4> reducedRow;
+      reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
+      double mean; 
+      double sigma;
+      MeanSigma(reducedRow,mean,sigma);
+      meanMtx(iCh,iSlice) = mean;
+      sigmaMtx(iCh,iSlice) = sigma;
+    }//eo loop on channels
+  }//eo loop on slices
+}
 //-------------------------------------------------------
+//Compute the mean of Raw spectra and also the mean/sigma of rebinned spectra 
+//Used like:
+//
 void meanOnCycles() throw(string) {
   list<string> listOfFiles;
@@ -166 +192 @@
   
   StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
-  TMatrix<r_4> sumOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+  TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
   uint_4 nSpectra=0;
   //Loop on files
@@ -193 +219 @@
       //How to see if the GetObject is ok?? Ask Reza
       nSpectra++;
-      sumOfSpectra+=aSpec;
+      meanOfSpectra+=aSpec;
     }//eo loop on spectra in a file
   }//eo loop on files
   
   //Normalisation
-  if(nSpectra>0)sumOfSpectra/=(r_4)(nSpectra);
+  if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
+
+  //Compute the reduced version of the mean and sigma
+  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
 
   {//Save the result
@@ -205 +236 @@
     string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
     cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
+    POutPersist fos(fileName);
+
     string tag = "mean";
+    fos << PPFNameTag(tag) << meanOfSpectra;
+    tag = "meanred";
+    fos << PPFNameTag(tag) << meanRedMtx;
+    tag = "sigmared";
+    fos << PPFNameTag(tag) << sigmaRedMtx;
+  }
+}
+//-------------------------------------------------------
+//Compute the median of Raw spectra and also the mean/sigma of rebinned spectra 
+//Used like:
+//
+void medianOnCycles() throw(string) {
+  list<string> listOfFiles;
+  string directoryName;
+  directoryName = para.inPath_ + "/" + para.sourceName_;
+
+  //Make the listing of the directory
+  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
+  
+  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
+  iFileEnd = listOfFiles.end();
+  
+
+  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough...
+
+  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
+  uint_4 nSpectra=0;
+  //Loop on files
+  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
+    if (para.debuglev_>90){
+      cout << "load file <" << *iFile << ">" << endl;
+    }
+    PInPersist fin(*iFile);
+    vector<string> vec = fin.GetNameTags();
+    list<string> listOfSpectra;
+    //Keep only required PPF objects
+    std::remove_copy_if(
+                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
+                        not1(match)
+                        );
+   
+    iSpecEnd = listOfSpectra.end();
+    listOfSpectra.sort(stringCompare);
+    //Loop of spectra matrix
+    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
+      if (para.debuglev_>90){
+        cout << " spactra <" << *iSpec << ">" << endl;
+      }
+      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      fin.GetObject(aSpec,*iSpec);
+
+      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;
+        tmp<<nSpectra;
+        throw "FATAL: SIZE ERROR: too many cycles:"+tmp.str();
+      }
+      nSpectra++;
+    }//eo loop on spectra in a file
+  }//eo loop on files
+  
+  //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());
+      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);
+    }
+  }
+
+
+  //Compute the reduced version of the mean and sigma
+  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
+  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
+
+  {//Save the result
+    stringstream tmp;
+    tmp << nSpectra;
+    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+    cout << "Save median based on " <<  (nSpectra+1) << " cycles " << endl;
     POutPersist fos(fileName);
-    fos << PPFNameTag(tag) << sumOfSpectra;
-  }
-}
+
+    string tag = "median";
+    fos << PPFNameTag(tag) << medianOfSpectra;
+    tag = "meanmedred";
+    fos << PPFNameTag(tag) << meanRedMtx;
+    tag = "sigmamedred";
+    fos << PPFNameTag(tag) << sigmaRedMtx;
+  }
+}
+
 //-------------------------------------------------------
-int main(int narg, char* arg[])
-{
+int main(int narg, char* arg[]) {
   
   int rc = 0; //return code
@@ -218 +357 @@
 
 
+
+  //default value for initial parameters (see Para structure on top of the file)
   string debuglev = "0";
-  string action = "default";
+  string action = "meanDiffOnOff";
   string inputPath = "."; 
   string outputPath = "."; 
   string sourceName = "Abell85";
   string ppfFile;
+  string nSliceInFreq = "32";
  
   
@@ -253 +395 @@
       ka+=2;
     }
+    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
+      nSliceInFreq=arg[ka+1];
+      ka+=2;
+    }
     else ka++;
   }//eo while
 
-  para.debuglev_ = atoi(debuglev.c_str());
-  para.inPath_   = inputPath;
-  para.outPath_  = outputPath;
+  para.debuglev_   = atoi(debuglev.c_str());
+  para.inPath_     = inputPath;
+  para.outPath_    = outputPath;
   para.sourceName_ = sourceName;
-  para.ppfFile_ = ppfFile;
+  para.ppfFile_    = ppfFile;
+  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
 
   cout << "Dump Initial parameters ............" << endl;
   cout << " action = " << action << "\n"
-       << " inputPath = " << inputPath << "\n" 
-       << " outputPath = " << outputPath << "\n"
-       << " sourceName = " << sourceName << "\n"
-       << " ppfFile = " << ppfFile << "\n"
-       << " debuglev = "  << debuglev  << "\n";
+       << " inputPath = " << para.inPath_  << "\n" 
+       << " outputPath = " <<  para.outPath_ << "\n"
+       << " sourceName = " << para.sourceName_ << "\n"
+       << " ppfFile = " <<  para.ppfFile_ << "\n"
+       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
+       << " debuglev = "  << para.debuglev_   << "\n";
   cout << "...................................." << endl;
 
@@ -284 +432 @@
 //     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
 
-    if ( action == "default" ) {
+    if ( action == "meanDiffOnOff" ) {
       meanOnCycles();
-    }
+    } else if (action == "medianDiffOnOff") {
+      medianOnCycles();
+    } else {
+      msg = "Unknown action " + action;
+      throw msg;
+    }
+
 
   }  catch (std::exception& sex) {