Changeset 591 for BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

Timestamp:

Oct 26, 2011, 11:36:56 AM (13 years ago)

Author:

campagne

Message:

continue improve analysis (jec)

File:

• BAORadio/AmasNancay/trunk/mergeAnaFiles.cc (modified) (7 diffs)

BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

@@ -111 +111 @@
 
 //-------------
+//JEC 25/10/11 Perform a median filtering with a sliding window of "halfwidth" half width
+//             It takes care of the edges and is based on the median function (above) 
+void medianFiltering(const TMatrix<r_4> mtx, 
+                     sa_size_t halfwidth,
+                     TMatrix<r_4>& vec) {
+  
+  sa_size_t nr = mtx.NRows();
+  sa_size_t nc = mtx.NCols();
+  sa_size_t chMin = 0;
+  sa_size_t chMax = nc-1;
+  
+  for (sa_size_t ir=0; ir<nr; ir++){
+    for (sa_size_t ic=0; ic<nc; ic++) {
+      sa_size_t chLow = ic-halfwidth;
+      chLow = (chLow >= chMin) ? chLow : chMin;
+      sa_size_t chHigh = ic+halfwidth;
+      chHigh = ( chHigh <= chMax ) ? chHigh : chMax;
+      TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
+      vector<r_4> val;
+      tmp.FillTo(val);
+      vec(ir,ic) = median(val.begin(),val.end());
+    }
+  }
+}
+//-------------
 void split(const string& str, const string& delimiters , vector<string>& tokens) {
     // Skip delimiters at beginning.
@@ -264 +289 @@
 
   //mean ON-OFF over the list of cycles
+  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0 
   TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0 
   TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
   TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
-  static const int NINFO=19;
+  static const int NINFO=21;
   char* onffTupleName[NINFO]={"cycle"
-                          ,"onoffRaw0","onoffRaw1"
-                          ,"onoffRun0","onoffRun1"
-                          ,"onoffCycle0","onoffCycle1"
-                          ,"onoffRaw01420","onoffRaw11420"
-                          ,"onoffRun01420","onoffRun11420"
-                          ,"onoffCycle01420","onoffCycle11420"
-                          ,"onoffRaw01420side","onoffRaw11420side"
-                          ,"onoffRun01420side","onoffRun11420side"
-                          ,"onoffCycle01420side","onoffCycle11420side"
+                              ,"onoffRaw0","onoffRaw1"
+                              ,"onoffRun0","onoffRun1"
+                              ,"onoffCycle0","onoffCycle1"
+                              ,"onoffRaw01420","onoffRaw11420"
+                              ,"onoffRun01420","onoffRun11420"
+                              ,"onoffCycle01420","onoffCycle11420"
+                              ,"onoffRaw01420side","onoffRaw11420side"
+                              ,"onoffRun01420side","onoffRun11420side"
+                              ,"onoffCycle01420side","onoffCycle11420side"
+                              ,"onoffRaw0f14001420","onoffRaw1f14001420"
   };
   NTuple onoffevolution(NINFO,onffTupleName);
@@ -294 +321 @@
   sa_size_t ch1420bLow  = freqToChan(1422);
   sa_size_t ch1420bHigh = freqToChan(1423);
-
+  
+  //25/10/11 follow 1400-1420Mhz
+  sa_size_t ch1420 = freqToChan(1420);
+  sa_size_t ch1400 = freqToChan(1400);
 
   if (para.debuglev_>0){
@@ -631 +661 @@
       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;
+
       TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
       meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
@@ -703 +742 @@
       xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
 
+
+      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
+      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
+      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
+      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
+      
+
       
       //store infos to Ntuple
@@ -717 +764 @@
   //Normalisation
   if(totalNumberCycles > 0){
+    //JEC 29/10 add ON-OFF/OFF
+    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
     meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
     meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
@@ -745 +794 @@
     string tag = "meanNoCalib";
     fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
+    
+    //JEC 29/10/11
+    tag = "meanOvOffNoCalib";
+    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
+
     tag = "meanPerRunCalib";
     fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;