Changeset 602

Timestamp:

Nov 14, 2011, 2:28:01 PM (13 years ago)

Author:

campagne

Message:

rewrite meanRawDiffOnOffCycles (jec)

Location:

BAORadio/AmasNancay/trunk

Files:

• makefile (modified) (1 diff)
• mergeAnaFiles.cc (modified) (4 diffs)

BAORadio/AmasNancay/trunk/makefile

@@ -7 +7 @@
 PKGOLIST =
 
-EXELIST = $(EXE)analyse $(EXE)mergeAnaFiles $(EXE)anaTools
+EXELIST = $(EXE)analyse $(EXE)mergeAnaFiles
 
-EXEOLIST = $(EXE)analyse.o $(EXE)mergeAnaFiles.o $(EXE)anaTools.o
+EXEOLIST = $(EXE)analyse.o $(EXE)mergeAnaFiles.o
 
 
-all : analyse mergeAnaFiles anaTools
+all : analyse mergeAnaFiles
 
 clean :

BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

@@ -38 +38 @@
 //Usage
 //
+//./Objs/mergeAnaFiles -act meanRawDiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src NGC4383 -ppfFile dataRaw -debug 1000 -mxcycles 10
 //./Objs/mergeAnaFiles -act meanCalibBAODiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src Abell85 -ppfFile dataRaw -debug 1 -mxcycles 500
 //./Objs/mergeAnaFiles -src Abell85 -act meanCalibBAODiffOnOff -mxcycles 500 -calibfreq 1410 -inPath /sps/baoradio/AmasNancay/JEC/ -ppfFile dataRaw -debug 1 >& mergeAna-500.log
@@ -644 +645 @@
       aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
       
+      //Load OFF phase
       ppftag = "specRawOff"+cycle.str();
       TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
@@ -823 +825 @@
   }//end of save
 }
+//JEC 14/11/11 New meanRawDiffOnOffCycles START
 //-------------------------------------------------------
-//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra 
+//Compute the mean of Diff ON-OFF/OFF from Raw spectra
 //Used like:
 //
@@ -837 +840 @@
   list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
   iFileEnd = listOfFiles.end();
-  
-  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
-  TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-  uint_4 nSpectra=0;
+
+  //mean ON-OFF over the list of cycles
+  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
+  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
+
+  int totalNumberCycles=0; //total number of cycles
+
   //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)
-                        );
+
+    vector<string> modeList;
+    modeList.push_back("On");
+    modeList.push_back("Off");
+    vector<string>::const_iterator iMode;
+    
+    map<string, list<int> > cycleModeCollect;
+    
+    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
+      list<string> listOfSpectra;
+      //Keep only required PPF objects
+      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
+      std::remove_copy_if(
+                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
+                          not1(StringMatch(matchstr))
+                          );
+      
+      listOfSpectra.sort(stringCompare);
+      iSpecEnd = listOfSpectra.end();
+      
+      matchstr = "[0-9]+";
+      //Loop of spectra matrix
+      list<int> listOfCycles;
+      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
+        int b,e;
+        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
+        if (para.debuglev_>90){
+          cout << " spactra <" << *iSpec << ">" << endl;
+          cout << " cycle " << iSpec->substr(b) << endl;
+        }
+        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
+      }//end loop spectra
+      cycleModeCollect[*iMode] = listOfCycles;
+    }//end of mode    
+
+    //Take the Intersection of the list Of cycles in mode Off and On    
+    list<int> commonCycles;
+    set_intersection(cycleModeCollect["On"].begin(),
+                     cycleModeCollect["On"].end(),
+                     cycleModeCollect["Off"].begin(),
+                     cycleModeCollect["Off"].end(),
+                     back_inserter(commonCycles)
+                     );
+    
+    if (para.debuglev_>90){
+      cout << "Liste of cycles common to On & Off: <";
+      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
+        cout << *i << " ";
+      }
+      cout << ">" << endl;
+    }
+
+    //Loop on cyles
+    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
+      
+      string ppftag;
+      //load ON phase
+      stringstream cycle;
+      cycle << *ic;
+      ppftag = "specRawOn"+cycle.str();
+      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      fin.GetObject(aSpecOn,ppftag);
+
+      //Load OFF phase
+      ppftag = "specRawOff"+cycle.str();
+      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      fin.GetObject(aSpecOff,ppftag);
+      
+      //Perform the difference ON-OFF
+      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
+      meanDiffONOFF_noCalib += diffOnOff_noCalib;
+      
+      //JEC 29/10/11 add ON-OFF/OFF
+      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
+      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
+      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
+      
+      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
+      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
+
+      totalNumberCycles++;
+
+      //Quit if enough
+      if (totalNumberCycles >= para.maxNumberCycles_) break;    
+    }//end of cycles
+
+    if (totalNumberCycles >= para.maxNumberCycles_) break;          
+
+  }//end files
+  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
+  //Normalisation
+  if(totalNumberCycles > 0){
+    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
+    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
+  }  
+
+  {//save results
+    stringstream tmp;
+    tmp << totalNumberCycles;
+    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+
+    POutPersist fos(fileName);
+    cout << "Save output in " << fileName << endl;
+
+    string tag = "meanNoCalib";
+    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
+    
+    tag = "meanOvOffNoCalib";
+    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
    
-    listOfSpectra.sort(stringCompare);
-    iSpecEnd = listOfSpectra.end();
-    //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);
-      //How to see if the GetObject is ok?? Ask Reza
-      nSpectra++;
-      meanOfSpectra+=aSpec;
-    }//eo loop on spectra in a file
-  }//eo loop on files
-  
-  //Normalisation
-  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
-    stringstream tmp;
-    tmp << nSpectra;
-    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;
-  }
-}
+  }//end save
+}
+//JEC 14/11/11 New meanRawDiffOnOffCycles END
+//-------------------------------------------------------
+//JEC 14/11/11 Obsolete START
+//-------------------------------------------------------
+//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra 
+//Used like:
+//
+// void meanRawDiffOnOffCycles() 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();
+  
+//   StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
+//   TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+//   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)
+//                      );
+   
+//     listOfSpectra.sort(stringCompare);
+//     iSpecEnd = listOfSpectra.end();
+//     //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);
+//       //How to see if the GetObject is ok?? Ask Reza
+//       nSpectra++;
+//       meanOfSpectra+=aSpec;
+//     }//eo loop on spectra in a file
+//   }//eo loop on files
+  
+//   //Normalisation
+//   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
+//     stringstream tmp;
+//     tmp << nSpectra;
+//     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;
+//   }
+// }
+//JEC 14/11/11 Obsolete END
 //-------------------------------------------------------
 //Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra