Changeset 524 for BAORadio

Timestamp:

Sep 23, 2011, 10:18:17 AM (13 years ago)

Author:

campagne

Message:

analyse ON-OFF sans calibration

File:

• BAORadio/AmasNancay/analyse.cc (modified) (14 diffs)

BAORadio/AmasNancay/analyse.cc

@@ -10 +10 @@
 #include <stdlib.h>
 #include <dirent.h>
-
+#include <matharr.h>
 
 // include standard c/c++
@@ -202 +202 @@
   virtual int processCmd() throw(string);
 };
+
+//JEC 22/9/11 Make ON-OFF analysis WO any calibration START
+//------------
+//Process ON/OFF Raw data
+//------------
+class ProcessONOFFRawData : public ProcessBase {
+
+public:
+  ProcessONOFFRawData(){}
+  virtual ~ProcessONOFFRawData(){}
+  
+  virtual int processCmd() throw(string);
+};
+//JEC 22/9/11 Make ON-OFF analysis WO any calibration END
+
 //------------
 //Process Gain
@@ -263 +278 @@
   //Init process types
   map<string,IProcess*> process;
+  //JEC 22/9/11 Make ON-OFF analysis WO any calibration START
+  process["rawOnOff"] = new ProcessONOFFRawData();
+  //JEC 22/9/11 Make ON-OFF analysis WO any calibration END
   process["dataOnOff"] = new ProcessONOFFData();
   process["gain"]      = new ProcessGain();
@@ -363 +381 @@
 
   //JEC 21/9/11 Give the input parameters START
-  cout << " action " << action << "\n"
+  cout << "Dump Iiitial parameters ............" << endl;
+  cout << " action = " << action << "\n"
        << " inputPath = " << inputPath << "\n" 
        << " outputPath = " << outputPath << "\n"
@@ -376 +395 @@
        << " numcycle = " << numcycle << "\n"
        << " calibrationOpt = " << calibrationOpt << endl;
+  cout << "...................................." << endl;
   //JEC 21/9/11 Give the input parameters END
 
@@ -388 +408 @@
     }
     
-
-
-
 
     //
@@ -450 +467 @@
       throw e.what();
     }
+
+    //JEC 22/9/11 Make ON-OFF analysis WO any calibration START
+//     try {
+//       ProcessONOFFRawData* procRawdata = dynamic_cast<ProcessONOFFRawData*>(process[action]);
+//     }
+//     catch(exception& e){
+//       throw e.what();
+//     }
+    //JEC 22/9/11 Make ON-OFF analysis WO any calibration END
+
 
     try {
@@ -543 +570 @@
 int Usage(bool flag) {
   cout << "Analyse.cc usage...." << endl;
-  cout << "analyse  -act <action_type>: dataOnOff, gain, calib\n"
+  cout << "analyse  -act <action_type>: dataOnOff, rawOnOff, gain, calib\n"
        << "         -inPath <path for input files: default='.'>\n"
        << "         -outPath <path for output files: default='.'>\n"
@@ -620 +647 @@
   return rc;
 }
+//JEC 22/9/11 Make ON-OFF analysis WO any calibration START
 //----------------------------------------------
-int ProcessONOFFData::processCmd() throw(string) {
+int ProcessONOFFRawData::processCmd() throw(string) {
   int rc = 0;
   try {
@@ -629 +657 @@
     throw s;
   }
-   if(debuglev_>0)cout << "Process Data" << endl;
+  if(debuglev_>0)cout << "Process Raw Data ON OFF" << endl;
   vector<string> modeList;
   modeList.push_back("On");
@@ -641 +669 @@
   //Process to get sucessively
   //Raw Spectra, 
-  //BAO Calibrated Spectra 
-  //and RT Calibrated Spectra
   //The pocesses are separated to allow intermediate save of results
 
@@ -665 +691 @@
     //
     for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
-      //TOBE FIXED      directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
       directoryName = "./" + mode + "/";
       stringstream sicycle;
@@ -697 +722 @@
     cout << "Save mean raw spectra" << endl;
     string fileName;
+    fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
+
+    POutPersist fos(fileName);
+    id=0;
+    iSpectreEnd = spectreCollect.end();
+    for (iSpectre = spectreCollect.begin();
+         iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
+      tag = "specRaw";
+      tag += (iSpectre->first).first;
+      stringstream sid;
+      sid << (iSpectre->first).second;
+      tag += sid.str();
+      if(debuglev_>9) {
+        cout << "save tag<" << tag << ">" << endl;
+      }
+      fos << PPFNameTag(tag) << iSpectre->second;
+    }
+  }//end of save fits
+  
+
+  //------------------------------------------
+  // Perform ON-OFF
+  //------------------------------------------
+  
+  map< sa_size_t, TMatrix<r_4> > diffCollect;
+  map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd;
+
+  TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero
+  r_4 nCycles = 0;
+  for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
+    nCycles++;
+    TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory 
+    TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory 
+    TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff;
+    diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) ));
+    diffMeanOnOff += diffOnOff;
+  }//end loop on cycle
+  if(nCycles>0) diffMeanOnOff/=nCycles;
+
+  //exctract channels and do the mean
+  TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0
+  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
+    meanOfChan += diffMeanOnOff.Row(iCh).Transpose();
+  }
+  meanOfChan /= (r_4)NUMBER_OF_CHANNELS;
+  
+
+
+  {//save diff ON-OFF on Raw data
+    if(debuglev_>0)cout << "save ON-OFF RAW spectra" << endl;
+    string fileName;
+    fileName = "./diffOnOffRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
+    POutPersist fos(fileName);
+    iDiffEnd = diffCollect.end();
+    id = 0;
+
+    //JEC 22/9/11 Mean & Sigma in 32-bins size START
+    sa_size_t nSliceFreq = 32; //TODO: put as an input parameter option ?
+    sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq;
+    //JEC 22/9/11 Mean & Sigma in 32-bins size END
+
+    for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
+      tag = "specONOFFRaw";
+      stringstream sid;
+      sid << iDiff->first;
+      tag += sid.str();
+      fos << PPFNameTag(tag) << iDiff->second;
+
+      //JEC 22/9/11 Mean & Sigma in 32-bins size START
+      if (debuglev_>9) {
+        cout << "Debug slicing: slice/delta " << nSliceFreq << " " << deltaFreq << endl;
+      }
+      TMatrix<r_4> reducedMeanDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
+      TMatrix<r_4> reducedSigmaDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
+      for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
+        sa_size_t freqLow= iSlice*deltaFreq;
+        sa_size_t freqHigh= freqLow + deltaFreq -1;
+        if (debuglev_>9) {
+          cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl;
+        }
+        for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+          if (debuglev_>9) {
+            cout << "Debug .......... Channel " << iCh;
+          }
+          TVector<r_4> reducedRow;
+          reducedRow = (iDiff->second).SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
+          double mean; 
+          double sigma;
+          MeanSigma(reducedRow,mean,sigma);
+          if (debuglev_>9) {
+            cout << "mean/signa " << mean << "/" << sigma << endl;
+          }
+          reducedMeanDiffOnOff(iCh,iSlice) = mean;
+          reducedSigmaDiffOnOff(iCh,iSlice) = sigma;
+        }//loop on Channel
+      }//loop on Freq. slice 
+      tag = "redMeanONOFFRaw";
+      tag += sid.str();
+      fos << PPFNameTag(tag) << reducedMeanDiffOnOff;
+      tag = "redSigmaONOFFRaw";
+      tag += sid.str();
+      fos << PPFNameTag(tag) << reducedSigmaDiffOnOff;
+      //JEC 22/9/11 END
+
+    }//loop on ON-OFF spectre
+    //save the mean also
+    fos << PPFNameTag("specONOFFRawMean") << diffMeanOnOff;
+
+    //JEC 22/9/11 START
+    TMatrix<r_4> reducedMeanDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
+    TMatrix<r_4> reducedSigmaDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
+    for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
+      sa_size_t freqLow= iSlice*deltaFreq;
+      sa_size_t freqHigh= freqLow + deltaFreq -1;
+      if (debuglev_>9) {
+        cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl;
+      }
+      for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
+        if (debuglev_>9) {
+          cout << "Debug .......... Channel " << iCh;
+        }
+        TVector<r_4> reducedRow;
+        reducedRow = diffMeanOnOff.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
+        double mean; 
+        double sigma;
+        MeanSigma(reducedRow,mean,sigma);
+        if (debuglev_>9) {
+          cout << "mean/signa " << mean << "/" << sigma << endl;
+        }
+        reducedMeanDiffOnOffAll(iCh,iSlice) = mean;
+        reducedSigmaDiffOnOffAll(iCh,iSlice) = sigma;
+      }//loop on Channel
+    }//loop on Freq. slice 
+    tag = "redMeanONOFFRawAll";
+    fos << PPFNameTag(tag) << reducedMeanDiffOnOffAll;
+    tag = "redSigmaONOFFRawAll";
+    fos << PPFNameTag(tag) << reducedSigmaDiffOnOffAll;
+    //JEC 22/9/11 END
+
+
+
+    fos << PPFNameTag("specONOFFRaw2ChanMean") << meanOfChan;
+  }//end of save fits
+
+  return rc;
+} //ProcessONOFFRawData::processCmd
+
+//JEC 22/9/11 Make ON-OFF analysis WO any calibration END
+//----------------------------------------------
+int ProcessONOFFData::processCmd() throw(string) {
+  int rc = 0;
+  try {
+    rc = ProcessBase::processCmd();
+  } 
+  catch (string s) {
+    throw s;
+  }
+   if(debuglev_>0)cout << "Process Data" << endl;
+  vector<string> modeList;
+  modeList.push_back("On");
+  modeList.push_back("Off");
+  vector<string>::const_iterator iMode;
+  
+  uint_4 id; 
+  string tag;
+
+  //
+  //Process to get sucessively
+  //Raw Spectra, 
+  //BAO Calibrated Spectra 
+  //and RT Calibrated Spectra
+  //The pocesses are separated to allow intermediate save of results
+
+  map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
+  map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
+  
+  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
+    string mode = *iMode;
+    if(debuglev_>0)cout << "Process RAW Mode " << mode << endl;
+
+    //------------------------------------------
+    //Produce Raw spectra per cycle
+    //------------------------------------------
+
+    string directoryName;
+    list<string> listOfSpecFiles;
+    list<string>::const_iterator iFile, iFileEnd;
+    
+        
+    //
+    //loop on cycles
+    //
+    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
+      //TOBE FIXED      directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
+      directoryName = "./" + mode + "/";
+      stringstream sicycle;
+      sicycle << icycle;
+      directoryName += spectraDirectory_ + sicycle.str() + "/";
+
+      //read directory
+      listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
+      
+
+      //compute mean of spectra created in a cycle
+      if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
+      TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
+      iFileEnd = listOfSpecFiles.end();
+      r_4 nSpectres  = 0;
+      for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
+        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
+        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+        aSpectrum >> spectre;
+        spectreMean += spectre;
+        nSpectres++;
+      }// end of for files 
+      if(nSpectres>0) spectreMean /= nSpectres;
+      
+      //save mean spectrum
+      spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
+    }//end of for cycles
+  }//end loop on mode for raw preocess
+
+  if(debuglev_>1) {//save mean spectra on file 
+    cout << "Save mean raw spectra" << endl;
+    string fileName;
     //TOBE FIXED    fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataRaw" + ".ppf";
     fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
@@ -734 +984 @@
     //Read BAO calibration files
     sa_size_t nr,nc; //values read 
-//JEC 20/9/11 use mean calibration coeff upon all cycles START
-//    bool first = true;
-//    TMatrix<r_4> calibBAOfactors;
-//     for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
-//       string calibFileName = inputPath_+ "/" 
-//      + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) 
-//      + "/calib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
-//      + mode + "_" + freqBAOCalibration_ + "MHz";
-//       stringstream channel;
-//       channel << iCh;
-//       calibFileName += "-Ch" + channel.str() + "Cycles.txt";
-//       if(debuglev_>0) cout << "Read file " << calibFileName << endl;
-//       ifstream ifs(calibFileName.c_str());
-//       if ( ! ifs.is_open() ) {
-//      rc = 999;
-//      throw calibFileName + " cannot be opened...";
-//       }      
-//       TVector<r_4> aCalib;
-//       if(debuglev_>9) cout << "Debug 1" << endl;
-//       aCalib.ReadASCII(ifs,nr,nc);
-//       if(debuglev_>9) cout << "Debug 2" << endl;
-//       if(first){
-//      first = false;
-//      if(debuglev_>9) cout << "Debug 3" << endl;
-//      calibBAOfactors.SetSize(NUMBER_OF_CHANNELS,nr);
-//      if(debuglev_>9) cout << "Debug 4" << endl;
-//       }
-//       calibBAOfactors( Range(iCh), Range::all() ) = aCalib.Transpose();
-//       if(debuglev_>9) cout << "Debug 5" << endl;
-//       ifs.close();
-//     }//end of loop on channels
-
     
     string calibFileName = inputPath_+ "/"