source: BAORadio/AmasNancay/trunk/mergeAnaFiles.cc @ 610

Last change on this file since 610 was 610, checked in by campagne, 13 years ago

cleaning
File size: 53.3 KB
Line 
1// Utilisation de SOPHYA pour faciliter les tests ...
2#include "sopnamsp.h"
3#include "machdefs.h"
4#include <dirent.h>
5#include <matharr.h>
6
7// include standard c/c++
8#include <regex.h>
9#include <stdio.h>
10#include <stdlib.h>
11#include <limits>
12#include <iostream>
13#include <fstream>
14#include <string>
15#include <vector>
16#include <map>
17#include <functional>
18#include <algorithm>
19#include <numeric>
20#include <list>
21#include <exception>
22
23// include sophya mesure ressource CPU/memoire ...
24#include "resusage.h"
25#include "ctimer.h"
26#include "timing.h"
27#include "timestamp.h"
28#include "strutilxx.h"
29#include "ntuple.h"
30#include "fioarr.h"
31#include "tarrinit.h"
32#include "histinit.h"
33#include "fitsioserver.h"
34#include "fiosinit.h"
35#include "ppersist.h"
36
37//-----------------------------------------------
38//Usage
39//
40//./Objs/mergeAnaFiles -act meanRawDiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src NGC4383 -ppfFile dataRaw -debug 1000 -mxcycles 10
41//./Objs/mergeAnaFiles -src Abell85 -act meanCalibBAODiffOnOff -mxcycles 500 -calibfreq 1410 -inPath /sps/baoradio/AmasNancay/JEC/ -ppfFile dataRaw -debug 1 >& mergeAna-500.log
42//
43//
44//-----------------------------------------------
45const sa_size_t NUMBER_OF_CHANNELS = 2;
46const sa_size_t NUMBER_OF_FREQ     = 8192;
47const r_4    LOWER_FREQUENCY       = 1250.0; //MHz
48const r_4    TOTAL_BANDWIDTH       = 250.0;  //MHz
49//Input parameters
50struct Param {
51  int debuglev_;      //debug
52  string inPath_;     //root directory of the input files
53  string outPath_;    //output files are located here
54  string sourceName_; //source name & subdirectory of the input files
55  string ppfFile_;    //generic name of the input files
56  int nSliceInFreq_;  //used by reduceSpectra() fnc
57  string calibFreq_;  //freq. value used for calibration
58  r_4 rcalibFreq_;    //float version
59  string calibBandFreq_;  //band of freq. used for calibration
60  r_4 rcalibBandFreq_;    //float version
61  int maxNumberCycles_;//maximum number of cycles to be treated
62} para;
63//--------------------------------------------------------------
64//Utility functions
65// template <class T>
66// void MeanSigma(TArray<T> const & a, T & mean, T & sig){
67//   if (a.NbDimensions() < 1)
68//     throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
69//   const T * pe;
70//   sa_size_t j,k;
71//   mean=0.;
72//   sig = 0.;
73//   T valok;
74//   if (a.AvgStep() > 0)   {  // regularly spaced elements
75//     sa_size_t step = a.AvgStep();
76//     sa_size_t maxx = a.Size()*step;
77//     pe = a.Data();
78//     for(k=0; k<maxx; k+=step )  {
79//       valok = pe[k];
80//       mean += valok;  sig += valok*valok;
81//     }
82//   }
83//   else {    // Non regular data spacing ...
84//     int_4 ka = a.MaxSizeKA();
85//     sa_size_t step = a.Step(ka);
86//     sa_size_t gpas = a.Size(ka)*step;
87//     sa_size_t naxa = a.Size()/a.Size(ka);
88//     for(j=0; j<naxa; j++)  {
89//       pe = a.DataBlock().Begin()+a.Offset(ka,j);
90//       for(k=0; k<gpas; k+=step) {
91//      valok = pe[k];
92//      mean += valok;  sig += valok*valok;
93//       }
94//     }
95//   }
96//   T dsz = (T)(a.Size());
97//   mean /= dsz;
98//   if (dsz > 1.5) {
99//     sig = sig/dsz - mean*mean;
100//     sig *= (dsz/(dsz-1));
101//     if (sig >= 0.) sig = sqrt(sig);
102//   }
103//   else sig = 0.;
104// }
105
106
107sa_size_t freqToChan(r_4 f){
108  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
109}
110//--------
111//COMPUTE the mean value on a freq. range for all channels
112//--------
113void meanInRange(const TMatrix<r_4> mtx, 
114                 sa_size_t chLow,
115                 sa_size_t chHigh, 
116                 TVector<r_4>& vec){
117 
118  sa_size_t nr = mtx.NRows();
119 
120  for (sa_size_t ir=0; ir<nr; ir++){
121    TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
122    r_8 mean, sigma;
123    MeanSigma(tmp,mean,sigma);
124    vec(ir) = mean;
125  }
126}
127//---------
128class median_of_empty_list_exception:public std::exception{
129    virtual const char* what() const throw() {
130    return "Attempt to take the median of an empty list of numbers.  "
131      "The median of an empty list is undefined.";
132  }
133};
134template<class RandAccessIter>
135double median(RandAccessIter begin, RandAccessIter end) 
136  throw(median_of_empty_list_exception){
137  if(begin == end){ throw median_of_empty_list_exception(); }
138  std::size_t size = end - begin;
139  std::size_t middleIdx = size/2;
140  RandAccessIter target = begin + middleIdx;
141  std::nth_element(begin, target, end);
142   
143  if(size % 2 != 0){ //Odd number of elements
144    return *target;
145  }else{            //Even number of elements
146    double a = *target;
147    RandAccessIter targetNeighbor= target-1;
148    std::nth_element(begin, targetNeighbor, end);
149    return (a+*targetNeighbor)/2.0;
150  }
151}
152
153//-------------
154//JEC 25/10/11 Perform a median filtering with a sliding window of "halfwidth" half width
155//             It takes care of the edges and is based on the median function (above)
156void medianFiltering(const TMatrix<r_4> mtx, 
157                     sa_size_t halfwidth,
158                     TMatrix<r_4>& vec) {
159 
160  sa_size_t nr = mtx.NRows();
161  sa_size_t nc = mtx.NCols();
162  sa_size_t chMin = 0;
163  sa_size_t chMax = nc-1;
164 
165  for (sa_size_t ir=0; ir<nr; ir++){
166    for (sa_size_t ic=0; ic<nc; ic++) {
167      sa_size_t chLow = ic-halfwidth;
168      chLow = (chLow >= chMin) ? chLow : chMin;
169      sa_size_t chHigh = ic+halfwidth;
170      chHigh = ( chHigh <= chMax ) ? chHigh : chMax;
171      TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
172      vector<r_4> val;
173      tmp.FillTo(val);
174      vec(ir,ic) = median(val.begin(),val.end());
175    }
176  }
177}
178//-------------
179void split(const string& str, const string& delimiters , vector<string>& tokens) {
180    // Skip delimiters at beginning.
181    string::size_type lastPos = str.find_first_not_of(delimiters, 0);
182    // Find first "non-delimiter".
183    string::size_type pos     = str.find_first_of(delimiters, lastPos);
184
185    while (string::npos != pos || string::npos != lastPos)
186    {
187        // Found a token, add it to the vector.
188        tokens.push_back(str.substr(lastPos, pos - lastPos));
189        // Skip delimiters.  Note the "not_of"
190        lastPos = str.find_first_not_of(delimiters, pos);
191        // Find next "non-delimiter"
192        pos = str.find_first_of(delimiters, lastPos);
193    }
194}
195//--------------------------------------------------------------
196char *regexp (const char *string, const char *patrn, int *begin, int *end) {   
197        int i, w=0, len;                 
198        char *word = NULL;
199        regex_t rgT;
200        regmatch_t match;
201        regcomp(&rgT,patrn,REG_EXTENDED);
202        if ((regexec(&rgT,string,1,&match,0)) == 0) {
203                *begin = (int)match.rm_so;
204                *end = (int)match.rm_eo;
205                len = *end-*begin;
206                word=(char*)malloc(len+1);
207                for (i=*begin; i<*end; i++) {
208                        word[w] = string[i];
209                        w++; }
210                word[w]=0;
211        }
212        regfree(&rgT);
213        return word;
214}
215//-------
216sa_size_t round_sa(r_4 r) {
217  return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
218}
219//-----
220string StringToLower(string strToConvert){
221  //change each element of the string to lower case
222  for(unsigned int i=0;i<strToConvert.length();i++) {
223    strToConvert[i] = tolower(strToConvert[i]);
224  }
225  return strToConvert;//return the converted string
226}
227//-----
228bool stringCompare( const string &left, const string &right ){
229   if( left.size() < right.size() )
230      return true;
231   for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
232      if( tolower( *lit ) < tolower( *rit ) )
233         return true;
234      else if( tolower( *lit ) > tolower( *rit ) )
235         return false;
236   return false;
237}
238//-----
239list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
240  list<string> theList;
241
242
243  DIR *dip;
244  struct dirent *dit;
245  string msg;  string fileName;
246  string fullFileName;
247  size_t found;
248
249  if ((dip=opendir(dir.c_str())) == NULL ) {
250    msg = "opendir failed on directory "+dir;
251    throw msg;
252  }
253  while ( (dit = readdir(dip)) != NULL ) {
254    fileName = dit->d_name;
255    found=fileName.find(filePettern);
256    if (found != string::npos) {
257      fullFileName = dir + "/";
258      fullFileName += fileName;
259      theList.push_back(fullFileName);
260    }
261  }//eo while
262  if (closedir(dip) == -1) {
263    msg = "closedir failed on directory "+dir;
264    throw msg;
265  }
266 
267  theList.sort(stringCompare);
268
269  return theList;
270
271}
272//
273class  StringMatch : public unary_function<string,bool> {
274public:
275  StringMatch(const string& pattern): pattern_(pattern){}
276  bool operator()(const string& aStr) const {
277
278
279    int b,e;
280    regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
281
282//     cout << "investigate " << aStr << " to find " << pattern_
283//       << "[" <<b<<","<<e<<"]"
284//       << endl;
285
286   
287    if (b != 0) return false;
288    if (e != (int)aStr.size()) return false;
289    return true;
290
291  }
292private:
293  string pattern_;
294};
295//-------------------------------------------------------
296//Rebin in frequence + compute mean and sigma
297void reduceSpectra(const TMatrix<r_4>& specMtxInPut, 
298                    TMatrix<r_4>& meanMtx, 
299                    TMatrix<r_4>& sigmaMtx) {
300  sa_size_t nSliceFreq = para.nSliceInFreq_;
301  sa_size_t deltaFreq =  NUMBER_OF_FREQ/nSliceFreq;
302  for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
303    sa_size_t freqLow= iSlice*deltaFreq;
304    sa_size_t freqHigh= freqLow + deltaFreq -1;
305    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
306      TVector<r_4> reducedRow;
307      reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
308      double mean; 
309      double sigma;
310      MeanSigma(reducedRow,mean,sigma);
311      meanMtx(iCh,iSlice) = mean;
312      sigmaMtx(iCh,iSlice) = sigma;
313    }//eo loop on channels
314  }//eo loop on slices
315}
316//-------------------------------------------------------
317//Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra
318//
319void meanCalibBAODiffOnOffCycles() throw(string) {
320
321  list<string> listOfFiles;
322  string directoryName;
323  directoryName = para.inPath_ + "/" + para.sourceName_;
324
325  //Make the listing of the directory
326  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
327 
328  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
329  iFileEnd = listOfFiles.end();
330
331  //mean ON-OFF over the list of cycles
332  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
333  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
334  TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
335  TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
336  static const int NINFO=25;
337  char* onffTupleName[NINFO]={"cycle"
338                              ,"onoffRaw0","onoffRaw1"
339                              ,"onoffRun0","onoffRun1"
340                              ,"onoffCycle0","onoffCycle1"
341                              ,"onoffRaw01420","onoffRaw11420"
342                              ,"onoffRun01420","onoffRun11420"
343                              ,"onoffCycle01420","onoffCycle11420"
344                              ,"onoffRaw01420side","onoffRaw11420side"
345                              ,"onoffRun01420side","onoffRun11420side"
346                              ,"onoffCycle01420side","onoffCycle11420side"
347                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
348                              ,"offRaw0f14101415","offRaw1f14101415"
349                              ,"onRaw0f14101415","onRaw1f14101415"
350  };
351  NTuple onoffevolution(NINFO,onffTupleName);
352  r_4 xnt[NINFO];
353
354  //Lower and Higher freq. arround the Calibration Freq. bin to perform mean follow up
355  sa_size_t chCalibLow  = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
356  sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
357  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
358  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
359  sa_size_t ch1420High = freqToChan(1420.4+0.2);
360
361  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
362  sa_size_t ch1420aLow  = freqToChan(1418);
363  sa_size_t ch1420aHigh = freqToChan(1419);
364  sa_size_t ch1420bLow  = freqToChan(1422);
365  sa_size_t ch1420bHigh = freqToChan(1423);
366 
367  //1400-1420Mhz
368  sa_size_t ch1400 = freqToChan(1400);
369  //  sa_size_t ch1405 = freqToChan(1400);
370  sa_size_t ch1410 = freqToChan(1410);
371  sa_size_t ch1415 = freqToChan(1415);
372  sa_size_t ch1420 = freqToChan(1420);
373
374  if (para.debuglev_>0){
375    cout << "freq. band for follow up [" <<  chCalibLow << ", "<< chCalibHigh << "]" << endl;
376    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
377    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
378    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
379  }
380 
381  //Loop on files/run
382
383  int totalNumberCycles=0; //total number of cycles for normalisation
384  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
385    if (para.debuglev_>90){
386      cout << "load file <" << *iFile << ">" << endl;
387    }
388
389    vector<string> tokens;
390    split(*iFile,"_",tokens);
391    string dateOfRun = tokens[1];
392    if (para.debuglev_>90){
393      cout << "date <" << dateOfRun << ">" << endl;
394    }
395    vector<string> tokens2;
396    split(tokens[2],".",tokens2);
397    string srcLower = tokens2[0];
398
399
400   
401    PInPersist fin(*iFile);
402    vector<string> vec = fin.GetNameTags();
403
404    vector<string> modeList;
405    modeList.push_back("On");
406    modeList.push_back("Off");
407    vector<string>::const_iterator iMode;
408   
409    map<string, list<int> > cycleModeCollect;
410   
411    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
412      list<string> listOfSpectra;
413      //Keep only required PPF objects
414      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
415      std::remove_copy_if(
416                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
417                          not1(StringMatch(matchstr))
418                          );
419     
420      listOfSpectra.sort(stringCompare);
421      iSpecEnd = listOfSpectra.end();
422     
423      matchstr = "[0-9]+";
424      //Loop of spectra matrix
425      list<int> listOfCycles;
426      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
427        int b,e;
428        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
429        if (para.debuglev_>90){
430          cout << " spactra <" << *iSpec << ">" << endl;
431          cout << " cycle " << iSpec->substr(b) << endl;
432        }
433        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
434      }//end loop spectra
435      cycleModeCollect[*iMode] = listOfCycles;
436    }//end of mode   
437
438    //Take the Intersection of the list Of cycles in mode Off and On   
439    list<int> commonCycles;
440    set_intersection(cycleModeCollect["On"].begin(),
441                     cycleModeCollect["On"].end(),
442                     cycleModeCollect["Off"].begin(),
443                     cycleModeCollect["Off"].end(),
444                     back_inserter(commonCycles)
445                     );
446   
447    if (para.debuglev_>90){
448      cout << "Liste of cycles common to On & Off: <";
449      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
450        cout << *i << " ";
451      }
452      cout << ">" << endl;
453    }
454   
455    //
456    //Load BAO Calibration factors "per Cycle and Channels"
457    //Compute the mean per Cycle to
458    //  fill factors "per Run and Channels" with the same cycle list
459    //
460    //
461    //TODO improve the code....
462
463    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
464    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
465    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
466    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
467   
468    string calibFileName;
469    ifstream ifs;
470    sa_size_t nr,nc; //values read
471
472    //OFF Cycle per Channel
473    calibFileName = directoryName + "/" 
474      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
475      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
476    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
477    ifs.open(calibFileName.c_str());
478    if ( ! ifs.is_open() ) {
479
480      throw calibFileName + " cannot be opened...";
481    }   
482    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
483    if(para.debuglev_>9){
484      cout << "(nr,nc): "<< nr << "," << nc << endl;
485      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
486      cout << endl;
487    }
488
489    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
490    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
491    {//Compute the mean
492      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
493      double mean,sigma;
494      MeanSigma(coef,mean,sigma);
495      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
496    }
497    if(para.debuglev_>9){
498      cout << "Fill calib. with mean value " << endl; 
499      calibBAOfactors_Off_Run_Ch0.Print(cout);
500      cout << endl;
501    }
502    ifs.close();
503
504    //
505    calibFileName = directoryName + "/" 
506      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
507      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
508    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
509    ifs.open(calibFileName.c_str());
510    if ( ! ifs.is_open() ) {
511
512      throw calibFileName + " cannot be opened...";
513    }   
514    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
515    if(para.debuglev_>9){
516      cout << "(nr,nc): "<< nr << "," << nc << endl;
517      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
518      cout << endl;
519    }
520    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
521    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
522    {//Compute the mean
523      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
524      double mean,sigma;
525      MeanSigma(coef,mean,sigma);
526      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
527      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
528    }
529    if(para.debuglev_>9){
530      cout << "Fill calib. with mean value " << endl; 
531      calibBAOfactors_Off_Run_Ch1.Print(cout);
532      cout << endl;
533    }
534    ifs.close();
535
536    //ON Cycle per Channel
537    calibFileName = directoryName + "/" 
538      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
539      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
540    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
541    ifs.open(calibFileName.c_str());
542    if ( ! ifs.is_open() ) {
543
544      throw calibFileName + " cannot be opened...";
545    }   
546    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
547    if(para.debuglev_>9){
548      cout << "(nr,nc): "<< nr << "," << nc << endl;
549      calibBAOfactors_On_Cycle_Ch0.Print(cout);
550      cout << endl;     
551    }
552
553    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
554    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
555    {//Compute the mean
556      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
557      double mean,sigma;
558      MeanSigma(coef,mean,sigma);
559      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
560      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
561    }
562    if(para.debuglev_>9){
563      cout << "Fill calib. with mean value " << endl; 
564      calibBAOfactors_On_Run_Ch0.Print(cout);
565      cout << endl;
566    }
567    ifs.close();
568
569   
570    calibFileName = directoryName + "/" 
571      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
572      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
573    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
574    ifs.open(calibFileName.c_str());
575    if ( ! ifs.is_open() ) {
576      throw calibFileName + " cannot be opened...";
577    }   
578    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
579    if(para.debuglev_>9){
580      cout << "(nr,nc): "<< nr << "," << nc << endl;
581      calibBAOfactors_On_Cycle_Ch1.Print(cout);
582      cout << endl;
583    }
584    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
585    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
586    {//Compute the mean
587      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
588      double mean,sigma;
589      MeanSigma(coef,mean,sigma);
590      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
591      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
592    }
593    if(para.debuglev_>9){
594      cout << "Fill calib. with mean value " << endl; 
595      calibBAOfactors_On_Run_Ch1.Print(cout);
596      cout << endl;
597    }
598
599    ifs.close();
600   
601    //link <cycle> - <calibration coefficient>
602    //We cannot rely on identical cycle list of the OFF and ON calibration
603    map<int,r_4> calibBAO_Off_Run_Ch0;
604    map<int,r_4> calibBAO_Off_Run_Ch1;
605    map<int,r_4> calibBAO_On_Run_Ch0;
606    map<int,r_4> calibBAO_On_Run_Ch1;
607
608    map<int,r_4> calibBAO_Off_Cycle_Ch0;
609    map<int,r_4> calibBAO_Off_Cycle_Ch1;
610    map<int,r_4> calibBAO_On_Cycle_Ch0;
611    map<int,r_4> calibBAO_On_Cycle_Ch1;
612
613    //per Run based BAO coefficients
614    nr = calibBAOfactors_Off_Run_Ch0.NRows();
615    for (sa_size_t ir=0; ir<nr; ++ir){
616      cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
617           << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
618
619      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
620        = calibBAOfactors_Off_Run_Ch0(ir,1);
621      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
622        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
623      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
624        = calibBAOfactors_Off_Run_Ch1(ir,1);
625      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
626        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
627    }//eo loop on coef Off
628   
629    nr = calibBAOfactors_On_Run_Ch0.NRows();
630    for (sa_size_t ir=0; ir<nr; ++ir){
631      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
632        = calibBAOfactors_On_Run_Ch0(ir,1);
633      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
634        = calibBAOfactors_On_Cycle_Ch0(ir,1);
635      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
636        = calibBAOfactors_On_Run_Ch1(ir,1);
637      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
638        = calibBAOfactors_On_Cycle_Ch1(ir,1);
639    }//eo loop on coef On
640     
641    //Loop on cyles
642    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
643
644      //Look if the cycle has been calibrated...
645      bool isCycleCalibrated = 
646        ( calibBAO_On_Run_Ch0.count(*ic)    *
647          calibBAO_On_Run_Ch1.count(*ic)    *
648          calibBAO_Off_Run_Ch0.count(*ic)   *
649          calibBAO_Off_Run_Ch1.count(*ic)   *
650          calibBAO_On_Cycle_Ch0.count(*ic)  *
651          calibBAO_On_Cycle_Ch1.count(*ic)  *
652          calibBAO_Off_Cycle_Ch0.count(*ic) *
653          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
654
655      if(para.debuglev_>9) {
656        cout << "Calibration coefficients for cycle "<<*ic << endl; 
657        if (isCycleCalibrated) {
658          cout << "Cycle calibrated " << endl;
659          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
660               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
661               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
662               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
663               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
664               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
665               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
666               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
667        } else {
668          cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
669        }
670      }//debug
671
672
673      if ( ! isCycleCalibrated ) continue;
674     
675      string ppftag;
676      //load ON phase
677      stringstream cycle;
678      cycle << *ic;
679     
680      ppftag = "specRawOn"+cycle.str();
681      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
682      fin.GetObject(aSpecOn,ppftag);
683
684      TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
685      aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
686      aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
687
688      TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
689      aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
690      aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
691     
692      //Load OFF phase
693      ppftag = "specRawOff"+cycle.str();
694      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
695      fin.GetObject(aSpecOff,ppftag);
696
697      TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
698      aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
699      aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
700
701      TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
702      aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
703      aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
704
705
706      //Perform the difference ON-OFF with the different calibration options
707      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
708      meanDiffONOFF_noCalib += diffOnOff_noCalib;
709     
710      //JEC 29/10/11 add ON-OFF/OFF
711      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
712      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
713      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
714      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
715     
716      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
717      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
718
719      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
720      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
721
722      TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
723      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
724
725      //
726      totalNumberCycles++;
727      //Fill NTuple
728      xnt[0] = totalNumberCycles;
729 
730      //Follow up arround the Calibration Freq.
731      TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
732      meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
733      xnt[1] = meanInRange_CalibFreq_noCalib(0);
734      xnt[2] = meanInRange_CalibFreq_noCalib(1);
735
736      TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
737      meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
738      xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
739      xnt[4] = meanInRange_CalibFreq_perRunCalib(1);
740
741      TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
742      meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
743      xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
744      xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);
745
746
747      //Follow up arround the 1420.4MHz Freq.
748      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
749      meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
750      xnt[7] = meanInRange_1420Freq_noCalib(0);
751      xnt[8] = meanInRange_1420Freq_noCalib(1);
752
753      TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
754      meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
755      xnt[9] = meanInRange_1420Freq_perRunCalib(0);
756      xnt[10] = meanInRange_1420Freq_perRunCalib(1);
757
758      TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
759      meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
760      xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
761      xnt[12] = meanInRange_1420Freq_perCycleCalib(1);
762
763
764      //Follow up below the 1420.4MHz Freq.
765      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
766      meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
767      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
768      meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
769
770      xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
771      xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
772
773
774      TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
775      meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
776      TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
777      meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);
778
779      xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
780      xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;
781
782
783      TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
784      meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
785      TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
786      meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);
787
788      xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
789      xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
790
791
792      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
793      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
794      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
795      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
796      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
797
798      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
799      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
800      aSpecOffovOff.Div(aSpecOffFitltered);
801
802      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
803      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
804      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
805
806      xnt[21] = meanInRange_1410a1415Freq_OFF_noCalib(0);
807      xnt[22] = meanInRange_1410a1415Freq_OFF_noCalib(1);
808
809      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
810      aSpecOnovOff.Div(aSpecOffFitltered);
811
812      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
813      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
814      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
815
816      xnt[23] = meanInRange_1410a1415Freq_ON_noCalib(0);
817      xnt[24] = meanInRange_1410a1415Freq_ON_noCalib(1);
818     
819      //store infos to Ntuple
820      onoffevolution.Fill(xnt);
821
822      //Quit if enough
823      if (totalNumberCycles >= para.maxNumberCycles_) break;   
824
825    }//eo loop on cycles
826    if (totalNumberCycles >= para.maxNumberCycles_) break;         
827 
828  }//eo loop on spectra in a file
829  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
830  //Normalisation
831  if(totalNumberCycles > 0){
832    //JEC 29/10 add ON-OFF/OFF
833    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
834    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
835    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
836    meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
837  } 
838 
839  //Compute the reduced version of the mean and sigma
840  TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
841  TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
842  reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
843
844  TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
845  TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
846  reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
847
848  TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
849  TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
850  reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
851
852  {//save the results
853    stringstream tmp;
854    tmp << totalNumberCycles;
855    string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
856
857    POutPersist fos(fileName);
858    cout << "Save output in " << fileName << endl;
859
860    string tag = "meanNoCalib";
861    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
862   
863    //JEC 29/10/11
864    tag = "meanOvOffNoCalib";
865    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
866
867    tag = "meanPerRunCalib";
868    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
869    tag = "meanPerCycleCalib";
870    fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
871
872    tag = "redmeanNoCalib";
873    fos << PPFNameTag(tag) << meanRedMtx_noCalib;
874    tag = "redsigmaNoCalib";
875    fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
876
877    tag = "redmeanPerRunCalib";
878    fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
879    tag = "redsigmaPerRunCalib";
880    fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
881
882    tag = "redmeanPerCycleCalib";
883    fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
884    tag = "redsigmaPerCycleCalib";
885    fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
886   
887    tag = "onoffevol";
888    fos << PPFNameTag(tag) << onoffevolution;   
889  }//end of save
890
891
892}
893//JEC 14/11/11 New meanRawDiffOnOffCycles START
894//-------------------------------------------------------
895//Compute the mean of Diff ON-OFF/OFF from Raw spectra
896//Used like:
897//
898void meanRawDiffOnOffCycles() throw(string) {
899  list<string> listOfFiles;
900  string directoryName;
901  directoryName = para.inPath_ + "/" + para.sourceName_;
902
903  //Make the listing of the directory
904  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
905 
906  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
907  iFileEnd = listOfFiles.end();
908
909  //mean ON-OFF over the list of cycles
910  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
911  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
912  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
913  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
914
915  //Tuple only for RAW things to follow
916  static const int NINFO=11;
917  char* onffTupleName[NINFO]={"cycle"
918                              ,"onoffRaw01420","onoffRaw11420"
919                              ,"onoffRaw01420side","onoffRaw11420side"
920                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
921                              ,"offRaw0f14101415","offRaw1f14101415"
922                              ,"onRaw0f14101415","onRaw1f14101415"
923  };
924  NTuple onoffevolution(NINFO,onffTupleName);
925  r_4 xnt[NINFO];
926
927  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
928  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
929  sa_size_t ch1420High = freqToChan(1420.4+0.2);
930
931  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
932  sa_size_t ch1420aLow  = freqToChan(1418);
933  sa_size_t ch1420aHigh = freqToChan(1419);
934  sa_size_t ch1420bLow  = freqToChan(1422);
935  sa_size_t ch1420bHigh = freqToChan(1423);
936 
937  //1400-1420Mhz
938  sa_size_t ch1400 = freqToChan(1400);
939  //  sa_size_t ch1405 = freqToChan(1400);
940  sa_size_t ch1410 = freqToChan(1410);
941  sa_size_t ch1415 = freqToChan(1415);
942  sa_size_t ch1420 = freqToChan(1420);
943
944  if (para.debuglev_>0){
945    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
946    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
947    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
948  }
949
950  int totalNumberCycles=0; //total number of cycles
951
952  //Loop on files
953  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
954    if (para.debuglev_>90){
955    }
956    PInPersist fin(*iFile);
957    vector<string> vec = fin.GetNameTags();
958
959    vector<string> modeList;
960    modeList.push_back("On");
961    modeList.push_back("Off");
962    vector<string>::const_iterator iMode;
963   
964    map<string, list<int> > cycleModeCollect;
965   
966    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
967      list<string> listOfSpectra;
968      //Keep only required PPF objects
969      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
970      std::remove_copy_if(
971                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
972                          not1(StringMatch(matchstr))
973                          );
974     
975      listOfSpectra.sort(stringCompare);
976      iSpecEnd = listOfSpectra.end();
977     
978      matchstr = "[0-9]+";
979      //Loop of spectra matrix
980      list<int> listOfCycles;
981      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
982        int b,e;
983        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
984        if (para.debuglev_>90){
985          cout << " spactra <" << *iSpec << ">" << endl;
986          cout << " cycle " << iSpec->substr(b) << endl;
987        }
988        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
989      }//end loop spectra
990      cycleModeCollect[*iMode] = listOfCycles;
991    }//end of mode   
992
993    //Take the Intersection of the list Of cycles in mode Off and On   
994    list<int> commonCycles;
995    set_intersection(cycleModeCollect["On"].begin(),
996                     cycleModeCollect["On"].end(),
997                     cycleModeCollect["Off"].begin(),
998                     cycleModeCollect["Off"].end(),
999                     back_inserter(commonCycles)
1000                     );
1001   
1002    if (para.debuglev_>90){
1003      cout << "Liste of cycles common to On & Off: <";
1004      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
1005        cout << *i << " ";
1006      }
1007      cout << ">" << endl;
1008    }
1009
1010    //Loop on cyles
1011    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1012     
1013      string ppftag;
1014      //load ON phase
1015      stringstream cycle;
1016      cycle << *ic;
1017      ppftag = "specRawOn"+cycle.str();
1018      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1019      fin.GetObject(aSpecOn,ppftag);
1020
1021      //Load OFF phase
1022      ppftag = "specRawOff"+cycle.str();
1023      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1024      fin.GetObject(aSpecOff,ppftag);
1025     
1026      //Perform the difference ON-OFF
1027      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1028      meanDiffONOFF_noCalib += diffOnOff_noCalib;
1029     
1030      //JEC 29/10/11 add ON-OFF/OFF
1031      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1032      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1033      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1034      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
1035     
1036      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
1037      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1038
1039
1040      //JEC 15/11/11 add ON/OFF and OFF/OFF
1041      TMatrix<r_4> onOvOff(aSpecOn,false);
1042      onOvOff.Div(aSpecOffFitltered);
1043      meanONovOFF_noCalib += onOvOff;
1044     
1045      TMatrix<r_4> offOvOff(aSpecOff,false);
1046      offOvOff.Div(aSpecOffFitltered);
1047      meanOFFovOFF_noCalib += offOvOff;
1048
1049      totalNumberCycles++;
1050
1051
1052      //Fill NTuple
1053      xnt[0] = totalNumberCycles;
1054 
1055
1056      //Follow up arround the 1420.4MHz Freq.
1057      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1058      meanInRange(diffOnOffOvOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1059      xnt[1] = meanInRange_1420Freq_noCalib(0);
1060      xnt[2] = meanInRange_1420Freq_noCalib(1);
1061
1062
1063      //Follow up below the 1420.4MHz Freq.
1064      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1065      meanInRange(diffOnOffOvOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1066      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1067      meanInRange(diffOnOffOvOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1068
1069      xnt[3] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1070      xnt[4] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1071
1072
1073      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1074      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
1075      meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
1076      xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
1077      xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
1078
1079      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
1080      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1081      aSpecOffovOff.Div(aSpecOffFitltered);
1082
1083      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1084      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1085
1086      xnt[7] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1087      xnt[8] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1088
1089      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1090      aSpecOnovOff.Div(aSpecOffFitltered);
1091
1092      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1093      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1094
1095      xnt[9] = meanInRange_1410a1415Freq_ON_noCalib(0);
1096      xnt[10] = meanInRange_1410a1415Freq_ON_noCalib(1);
1097     
1098      //store infos to Ntuple
1099      onoffevolution.Fill(xnt);
1100
1101      //Quit if enough
1102      if (totalNumberCycles >= para.maxNumberCycles_) break;   
1103    }//end of cycles
1104
1105    if (totalNumberCycles >= para.maxNumberCycles_) break;         
1106
1107  }//end files
1108  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1109  //Normalisation
1110  if(totalNumberCycles > 0){
1111    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
1112    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
1113    meanONovOFF_noCalib         /= (r_4)totalNumberCycles;
1114    meanOFFovOFF_noCalib        /= (r_4)totalNumberCycles;
1115  } 
1116
1117  {//save results
1118    stringstream tmp;
1119    tmp << totalNumberCycles;
1120    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1121
1122    POutPersist fos(fileName);
1123    cout << "Save output in " << fileName << endl;
1124
1125    string tag = "meanNoCalib";
1126    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1127   
1128    tag = "meanOvOffNoCalib";
1129    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1130   
1131    tag = "meanOnovOffNoCalib";
1132    fos << PPFNameTag(tag) << meanONovOFF_noCalib;
1133
1134    tag = "meanOffovOffNoCalib";
1135    fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;
1136
1137    tag = "onoffevol";
1138    fos << PPFNameTag(tag) << onoffevolution;   
1139
1140  }//end save
1141
1142
1143}
1144//JEC 14/11/11 New meanRawDiffOnOffCycles END
1145//-------------------------------------------------------
1146//JEC 14/11/11 Obsolete START
1147//-------------------------------------------------------
1148//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1149//Used like:
1150//
1151// void meanRawDiffOnOffCycles() throw(string) {
1152//   list<string> listOfFiles;
1153//   string directoryName;
1154//   directoryName = para.inPath_ + "/" + para.sourceName_;
1155
1156//   //Make the listing of the directory
1157//   listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1158 
1159//   list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1160//   iFileEnd = listOfFiles.end();
1161 
1162//   StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1163//   TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1164//   uint_4 nSpectra=0;
1165//   //Loop on files
1166//   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1167//     if (para.debuglev_>90){
1168//       cout << "load file <" << *iFile << ">" << endl;
1169//     }
1170//     PInPersist fin(*iFile);
1171//     vector<string> vec = fin.GetNameTags();
1172//     list<string> listOfSpectra;
1173//     //Keep only required PPF objects
1174//     std::remove_copy_if(
1175//                      vec.begin(), vec.end(), back_inserter(listOfSpectra),
1176//                      not1(match)
1177//                      );
1178   
1179//     listOfSpectra.sort(stringCompare);
1180//     iSpecEnd = listOfSpectra.end();
1181//     //Loop of spectra matrix
1182//     for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
1183//       if (para.debuglev_>90){
1184//      cout << " spactra <" << *iSpec << ">" << endl;
1185//       }
1186//       TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1187//       fin.GetObject(aSpec,*iSpec);
1188//       //How to see if the GetObject is ok?? Ask Reza
1189//       nSpectra++;
1190//       meanOfSpectra+=aSpec;
1191//     }//eo loop on spectra in a file
1192//   }//eo loop on files
1193 
1194//   //Normalisation
1195//   if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
1196
1197//   //Compute the reduced version of the mean and sigma
1198//   TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1199//   TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1200//   reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
1201
1202//   {//Save the result
1203//     stringstream tmp;
1204//     tmp << nSpectra;
1205//     string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1206//     cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
1207//     POutPersist fos(fileName);
1208
1209//     string tag = "mean";
1210//     fos << PPFNameTag(tag) << meanOfSpectra;
1211//     tag = "meanred";
1212//     fos << PPFNameTag(tag) << meanRedMtx;
1213//     tag = "sigmared";
1214//     fos << PPFNameTag(tag) << sigmaRedMtx;
1215//   }
1216// }
1217//JEC 14/11/11 Obsolete END
1218//-------------------------------------------------------
1219//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1220//Used like:
1221//
1222void medianRawDiffOnOffCycles() throw(string) {
1223  list<string> listOfFiles;
1224  string directoryName;
1225  directoryName = para.inPath_ + "/" + para.sourceName_;
1226
1227  //Make the listing of the directory
1228  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1229 
1230  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1231  iFileEnd = listOfFiles.end();
1232 
1233
1234  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
1235
1236  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1237  uint_4 nSpectra=0;
1238  //Loop on files
1239  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1240    if (para.debuglev_>90){
1241      cout << "load file <" << *iFile << ">" << endl;
1242    }
1243    PInPersist fin(*iFile);
1244    vector<string> vec = fin.GetNameTags();
1245    list<string> listOfSpectra;
1246    //Keep only required PPF objects
1247    std::remove_copy_if(
1248                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
1249                        not1(match)
1250                        );
1251   
1252    listOfSpectra.sort(stringCompare);
1253    iSpecEnd = listOfSpectra.end();
1254    //Loop of spectra matrix
1255    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ;  ++iSpec){
1256      if (para.debuglev_>90){
1257        cout << " spactra <" << *iSpec << ">" << endl;
1258      }
1259      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1260      fin.GetObject(aSpec,*iSpec);
1261
1262      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
1263
1264      nSpectra++;
1265    }//eo loop on spectra in a file
1266  }//eo loop on files
1267 
1268
1269 
1270  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1271  //Compute the median for each freq. and Channel
1272  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1273    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
1274      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
1275      vector<r_4> tmp;
1276      tmp0.FillTo(tmp);
1277      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
1278    }
1279  }
1280
1281
1282  //Compute the reduced version of the mean and sigma
1283  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1284  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1285  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
1286
1287
1288  sa_size_t f1320=freqToChan(1320.);
1289  sa_size_t f1345=freqToChan(1345.);
1290  sa_size_t f1355=freqToChan(1355.);
1291  sa_size_t f1380=freqToChan(1380.);
1292  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
1293  if (para.debuglev_>9){
1294    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
1295  }
1296  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
1297  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1298    double meanMed1;
1299    double sigmaMed1;
1300    TVector<r_4> band1;
1301    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
1302    MeanSigma(band1,meanMed1,sigmaMed1);
1303    double meanMed2;
1304    double sigmaMed2;
1305    TVector<r_4> band2;
1306    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
1307    MeanSigma(band2,meanMed2,sigmaMed2);
1308    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
1309  } 
1310  meanMed.Print(cout);
1311  cout << endl;
1312
1313 
1314  //Compute the sigma in the range 1320MHz-1380MHz
1315  if (para.debuglev_>9){
1316    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
1317  }
1318  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
1319  sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1320  sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1321
1322  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1323    double meanSigma;
1324    double sigmaSigma;
1325    TVector<r_4> band;
1326    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
1327    MeanSigma(band,meanSigma,sigmaSigma);
1328    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
1329    sigmaMed(iCh) = meanSigma;
1330  }
1331  sigmaMed.Print(cout);
1332  cout << endl;
1333
1334 
1335 
1336  if (para.debuglev_>9){
1337    cout << "Compute medianOfSpectraNorm" << endl;
1338  }
1339  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
1340  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1341    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
1342    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
1343  }
1344
1345 
1346
1347  {//Save the result
1348    stringstream tmp;
1349    tmp << nSpectra;
1350    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1351    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
1352    POutPersist fos(fileName);
1353
1354    string tag = "median";
1355    fos << PPFNameTag(tag) << medianOfSpectra;
1356
1357    tag = "medianNorm";
1358    fos << PPFNameTag(tag) << medianOfSpectraNorm;
1359   
1360
1361    tag = "meanmedred";
1362    fos << PPFNameTag(tag) << meanRedMtx;
1363    tag = "sigmamedred";
1364    fos << PPFNameTag(tag) << sigmaRedMtx;
1365    tag = "cycleVsfreq";
1366   
1367    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1368    fos << PPFNameTag(tag) << tarr;
1369  }
1370}
1371
1372//-------------------------------------------------------
1373int main(int narg, char* arg[]) {
1374 
1375  int rc = 0; //return code
1376  string msg; //message used in Exceptions
1377
1378
1379
1380  //default value for initial parameters (see Para structure on top of the file)
1381  string debuglev = "0";
1382  string action = "meanDiffOnOff";
1383  string inputPath = "."; 
1384  string outputPath = "."; 
1385  string sourceName = "Abell85";
1386  string ppfFile;
1387  string nSliceInFreq = "32";
1388  string typeOfCalib="perRun";
1389  string calibFreq = "1346";
1390  string calibBandFreq="6.25";
1391  string mxcycles;
1392
1393  //  bool okarg=false;
1394  int ka=1;
1395  while (ka<narg) {
1396    if (strcmp(arg[ka],"-h")==0) {
1397      cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1398           << " -mxcycles <number> (max. number of cycles to be treated)\n"
1399           << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1400           << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1401           << " -src [Abell85]\n" 
1402           << " -inPath [.]|<top_root_dir of the ppf file>\n" 
1403           << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
1404           << " -outPath [.]|<dir of the output> \n"
1405           << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1406           << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1407           << " -debug <level> "
1408           << endl;
1409      return 0;
1410    }
1411    else if (strcmp(arg[ka],"-debug")==0) {
1412      debuglev=arg[ka+1];
1413      ka+=2;
1414    }
1415    else if (strcmp(arg[ka],"-act")==0) {
1416      action=arg[ka+1];
1417      ka+=2;
1418    }
1419    else if (strcmp(arg[ka],"-calibfreq")==0) {
1420      calibFreq=arg[ka+1];
1421      ka+=2;
1422    }   
1423    else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1424      calibBandFreq=arg[ka+1];
1425      ka+=2;
1426    }   
1427    else if (strcmp(arg[ka],"-mxcycles")==0) {
1428      mxcycles=arg[ka+1];
1429      ka+=2;
1430    }   
1431    else if (strcmp(arg[ka],"-inPath")==0) {
1432      inputPath=arg[ka+1];
1433      ka+=2;
1434    }
1435    else if (strcmp(arg[ka],"-outPath")==0) {
1436      outputPath=arg[ka+1];
1437      ka+=2;
1438    }
1439    else if (strcmp(arg[ka],"-src")==0) {
1440      sourceName=arg[ka+1];
1441      ka+=2;
1442    }
1443    else if (strcmp(arg[ka],"-ppfFile")==0) {
1444      ppfFile=arg[ka+1];
1445      ka+=2;
1446    }
1447    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1448      nSliceInFreq=arg[ka+1];
1449      ka+=2;
1450    }
1451    else ka++;
1452  }//eo while
1453
1454  para.debuglev_   = atoi(debuglev.c_str());
1455  para.inPath_     = inputPath;
1456  para.outPath_    = outputPath;
1457  para.sourceName_ = sourceName;
1458  para.ppfFile_    = ppfFile;
1459  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1460  para.calibFreq_   = calibFreq;
1461  para.calibBandFreq_ = calibBandFreq;
1462  para.rcalibFreq_   = atof(calibFreq.c_str());
1463  para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1464  if (mxcycles != "") {
1465    para.maxNumberCycles_ = atoi(mxcycles.c_str());
1466  } else {
1467    para.maxNumberCycles_ = std::numeric_limits<int>::max();
1468  }
1469
1470  cout << "Dump Initial parameters ............" << endl;
1471  cout << " action = " << action << "\n"
1472       << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1473       << " inputPath = " << para.inPath_  << "\n" 
1474       << " outputPath = " <<  para.outPath_ << "\n"
1475       << " sourceName = " << para.sourceName_ << "\n"
1476       << " ppfFile = " <<  para.ppfFile_ << "\n"
1477       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
1478       << " calibFreq = " <<  para.calibFreq_ << "\n"
1479       << " calibBandFreq = " <<  para.calibBandFreq_ << "\n"
1480       << " debuglev = "  << para.debuglev_   << "\n";
1481  cout << "...................................." << endl;
1482
1483  if ( "" == ppfFile ) {
1484    cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1485         << endl;
1486    return 999;
1487  }
1488
1489
1490  try {
1491
1492//     int b,e;
1493//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1494//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1495
1496    if ( action == "meanRawDiffOnOff" ) {
1497      meanRawDiffOnOffCycles();
1498    } else if (action == "medianRawDiffOnOff") {
1499      medianRawDiffOnOffCycles();
1500    } else if (action == "meanCalibBAODiffOnOff") {
1501      meanCalibBAODiffOnOffCycles();
1502    } else {
1503      msg = "Unknown action " + action;
1504      throw msg;
1505    }
1506
1507
1508  }  catch (std::exception& sex) {
1509    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
1510         << "\n msg= " << sex.what() << endl;
1511    rc = 78;
1512  }
1513  catch ( string str ) {
1514    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1515  }
1516  catch (...) {
1517    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
1518    rc = 79; 
1519  } 
1520
1521  return 0;
1522
1523}
Note: See TracBrowser for help on using the repository browser.