Index: /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc
===================================================================
--- /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc	(revision 615)
+++ /BAORadio/AmasNancay/trunk/mergeAnaFiles.cc	(revision 616)
@@ -63,4 +63,45 @@
 //--------------------------------------------------------------
 //Utility functions
+// template <class T>
+// void MeanSigma(TArray<T> const & a, T & mean, T & sig){
+//   if (a.NbDimensions() < 1) 
+//     throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
+//   const T * pe;
+//   sa_size_t j,k;
+//   mean=0.;
+//   sig = 0.;
+//   T valok;
+//   if (a.AvgStep() > 0)   {  // regularly spaced elements
+//     sa_size_t step = a.AvgStep(); 
+//     sa_size_t maxx = a.Size()*step;
+//     pe = a.Data();
+//     for(k=0; k<maxx; k+=step )  { 
+//       valok = pe[k]; 
+//       mean += valok;  sig += valok*valok;
+//     }
+//   }
+//   else {    // Non regular data spacing ...
+//     int_4 ka = a.MaxSizeKA();
+//     sa_size_t step = a.Step(ka);
+//     sa_size_t gpas = a.Size(ka)*step;
+//     sa_size_t naxa = a.Size()/a.Size(ka);
+//     for(j=0; j<naxa; j++)  {
+//       pe = a.DataBlock().Begin()+a.Offset(ka,j);
+//       for(k=0; k<gpas; k+=step) { 
+// 	valok = pe[k]; 
+// 	mean += valok;  sig += valok*valok;
+//       }
+//     }
+//   }
+//   T dsz = (T)(a.Size());
+//   mean /= dsz;
+//   if (dsz > 1.5) {
+//     sig = sig/dsz - mean*mean;
+//     sig *= (dsz/(dsz-1));
+//     if (sig >= 0.) sig = sqrt(sig);
+//   }
+//   else sig = 0.;
+// }
+
 
 sa_size_t freqToChan(r_4 f){
@@ -79,5 +120,5 @@
   for (sa_size_t ir=0; ir<nr; ir++){
     TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
-    double mean, sigma;
+    r_8 mean, sigma;
     MeanSigma(tmp,mean,sigma);
     vec(ir) = mean;
@@ -273,4 +314,5 @@
   }//eo loop on slices
 }
+
 //AST 18/11/11 
 //-------------------------------------------------------
@@ -449,5 +491,4 @@
     }
     ifs.close();
-
     //
     calibFileName = directoryName + "/" 
@@ -621,5 +662,4 @@
     }
     ifs.close();
-
     //link <cycle> - <calibration coefficient>
     //We cannot rely on identical cycle list of the OFF and ON calibration
@@ -770,4 +810,5 @@
 
 }
+
 //-------------------------------------------------------
 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra 
@@ -812,19 +853,7 @@
   sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
   //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
-  sa_size_t ch1420Low;
-  sa_size_t ch1420High;
-  if (para.sourceName_ == "Abell85") {
-    ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
-    ch1420High = freqToChan(1420.4+0.2);
-  } else if (para.sourceName_ == "Abell1205") {
-    ch1420Low  = freqToChan(1420.4-0.3);  // Abell1205
-    ch1420High = freqToChan(1420.4+0.2);
-  } else if (para.sourceName_ == "Abell2440") {
-    ch1420Low  = freqToChan(1420.4);      // Abell2440
-    ch1420High = freqToChan(1420.4+0.3);
-  } else {
-    ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
-    ch1420High = freqToChan(1420.4+0.2);
-  }
+  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
+  sa_size_t ch1420High = freqToChan(1420.4+0.2);
+
   //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
   sa_size_t ch1420aLow  = freqToChan(1418);
@@ -865,6 +894,4 @@
     string srcLower = tokens2[0];
 
-
-    
     PInPersist fin(*iFile);
     vector<string> vec = fin.GetNameTags();
@@ -1082,6 +1109,6 @@
     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;
+      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)]
@@ -1134,5 +1161,5 @@
 	       << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
 	} else {
-	  cout << "Cycle NOT calibrated " << endl;
+	  cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
 	}
       }//debug
@@ -1178,9 +1205,9 @@
       //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);
+      TMatrix<r_4> aSpecOffFiltered(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);
+      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
       
-      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
+      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
       meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
 
@@ -1191,6 +1218,6 @@
       meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
 
-      //
       totalNumberCycles++;
+
       //Fill NTuple
       xnt[0] = totalNumberCycles;
@@ -1266,7 +1293,8 @@
       //JEC 18/11/11 follow up the 1400-1420MHz OFF only
       TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
-      aSpecOffovOff.Div(aSpecOffFitltered);
+      aSpecOffovOff.Div(aSpecOffFiltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
+      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
       meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
 
@@ -1275,7 +1303,8 @@
 
       TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
-      aSpecOnovOff.Div(aSpecOffFitltered);
+      aSpecOnovOff.Div(aSpecOffFiltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
+      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
       meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
 
@@ -1352,6 +1381,9 @@
     
     tag = "onoffevol";
-    fos << PPFNameTag(tag) << onoffevolution;    
+    fos << PPFNameTag(tag) << onoffevolution;
+ 
   }//end of save
+
+
 }
 //JEC 14/11/11 New meanRawDiffOnOffCycles START
@@ -1376,5 +1408,4 @@
   TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
   TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
-
   //Tuple only for RAW things to follow 
   static const int NINFO=11;
@@ -1417,5 +1448,17 @@
   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
     if (para.debuglev_>90){
-    }
+      cout << "load file <" << *iFile << ">" << endl;
+    }
+
+    vector<string> tokens;
+    split(*iFile,"_",tokens);
+    string dateOfRun = tokens[1];
+    if (para.debuglev_>90){
+      cout << "date <" << dateOfRun << ">" << endl;
+    }
+    vector<string> tokens2;
+    split(tokens[2],".",tokens2);
+    string srcLower = tokens2[0];
+
     PInPersist fin(*iFile);
     vector<string> vec = fin.GetNameTags();
@@ -1447,5 +1490,5 @@
 	regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
 	if (para.debuglev_>90){
-	  cout << " spactra <" << *iSpec << ">" << endl;
+	  cout << " spectra <" << *iSpec << ">" << endl;
 	  cout << " cycle " << iSpec->substr(b) << endl;
 	}
@@ -1475,4 +1518,22 @@
     for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
       
+      // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440
+      if ( *ic == 1 && srcLower == "abell1205" ) {
+	  if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {
+	    cout << "Skipping non-calibrated cycle " << *ic << endl;
+	    continue;
+	  }
+      } else if ( *ic == 1 && srcLower == "abell2440" ) {
+	  if ( dateOfRun == "20110516" ) {
+	    cout << "Skipping non-calibrated cycle " << *ic << endl;
+	    continue;
+	  }
+      } else if ( *ic == 3 && srcLower == "abell1205" ) {
+	  if ( dateOfRun == "20110810" ) {
+	    cout << "Skipping non-calibrated cycle " << *ic << endl;
+	    continue;
+	  }
+      }
+
       string ppftag;
       //load ON phase
@@ -1490,32 +1551,30 @@
       //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);
+      TMatrix<r_4> aSpecOffFiltered(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);
+      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
       
-      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
+      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
       meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
-
 
       //JEC 15/11/11 add ON/OFF and OFF/OFF
       TMatrix<r_4> onOvOff(aSpecOn,false);
-      onOvOff.Div(aSpecOffFitltered);
+      onOvOff.Div(aSpecOffFiltered);
       meanONovOFF_noCalib += onOvOff;
       
       TMatrix<r_4> offOvOff(aSpecOff,false);
-      offOvOff.Div(aSpecOffFitltered);
+      offOvOff.Div(aSpecOffFiltered);
       meanOFFovOFF_noCalib += offOvOff;
 
       totalNumberCycles++;
-
 
       //Fill NTuple
       xnt[0] = totalNumberCycles;
  
-
       //Follow up arround the 1420.4MHz Freq.
       TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
@@ -1535,13 +1594,13 @@
 
 
-      //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(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
-      xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
-      xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
-
-      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
+      //JEC 25/10/11 follow 1410-1415 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
+      TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib);
+      xnt[5] = meanInRange_1410a1415Freq_noCalib(0);
+      xnt[6] = meanInRange_1410a1415Freq_noCalib(1);
+
+      //JEC 18/11/11 follow up the 1410-1415MHz OFF only
       TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
-      aSpecOffovOff.Div(aSpecOffFitltered);
+      aSpecOffovOff.Div(aSpecOffFiltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
@@ -1552,5 +1611,5 @@
 
       TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
-      aSpecOnovOff.Div(aSpecOffFitltered);
+      aSpecOnovOff.Div(aSpecOffFiltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
@@ -1600,7 +1659,9 @@
 
     tag = "onoffevol";
-    fos << PPFNameTag(tag) << onoffevolution;    
+    fos << PPFNameTag(tag) << onoffevolution;
 
   }//end save
+
+
 }
 //JEC 14/11/11 New meanRawDiffOnOffCycles END
@@ -1962,6 +2023,4 @@
     } else if (action == "meanCalibBAODiffOnOff") {
       meanCalibBAODiffOnOffCycles();
-    } else if (action == "calibCoeffNtp") {
-      calibCoeffNtp();
     } else {
       msg = "Unknown action " + action;
