Changeset 563

Timestamp:

Oct 10, 2011, 4:31:37 PM (13 years ago)

Author:

campagne

Message:

improve mean ON-OFF with Calib BAO (jec)

Location:

BAORadio/AmasNancay/trunk

Files:

• etude_mergeAna.pic (added)
• mergeAnaFiles.cc (modified) (15 diffs)
• plotDiffOnOff.pic (modified) (1 diff)

BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

@@ -36 +36 @@
 
 //-----------------------------------------------
+//Usage
+//
+// ./Objs/mergeAnaFiles -act meanCalibBAODiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src Abell85 -ppfFile dataRaw -debug 1 -mxcycles 500
+//
+//
+//
+//-----------------------------------------------
 const sa_size_t NUMBER_OF_CHANNELS = 2;
 const sa_size_t NUMBER_OF_FREQ     = 8192;
 const r_4    LOWER_FREQUENCY       = 1250.0; //MHz
 const r_4    TOTAL_BANDWIDTH       = 250.0;  //MHz
-//-----------------------------------------------
 //Input parameters
 struct Param {
@@ -68 +74 @@
                  sa_size_t chHigh, 
                  TVector<r_4>& vec){
+  
   sa_size_t nr = mtx.NRows();
+  
   for (sa_size_t ir=0; ir<nr; ir++){
     TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
@@ -242 +250 @@
 //-------------------------------------------------------
 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra 
-//Used like: 
 //
-// void meanCalibBAODiffOnOffCycles() 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, iSpecOff, iSpecOffEnd, iSpecOn, iSpecOnEnd;
-//   iFileEnd = listOfFiles.end();
-  
-
-//   //Loop on files
-//   uint_4 nRuns=0;
-//   TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
-
-
-//   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
-//     if (para.debuglev_>90){
-//       cout << "load file <" << *iFile << ">" << endl;
-//     }
-
-//     vector<string> tokens;
-//     split(*File,"_",tokens);
-//     string dateOfRun = tokens[1];
-//     string srcLower = tokens[2];
-//     if (para.debuglev_>90){
-//       cout << "date <" << dateOfRun << ">" << endl;
-//     }
-    
-//     PInPersist fin(*iFile);
-//     vector<string> vec = fin.GetNameTags();
-
-//     if (para.typeOfCalib_ == "perRun") {
-//       ///////////////////
-//       //make the calibration of the mean of all Off and On of the run and perform the difference
-
-//       vector<string> modeList;
-//       modeList.push_back("On");
-//       modeList.push_back("Off");
-//       vector<string>::const_iterator iMode;
-//       map<string, TMatrix<r_4> > spectreModeCollect;
-
-//       for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
-//      ///////////////////
-//      //
-//      //Compute the mean of the mode
-//      //
-//      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();
-//      TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
-//      uint_4 nSpectra=0;
-//      //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;
-//      }//end loop Off
-//      //Mean
-//      if(nSpectra>0)meanOfSpectra=(r_4)(nSpectra);
-
-//      //BAO Calibrator
-//      string calibFileName = directoryName + "/" 
-//        + "calib_" + dateOfRun + "_" + srcLower + "_"+(*iMode)+"_" 
-//        + para.calibFreq_ +"MHz-All.txt";
-//      if(debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
-//      ifstream ifs(calibFileName.c_str());
-//      if ( ! ifs.is_open() ) {
-//        rc = 999;
-//        throw calibFileName + " cannot be opened...";
-//      }       
-//      TVector<r_4> calibBAOfactors;
-//      sa_size_t nr,nc; //values read 
-//      calibBAOfactorsOff.ReadASCII(ifs,nr,nc);
-//      if(debuglev_>9){
-//        cout << "(nr,nc): "<< nr << "," << nc << endl;
-//        calibBAOfactors.Print(cout);
-//      }
-        
-//      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
-//        meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);
-//      }
-//      //storage
-//      spectreModeCollect.insert(pair<string, TMatrix<r_4> >(*iMode,TMatrix<r_4>(meanOfSpectra,false))); //do not share data (cf. SOPHYA) 
-//       }//end of mode
-
-//       //Take the difference ON-OFF in current run
-//       TMatrix<r_4>diffOnOff = spectreModeCollect["On"]-spectreModeCollect["Off"];
-
-
-//     } else if (para.typeOfCalib_ == "perCycle") {
-//       //perform the calibration of the OFF and ON per cycle, then make the mean and take the diff
-//     } else {
-//       string msg="FATAL (meanCalibBAODiffOnOffCycles); unknown calibration mode " 
-//      + para.typeOfCalib_ ;
-//       throw(msg);
-//     }  
-
-
-//     nRuns++;
-    
-//   }//eo loop on spectra in a file
-// }
 void meanCalibBAODiffOnOffCycles() throw(string) {
 
@@ -378 +268 @@
   TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
   char* onffTupleName[7]={"cycle",
-                          "onoffRaw0","onoffRun0","onoffCycle0",
-                          "onoffRaw1","onoffRun1","onoffCycle1"};
+                          "onoffRaw0","onoffRaw1",
+                          "onoffRun0","onoffRun1",
+                          "onoffCycle0","onoffCycle1"};
   NTuple onoffevolution(7,onffTupleName);
   r_4 xnt[7];
 
   //Lower and Higher freq. bin to perform mean follow up
-  sa_size_t chLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
+  sa_size_t chLow  = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
   sa_size_t chHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
+  if (para.debuglev_>90){
+    cout << "freq. band for follow up [" <<  chLow << ", "<< chHigh << "]" << endl;
+  }
   
   //Loop on files/run
@@ -501 +395 @@
       double mean,sigma;
       MeanSigma(coef,mean,sigma);
-      cout << "Mean: " << mean << " sigma " << sigma << endl;
       calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
     }
@@ -533 +426 @@
       double mean,sigma;
       MeanSigma(coef,mean,sigma);
-      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
       calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
     }
@@ -566 +459 @@
       double mean,sigma;
       MeanSigma(coef,mean,sigma);
-      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
       calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
     }
@@ -597 +490 @@
       double mean,sigma;
       MeanSigma(coef,mean,sigma);
-      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
       calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
     }
@@ -623 +516 @@
     nr = calibBAOfactors_Off_Run_Ch0.NRows();
     for (sa_size_t ir=0; ir<nr; ++ir){
+//       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
+//         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
+
       calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
         = calibBAOfactors_Off_Run_Ch0(ir,1);
@@ -631 +527 @@
       calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
         = calibBAOfactors_Off_Cycle_Ch1(ir,1);
-    }
+    }//eo loop on coef Off
     
     nr = calibBAOfactors_On_Run_Ch0.NRows();
@@ -643 +539 @@
       calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
         = calibBAOfactors_On_Cycle_Ch1(ir,1);
-    }
+    }//eo loop on coef On
       
-
-    
-//      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
-//        meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);
-//      }
-
     //Loop on cyles
     for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
 
-      if(para.debuglev_>9){
+      //Look if the cycle has been calibrated...
+      bool isCycleCalibrated = 
+        ( calibBAO_On_Run_Ch0.count(*ic)    *
+          calibBAO_On_Run_Ch1.count(*ic)    *
+          calibBAO_Off_Run_Ch0.count(*ic)   *
+          calibBAO_Off_Run_Ch1.count(*ic)   *
+          calibBAO_On_Cycle_Ch0.count(*ic)  *
+          calibBAO_On_Cycle_Ch1.count(*ic)  *
+          calibBAO_Off_Cycle_Ch0.count(*ic) *
+          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
+
+      if(para.debuglev_>9) {
         cout << "Calibration coefficients for cycle "<<*ic << endl; 
-        cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
-             << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
-             << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
-             << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
-             << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
-             << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
-             << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
-             << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
-      }
+        if (isCycleCalibrated) {
+          cout << "Cycle calibrated " << endl;
+          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
+               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
+               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
+               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
+        } else {
+          cout << "Cycle NOT calibrated " << endl;
+        }
+      }//debug
+
+
+      if ( ! isCycleCalibrated ) continue;
       
       string ppftag;
@@ -706 +615 @@
       meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
 
+      //
+      totalNumberCycles++;
       //Fill NTuple
-      xnt[0] = *ic;
+      xnt[0] = totalNumberCycles;
       
       TVector<r_4> meanInRange_noCalib(NUMBER_OF_CHANNELS);
@@ -721 +632 @@
       TVector<r_4> meanInRange_perCycleCalib(NUMBER_OF_CHANNELS);
       meanInRange(diffOnOff_perCycleCalib,chLow,chHigh,meanInRange_perCycleCalib);
-      xnt[3] = meanInRange_perCycleCalib(0);
-      xnt[4] = meanInRange_perCycleCalib(1);
+      xnt[5] = meanInRange_perCycleCalib(0);
+      xnt[6] = meanInRange_perCycleCalib(1);
       
       onoffevolution.Fill(xnt);
 
-      totalNumberCycles++;
+      //Quit if enough
       if (totalNumberCycles >= para.maxNumberCycles_) break;    
 
@@ -735 +646 @@
   cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
   //Normalisation
-  if(totalNumberCycles>0){
+  if(totalNumberCycles > 0){
     meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
     meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
@@ -758 +669 @@
     tmp << totalNumberCycles;
     string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+
     POutPersist fos(fileName);
+    cout << "Save output in " << fileName << endl;
 
     string tag = "meanNoCalib";

BAORadio/AmasNancay/trunk/plotDiffOnOff.pic

@@ -8 +8 @@
 set action $4
 
-set defatt "font=helvetica,bold,20"
+set defatt "font=helvetica,bold,20 fixedfontsize"