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

Last change on this file since 611 was 611, checked in by torrento, 13 years ago

Add pic analysis files and emptydir function for proc_specmfib
File size: 54.2 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    PInPersist fin(*iFile);
400    vector<string> vec = fin.GetNameTags();
401
402    vector<string> modeList;
403    modeList.push_back("On");
404    modeList.push_back("Off");
405    vector<string>::const_iterator iMode;
406   
407    map<string, list<int> > cycleModeCollect;
408   
409    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
410      list<string> listOfSpectra;
411      //Keep only required PPF objects
412      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
413      std::remove_copy_if(
414                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
415                          not1(StringMatch(matchstr))
416                          );
417     
418      listOfSpectra.sort(stringCompare);
419      iSpecEnd = listOfSpectra.end();
420     
421      matchstr = "[0-9]+";
422      //Loop of spectra matrix
423      list<int> listOfCycles;
424      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
425        int b,e;
426        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
427        if (para.debuglev_>90){
428          cout << " spactra <" << *iSpec << ">" << endl;
429          cout << " cycle " << iSpec->substr(b) << endl;
430        }
431        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
432      }//end loop spectra
433      cycleModeCollect[*iMode] = listOfCycles;
434    }//end of mode   
435
436    //Take the Intersection of the list Of cycles in mode Off and On   
437    list<int> commonCycles;
438    set_intersection(cycleModeCollect["On"].begin(),
439                     cycleModeCollect["On"].end(),
440                     cycleModeCollect["Off"].begin(),
441                     cycleModeCollect["Off"].end(),
442                     back_inserter(commonCycles)
443                     );
444   
445    if (para.debuglev_>90){
446      cout << "Liste of cycles common to On & Off: <";
447      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
448        cout << *i << " ";
449      }
450      cout << ">" << endl;
451    }
452   
453    //
454    //Load BAO Calibration factors "per Cycle and Channels"
455    //Compute the mean per Cycle to
456    //  fill factors "per Run and Channels" with the same cycle list
457    //
458    //
459    //TODO improve the code....
460
461    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
462    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
463    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
464    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
465   
466    string calibFileName;
467    ifstream ifs;
468    sa_size_t nr,nc; //values read
469
470    //OFF Cycle per Channel
471    calibFileName = directoryName + "/" 
472      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
473      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
474    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
475    ifs.open(calibFileName.c_str());
476    if ( ! ifs.is_open() ) {
477
478      throw calibFileName + " cannot be opened...";
479    }   
480    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
481    if(para.debuglev_>9){
482      cout << "(nr,nc): "<< nr << "," << nc << endl;
483      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
484      cout << endl;
485    }
486
487    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
488    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
489    {//Compute the mean
490      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
491      double mean,sigma;
492      MeanSigma(coef,mean,sigma);
493      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
494    }
495    if(para.debuglev_>9){
496      cout << "Fill calib. with mean value " << endl; 
497      calibBAOfactors_Off_Run_Ch0.Print(cout);
498      cout << endl;
499    }
500    ifs.close();
501
502    //
503    calibFileName = directoryName + "/" 
504      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
505      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
506    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
507    ifs.open(calibFileName.c_str());
508    if ( ! ifs.is_open() ) {
509
510      throw calibFileName + " cannot be opened...";
511    }   
512    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
513    if(para.debuglev_>9){
514      cout << "(nr,nc): "<< nr << "," << nc << endl;
515      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
516      cout << endl;
517    }
518    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
519    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
520    {//Compute the mean
521      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
522      double mean,sigma;
523      MeanSigma(coef,mean,sigma);
524      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
525      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
526    }
527    if(para.debuglev_>9){
528      cout << "Fill calib. with mean value " << endl; 
529      calibBAOfactors_Off_Run_Ch1.Print(cout);
530      cout << endl;
531    }
532    ifs.close();
533
534    //ON Cycle per Channel
535    calibFileName = directoryName + "/" 
536      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
537      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
538    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
539    ifs.open(calibFileName.c_str());
540    if ( ! ifs.is_open() ) {
541
542      throw calibFileName + " cannot be opened...";
543    }   
544    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
545    if(para.debuglev_>9){
546      cout << "(nr,nc): "<< nr << "," << nc << endl;
547      calibBAOfactors_On_Cycle_Ch0.Print(cout);
548      cout << endl;     
549    }
550
551    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
552    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
553    {//Compute the mean
554      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
555      double mean,sigma;
556      MeanSigma(coef,mean,sigma);
557      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
558      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
559    }
560    if(para.debuglev_>9){
561      cout << "Fill calib. with mean value " << endl; 
562      calibBAOfactors_On_Run_Ch0.Print(cout);
563      cout << endl;
564    }
565    ifs.close();
566
567   
568    calibFileName = directoryName + "/" 
569      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
570      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
571    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
572    ifs.open(calibFileName.c_str());
573    if ( ! ifs.is_open() ) {
574      throw calibFileName + " cannot be opened...";
575    }   
576    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
577    if(para.debuglev_>9){
578      cout << "(nr,nc): "<< nr << "," << nc << endl;
579      calibBAOfactors_On_Cycle_Ch1.Print(cout);
580      cout << endl;
581    }
582    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
583    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
584    {//Compute the mean
585      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
586      double mean,sigma;
587      MeanSigma(coef,mean,sigma);
588      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
589      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
590    }
591    if(para.debuglev_>9){
592      cout << "Fill calib. with mean value " << endl; 
593      calibBAOfactors_On_Run_Ch1.Print(cout);
594      cout << endl;
595    }
596
597    ifs.close();
598   
599    //link <cycle> - <calibration coefficient>
600    //We cannot rely on identical cycle list of the OFF and ON calibration
601    map<int,r_4> calibBAO_Off_Run_Ch0;
602    map<int,r_4> calibBAO_Off_Run_Ch1;
603    map<int,r_4> calibBAO_On_Run_Ch0;
604    map<int,r_4> calibBAO_On_Run_Ch1;
605
606    map<int,r_4> calibBAO_Off_Cycle_Ch0;
607    map<int,r_4> calibBAO_Off_Cycle_Ch1;
608    map<int,r_4> calibBAO_On_Cycle_Ch0;
609    map<int,r_4> calibBAO_On_Cycle_Ch1;
610
611    //per Run based BAO coefficients
612    nr = calibBAOfactors_Off_Run_Ch0.NRows();
613    for (sa_size_t ir=0; ir<nr; ++ir){
614      cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
615           << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
616
617      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
618        = calibBAOfactors_Off_Run_Ch0(ir,1);
619      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
620        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
621      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
622        = calibBAOfactors_Off_Run_Ch1(ir,1);
623      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
624        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
625    }//eo loop on coef Off
626   
627    nr = calibBAOfactors_On_Run_Ch0.NRows();
628    for (sa_size_t ir=0; ir<nr; ++ir){
629      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
630        = calibBAOfactors_On_Run_Ch0(ir,1);
631      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
632        = calibBAOfactors_On_Cycle_Ch0(ir,1);
633      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
634        = calibBAOfactors_On_Run_Ch1(ir,1);
635      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
636        = calibBAOfactors_On_Cycle_Ch1(ir,1);
637    }//eo loop on coef On
638     
639    //Loop on cyles
640    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
641
642      //Look if the cycle has been calibrated...
643      bool isCycleCalibrated = 
644        ( calibBAO_On_Run_Ch0.count(*ic)    *
645          calibBAO_On_Run_Ch1.count(*ic)    *
646          calibBAO_Off_Run_Ch0.count(*ic)   *
647          calibBAO_Off_Run_Ch1.count(*ic)   *
648          calibBAO_On_Cycle_Ch0.count(*ic)  *
649          calibBAO_On_Cycle_Ch1.count(*ic)  *
650          calibBAO_Off_Cycle_Ch0.count(*ic) *
651          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
652
653      if(para.debuglev_>9) {
654        cout << "Calibration coefficients for cycle "<<*ic << endl; 
655        if (isCycleCalibrated) {
656          cout << "Cycle calibrated " << endl;
657          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
658               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
659               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
660               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
661               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
662               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
663               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
664               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
665        } else {
666          cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
667        }
668      }//debug
669
670
671      if ( ! isCycleCalibrated ) continue;
672     
673      string ppftag;
674      //load ON phase
675      stringstream cycle;
676      cycle << *ic;
677     
678      ppftag = "specRawOn"+cycle.str();
679      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
680      fin.GetObject(aSpecOn,ppftag);
681
682      TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
683      aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
684      aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
685
686      TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
687      aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
688      aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
689     
690      //Load OFF phase
691      ppftag = "specRawOff"+cycle.str();
692      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
693      fin.GetObject(aSpecOff,ppftag);
694
695      TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
696      aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
697      aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
698
699      TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
700      aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
701      aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
702
703
704      //Perform the difference ON-OFF with the different calibration options
705      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
706      meanDiffONOFF_noCalib += diffOnOff_noCalib;
707     
708      //JEC 29/10/11 add ON-OFF/OFF
709      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
710      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
711      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
712      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
713     
714      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
715      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
716
717      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
718      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
719
720      TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
721      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
722
723      totalNumberCycles++;
724
725      //Fill NTuple
726      xnt[0] = totalNumberCycles;
727 
728      //Follow up arround the Calibration Freq.
729      TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
730      meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
731      xnt[1] = meanInRange_CalibFreq_noCalib(0);
732      xnt[2] = meanInRange_CalibFreq_noCalib(1);
733
734      TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
735      meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
736      xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
737      xnt[4] = meanInRange_CalibFreq_perRunCalib(1);
738
739      TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
740      meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
741      xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
742      xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);
743
744
745      //Follow up arround the 1420.4MHz Freq.
746      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
747      meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
748      xnt[7] = meanInRange_1420Freq_noCalib(0);
749      xnt[8] = meanInRange_1420Freq_noCalib(1);
750
751      TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
752      meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
753      xnt[9] = meanInRange_1420Freq_perRunCalib(0);
754      xnt[10] = meanInRange_1420Freq_perRunCalib(1);
755
756      TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
757      meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
758      xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
759      xnt[12] = meanInRange_1420Freq_perCycleCalib(1);
760
761
762      //Follow up below the 1420.4MHz Freq.
763      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
764      meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
765      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
766      meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
767
768      xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
769      xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
770
771
772      TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
773      meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
774      TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
775      meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);
776
777      xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
778      xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;
779
780
781      TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
782      meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
783      TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
784      meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);
785
786      xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
787      xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
788
789
790      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
791      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
792      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
793      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
794      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
795
796      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
797      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
798      aSpecOffovOff.Div(aSpecOffFiltered);
799
800      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
801      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
802      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
803
804      xnt[21] = meanInRange_1410a1415Freq_OFF_noCalib(0);
805      xnt[22] = meanInRange_1410a1415Freq_OFF_noCalib(1);
806
807      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
808      aSpecOnovOff.Div(aSpecOffFiltered);
809
810      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
811      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
812      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
813
814      xnt[23] = meanInRange_1410a1415Freq_ON_noCalib(0);
815      xnt[24] = meanInRange_1410a1415Freq_ON_noCalib(1);
816     
817      //store infos to Ntuple
818      onoffevolution.Fill(xnt);
819
820      //Quit if enough
821      if (totalNumberCycles >= para.maxNumberCycles_) break;   
822
823    }//eo loop on cycles
824    if (totalNumberCycles >= para.maxNumberCycles_) break;         
825 
826  }//eo loop on spectra in a file
827  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
828  //Normalisation
829  if(totalNumberCycles > 0){
830    //JEC 29/10 add ON-OFF/OFF
831    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
832    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
833    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
834    meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
835  } 
836 
837  //Compute the reduced version of the mean and sigma
838  TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
839  TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
840  reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
841
842  TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
843  TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
844  reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
845
846  TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
847  TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
848  reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
849
850  {//save the results
851    stringstream tmp;
852    tmp << totalNumberCycles;
853    string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
854
855    POutPersist fos(fileName);
856    cout << "Save output in " << fileName << endl;
857
858    string tag = "meanNoCalib";
859    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
860   
861    //JEC 29/10/11
862    tag = "meanOvOffNoCalib";
863    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
864
865    tag = "meanPerRunCalib";
866    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
867    tag = "meanPerCycleCalib";
868    fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
869
870    tag = "redmeanNoCalib";
871    fos << PPFNameTag(tag) << meanRedMtx_noCalib;
872    tag = "redsigmaNoCalib";
873    fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
874
875    tag = "redmeanPerRunCalib";
876    fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
877    tag = "redsigmaPerRunCalib";
878    fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
879
880    tag = "redmeanPerCycleCalib";
881    fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
882    tag = "redsigmaPerCycleCalib";
883    fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
884   
885    tag = "onoffevol";
886    fos << PPFNameTag(tag) << onoffevolution;
887 
888  }//end of save
889
890
891}
892//JEC 14/11/11 New meanRawDiffOnOffCycles START
893//-------------------------------------------------------
894//Compute the mean of Diff ON-OFF/OFF from Raw spectra
895//Used like:
896//
897void meanRawDiffOnOffCycles() throw(string) {
898  list<string> listOfFiles;
899  string directoryName;
900  directoryName = para.inPath_ + "/" + para.sourceName_;
901
902  //Make the listing of the directory
903  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
904 
905  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
906  iFileEnd = listOfFiles.end();
907
908  //mean ON-OFF over the list of cycles
909  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
910  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
911  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
912  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
913  //Tuple only for RAW things to follow
914  static const int NINFO=11;
915  char* onffTupleName[NINFO]={"cycle"
916                              ,"onoffRaw01420","onoffRaw11420"
917                              ,"onoffRaw01420side","onoffRaw11420side"
918                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
919                              ,"offRaw0f14101415","offRaw1f14101415"
920                              ,"onRaw0f14101415","onRaw1f14101415"
921  };
922  NTuple onoffevolution(NINFO,onffTupleName);
923  r_4 xnt[NINFO];
924
925  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
926  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
927  sa_size_t ch1420High = freqToChan(1420.4+0.2);
928
929  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
930  sa_size_t ch1420aLow  = freqToChan(1418);
931  sa_size_t ch1420aHigh = freqToChan(1419);
932  sa_size_t ch1420bLow  = freqToChan(1422);
933  sa_size_t ch1420bHigh = freqToChan(1423);
934 
935  //1400-1420Mhz
936  sa_size_t ch1400 = freqToChan(1400);
937  //  sa_size_t ch1405 = freqToChan(1400);
938  sa_size_t ch1410 = freqToChan(1410);
939  sa_size_t ch1415 = freqToChan(1415);
940  sa_size_t ch1420 = freqToChan(1420);
941
942  if (para.debuglev_>0){
943    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
944    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
945    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
946  }
947
948  int totalNumberCycles=0; //total number of cycles
949
950  //Loop on files
951  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
952    if (para.debuglev_>90){
953      cout << "load file <" << *iFile << ">" << endl;
954    }
955
956    vector<string> tokens;
957    split(*iFile,"_",tokens);
958    string dateOfRun = tokens[1];
959    if (para.debuglev_>90){
960      cout << "date <" << dateOfRun << ">" << endl;
961    }
962    vector<string> tokens2;
963    split(tokens[2],".",tokens2);
964    string srcLower = tokens2[0];
965
966    PInPersist fin(*iFile);
967    vector<string> vec = fin.GetNameTags();
968
969    vector<string> modeList;
970    modeList.push_back("On");
971    modeList.push_back("Off");
972    vector<string>::const_iterator iMode;
973   
974    map<string, list<int> > cycleModeCollect;
975   
976    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
977      list<string> listOfSpectra;
978      //Keep only required PPF objects
979      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
980      std::remove_copy_if(
981                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
982                          not1(StringMatch(matchstr))
983                          );
984     
985      listOfSpectra.sort(stringCompare);
986      iSpecEnd = listOfSpectra.end();
987     
988      matchstr = "[0-9]+";
989      //Loop of spectra matrix
990      list<int> listOfCycles;
991      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
992        int b,e;
993        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
994        if (para.debuglev_>90){
995          cout << " spectra <" << *iSpec << ">" << endl;
996          cout << " cycle " << iSpec->substr(b) << endl;
997        }
998        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
999      }//end loop spectra
1000      cycleModeCollect[*iMode] = listOfCycles;
1001    }//end of mode   
1002
1003    //Take the Intersection of the list Of cycles in mode Off and On   
1004    list<int> commonCycles;
1005    set_intersection(cycleModeCollect["On"].begin(),
1006                     cycleModeCollect["On"].end(),
1007                     cycleModeCollect["Off"].begin(),
1008                     cycleModeCollect["Off"].end(),
1009                     back_inserter(commonCycles)
1010                     );
1011   
1012    if (para.debuglev_>90){
1013      cout << "Liste of cycles common to On & Off: <";
1014      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
1015        cout << *i << " ";
1016      }
1017      cout << ">" << endl;
1018    }
1019
1020    //Loop on cyles
1021    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1022     
1023      // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440
1024      if ( *ic == 1 && srcLower == "abell1205" ) {
1025          if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {
1026            cout << "Skipping non-calibrated cycle " << *ic << endl;
1027            continue;
1028          }
1029      } else if ( *ic == 1 && srcLower == "abell2440" ) {
1030          if ( dateOfRun == "20110516" ) {
1031            cout << "Skipping non-calibrated cycle " << *ic << endl;
1032            continue;
1033          }
1034      } else if ( *ic == 3 && srcLower == "abell1205" ) {
1035          if ( dateOfRun == "20110810" ) {
1036            cout << "Skipping non-calibrated cycle " << *ic << endl;
1037            continue;
1038          }
1039      }
1040
1041      string ppftag;
1042      //load ON phase
1043      stringstream cycle;
1044      cycle << *ic;
1045      ppftag = "specRawOn"+cycle.str();
1046      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1047      fin.GetObject(aSpecOn,ppftag);
1048
1049      //Load OFF phase
1050      ppftag = "specRawOff"+cycle.str();
1051      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1052      fin.GetObject(aSpecOff,ppftag);
1053     
1054      //Perform the difference ON-OFF
1055      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1056
1057      meanDiffONOFF_noCalib += diffOnOff_noCalib;
1058     
1059      //JEC 29/10/11 add ON-OFF/OFF
1060      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1061      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1062      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1063      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
1064     
1065      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
1066      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1067
1068      //JEC 15/11/11 add ON/OFF and OFF/OFF
1069      TMatrix<r_4> onOvOff(aSpecOn,false);
1070      onOvOff.Div(aSpecOffFiltered);
1071      meanONovOFF_noCalib += onOvOff;
1072     
1073      TMatrix<r_4> offOvOff(aSpecOff,false);
1074      offOvOff.Div(aSpecOffFiltered);
1075      meanOFFovOFF_noCalib += offOvOff;
1076
1077      totalNumberCycles++;
1078
1079      //Fill NTuple
1080      xnt[0] = totalNumberCycles;
1081 
1082      //Follow up arround the 1420.4MHz Freq.
1083      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1084      meanInRange(diffOnOffOvOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1085      xnt[1] = meanInRange_1420Freq_noCalib(0);
1086      xnt[2] = meanInRange_1420Freq_noCalib(1);
1087
1088
1089      //Follow up below the 1420.4MHz Freq.
1090      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1091      meanInRange(diffOnOffOvOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1092      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1093      meanInRange(diffOnOffOvOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1094
1095      xnt[3] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1096      xnt[4] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1097
1098
1099      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1100      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
1101      meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
1102      xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
1103      xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
1104
1105      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
1106      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1107      aSpecOffovOff.Div(aSpecOffFiltered);
1108
1109      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1110      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1111
1112      xnt[7] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1113      xnt[8] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1114
1115      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1116      aSpecOnovOff.Div(aSpecOffFiltered);
1117
1118      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1119      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1120
1121      xnt[9] = meanInRange_1410a1415Freq_ON_noCalib(0);
1122      xnt[10] = meanInRange_1410a1415Freq_ON_noCalib(1);
1123     
1124      //store infos to Ntuple
1125      onoffevolution.Fill(xnt);
1126
1127      //Quit if enough
1128      if (totalNumberCycles >= para.maxNumberCycles_) break;   
1129    }//end of cycles
1130
1131    if (totalNumberCycles >= para.maxNumberCycles_) break;         
1132
1133  }//end files
1134  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1135  //Normalisation
1136  if(totalNumberCycles > 0){
1137    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
1138    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
1139    meanONovOFF_noCalib         /= (r_4)totalNumberCycles;
1140    meanOFFovOFF_noCalib        /= (r_4)totalNumberCycles;
1141  } 
1142
1143  {//save results
1144    stringstream tmp;
1145    tmp << totalNumberCycles;
1146    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1147
1148    POutPersist fos(fileName);
1149    cout << "Save output in " << fileName << endl;
1150
1151    string tag = "meanNoCalib";
1152    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1153   
1154    tag = "meanOvOffNoCalib";
1155    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1156   
1157    tag = "meanOnovOffNoCalib";
1158    fos << PPFNameTag(tag) << meanONovOFF_noCalib;
1159
1160    tag = "meanOffovOffNoCalib";
1161    fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;
1162
1163    tag = "onoffevol";
1164    fos << PPFNameTag(tag) << onoffevolution;
1165
1166  }//end save
1167
1168
1169}
1170//JEC 14/11/11 New meanRawDiffOnOffCycles END
1171//-------------------------------------------------------
1172//JEC 14/11/11 Obsolete START
1173//-------------------------------------------------------
1174//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1175//Used like:
1176//
1177// void meanRawDiffOnOffCycles() throw(string) {
1178//   list<string> listOfFiles;
1179//   string directoryName;
1180//   directoryName = para.inPath_ + "/" + para.sourceName_;
1181
1182//   //Make the listing of the directory
1183//   listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1184 
1185//   list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1186//   iFileEnd = listOfFiles.end();
1187 
1188//   StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1189//   TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1190//   uint_4 nSpectra=0;
1191//   //Loop on files
1192//   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1193//     if (para.debuglev_>90){
1194//       cout << "load file <" << *iFile << ">" << endl;
1195//     }
1196//     PInPersist fin(*iFile);
1197//     vector<string> vec = fin.GetNameTags();
1198//     list<string> listOfSpectra;
1199//     //Keep only required PPF objects
1200//     std::remove_copy_if(
1201//                      vec.begin(), vec.end(), back_inserter(listOfSpectra),
1202//                      not1(match)
1203//                      );
1204   
1205//     listOfSpectra.sort(stringCompare);
1206//     iSpecEnd = listOfSpectra.end();
1207//     //Loop of spectra matrix
1208//     for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
1209//       if (para.debuglev_>90){
1210//      cout << " spactra <" << *iSpec << ">" << endl;
1211//       }
1212//       TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1213//       fin.GetObject(aSpec,*iSpec);
1214//       //How to see if the GetObject is ok?? Ask Reza
1215//       nSpectra++;
1216//       meanOfSpectra+=aSpec;
1217//     }//eo loop on spectra in a file
1218//   }//eo loop on files
1219 
1220//   //Normalisation
1221//   if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
1222
1223//   //Compute the reduced version of the mean and sigma
1224//   TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1225//   TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1226//   reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
1227
1228//   {//Save the result
1229//     stringstream tmp;
1230//     tmp << nSpectra;
1231//     string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1232//     cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
1233//     POutPersist fos(fileName);
1234
1235//     string tag = "mean";
1236//     fos << PPFNameTag(tag) << meanOfSpectra;
1237//     tag = "meanred";
1238//     fos << PPFNameTag(tag) << meanRedMtx;
1239//     tag = "sigmared";
1240//     fos << PPFNameTag(tag) << sigmaRedMtx;
1241//   }
1242// }
1243//JEC 14/11/11 Obsolete END
1244//-------------------------------------------------------
1245//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1246//Used like:
1247//
1248void medianRawDiffOnOffCycles() throw(string) {
1249  list<string> listOfFiles;
1250  string directoryName;
1251  directoryName = para.inPath_ + "/" + para.sourceName_;
1252
1253  //Make the listing of the directory
1254  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1255 
1256  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1257  iFileEnd = listOfFiles.end();
1258 
1259
1260  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
1261
1262  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1263  uint_4 nSpectra=0;
1264  //Loop on files
1265  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1266    if (para.debuglev_>90){
1267      cout << "load file <" << *iFile << ">" << endl;
1268    }
1269    PInPersist fin(*iFile);
1270    vector<string> vec = fin.GetNameTags();
1271    list<string> listOfSpectra;
1272    //Keep only required PPF objects
1273    std::remove_copy_if(
1274                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
1275                        not1(match)
1276                        );
1277   
1278    listOfSpectra.sort(stringCompare);
1279    iSpecEnd = listOfSpectra.end();
1280    //Loop of spectra matrix
1281    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ;  ++iSpec){
1282      if (para.debuglev_>90){
1283        cout << " spactra <" << *iSpec << ">" << endl;
1284      }
1285      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1286      fin.GetObject(aSpec,*iSpec);
1287
1288      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
1289
1290      nSpectra++;
1291    }//eo loop on spectra in a file
1292  }//eo loop on files
1293 
1294
1295 
1296  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1297  //Compute the median for each freq. and Channel
1298  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1299    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
1300      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
1301      vector<r_4> tmp;
1302      tmp0.FillTo(tmp);
1303      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
1304    }
1305  }
1306
1307
1308  //Compute the reduced version of the mean and sigma
1309  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1310  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1311  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
1312
1313
1314  sa_size_t f1320=freqToChan(1320.);
1315  sa_size_t f1345=freqToChan(1345.);
1316  sa_size_t f1355=freqToChan(1355.);
1317  sa_size_t f1380=freqToChan(1380.);
1318  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
1319  if (para.debuglev_>9){
1320    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
1321  }
1322  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
1323  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1324    double meanMed1;
1325    double sigmaMed1;
1326    TVector<r_4> band1;
1327    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
1328    MeanSigma(band1,meanMed1,sigmaMed1);
1329    double meanMed2;
1330    double sigmaMed2;
1331    TVector<r_4> band2;
1332    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
1333    MeanSigma(band2,meanMed2,sigmaMed2);
1334    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
1335  } 
1336  meanMed.Print(cout);
1337  cout << endl;
1338
1339 
1340  //Compute the sigma in the range 1320MHz-1380MHz
1341  if (para.debuglev_>9){
1342    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
1343  }
1344  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
1345  sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1346  sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1347
1348  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1349    double meanSigma;
1350    double sigmaSigma;
1351    TVector<r_4> band;
1352    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
1353    MeanSigma(band,meanSigma,sigmaSigma);
1354    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
1355    sigmaMed(iCh) = meanSigma;
1356  }
1357  sigmaMed.Print(cout);
1358  cout << endl;
1359
1360 
1361 
1362  if (para.debuglev_>9){
1363    cout << "Compute medianOfSpectraNorm" << endl;
1364  }
1365  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
1366  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1367    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
1368    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
1369  }
1370
1371 
1372
1373  {//Save the result
1374    stringstream tmp;
1375    tmp << nSpectra;
1376    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1377    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
1378    POutPersist fos(fileName);
1379
1380    string tag = "median";
1381    fos << PPFNameTag(tag) << medianOfSpectra;
1382
1383    tag = "medianNorm";
1384    fos << PPFNameTag(tag) << medianOfSpectraNorm;
1385   
1386
1387    tag = "meanmedred";
1388    fos << PPFNameTag(tag) << meanRedMtx;
1389    tag = "sigmamedred";
1390    fos << PPFNameTag(tag) << sigmaRedMtx;
1391    tag = "cycleVsfreq";
1392   
1393    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1394    fos << PPFNameTag(tag) << tarr;
1395  }
1396}
1397
1398//-------------------------------------------------------
1399int main(int narg, char* arg[]) {
1400 
1401  int rc = 0; //return code
1402  string msg; //message used in Exceptions
1403
1404
1405
1406  //default value for initial parameters (see Para structure on top of the file)
1407  string debuglev = "0";
1408  string action = "meanDiffOnOff";
1409  string inputPath = "."; 
1410  string outputPath = "."; 
1411  string sourceName = "Abell85";
1412  string ppfFile;
1413  string nSliceInFreq = "32";
1414  string typeOfCalib="perRun";
1415  string calibFreq = "1346";
1416  string calibBandFreq="6.25";
1417  string mxcycles;
1418
1419  //  bool okarg=false;
1420  int ka=1;
1421  while (ka<narg) {
1422    if (strcmp(arg[ka],"-h")==0) {
1423      cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1424           << " -mxcycles <number> (max. number of cycles to be treated)\n"
1425           << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1426           << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1427           << " -src [Abell85]\n" 
1428           << " -inPath [.]|<top_root_dir of the ppf file>\n" 
1429           << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
1430           << " -outPath [.]|<dir of the output> \n"
1431           << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1432           << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1433           << " -debug <level> "
1434           << endl;
1435      return 0;
1436    }
1437    else if (strcmp(arg[ka],"-debug")==0) {
1438      debuglev=arg[ka+1];
1439      ka+=2;
1440    }
1441    else if (strcmp(arg[ka],"-act")==0) {
1442      action=arg[ka+1];
1443      ka+=2;
1444    }
1445    else if (strcmp(arg[ka],"-calibfreq")==0) {
1446      calibFreq=arg[ka+1];
1447      ka+=2;
1448    }   
1449    else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1450      calibBandFreq=arg[ka+1];
1451      ka+=2;
1452    }   
1453    else if (strcmp(arg[ka],"-mxcycles")==0) {
1454      mxcycles=arg[ka+1];
1455      ka+=2;
1456    }   
1457    else if (strcmp(arg[ka],"-inPath")==0) {
1458      inputPath=arg[ka+1];
1459      ka+=2;
1460    }
1461    else if (strcmp(arg[ka],"-outPath")==0) {
1462      outputPath=arg[ka+1];
1463      ka+=2;
1464    }
1465    else if (strcmp(arg[ka],"-src")==0) {
1466      sourceName=arg[ka+1];
1467      ka+=2;
1468    }
1469    else if (strcmp(arg[ka],"-ppfFile")==0) {
1470      ppfFile=arg[ka+1];
1471      ka+=2;
1472    }
1473    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1474      nSliceInFreq=arg[ka+1];
1475      ka+=2;
1476    }
1477    else ka++;
1478  }//eo while
1479
1480  para.debuglev_   = atoi(debuglev.c_str());
1481  para.inPath_     = inputPath;
1482  para.outPath_    = outputPath;
1483  para.sourceName_ = sourceName;
1484  para.ppfFile_    = ppfFile;
1485  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1486  para.calibFreq_   = calibFreq;
1487  para.calibBandFreq_ = calibBandFreq;
1488  para.rcalibFreq_   = atof(calibFreq.c_str());
1489  para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1490  if (mxcycles != "") {
1491    para.maxNumberCycles_ = atoi(mxcycles.c_str());
1492  } else {
1493    para.maxNumberCycles_ = std::numeric_limits<int>::max();
1494  }
1495
1496  cout << "Dump Initial parameters ............" << endl;
1497  cout << " action = " << action << "\n"
1498       << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1499       << " inputPath = " << para.inPath_  << "\n" 
1500       << " outputPath = " <<  para.outPath_ << "\n"
1501       << " sourceName = " << para.sourceName_ << "\n"
1502       << " ppfFile = " <<  para.ppfFile_ << "\n"
1503       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
1504       << " calibFreq = " <<  para.calibFreq_ << "\n"
1505       << " calibBandFreq = " <<  para.calibBandFreq_ << "\n"
1506       << " debuglev = "  << para.debuglev_   << "\n";
1507  cout << "...................................." << endl;
1508
1509  if ( "" == ppfFile ) {
1510    cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1511         << endl;
1512    return 999;
1513  }
1514
1515
1516  try {
1517
1518//     int b,e;
1519//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1520//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1521
1522    if ( action == "meanRawDiffOnOff" ) {
1523      meanRawDiffOnOffCycles();
1524    } else if (action == "medianRawDiffOnOff") {
1525      medianRawDiffOnOffCycles();
1526    } else if (action == "meanCalibBAODiffOnOff") {
1527      meanCalibBAODiffOnOffCycles();
1528    } else {
1529      msg = "Unknown action " + action;
1530      throw msg;
1531    }
1532
1533
1534  }  catch (std::exception& sex) {
1535    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
1536         << "\n msg= " << sex.what() << endl;
1537    rc = 78;
1538  }
1539  catch ( string str ) {
1540    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1541  }
1542  catch (...) {
1543    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
1544    rc = 79; 
1545  } 
1546
1547  return 0;
1548
1549}
Note: See TracBrowser for help on using the repository browser.