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

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

continue w/wo calibration stuff (jec)
File size: 38.6 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//-----------------------------------------------
38const sa_size_t NUMBER_OF_CHANNELS = 2;
39const sa_size_t NUMBER_OF_FREQ     = 8192;
40const r_4    LOWER_FREQUENCY       = 1250.0; //MHz
41const r_4    TOTAL_BANDWIDTH       = 250.0;  //MHz
42//-----------------------------------------------
43//Input parameters
44struct Param {
45  int debuglev_;      //debug
46  string inPath_;     //root directory of the input files
47  string outPath_;    //output files are located here
48  string sourceName_; //source name & subdirectory of the input files
49  string ppfFile_;    //generic name of the input files
50  int nSliceInFreq_;  //used by reduceSpectra() fnc
51  string calibFreq_;  //freq. value used for calibration
52  r_4 rcalibFreq_;    //float version
53  string calibBandFreq_;  //band of freq. used for calibration
54  r_4 rcalibBandFreq_;    //float version
55  int maxNumberCycles_;//maximum number of cycles to be treated
56} para;
57//--------------------------------------------------------------
58//Utility functions
59
60sa_size_t freqToChan(r_4 f){
61  return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
62}
63//--------
64//COMPUTE the mean value on a freq. range for all channels
65//--------
66void meanInRange(const TMatrix<r_4> mtx, 
67                 sa_size_t chLow,
68                 sa_size_t chHigh, 
69                 TVector<r_4>& vec){
70  sa_size_t nr = mtx.NRows();
71  for (sa_size_t ir=0; ir<nr; ir++){
72    TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
73    double mean, sigma;
74    MeanSigma(tmp,mean,sigma);
75    vec(ir) = mean;
76  }
77}
78//---------
79class median_of_empty_list_exception:public std::exception{
80    virtual const char* what() const throw() {
81    return "Attempt to take the median of an empty list of numbers.  "
82      "The median of an empty list is undefined.";
83  }
84  };
85  template<class RandAccessIter>
86    double median(RandAccessIter begin, RandAccessIter end) 
87    throw(median_of_empty_list_exception){
88    if(begin == end){ throw median_of_empty_list_exception(); }
89    std::size_t size = end - begin;
90    std::size_t middleIdx = size/2;
91    RandAccessIter target = begin + middleIdx;
92    std::nth_element(begin, target, end);
93   
94    if(size % 2 != 0){ //Odd number of elements
95      return *target;
96    }else{            //Even number of elements
97      double a = *target;
98      RandAccessIter targetNeighbor= target-1;
99      std::nth_element(begin, targetNeighbor, end);
100      return (a+*targetNeighbor)/2.0;
101    }
102  }
103
104//-------------
105void split(const string& str, const string& delimiters , vector<string>& tokens) {
106    // Skip delimiters at beginning.
107    string::size_type lastPos = str.find_first_not_of(delimiters, 0);
108    // Find first "non-delimiter".
109    string::size_type pos     = str.find_first_of(delimiters, lastPos);
110
111    while (string::npos != pos || string::npos != lastPos)
112    {
113        // Found a token, add it to the vector.
114        tokens.push_back(str.substr(lastPos, pos - lastPos));
115        // Skip delimiters.  Note the "not_of"
116        lastPos = str.find_first_not_of(delimiters, pos);
117        // Find next "non-delimiter"
118        pos = str.find_first_of(delimiters, lastPos);
119    }
120}
121//--------------------------------------------------------------
122char *regexp (const char *string, const char *patrn, int *begin, int *end) {   
123        int i, w=0, len;                 
124        char *word = NULL;
125        regex_t rgT;
126        regmatch_t match;
127        regcomp(&rgT,patrn,REG_EXTENDED);
128        if ((regexec(&rgT,string,1,&match,0)) == 0) {
129                *begin = (int)match.rm_so;
130                *end = (int)match.rm_eo;
131                len = *end-*begin;
132                word=(char*)malloc(len+1);
133                for (i=*begin; i<*end; i++) {
134                        word[w] = string[i];
135                        w++; }
136                word[w]=0;
137        }
138        regfree(&rgT);
139        return word;
140}
141//-------
142sa_size_t round_sa(r_4 r) {
143  return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
144}
145//-----
146string StringToLower(string strToConvert){
147  //change each element of the string to lower case
148  for(unsigned int i=0;i<strToConvert.length();i++) {
149    strToConvert[i] = tolower(strToConvert[i]);
150  }
151  return strToConvert;//return the converted string
152}
153//-----
154bool stringCompare( const string &left, const string &right ){
155   if( left.size() < right.size() )
156      return true;
157   for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
158      if( tolower( *lit ) < tolower( *rit ) )
159         return true;
160      else if( tolower( *lit ) > tolower( *rit ) )
161         return false;
162   return false;
163}
164//-----
165list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
166  list<string> theList;
167
168
169  DIR *dip;
170  struct dirent *dit;
171  string msg;  string fileName;
172  string fullFileName;
173  size_t found;
174
175  if ((dip=opendir(dir.c_str())) == NULL ) {
176    msg = "opendir failed on directory "+dir;
177    throw msg;
178  }
179  while ( (dit = readdir(dip)) != NULL ) {
180    fileName = dit->d_name;
181    found=fileName.find(filePettern);
182    if (found != string::npos) {
183      fullFileName = dir + "/";
184      fullFileName += fileName;
185      theList.push_back(fullFileName);
186    }
187  }//eo while
188  if (closedir(dip) == -1) {
189    msg = "closedir failed on directory "+dir;
190    throw msg;
191  }
192 
193  theList.sort(stringCompare);
194
195  return theList;
196
197}
198//
199class  StringMatch : public unary_function<string,bool> {
200public:
201  StringMatch(const string& pattern): pattern_(pattern){}
202  bool operator()(const string& aStr) const {
203
204
205    int b,e;
206    regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
207
208//     cout << "investigate " << aStr << " to find " << pattern_
209//       << "[" <<b<<","<<e<<"]"
210//       << endl;
211
212   
213    if (b != 0) return false;
214    if (e != (int)aStr.size()) return false;
215    return true;
216
217  }
218private:
219  string pattern_;
220};
221//-------------------------------------------------------
222//Rebin in frequence + compute mean and sigma
223void reduceSpectra(const TMatrix<r_4>& specMtxInPut, 
224                    TMatrix<r_4>& meanMtx, 
225                    TMatrix<r_4>& sigmaMtx) {
226  sa_size_t nSliceFreq = para.nSliceInFreq_;
227  sa_size_t deltaFreq =  NUMBER_OF_FREQ/nSliceFreq;
228  for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
229    sa_size_t freqLow= iSlice*deltaFreq;
230    sa_size_t freqHigh= freqLow + deltaFreq -1;
231    for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
232      TVector<r_4> reducedRow;
233      reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
234      double mean; 
235      double sigma;
236      MeanSigma(reducedRow,mean,sigma);
237      meanMtx(iCh,iSlice) = mean;
238      sigmaMtx(iCh,iSlice) = sigma;
239    }//eo loop on channels
240  }//eo loop on slices
241}
242//-------------------------------------------------------
243//Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra
244//Used like:
245//
246// void meanCalibBAODiffOnOffCycles() throw(string) {
247
248//   list<string> listOfFiles;
249//   string directoryName;
250//   directoryName = para.inPath_ + "/" + para.sourceName_;
251
252//   //Make the listing of the directory
253//   listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
254 
255//   list<string>::const_iterator iFile, iFileEnd, iSpecOff, iSpecOffEnd, iSpecOn, iSpecOnEnd;
256//   iFileEnd = listOfFiles.end();
257 
258
259//   //Loop on files
260//   uint_4 nRuns=0;
261//   TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
262
263
264//   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
265//     if (para.debuglev_>90){
266//       cout << "load file <" << *iFile << ">" << endl;
267//     }
268
269//     vector<string> tokens;
270//     split(*File,"_",tokens);
271//     string dateOfRun = tokens[1];
272//     string srcLower = tokens[2];
273//     if (para.debuglev_>90){
274//       cout << "date <" << dateOfRun << ">" << endl;
275//     }
276   
277//     PInPersist fin(*iFile);
278//     vector<string> vec = fin.GetNameTags();
279
280//     if (para.typeOfCalib_ == "perRun") {
281//       ///////////////////
282//       //make the calibration of the mean of all Off and On of the run and perform the difference
283
284//       vector<string> modeList;
285//       modeList.push_back("On");
286//       modeList.push_back("Off");
287//       vector<string>::const_iterator iMode;
288//       map<string, TMatrix<r_4> > spectreModeCollect;
289
290//       for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
291//      ///////////////////
292//      //
293//      //Compute the mean of the mode
294//      //
295//      list<string> listOfSpectra;
296//      //Keep only required PPF objects
297//      string matchstr = "specRaw"+(*iMode)+"[0-9]+";
298//      std::remove_copy_if(
299//                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
300//                          not1(StringMatch(matchstr))
301//                          );
302       
303//      listOfSpectra.sort(stringCompare);
304//      iSpecEnd = listOfSpectra.end();
305//      TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
306//      uint_4 nSpectra=0;
307//      //Loop of spectra matrix
308//      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
309//        if (para.debuglev_>90){
310//          cout << " spactra <" << *iSpec << ">" << endl;
311//        }
312//        TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
313//        fin.GetObject(aSpec,*iSpec);
314//        //How to see if the GetObject is ok?? Ask Reza
315//        nSpectra++;
316//        meanOfSpectra+=aSpec;
317//      }//end loop Off
318//      //Mean
319//      if(nSpectra>0)meanOfSpectra=(r_4)(nSpectra);
320
321//      //BAO Calibrator
322//      string calibFileName = directoryName + "/"
323//        + "calib_" + dateOfRun + "_" + srcLower + "_"+(*iMode)+"_"
324//        + para.calibFreq_ +"MHz-All.txt";
325//      if(debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
326//      ifstream ifs(calibFileName.c_str());
327//      if ( ! ifs.is_open() ) {
328//        rc = 999;
329//        throw calibFileName + " cannot be opened...";
330//      }       
331//      TVector<r_4> calibBAOfactors;
332//      sa_size_t nr,nc; //values read
333//      calibBAOfactorsOff.ReadASCII(ifs,nr,nc);
334//      if(debuglev_>9){
335//        cout << "(nr,nc): "<< nr << "," << nc << endl;
336//        calibBAOfactors.Print(cout);
337//      }
338       
339//      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
340//        meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);
341//      }
342//      //storage
343//      spectreModeCollect.insert(pair<string, TMatrix<r_4> >(*iMode,TMatrix<r_4>(meanOfSpectra,false))); //do not share data (cf. SOPHYA)
344//       }//end of mode
345
346//       //Take the difference ON-OFF in current run
347//       TMatrix<r_4>diffOnOff = spectreModeCollect["On"]-spectreModeCollect["Off"];
348
349
350//     } else if (para.typeOfCalib_ == "perCycle") {
351//       //perform the calibration of the OFF and ON per cycle, then make the mean and take the diff
352//     } else {
353//       string msg="FATAL (meanCalibBAODiffOnOffCycles); unknown calibration mode "
354//      + para.typeOfCalib_ ;
355//       throw(msg);
356//     } 
357
358
359//     nRuns++;
360   
361//   }//eo loop on spectra in a file
362// }
363void meanCalibBAODiffOnOffCycles() throw(string) {
364
365  list<string> listOfFiles;
366  string directoryName;
367  directoryName = para.inPath_ + "/" + para.sourceName_;
368
369  //Make the listing of the directory
370  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
371 
372  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
373  iFileEnd = listOfFiles.end();
374
375  //mean ON-OFF over the list of cycles
376  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
377  TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
378  TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
379  char* onffTupleName[7]={"cycle",
380                          "onoffRaw0","onoffRun0","onoffCycle0",
381                          "onoffRaw1","onoffRun1","onoffCycle1"};
382  NTuple onoffevolution(7,onffTupleName);
383  r_4 xnt[7];
384
385  //Lower and Higher freq. bin to perform mean follow up
386  sa_size_t chLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
387  sa_size_t chHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
388 
389  //Loop on files/run
390
391  int totalNumberCycles=0; //total number of cycles for normalisation
392  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
393    if (para.debuglev_>90){
394      cout << "load file <" << *iFile << ">" << endl;
395    }
396
397    vector<string> tokens;
398    split(*iFile,"_",tokens);
399    string dateOfRun = tokens[1];
400    if (para.debuglev_>90){
401      cout << "date <" << dateOfRun << ">" << endl;
402    }
403    vector<string> tokens2;
404    split(tokens[2],".",tokens2);
405    string srcLower = tokens2[0];
406
407
408   
409    PInPersist fin(*iFile);
410    vector<string> vec = fin.GetNameTags();
411
412    vector<string> modeList;
413    modeList.push_back("On");
414    modeList.push_back("Off");
415    vector<string>::const_iterator iMode;
416   
417    map<string, list<int> > cycleModeCollect;
418   
419    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
420      list<string> listOfSpectra;
421      //Keep only required PPF objects
422      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
423      std::remove_copy_if(
424                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
425                          not1(StringMatch(matchstr))
426                          );
427     
428      listOfSpectra.sort(stringCompare);
429      iSpecEnd = listOfSpectra.end();
430     
431      matchstr = "[0-9]+";
432      //Loop of spectra matrix
433      list<int> listOfCycles;
434      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
435        int b,e;
436        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
437        if (para.debuglev_>90){
438          cout << " spactra <" << *iSpec << ">" << endl;
439          cout << " cycle " << iSpec->substr(b) << endl;
440        }
441        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
442      }//end loop spectra
443      cycleModeCollect[*iMode] = listOfCycles;
444    }//end of mode   
445
446    //Take the Intersection of the list Of cycles in mode Off and On   
447    list<int> commonCycles;
448    set_intersection(cycleModeCollect["On"].begin(),
449                     cycleModeCollect["On"].end(),
450                     cycleModeCollect["Off"].begin(),
451                     cycleModeCollect["Off"].end(),
452                     back_inserter(commonCycles)
453                     );
454   
455    if (para.debuglev_>90){
456      cout << "Liste of cycles common to On & Off: <";
457      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
458        cout << *i << " ";
459      }
460      cout << ">" << endl;
461    }
462   
463    //
464    //Load BAO Calibration factors "per Cycle and Channels"
465    //Compute the mean per Cycle to
466    //  fill factors "per Run and Channels" with the same cycle list
467    //
468    //
469    //TODO improve the code....
470
471    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
472    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
473    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
474    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
475   
476    string calibFileName;
477    ifstream ifs;
478    sa_size_t nr,nc; //values read
479
480    //OFF Cycle per Channel
481    calibFileName = directoryName + "/" 
482      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
483      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
484    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
485    ifs.open(calibFileName.c_str());
486    if ( ! ifs.is_open() ) {
487
488      throw calibFileName + " cannot be opened...";
489    }   
490    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
491    if(para.debuglev_>9){
492      cout << "(nr,nc): "<< nr << "," << nc << endl;
493      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
494      cout << endl;
495    }
496
497    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
498    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
499    {//Compute the mean
500      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
501      double mean,sigma;
502      MeanSigma(coef,mean,sigma);
503      cout << "Mean: " << mean << " sigma " << sigma << endl;
504      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
505    }
506    if(para.debuglev_>9){
507      cout << "Fill calib. with mean value " << endl; 
508      calibBAOfactors_Off_Run_Ch0.Print(cout);
509      cout << endl;
510    }
511    ifs.close();
512
513    //
514    calibFileName = directoryName + "/" 
515      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
516      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
517    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
518    ifs.open(calibFileName.c_str());
519    if ( ! ifs.is_open() ) {
520
521      throw calibFileName + " cannot be opened...";
522    }   
523    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
524    if(para.debuglev_>9){
525      cout << "(nr,nc): "<< nr << "," << nc << endl;
526      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
527      cout << endl;
528    }
529    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
530    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
531    {//Compute the mean
532      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
533      double mean,sigma;
534      MeanSigma(coef,mean,sigma);
535      cout << "Mean: " << mean << " sigma " << sigma << endl;
536      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
537    }
538    if(para.debuglev_>9){
539      cout << "Fill calib. with mean value " << endl; 
540      calibBAOfactors_Off_Run_Ch1.Print(cout);
541      cout << endl;
542    }
543    ifs.close();
544
545    //ON Cycle per Channel
546    calibFileName = directoryName + "/" 
547      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
548      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
549    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
550    ifs.open(calibFileName.c_str());
551    if ( ! ifs.is_open() ) {
552
553      throw calibFileName + " cannot be opened...";
554    }   
555    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
556    if(para.debuglev_>9){
557      cout << "(nr,nc): "<< nr << "," << nc << endl;
558      calibBAOfactors_On_Cycle_Ch0.Print(cout);
559      cout << endl;     
560    }
561
562    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
563    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
564    {//Compute the mean
565      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
566      double mean,sigma;
567      MeanSigma(coef,mean,sigma);
568      cout << "Mean: " << mean << " sigma " << sigma << endl;
569      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
570    }
571    if(para.debuglev_>9){
572      cout << "Fill calib. with mean value " << endl; 
573      calibBAOfactors_On_Run_Ch0.Print(cout);
574      cout << endl;
575    }
576    ifs.close();
577
578   
579    calibFileName = directoryName + "/" 
580      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
581      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
582    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
583    ifs.open(calibFileName.c_str());
584    if ( ! ifs.is_open() ) {
585      throw calibFileName + " cannot be opened...";
586    }   
587    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
588    if(para.debuglev_>9){
589      cout << "(nr,nc): "<< nr << "," << nc << endl;
590      calibBAOfactors_On_Cycle_Ch1.Print(cout);
591      cout << endl;
592    }
593    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
594    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
595    {//Compute the mean
596      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
597      double mean,sigma;
598      MeanSigma(coef,mean,sigma);
599      cout << "Mean: " << mean << " sigma " << sigma << endl;
600      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
601    }
602    if(para.debuglev_>9){
603      cout << "Fill calib. with mean value " << endl; 
604      calibBAOfactors_On_Run_Ch1.Print(cout);
605      cout << endl;
606    }
607
608    ifs.close();
609   
610    //link <cycle> - <calibration coefficient>
611    //We cannot rely on identical cycle list of the OFF and ON calibration
612    map<int,r_4> calibBAO_Off_Run_Ch0;
613    map<int,r_4> calibBAO_Off_Run_Ch1;
614    map<int,r_4> calibBAO_On_Run_Ch0;
615    map<int,r_4> calibBAO_On_Run_Ch1;
616
617    map<int,r_4> calibBAO_Off_Cycle_Ch0;
618    map<int,r_4> calibBAO_Off_Cycle_Ch1;
619    map<int,r_4> calibBAO_On_Cycle_Ch0;
620    map<int,r_4> calibBAO_On_Cycle_Ch1;
621
622    //per Run based BAO coefficients
623    nr = calibBAOfactors_Off_Run_Ch0.NRows();
624    for (sa_size_t ir=0; ir<nr; ++ir){
625      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
626        = calibBAOfactors_Off_Run_Ch0(ir,1);
627      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
628        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
629      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
630        = calibBAOfactors_Off_Run_Ch1(ir,1);
631      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
632        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
633    }
634   
635    nr = calibBAOfactors_On_Run_Ch0.NRows();
636    for (sa_size_t ir=0; ir<nr; ++ir){
637      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
638        = calibBAOfactors_On_Run_Ch0(ir,1);
639      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
640        = calibBAOfactors_On_Cycle_Ch0(ir,1);
641      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
642        = calibBAOfactors_On_Run_Ch1(ir,1);
643      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
644        = calibBAOfactors_On_Cycle_Ch1(ir,1);
645    }
646     
647
648   
649//      for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
650//        meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);
651//      }
652
653    //Loop on cyles
654    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
655
656      if(para.debuglev_>9){
657        cout << "Calibration coefficients for cycle "<<*ic << endl; 
658        cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
659             << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
660             << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
661             << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
662             << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
663             << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
664             << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
665             << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
666      }
667     
668      string ppftag;
669      //load ON phase
670      stringstream cycle;
671      cycle << *ic;
672     
673      ppftag = "specRawOn"+cycle.str();
674      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
675      fin.GetObject(aSpecOn,ppftag);
676
677      TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
678      aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
679      aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
680
681      TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
682      aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
683      aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
684     
685      ppftag = "specRawOff"+cycle.str();
686      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
687      fin.GetObject(aSpecOff,ppftag);
688
689      TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
690      aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
691      aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
692
693      TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
694      aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
695      aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
696
697
698      //Perform the difference ON-OFF with the different calibration options
699      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
700      meanDiffONOFF_noCalib += diffOnOff_noCalib;
701     
702      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
703      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
704
705      TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
706      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
707
708      //Fill NTuple
709      xnt[0] = *ic;
710     
711      TVector<r_4> meanInRange_noCalib(NUMBER_OF_CHANNELS);
712      meanInRange(diffOnOff_noCalib,chLow,chHigh,meanInRange_noCalib);
713      xnt[1] = meanInRange_noCalib(0);
714      xnt[2] = meanInRange_noCalib(1);
715
716      TVector<r_4> meanInRange_perRunCalib(NUMBER_OF_CHANNELS);
717      meanInRange(diffOnOff_perRunCalib,chLow,chHigh,meanInRange_perRunCalib);
718      xnt[3] = meanInRange_perRunCalib(0);
719      xnt[4] = meanInRange_perRunCalib(1);
720
721      TVector<r_4> meanInRange_perCycleCalib(NUMBER_OF_CHANNELS);
722      meanInRange(diffOnOff_perCycleCalib,chLow,chHigh,meanInRange_perCycleCalib);
723      xnt[3] = meanInRange_perCycleCalib(0);
724      xnt[4] = meanInRange_perCycleCalib(1);
725     
726      onoffevolution.Fill(xnt);
727
728      totalNumberCycles++;
729      if (totalNumberCycles >= para.maxNumberCycles_) break;   
730
731    }//eo loop on cycles
732    if (totalNumberCycles >= para.maxNumberCycles_) break;         
733 
734  }//eo loop on spectra in a file
735  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
736  //Normalisation
737  if(totalNumberCycles>0){
738    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
739    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
740    meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
741  } 
742 
743  //Compute the reduced version of the mean and sigma
744  TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
745  TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
746  reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
747
748  TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
749  TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
750  reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
751
752  TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
753  TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
754  reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
755
756  {//save the results
757    stringstream tmp;
758    tmp << totalNumberCycles;
759    string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
760    POutPersist fos(fileName);
761
762    string tag = "meanNoCalib";
763    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
764    tag = "meanPerRunCalib";
765    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
766    tag = "meanPerCycleCalib";
767    fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
768
769    tag = "redmeanNoCalib";
770    fos << PPFNameTag(tag) << meanRedMtx_noCalib;
771    tag = "redsigmaNoCalib";
772    fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
773
774    tag = "redmeanPerRunCalib";
775    fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
776    tag = "redsigmaPerRunCalib";
777    fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
778
779    tag = "redmeanPerCycleCalib";
780    fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
781    tag = "redsigmaPerCycleCalib";
782    fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
783   
784    tag = "onoffevol";
785    fos << PPFNameTag(tag) << onoffevolution;   
786  }//end of save
787}
788//-------------------------------------------------------
789//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
790//Used like:
791//
792void meanRawDiffOnOffCycles() throw(string) {
793  list<string> listOfFiles;
794  string directoryName;
795  directoryName = para.inPath_ + "/" + para.sourceName_;
796
797  //Make the listing of the directory
798  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
799 
800  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
801  iFileEnd = listOfFiles.end();
802 
803  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
804  TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
805  uint_4 nSpectra=0;
806  //Loop on files
807  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
808    if (para.debuglev_>90){
809      cout << "load file <" << *iFile << ">" << endl;
810    }
811    PInPersist fin(*iFile);
812    vector<string> vec = fin.GetNameTags();
813    list<string> listOfSpectra;
814    //Keep only required PPF objects
815    std::remove_copy_if(
816                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
817                        not1(match)
818                        );
819   
820    listOfSpectra.sort(stringCompare);
821    iSpecEnd = listOfSpectra.end();
822    //Loop of spectra matrix
823    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
824      if (para.debuglev_>90){
825        cout << " spactra <" << *iSpec << ">" << endl;
826      }
827      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
828      fin.GetObject(aSpec,*iSpec);
829      //How to see if the GetObject is ok?? Ask Reza
830      nSpectra++;
831      meanOfSpectra+=aSpec;
832    }//eo loop on spectra in a file
833  }//eo loop on files
834 
835  //Normalisation
836  if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
837
838  //Compute the reduced version of the mean and sigma
839  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
840  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
841  reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
842
843  {//Save the result
844    stringstream tmp;
845    tmp << nSpectra;
846    string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
847    cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
848    POutPersist fos(fileName);
849
850    string tag = "mean";
851    fos << PPFNameTag(tag) << meanOfSpectra;
852    tag = "meanred";
853    fos << PPFNameTag(tag) << meanRedMtx;
854    tag = "sigmared";
855    fos << PPFNameTag(tag) << sigmaRedMtx;
856  }
857}
858//-------------------------------------------------------
859//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
860//Used like:
861//
862void medianRawDiffOnOffCycles() throw(string) {
863  list<string> listOfFiles;
864  string directoryName;
865  directoryName = para.inPath_ + "/" + para.sourceName_;
866
867  //Make the listing of the directory
868  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
869 
870  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
871  iFileEnd = listOfFiles.end();
872 
873
874  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
875
876  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
877  uint_4 nSpectra=0;
878  //Loop on files
879  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
880    if (para.debuglev_>90){
881      cout << "load file <" << *iFile << ">" << endl;
882    }
883    PInPersist fin(*iFile);
884    vector<string> vec = fin.GetNameTags();
885    list<string> listOfSpectra;
886    //Keep only required PPF objects
887    std::remove_copy_if(
888                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
889                        not1(match)
890                        );
891   
892    listOfSpectra.sort(stringCompare);
893    iSpecEnd = listOfSpectra.end();
894    //Loop of spectra matrix
895    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ;  ++iSpec){
896      if (para.debuglev_>90){
897        cout << " spactra <" << *iSpec << ">" << endl;
898      }
899      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
900      fin.GetObject(aSpec,*iSpec);
901
902      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
903
904      nSpectra++;
905    }//eo loop on spectra in a file
906  }//eo loop on files
907 
908
909 
910  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
911  //Compute the median for each freq. and Channel
912  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
913    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
914      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
915      vector<r_4> tmp;
916      tmp0.FillTo(tmp);
917      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
918    }
919  }
920
921
922  //Compute the reduced version of the mean and sigma
923  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
924  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
925  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
926
927
928  sa_size_t f1320=freqToChan(1320.);
929  sa_size_t f1345=freqToChan(1345.);
930  sa_size_t f1355=freqToChan(1355.);
931  sa_size_t f1380=freqToChan(1380.);
932  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
933  if (para.debuglev_>9){
934    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
935  }
936  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
937  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
938    double meanMed1;
939    double sigmaMed1;
940    TVector<r_4> band1;
941    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
942    MeanSigma(band1,meanMed1,sigmaMed1);
943    double meanMed2;
944    double sigmaMed2;
945    TVector<r_4> band2;
946    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
947    MeanSigma(band2,meanMed2,sigmaMed2);
948    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
949  } 
950  meanMed.Print(cout);
951  cout << endl;
952
953 
954  //Compute the sigma in the range 1320MHz-1380MHz
955  if (para.debuglev_>9){
956    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
957  }
958  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
959  sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
960  sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
961
962  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
963    double meanSigma;
964    double sigmaSigma;
965    TVector<r_4> band;
966    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
967    MeanSigma(band,meanSigma,sigmaSigma);
968    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
969    sigmaMed(iCh) = meanSigma;
970  }
971  sigmaMed.Print(cout);
972  cout << endl;
973
974 
975 
976  if (para.debuglev_>9){
977    cout << "Compute medianOfSpectraNorm" << endl;
978  }
979  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
980  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
981    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
982    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
983  }
984
985 
986
987  {//Save the result
988    stringstream tmp;
989    tmp << nSpectra;
990    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
991    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
992    POutPersist fos(fileName);
993
994    string tag = "median";
995    fos << PPFNameTag(tag) << medianOfSpectra;
996
997    tag = "medianNorm";
998    fos << PPFNameTag(tag) << medianOfSpectraNorm;
999   
1000
1001    tag = "meanmedred";
1002    fos << PPFNameTag(tag) << meanRedMtx;
1003    tag = "sigmamedred";
1004    fos << PPFNameTag(tag) << sigmaRedMtx;
1005    tag = "cycleVsfreq";
1006   
1007    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1008    fos << PPFNameTag(tag) << tarr;
1009  }
1010}
1011
1012//-------------------------------------------------------
1013int main(int narg, char* arg[]) {
1014 
1015  int rc = 0; //return code
1016  string msg; //message used in Exceptions
1017
1018
1019
1020  //default value for initial parameters (see Para structure on top of the file)
1021  string debuglev = "0";
1022  string action = "meanDiffOnOff";
1023  string inputPath = "."; 
1024  string outputPath = "."; 
1025  string sourceName = "Abell85";
1026  string ppfFile;
1027  string nSliceInFreq = "32";
1028  string typeOfCalib="perRun";
1029  string calibFreq = "1346";
1030  string calibBandFreq="6.25";
1031  string mxcycles;
1032
1033  //  bool okarg=false;
1034  int ka=1;
1035  while (ka<narg) {
1036    if (strcmp(arg[ka],"-h")==0) {
1037      cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1038           << " -mxcycles <number> (max. number of cycles to be treated)\n"
1039           << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1040           << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1041           << " -src [Abell85]\n -inPath [.]|<top_root_dir of the ppf file>\n" 
1042           << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
1043           << " -outPath [.]|<dir of the output> \n"
1044           << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1045           << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1046           << " -debug <level> "
1047           << endl;
1048      return 0;
1049    }
1050    else if (strcmp(arg[ka],"-debug")==0) {
1051      debuglev=arg[ka+1];
1052      ka+=2;
1053    }
1054    else if (strcmp(arg[ka],"-act")==0) {
1055      action=arg[ka+1];
1056      ka+=2;
1057    }
1058    else if (strcmp(arg[ka],"-calibfreq")==0) {
1059      calibFreq=arg[ka+1];
1060      ka+=2;
1061    }   
1062    else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1063      calibBandFreq=arg[ka+1];
1064      ka+=2;
1065    }   
1066    else if (strcmp(arg[ka],"-mxcycles")==0) {
1067      mxcycles=arg[ka+1];
1068      ka+=2;
1069    }   
1070    else if (strcmp(arg[ka],"-inPath")==0) {
1071      inputPath=arg[ka+1];
1072      ka+=2;
1073    }
1074    else if (strcmp(arg[ka],"-outPath")==0) {
1075      outputPath=arg[ka+1];
1076      ka+=2;
1077    }
1078    else if (strcmp(arg[ka],"-src")==0) {
1079      sourceName=arg[ka+1];
1080      ka+=2;
1081    }
1082    else if (strcmp(arg[ka],"-ppfFile")==0) {
1083      ppfFile=arg[ka+1];
1084      ka+=2;
1085    }
1086    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1087      nSliceInFreq=arg[ka+1];
1088      ka+=2;
1089    }
1090    else ka++;
1091  }//eo while
1092
1093  para.debuglev_   = atoi(debuglev.c_str());
1094  para.inPath_     = inputPath;
1095  para.outPath_    = outputPath;
1096  para.sourceName_ = sourceName;
1097  para.ppfFile_    = ppfFile;
1098  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1099  para.calibFreq_   = calibFreq;
1100  para.calibBandFreq_ = calibBandFreq;
1101  para.rcalibFreq_   = atof(calibFreq.c_str());
1102  para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1103  if (mxcycles != "") {
1104    para.maxNumberCycles_ = atoi(mxcycles.c_str());
1105  } else {
1106    para.maxNumberCycles_ = std::numeric_limits<int>::max();
1107  }
1108
1109  cout << "Dump Initial parameters ............" << endl;
1110  cout << " action = " << action << "\n"
1111       << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1112       << " inputPath = " << para.inPath_  << "\n" 
1113       << " outputPath = " <<  para.outPath_ << "\n"
1114       << " sourceName = " << para.sourceName_ << "\n"
1115       << " ppfFile = " <<  para.ppfFile_ << "\n"
1116       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
1117       << " calibFreq = " <<  para.calibFreq_ << "\n"
1118       << " calibBandFreq = " <<  para.calibBandFreq_ << "\n"
1119       << " debuglev = "  << para.debuglev_   << "\n";
1120  cout << "...................................." << endl;
1121
1122  if ( "" == ppfFile ) {
1123    cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1124         << endl;
1125    return 999;
1126  }
1127
1128
1129  try {
1130
1131//     int b,e;
1132//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1133//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1134
1135    if ( action == "meanRawDiffOnOff" ) {
1136      meanRawDiffOnOffCycles();
1137    } else if (action == "medianRawDiffOnOff") {
1138      medianRawDiffOnOffCycles();
1139    } else if (action == "meanCalibBAODiffOnOff") {
1140      meanCalibBAODiffOnOffCycles();
1141    } else {
1142      msg = "Unknown action " + action;
1143      throw msg;
1144    }
1145
1146
1147  }  catch (std::exception& sex) {
1148    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
1149         << "\n msg= " << sex.what() << endl;
1150    rc = 78;
1151  }
1152  catch ( string str ) {
1153    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1154  }
1155  catch (...) {
1156    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
1157    rc = 79; 
1158  } 
1159
1160  return 0;
1161
1162}
Note: See TracBrowser for help on using the repository browser.