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

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

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