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

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

Added calibcoeffNtp function (this is the good one)
File size: 72.8 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//AST 18/11/11
318//-------------------------------------------------------
319//Make n-tuple with calibration coefficients
320//
321void calibCoeffNtp() throw(string) {
322  list<string> listOfFiles;
323  string directoryName;
324  directoryName = para.inPath_ + "/" + para.sourceName_;
325
326  //Make the listing of the directory
327  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
328 
329  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
330  iFileEnd = listOfFiles.end();
331
332  static const int NINFO=19;
333  char* calibTupleName[NINFO]={"date","run","cycle"
334                              ,"meancoeffoff0","meancoeffoff1"
335                              ,"meancoeffon0","meancoeffon1"
336                              ,"coeffoff0","coeffoff1"                             
337                              ,"coeffon0","coeffon1"
338                              ,"invmeancoeffoff0","invmeancoeffoff1"
339                              ,"invmeancoeffon0","invmeancoeffon1"
340                              ,"invcoeffoff0","invcoeffoff1"                             
341                              ,"invcoeffon0","invcoeffon1"
342  };
343  NTuple calibcoeffs(NINFO,calibTupleName); 
344  r_4 xnt[NINFO];
345 
346  int totalNumberRuns=0;   //total number of runs
347  int totalNumberCycles=0; //total number of cycles for normalisation
348  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
349    if (para.debuglev_>90){
350      cout << "load file <" << *iFile << ">" << endl;
351    }
352
353    vector<string> tokens;
354    split(*iFile,"_",tokens);
355    string dateOfRun = tokens[1];
356    r_4 rdateOfRun = atoi(dateOfRun.c_str()) - 2011*1.e4;
357
358    if (para.debuglev_>90){
359      cout << "date <" << dateOfRun << ">" << endl;
360      cout << "rdate <" << rdateOfRun << ">" << endl;
361    }
362    vector<string> tokens2;
363    split(tokens[2],".",tokens2);
364    string srcLower = tokens2[0];
365   
366    PInPersist fin(*iFile);
367    vector<string> vec = fin.GetNameTags();
368
369    vector<string> modeList;
370    modeList.push_back("On");
371    modeList.push_back("Off");
372    vector<string>::const_iterator iMode;
373   
374    map<string, list<int> > cycleModeCollect;
375   
376    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
377      list<string> listOfSpectra;
378      //Keep only required PPF objects
379      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
380      std::remove_copy_if(
381                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
382                          not1(StringMatch(matchstr))
383                          );
384     
385      listOfSpectra.sort(stringCompare);
386      iSpecEnd = listOfSpectra.end();
387     
388      matchstr = "[0-9]+";
389      //Loop of spectra matrix
390      list<int> listOfCycles;
391      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
392        int b,e;
393        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
394        if (para.debuglev_>90){
395          cout << " spectra <" << *iSpec << ">" << endl;
396          cout << " cycle " << iSpec->substr(b) << endl;
397        }
398        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
399      }//end loop spectra
400      cycleModeCollect[*iMode] = listOfCycles;
401    }//end of mode   
402
403    //Take the Intersection of the list Of cycles in mode Off and On   
404    list<int> commonCycles;
405    set_intersection(cycleModeCollect["On"].begin(),
406                     cycleModeCollect["On"].end(),
407                     cycleModeCollect["Off"].begin(),
408                     cycleModeCollect["Off"].end(),
409                     back_inserter(commonCycles)
410                     );
411   
412    if (para.debuglev_>90){
413      cout << "Liste of cycles common to On & Off: <";
414      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
415        cout << *i << " ";
416      }
417      cout << ">" << endl;
418    }
419   
420    //
421    //Load BAO Calibration factors "per Cycle and Channels"
422    //Compute the mean per Cycle to
423    //  fill factors "per Run and Channels" with the same cycle list
424    //
425    //
426    //TODO improve the code....
427
428    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
429    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
430    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
431    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
432   
433    string calibFileName;
434    ifstream ifs;
435    sa_size_t nr,nc; //values read
436
437    //OFF Cycle per Channel
438    calibFileName = directoryName + "/" 
439      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
440      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
441    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
442    ifs.open(calibFileName.c_str());
443    if ( ! ifs.is_open() ) {
444      throw calibFileName + " cannot be opened...";
445    }   
446    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
447    if(para.debuglev_>9){
448      cout << "(nr,nc): "<< nr << "," << nc << endl;
449      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
450      cout << endl;
451    }
452
453    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
454    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
455    {//Compute the mean
456      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
457      double mean,sigma;
458      MeanSigma(coef,mean,sigma);
459      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
460    }
461    if(para.debuglev_>9){
462      cout << "Fill calib. with mean value " << endl; 
463      calibBAOfactors_Off_Run_Ch0.Print(cout);
464      cout << endl;
465    }
466
467    TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch0(nr,nc);
468    invcalibBAOfactors_Off_Cycle_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
469    {//Compute the inverse value
470      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
471      invcalibBAOfactors_Off_Cycle_Ch0.Column(1).SetCst(1);
472      invcalibBAOfactors_Off_Cycle_Ch0.Column(1).Div(coef);
473      if(para.debuglev_>9){
474        cout << "Fill calib. with inverse value " << endl;
475        invcalibBAOfactors_Off_Cycle_Ch0.Print(cout);
476        cout << endl;
477      }
478    }
479    TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch0(nr,nc);
480    invcalibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
481    {//Compute the inverse mean
482      TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
483      double mean,sigma;
484      MeanSigma(coef,mean,sigma);
485      invcalibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
486    }
487    if(para.debuglev_>9){
488      cout << "Fill calib. with inverse mean value " << endl; 
489      invcalibBAOfactors_Off_Run_Ch0.Print(cout);
490      cout << endl;
491    }
492    ifs.close();
493    //
494    calibFileName = directoryName + "/" 
495      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
496      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
497    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
498    ifs.open(calibFileName.c_str());
499    if ( ! ifs.is_open() ) {
500
501      throw calibFileName + " cannot be opened...";
502    }   
503    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
504    if(para.debuglev_>9){
505      cout << "(nr,nc): "<< nr << "," << nc << endl;
506      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
507      cout << endl;
508    }
509    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
510    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
511    {//Compute the mean
512      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
513      double mean,sigma;
514      MeanSigma(coef,mean,sigma);
515      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
516      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
517    }
518    if(para.debuglev_>9){
519      cout << "Fill calib. with mean value " << endl; 
520      calibBAOfactors_Off_Run_Ch1.Print(cout);
521      cout << endl;
522    }
523    TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch1(nr,nc);
524    invcalibBAOfactors_Off_Cycle_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
525    {//Compute the inverse value
526      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
527      invcalibBAOfactors_Off_Cycle_Ch1.Column(1).SetCst(1);
528      invcalibBAOfactors_Off_Cycle_Ch1.Column(1).Div(coef);
529      if(para.debuglev_>9){
530        cout << "Fill calib. with inverse value " << endl;
531        invcalibBAOfactors_Off_Cycle_Ch1.Print(cout);
532        cout << endl;
533      }
534    }
535    TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch1(nr,nc);
536    invcalibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
537    {//Compute the inverse mean
538      TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
539      double mean,sigma;
540      MeanSigma(coef,mean,sigma);
541      invcalibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
542    }
543    if(para.debuglev_>9){
544      cout << "Fill calib. with inverse mean value " << endl; 
545      invcalibBAOfactors_Off_Run_Ch1.Print(cout);
546      cout << endl;
547    }
548    ifs.close();
549
550    //ON Cycle per Channel
551    calibFileName = directoryName + "/" 
552      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
553      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
554    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
555    ifs.open(calibFileName.c_str());
556    if ( ! ifs.is_open() ) {
557
558      throw calibFileName + " cannot be opened...";
559    }   
560    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
561    if(para.debuglev_>9){
562      cout << "(nr,nc): "<< nr << "," << nc << endl;
563      calibBAOfactors_On_Cycle_Ch0.Print(cout);
564      cout << endl;     
565    }
566
567    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
568    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
569    {//Compute the mean
570      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
571      double mean,sigma;
572      MeanSigma(coef,mean,sigma);
573      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
574      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
575    }
576    if(para.debuglev_>9){
577      cout << "Fill calib. with mean value " << endl; 
578      calibBAOfactors_On_Run_Ch0.Print(cout);
579      cout << endl;
580    }
581
582    TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch0(nr,nc);
583    invcalibBAOfactors_On_Cycle_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
584    {//Compute the inverse value
585      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
586      invcalibBAOfactors_On_Cycle_Ch0.Column(1).SetCst(1);
587      invcalibBAOfactors_On_Cycle_Ch0.Column(1).Div(coef);
588      if(para.debuglev_>9){
589        cout << "Fill calib. with inverse value " << endl;
590        invcalibBAOfactors_On_Cycle_Ch0.Print(cout);
591        cout << endl;
592      }
593    }
594    TMatrix<r_4> invcalibBAOfactors_On_Run_Ch0(nr,nc);
595    invcalibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
596    {//Compute the inverse mean
597      TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
598      double mean,sigma;
599      MeanSigma(coef,mean,sigma);
600      invcalibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
601    }
602    if(para.debuglev_>9){
603      cout << "Fill calib. with inverse mean value " << endl; 
604      invcalibBAOfactors_On_Run_Ch0.Print(cout);
605      cout << endl;
606    }
607    ifs.close();
608   
609    calibFileName = directoryName + "/" 
610      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
611      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
612    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
613    ifs.open(calibFileName.c_str());
614    if ( ! ifs.is_open() ) {
615      throw calibFileName + " cannot be opened...";
616    }   
617    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
618    if(para.debuglev_>9){
619      cout << "(nr,nc): "<< nr << "," << nc << endl;
620      calibBAOfactors_On_Cycle_Ch1.Print(cout);
621      cout << endl;
622    }
623    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
624    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
625    {//Compute the mean
626      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
627      double mean,sigma;
628      MeanSigma(coef,mean,sigma);
629      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
630      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
631    }
632    if(para.debuglev_>9){
633      cout << "Fill calib. with mean value " << endl; 
634      calibBAOfactors_On_Run_Ch1.Print(cout);
635      cout << endl;
636    }
637
638    TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch1(nr,nc);
639    invcalibBAOfactors_On_Cycle_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
640    {//Compute the inverse value
641      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
642      invcalibBAOfactors_On_Cycle_Ch1.Column(1).SetCst(1);
643      invcalibBAOfactors_On_Cycle_Ch1.Column(1).Div(coef);
644      if(para.debuglev_>9){
645        cout << "Fill calib. with inverse value " << endl;
646        invcalibBAOfactors_On_Cycle_Ch1.Print(cout);
647        cout << endl;
648      }
649    }
650    TMatrix<r_4> invcalibBAOfactors_On_Run_Ch1(nr,nc);
651    invcalibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
652    {//Compute the inverse mean
653      TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
654      double mean,sigma;
655      MeanSigma(coef,mean,sigma);
656      invcalibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
657    }
658    if(para.debuglev_>9){
659      cout << "Fill calib. with inverse mean value " << endl; 
660      invcalibBAOfactors_On_Run_Ch1.Print(cout);
661      cout << endl;
662    }
663    ifs.close();
664    //link <cycle> - <calibration coefficient>
665    //We cannot rely on identical cycle list of the OFF and ON calibration
666    map<int,r_4> calibBAO_Off_Run_Ch0;
667    map<int,r_4> calibBAO_Off_Run_Ch1;
668    map<int,r_4> calibBAO_On_Run_Ch0;
669    map<int,r_4> calibBAO_On_Run_Ch1;
670
671    map<int,r_4> calibBAO_Off_Cycle_Ch0;
672    map<int,r_4> calibBAO_Off_Cycle_Ch1;
673    map<int,r_4> calibBAO_On_Cycle_Ch0;
674    map<int,r_4> calibBAO_On_Cycle_Ch1;
675
676    map<int,r_4> invcalibBAO_Off_Run_Ch0;
677    map<int,r_4> invcalibBAO_Off_Run_Ch1;
678    map<int,r_4> invcalibBAO_On_Run_Ch0;
679    map<int,r_4> invcalibBAO_On_Run_Ch1;
680
681    map<int,r_4> invcalibBAO_Off_Cycle_Ch0;
682    map<int,r_4> invcalibBAO_Off_Cycle_Ch1;
683    map<int,r_4> invcalibBAO_On_Cycle_Ch0;
684    map<int,r_4> invcalibBAO_On_Cycle_Ch1;
685
686    //per Run based BAO coefficients
687    nr = calibBAOfactors_Off_Run_Ch0.NRows();
688    for (sa_size_t ir=0; ir<nr; ++ir){
689//       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
690//         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
691
692      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
693        = calibBAOfactors_Off_Run_Ch0(ir,1);
694      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
695        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
696      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
697        = calibBAOfactors_Off_Run_Ch1(ir,1);
698      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
699        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
700
701      invcalibBAO_Off_Run_Ch0[(int)invcalibBAOfactors_Off_Run_Ch0(ir,0)]
702        = invcalibBAOfactors_Off_Run_Ch0(ir,1);
703      invcalibBAO_Off_Cycle_Ch0[(int)invcalibBAOfactors_Off_Cycle_Ch0(ir,0)]
704        = invcalibBAOfactors_Off_Cycle_Ch0(ir,1);
705      invcalibBAO_Off_Run_Ch1[(int)invcalibBAOfactors_Off_Run_Ch1(ir,0)]
706        = invcalibBAOfactors_Off_Run_Ch1(ir,1);
707      invcalibBAO_Off_Cycle_Ch1[(int)invcalibBAOfactors_Off_Cycle_Ch1(ir,0)]
708        = invcalibBAOfactors_Off_Cycle_Ch1(ir,1);
709    }//eo loop on coef Off
710   
711    nr = calibBAOfactors_On_Run_Ch0.NRows();
712    for (sa_size_t ir=0; ir<nr; ++ir){
713      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
714        = calibBAOfactors_On_Run_Ch0(ir,1);
715      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
716        = calibBAOfactors_On_Cycle_Ch0(ir,1);
717      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
718        = calibBAOfactors_On_Run_Ch1(ir,1);
719      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
720        = calibBAOfactors_On_Cycle_Ch1(ir,1);
721
722      invcalibBAO_On_Run_Ch0[(int)invcalibBAOfactors_On_Run_Ch0(ir,0)]
723        = invcalibBAOfactors_On_Run_Ch0(ir,1);
724      invcalibBAO_On_Cycle_Ch0[(int)invcalibBAOfactors_On_Cycle_Ch0(ir,0)]
725        = invcalibBAOfactors_On_Cycle_Ch0(ir,1);
726      invcalibBAO_On_Run_Ch1[(int)invcalibBAOfactors_On_Run_Ch1(ir,0)]
727        = invcalibBAOfactors_On_Run_Ch1(ir,1);
728      invcalibBAO_On_Cycle_Ch1[(int)invcalibBAOfactors_On_Cycle_Ch1(ir,0)]
729        = invcalibBAOfactors_On_Cycle_Ch1(ir,1);
730    }//eo loop on coef On
731     
732    //Loop on cyles
733    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
734
735      //Look if the cycle has been calibrated...
736      bool isCycleCalibrated = 
737        ( calibBAO_On_Run_Ch0.count(*ic)    *
738          calibBAO_On_Run_Ch1.count(*ic)    *
739          calibBAO_Off_Run_Ch0.count(*ic)   *
740          calibBAO_Off_Run_Ch1.count(*ic)   *
741          calibBAO_On_Cycle_Ch0.count(*ic)  *
742          calibBAO_On_Cycle_Ch1.count(*ic)  *
743          calibBAO_Off_Cycle_Ch0.count(*ic) *
744          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
745
746      if(para.debuglev_>9) {
747        cout << "Calibration coefficients for cycle "<<*ic << endl; 
748        if (isCycleCalibrated) {
749          cout << "Cycle calibrated " << endl;
750          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
751               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
752               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
753               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
754               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
755               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
756               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
757               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
758        } else {
759          cout << "Cycle NOT calibrated " << endl;
760        }
761      }//debug
762
763      if ( ! isCycleCalibrated ) continue;
764
765      totalNumberCycles++;
766      //Fill NTuple
767      xnt[0] = rdateOfRun;
768      xnt[1] = totalNumberRuns;
769      xnt[2] = totalNumberCycles;
770      xnt[3] = calibBAO_Off_Run_Ch0[*ic];
771      xnt[4] = calibBAO_Off_Run_Ch1[*ic];
772      xnt[5] = calibBAO_On_Run_Ch0[*ic];
773      xnt[6] = calibBAO_On_Run_Ch1[*ic];
774      xnt[7] = calibBAO_Off_Cycle_Ch0[*ic];
775      xnt[8] = calibBAO_Off_Cycle_Ch1[*ic];
776      xnt[9] = calibBAO_On_Cycle_Ch0[*ic];
777      xnt[10] = calibBAO_On_Cycle_Ch1[*ic];
778      xnt[11] = invcalibBAO_Off_Run_Ch0[*ic];
779      xnt[12] = invcalibBAO_Off_Run_Ch1[*ic];
780      xnt[13] = invcalibBAO_On_Run_Ch0[*ic];
781      xnt[14] = invcalibBAO_On_Run_Ch1[*ic];
782      xnt[15] = invcalibBAO_Off_Cycle_Ch0[*ic];
783      xnt[16] = invcalibBAO_Off_Cycle_Ch1[*ic];
784      xnt[17] = invcalibBAO_On_Cycle_Ch0[*ic];
785      xnt[18] = invcalibBAO_On_Cycle_Ch1[*ic];
786
787      calibcoeffs.Fill(xnt);
788
789      //Quit if enough
790      if (totalNumberCycles >= para.maxNumberCycles_) break;   
791
792    }//eo loop on cycles
793    if (totalNumberCycles >= para.maxNumberCycles_) break;         
794    totalNumberRuns++;
795  }//eo loop on spectra in a file
796  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
797
798  {//save the results
799    stringstream tmp;
800    tmp << totalNumberCycles;
801    string fileName = para.outPath_+"/calibcoeffTuple_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
802
803    POutPersist fos(fileName);
804    cout << "Save output in " << fileName << endl;
805
806    string tag = "calib";
807    fos << PPFNameTag(tag) << calibcoeffs;
808
809  }//end of save
810}
811
812//-------------------------------------------------------
813//Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra
814//
815void meanCalibBAODiffOnOffCycles() throw(string) {
816
817  list<string> listOfFiles;
818  string directoryName;
819  directoryName = para.inPath_ + "/" + para.sourceName_;
820
821  //Make the listing of the directory
822  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
823 
824  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
825  iFileEnd = listOfFiles.end();
826
827  //mean ON-OFF over the list of cycles
828  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);     //set to 0
829  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
830  TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
831  TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
832  static const int NINFO=25;
833  char* onffTupleName[NINFO]={"cycle"
834                              ,"onoffRaw0","onoffRaw1"
835                              ,"onoffRun0","onoffRun1"
836                              ,"onoffCycle0","onoffCycle1"
837                              ,"onoffRaw01420","onoffRaw11420"
838                              ,"onoffRun01420","onoffRun11420"
839                              ,"onoffCycle01420","onoffCycle11420"
840                              ,"onoffRaw01420side","onoffRaw11420side"
841                              ,"onoffRun01420side","onoffRun11420side"
842                              ,"onoffCycle01420side","onoffCycle11420side"
843                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
844                              ,"offRaw0f14101415","offRaw1f14101415"
845                              ,"onRaw0f14101415","onRaw1f14101415"
846  };
847  NTuple onoffevolution(NINFO,onffTupleName);
848  r_4 xnt[NINFO];
849
850  //Lower and Higher freq. arround the Calibration Freq. bin to perform mean follow up
851  sa_size_t chCalibLow  = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
852  sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
853  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
854  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
855  sa_size_t ch1420High = freqToChan(1420.4+0.2);
856
857  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
858  sa_size_t ch1420aLow  = freqToChan(1418);
859  sa_size_t ch1420aHigh = freqToChan(1419);
860  sa_size_t ch1420bLow  = freqToChan(1422);
861  sa_size_t ch1420bHigh = freqToChan(1423);
862 
863  //1400-1420Mhz
864  sa_size_t ch1400 = freqToChan(1400);
865  //  sa_size_t ch1405 = freqToChan(1400);
866  sa_size_t ch1410 = freqToChan(1410);
867  sa_size_t ch1415 = freqToChan(1415);
868  sa_size_t ch1420 = freqToChan(1420);
869
870  if (para.debuglev_>0){
871    cout << "freq. band for follow up [" <<  chCalibLow << ", "<< chCalibHigh << "]" << endl;
872    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
873    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
874    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
875  }
876 
877  //Loop on files/run
878
879  int totalNumberCycles=0; //total number of cycles for normalisation
880  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
881    if (para.debuglev_>90){
882      cout << "load file <" << *iFile << ">" << endl;
883    }
884
885    vector<string> tokens;
886    split(*iFile,"_",tokens);
887    string dateOfRun = tokens[1];
888    if (para.debuglev_>90){
889      cout << "date <" << dateOfRun << ">" << endl;
890    }
891    vector<string> tokens2;
892    split(tokens[2],".",tokens2);
893    string srcLower = tokens2[0];
894
895    PInPersist fin(*iFile);
896    vector<string> vec = fin.GetNameTags();
897
898    vector<string> modeList;
899    modeList.push_back("On");
900    modeList.push_back("Off");
901    vector<string>::const_iterator iMode;
902   
903    map<string, list<int> > cycleModeCollect;
904   
905    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
906      list<string> listOfSpectra;
907      //Keep only required PPF objects
908      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
909      std::remove_copy_if(
910                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
911                          not1(StringMatch(matchstr))
912                          );
913     
914      listOfSpectra.sort(stringCompare);
915      iSpecEnd = listOfSpectra.end();
916     
917      matchstr = "[0-9]+";
918      //Loop of spectra matrix
919      list<int> listOfCycles;
920      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
921        int b,e;
922        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
923        if (para.debuglev_>90){
924          cout << " spactra <" << *iSpec << ">" << endl;
925          cout << " cycle " << iSpec->substr(b) << endl;
926        }
927        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
928      }//end loop spectra
929      cycleModeCollect[*iMode] = listOfCycles;
930    }//end of mode   
931
932    //Take the Intersection of the list Of cycles in mode Off and On   
933    list<int> commonCycles;
934    set_intersection(cycleModeCollect["On"].begin(),
935                     cycleModeCollect["On"].end(),
936                     cycleModeCollect["Off"].begin(),
937                     cycleModeCollect["Off"].end(),
938                     back_inserter(commonCycles)
939                     );
940   
941    if (para.debuglev_>90){
942      cout << "Liste of cycles common to On & Off: <";
943      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
944        cout << *i << " ";
945      }
946      cout << ">" << endl;
947    }
948   
949    //
950    //Load BAO Calibration factors "per Cycle and Channels"
951    //Compute the mean per Cycle to
952    //  fill factors "per Run and Channels" with the same cycle list
953    //
954    //
955    //TODO improve the code....
956
957    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
958    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
959    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
960    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
961   
962    string calibFileName;
963    ifstream ifs;
964    sa_size_t nr,nc; //values read
965
966    //OFF Cycle per Channel
967    calibFileName = directoryName + "/" 
968      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
969      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
970    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
971    ifs.open(calibFileName.c_str());
972    if ( ! ifs.is_open() ) {
973
974      throw calibFileName + " cannot be opened...";
975    }   
976    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
977    if(para.debuglev_>9){
978      cout << "(nr,nc): "<< nr << "," << nc << endl;
979      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
980      cout << endl;
981    }
982
983    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
984    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
985    {//Compute the mean
986      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
987      double mean,sigma;
988      MeanSigma(coef,mean,sigma);
989      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
990    }
991    if(para.debuglev_>9){
992      cout << "Fill calib. with mean value " << endl; 
993      calibBAOfactors_Off_Run_Ch0.Print(cout);
994      cout << endl;
995    }
996    ifs.close();
997
998    //
999    calibFileName = directoryName + "/" 
1000      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
1001      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
1002    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
1003    ifs.open(calibFileName.c_str());
1004    if ( ! ifs.is_open() ) {
1005
1006      throw calibFileName + " cannot be opened...";
1007    }   
1008    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
1009    if(para.debuglev_>9){
1010      cout << "(nr,nc): "<< nr << "," << nc << endl;
1011      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
1012      cout << endl;
1013    }
1014    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
1015    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
1016    {//Compute the mean
1017      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
1018      double mean,sigma;
1019      MeanSigma(coef,mean,sigma);
1020      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
1021      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
1022    }
1023    if(para.debuglev_>9){
1024      cout << "Fill calib. with mean value " << endl; 
1025      calibBAOfactors_Off_Run_Ch1.Print(cout);
1026      cout << endl;
1027    }
1028    ifs.close();
1029
1030    //ON Cycle per Channel
1031    calibFileName = directoryName + "/" 
1032      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
1033      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
1034    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
1035    ifs.open(calibFileName.c_str());
1036    if ( ! ifs.is_open() ) {
1037
1038      throw calibFileName + " cannot be opened...";
1039    }   
1040    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
1041    if(para.debuglev_>9){
1042      cout << "(nr,nc): "<< nr << "," << nc << endl;
1043      calibBAOfactors_On_Cycle_Ch0.Print(cout);
1044      cout << endl;     
1045    }
1046
1047    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
1048    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
1049    {//Compute the mean
1050      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
1051      double mean,sigma;
1052      MeanSigma(coef,mean,sigma);
1053      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
1054      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
1055    }
1056    if(para.debuglev_>9){
1057      cout << "Fill calib. with mean value " << endl; 
1058      calibBAOfactors_On_Run_Ch0.Print(cout);
1059      cout << endl;
1060    }
1061    ifs.close();
1062
1063   
1064    calibFileName = directoryName + "/" 
1065      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
1066      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
1067    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
1068    ifs.open(calibFileName.c_str());
1069    if ( ! ifs.is_open() ) {
1070      throw calibFileName + " cannot be opened...";
1071    }   
1072    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
1073    if(para.debuglev_>9){
1074      cout << "(nr,nc): "<< nr << "," << nc << endl;
1075      calibBAOfactors_On_Cycle_Ch1.Print(cout);
1076      cout << endl;
1077    }
1078    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
1079    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
1080    {//Compute the mean
1081      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
1082      double mean,sigma;
1083      MeanSigma(coef,mean,sigma);
1084      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
1085      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
1086    }
1087    if(para.debuglev_>9){
1088      cout << "Fill calib. with mean value " << endl; 
1089      calibBAOfactors_On_Run_Ch1.Print(cout);
1090      cout << endl;
1091    }
1092
1093    ifs.close();
1094   
1095    //link <cycle> - <calibration coefficient>
1096    //We cannot rely on identical cycle list of the OFF and ON calibration
1097    map<int,r_4> calibBAO_Off_Run_Ch0;
1098    map<int,r_4> calibBAO_Off_Run_Ch1;
1099    map<int,r_4> calibBAO_On_Run_Ch0;
1100    map<int,r_4> calibBAO_On_Run_Ch1;
1101
1102    map<int,r_4> calibBAO_Off_Cycle_Ch0;
1103    map<int,r_4> calibBAO_Off_Cycle_Ch1;
1104    map<int,r_4> calibBAO_On_Cycle_Ch0;
1105    map<int,r_4> calibBAO_On_Cycle_Ch1;
1106
1107    //per Run based BAO coefficients
1108    nr = calibBAOfactors_Off_Run_Ch0.NRows();
1109    for (sa_size_t ir=0; ir<nr; ++ir){
1110      cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
1111           << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
1112
1113      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
1114        = calibBAOfactors_Off_Run_Ch0(ir,1);
1115      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
1116        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
1117      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
1118        = calibBAOfactors_Off_Run_Ch1(ir,1);
1119      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
1120        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
1121    }//eo loop on coef Off
1122   
1123    nr = calibBAOfactors_On_Run_Ch0.NRows();
1124    for (sa_size_t ir=0; ir<nr; ++ir){
1125      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
1126        = calibBAOfactors_On_Run_Ch0(ir,1);
1127      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
1128        = calibBAOfactors_On_Cycle_Ch0(ir,1);
1129      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
1130        = calibBAOfactors_On_Run_Ch1(ir,1);
1131      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
1132        = calibBAOfactors_On_Cycle_Ch1(ir,1);
1133    }//eo loop on coef On
1134     
1135    //Loop on cyles
1136    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1137
1138      //Look if the cycle has been calibrated...
1139      bool isCycleCalibrated = 
1140        ( calibBAO_On_Run_Ch0.count(*ic)    *
1141          calibBAO_On_Run_Ch1.count(*ic)    *
1142          calibBAO_Off_Run_Ch0.count(*ic)   *
1143          calibBAO_Off_Run_Ch1.count(*ic)   *
1144          calibBAO_On_Cycle_Ch0.count(*ic)  *
1145          calibBAO_On_Cycle_Ch1.count(*ic)  *
1146          calibBAO_Off_Cycle_Ch0.count(*ic) *
1147          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
1148
1149      if(para.debuglev_>9) {
1150        cout << "Calibration coefficients for cycle "<<*ic << endl; 
1151        if (isCycleCalibrated) {
1152          cout << "Cycle calibrated " << endl;
1153          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
1154               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
1155               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
1156               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
1157               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
1158               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
1159               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
1160               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
1161        } else {
1162          cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
1163        }
1164      }//debug
1165
1166
1167      if ( ! isCycleCalibrated ) continue;
1168     
1169      string ppftag;
1170      //load ON phase
1171      stringstream cycle;
1172      cycle << *ic;
1173     
1174      ppftag = "specRawOn"+cycle.str();
1175      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1176      fin.GetObject(aSpecOn,ppftag);
1177
1178      TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
1179      aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
1180      aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
1181
1182      TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
1183      aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
1184      aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
1185     
1186      //Load OFF phase
1187      ppftag = "specRawOff"+cycle.str();
1188      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1189      fin.GetObject(aSpecOff,ppftag);
1190
1191      TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
1192      aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
1193      aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
1194
1195      TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
1196      aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
1197      aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
1198
1199
1200      //Perform the difference ON-OFF with the different calibration options
1201      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1202      meanDiffONOFF_noCalib += diffOnOff_noCalib;
1203     
1204      //JEC 29/10/11 add ON-OFF/OFF
1205      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1206      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1207      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1208      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
1209     
1210      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
1211      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1212
1213      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
1214      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
1215
1216      TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
1217      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
1218
1219      totalNumberCycles++;
1220
1221      //Fill NTuple
1222      xnt[0] = totalNumberCycles;
1223 
1224      //Follow up arround the Calibration Freq.
1225      TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
1226      meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
1227      xnt[1] = meanInRange_CalibFreq_noCalib(0);
1228      xnt[2] = meanInRange_CalibFreq_noCalib(1);
1229
1230      TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
1231      meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
1232      xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
1233      xnt[4] = meanInRange_CalibFreq_perRunCalib(1);
1234
1235      TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
1236      meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
1237      xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
1238      xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);
1239
1240
1241      //Follow up arround the 1420.4MHz Freq.
1242      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1243      meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1244      xnt[7] = meanInRange_1420Freq_noCalib(0);
1245      xnt[8] = meanInRange_1420Freq_noCalib(1);
1246
1247      TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
1248      meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
1249      xnt[9] = meanInRange_1420Freq_perRunCalib(0);
1250      xnt[10] = meanInRange_1420Freq_perRunCalib(1);
1251
1252      TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
1253      meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
1254      xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
1255      xnt[12] = meanInRange_1420Freq_perCycleCalib(1);
1256
1257
1258      //Follow up below the 1420.4MHz Freq.
1259      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1260      meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1261      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1262      meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1263
1264      xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1265      xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1266
1267
1268      TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
1269      meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
1270      TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
1271      meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);
1272
1273      xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
1274      xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;
1275
1276
1277      TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
1278      meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
1279      TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
1280      meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);
1281
1282      xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
1283      xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
1284
1285
1286      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1287      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
1288      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
1289      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
1290      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
1291
1292      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
1293      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1294      aSpecOffovOff.Div(aSpecOffFiltered);
1295
1296      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1297      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1298      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1299
1300      xnt[21] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1301      xnt[22] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1302
1303      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1304      aSpecOnovOff.Div(aSpecOffFiltered);
1305
1306      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1307      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1308      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1309
1310      xnt[23] = meanInRange_1410a1415Freq_ON_noCalib(0);
1311      xnt[24] = meanInRange_1410a1415Freq_ON_noCalib(1);
1312     
1313      //store infos to Ntuple
1314      onoffevolution.Fill(xnt);
1315
1316      //Quit if enough
1317      if (totalNumberCycles >= para.maxNumberCycles_) break;   
1318
1319    }//eo loop on cycles
1320    if (totalNumberCycles >= para.maxNumberCycles_) break;         
1321 
1322  }//eo loop on spectra in a file
1323  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1324  //Normalisation
1325  if(totalNumberCycles > 0){
1326    //JEC 29/10 add ON-OFF/OFF
1327    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
1328    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
1329    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
1330    meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
1331  } 
1332 
1333  //Compute the reduced version of the mean and sigma
1334  TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1335  TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1336  reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
1337
1338  TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1339  TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1340  reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
1341
1342  TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1343  TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1344  reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
1345
1346  {//save the results
1347    stringstream tmp;
1348    tmp << totalNumberCycles;
1349    string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1350
1351    POutPersist fos(fileName);
1352    cout << "Save output in " << fileName << endl;
1353
1354    string tag = "meanNoCalib";
1355    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1356   
1357    //JEC 29/10/11
1358    tag = "meanOvOffNoCalib";
1359    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1360
1361    tag = "meanPerRunCalib";
1362    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
1363    tag = "meanPerCycleCalib";
1364    fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
1365
1366    tag = "redmeanNoCalib";
1367    fos << PPFNameTag(tag) << meanRedMtx_noCalib;
1368    tag = "redsigmaNoCalib";
1369    fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
1370
1371    tag = "redmeanPerRunCalib";
1372    fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
1373    tag = "redsigmaPerRunCalib";
1374    fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
1375
1376    tag = "redmeanPerCycleCalib";
1377    fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
1378    tag = "redsigmaPerCycleCalib";
1379    fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
1380   
1381    tag = "onoffevol";
1382    fos << PPFNameTag(tag) << onoffevolution;
1383 
1384  }//end of save
1385
1386
1387}
1388//JEC 14/11/11 New meanRawDiffOnOffCycles START
1389//-------------------------------------------------------
1390//Compute the mean of Diff ON-OFF/OFF from Raw spectra
1391//Used like:
1392//
1393void meanRawDiffOnOffCycles() throw(string) {
1394  list<string> listOfFiles;
1395  string directoryName;
1396  directoryName = para.inPath_ + "/" + para.sourceName_;
1397
1398  //Make the listing of the directory
1399  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1400 
1401  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1402  iFileEnd = listOfFiles.end();
1403
1404  //mean ON-OFF over the list of cycles
1405  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
1406  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
1407  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
1408  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
1409  //Tuple only for RAW things to follow
1410  static const int NINFO=11;
1411  char* onffTupleName[NINFO]={"cycle"
1412                              ,"onoffRaw01420","onoffRaw11420"
1413                              ,"onoffRaw01420side","onoffRaw11420side"
1414                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
1415                              ,"offRaw0f14101415","offRaw1f14101415"
1416                              ,"onRaw0f14101415","onRaw1f14101415"
1417  };
1418  NTuple onoffevolution(NINFO,onffTupleName);
1419  r_4 xnt[NINFO];
1420
1421  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
1422  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
1423  sa_size_t ch1420High = freqToChan(1420.4+0.2);
1424
1425  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
1426  sa_size_t ch1420aLow  = freqToChan(1418);
1427  sa_size_t ch1420aHigh = freqToChan(1419);
1428  sa_size_t ch1420bLow  = freqToChan(1422);
1429  sa_size_t ch1420bHigh = freqToChan(1423);
1430 
1431  //1400-1420Mhz
1432  sa_size_t ch1400 = freqToChan(1400);
1433  //  sa_size_t ch1405 = freqToChan(1400);
1434  sa_size_t ch1410 = freqToChan(1410);
1435  sa_size_t ch1415 = freqToChan(1415);
1436  sa_size_t ch1420 = freqToChan(1420);
1437
1438  if (para.debuglev_>0){
1439    cout << "freq. band for follow up [" <<  ch1420Low << ", "<< ch1420High << "]" << endl;
1440    cout << "freq. band for follow up [" <<  ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
1441    cout << "freq. band for follow up [" <<  ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
1442  }
1443
1444  int totalNumberCycles=0; //total number of cycles
1445
1446  //Loop on files
1447  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1448    if (para.debuglev_>90){
1449      cout << "load file <" << *iFile << ">" << endl;
1450    }
1451
1452    vector<string> tokens;
1453    split(*iFile,"_",tokens);
1454    string dateOfRun = tokens[1];
1455    if (para.debuglev_>90){
1456      cout << "date <" << dateOfRun << ">" << endl;
1457    }
1458    vector<string> tokens2;
1459    split(tokens[2],".",tokens2);
1460    string srcLower = tokens2[0];
1461
1462    PInPersist fin(*iFile);
1463    vector<string> vec = fin.GetNameTags();
1464
1465    vector<string> modeList;
1466    modeList.push_back("On");
1467    modeList.push_back("Off");
1468    vector<string>::const_iterator iMode;
1469   
1470    map<string, list<int> > cycleModeCollect;
1471   
1472    for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
1473      list<string> listOfSpectra;
1474      //Keep only required PPF objects
1475      string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
1476      std::remove_copy_if(
1477                          vec.begin(), vec.end(), back_inserter(listOfSpectra),
1478                          not1(StringMatch(matchstr))
1479                          );
1480     
1481      listOfSpectra.sort(stringCompare);
1482      iSpecEnd = listOfSpectra.end();
1483     
1484      matchstr = "[0-9]+";
1485      //Loop of spectra matrix
1486      list<int> listOfCycles;
1487      for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
1488        int b,e;
1489        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
1490        if (para.debuglev_>90){
1491          cout << " spectra <" << *iSpec << ">" << endl;
1492          cout << " cycle " << iSpec->substr(b) << endl;
1493        }
1494        listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
1495      }//end loop spectra
1496      cycleModeCollect[*iMode] = listOfCycles;
1497    }//end of mode   
1498
1499    //Take the Intersection of the list Of cycles in mode Off and On   
1500    list<int> commonCycles;
1501    set_intersection(cycleModeCollect["On"].begin(),
1502                     cycleModeCollect["On"].end(),
1503                     cycleModeCollect["Off"].begin(),
1504                     cycleModeCollect["Off"].end(),
1505                     back_inserter(commonCycles)
1506                     );
1507   
1508    if (para.debuglev_>90){
1509      cout << "Liste of cycles common to On & Off: <";
1510      for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
1511        cout << *i << " ";
1512      }
1513      cout << ">" << endl;
1514    }
1515
1516    //Loop on cyles
1517    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1518     
1519      // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440
1520      if ( *ic == 1 && srcLower == "abell1205" ) {
1521          if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {
1522            cout << "Skipping non-calibrated cycle " << *ic << endl;
1523            continue;
1524          }
1525      } else if ( *ic == 1 && srcLower == "abell2440" ) {
1526          if ( dateOfRun == "20110516" ) {
1527            cout << "Skipping non-calibrated cycle " << *ic << endl;
1528            continue;
1529          }
1530      } else if ( *ic == 3 && srcLower == "abell1205" ) {
1531          if ( dateOfRun == "20110810" ) {
1532            cout << "Skipping non-calibrated cycle " << *ic << endl;
1533            continue;
1534          }
1535      }
1536
1537      string ppftag;
1538      //load ON phase
1539      stringstream cycle;
1540      cycle << *ic;
1541      ppftag = "specRawOn"+cycle.str();
1542      TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1543      fin.GetObject(aSpecOn,ppftag);
1544
1545      //Load OFF phase
1546      ppftag = "specRawOff"+cycle.str();
1547      TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1548      fin.GetObject(aSpecOff,ppftag);
1549     
1550      //Perform the difference ON-OFF
1551      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1552
1553      meanDiffONOFF_noCalib += diffOnOff_noCalib;
1554     
1555      //JEC 29/10/11 add ON-OFF/OFF
1556      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1557      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1558      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1559      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
1560     
1561      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
1562      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1563
1564      //JEC 15/11/11 add ON/OFF and OFF/OFF
1565      TMatrix<r_4> onOvOff(aSpecOn,false);
1566      onOvOff.Div(aSpecOffFiltered);
1567      meanONovOFF_noCalib += onOvOff;
1568     
1569      TMatrix<r_4> offOvOff(aSpecOff,false);
1570      offOvOff.Div(aSpecOffFiltered);
1571      meanOFFovOFF_noCalib += offOvOff;
1572
1573      totalNumberCycles++;
1574
1575      //Fill NTuple
1576      xnt[0] = totalNumberCycles;
1577 
1578      //Follow up arround the 1420.4MHz Freq.
1579      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1580      meanInRange(diffOnOffOvOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1581      xnt[1] = meanInRange_1420Freq_noCalib(0);
1582      xnt[2] = meanInRange_1420Freq_noCalib(1);
1583
1584
1585      //Follow up below the 1420.4MHz Freq.
1586      TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1587      meanInRange(diffOnOffOvOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1588      TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1589      meanInRange(diffOnOffOvOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1590
1591      xnt[3] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1592      xnt[4] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1593
1594
1595      //JEC 25/10/11 follow 1410-1415 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1596      TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS);
1597      meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib);
1598      xnt[5] = meanInRange_1410a1415Freq_noCalib(0);
1599      xnt[6] = meanInRange_1410a1415Freq_noCalib(1);
1600
1601      //JEC 18/11/11 follow up the 1410-1415MHz OFF only
1602      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1603      aSpecOffovOff.Div(aSpecOffFiltered);
1604
1605      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1606      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1607
1608      xnt[7] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1609      xnt[8] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1610
1611      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1612      aSpecOnovOff.Div(aSpecOffFiltered);
1613
1614      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1615      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1616
1617      xnt[9] = meanInRange_1410a1415Freq_ON_noCalib(0);
1618      xnt[10] = meanInRange_1410a1415Freq_ON_noCalib(1);
1619     
1620      //store infos to Ntuple
1621      onoffevolution.Fill(xnt);
1622
1623      //Quit if enough
1624      if (totalNumberCycles >= para.maxNumberCycles_) break;   
1625    }//end of cycles
1626
1627    if (totalNumberCycles >= para.maxNumberCycles_) break;         
1628
1629  }//end files
1630  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1631  //Normalisation
1632  if(totalNumberCycles > 0){
1633    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
1634    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
1635    meanONovOFF_noCalib         /= (r_4)totalNumberCycles;
1636    meanOFFovOFF_noCalib        /= (r_4)totalNumberCycles;
1637  } 
1638
1639  {//save results
1640    stringstream tmp;
1641    tmp << totalNumberCycles;
1642    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1643
1644    POutPersist fos(fileName);
1645    cout << "Save output in " << fileName << endl;
1646
1647    string tag = "meanNoCalib";
1648    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1649   
1650    tag = "meanOvOffNoCalib";
1651    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1652   
1653    tag = "meanOnovOffNoCalib";
1654    fos << PPFNameTag(tag) << meanONovOFF_noCalib;
1655
1656    tag = "meanOffovOffNoCalib";
1657    fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;
1658
1659    tag = "onoffevol";
1660    fos << PPFNameTag(tag) << onoffevolution;
1661
1662  }//end save
1663
1664
1665}
1666//JEC 14/11/11 New meanRawDiffOnOffCycles END
1667//-------------------------------------------------------
1668//JEC 14/11/11 Obsolete START
1669//-------------------------------------------------------
1670//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1671//Used like:
1672//
1673// void meanRawDiffOnOffCycles() throw(string) {
1674//   list<string> listOfFiles;
1675//   string directoryName;
1676//   directoryName = para.inPath_ + "/" + para.sourceName_;
1677
1678//   //Make the listing of the directory
1679//   listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1680 
1681//   list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1682//   iFileEnd = listOfFiles.end();
1683 
1684//   StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1685//   TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1686//   uint_4 nSpectra=0;
1687//   //Loop on files
1688//   for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1689//     if (para.debuglev_>90){
1690//       cout << "load file <" << *iFile << ">" << endl;
1691//     }
1692//     PInPersist fin(*iFile);
1693//     vector<string> vec = fin.GetNameTags();
1694//     list<string> listOfSpectra;
1695//     //Keep only required PPF objects
1696//     std::remove_copy_if(
1697//                      vec.begin(), vec.end(), back_inserter(listOfSpectra),
1698//                      not1(match)
1699//                      );
1700   
1701//     listOfSpectra.sort(stringCompare);
1702//     iSpecEnd = listOfSpectra.end();
1703//     //Loop of spectra matrix
1704//     for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd;  ++iSpec){
1705//       if (para.debuglev_>90){
1706//      cout << " spactra <" << *iSpec << ">" << endl;
1707//       }
1708//       TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1709//       fin.GetObject(aSpec,*iSpec);
1710//       //How to see if the GetObject is ok?? Ask Reza
1711//       nSpectra++;
1712//       meanOfSpectra+=aSpec;
1713//     }//eo loop on spectra in a file
1714//   }//eo loop on files
1715 
1716//   //Normalisation
1717//   if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
1718
1719//   //Compute the reduced version of the mean and sigma
1720//   TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1721//   TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1722//   reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
1723
1724//   {//Save the result
1725//     stringstream tmp;
1726//     tmp << nSpectra;
1727//     string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1728//     cout << "Save mean based on " <<  nSpectra << " cycles " << endl;
1729//     POutPersist fos(fileName);
1730
1731//     string tag = "mean";
1732//     fos << PPFNameTag(tag) << meanOfSpectra;
1733//     tag = "meanred";
1734//     fos << PPFNameTag(tag) << meanRedMtx;
1735//     tag = "sigmared";
1736//     fos << PPFNameTag(tag) << sigmaRedMtx;
1737//   }
1738// }
1739//JEC 14/11/11 Obsolete END
1740//-------------------------------------------------------
1741//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1742//Used like:
1743//
1744void medianRawDiffOnOffCycles() throw(string) {
1745  list<string> listOfFiles;
1746  string directoryName;
1747  directoryName = para.inPath_ + "/" + para.sourceName_;
1748
1749  //Make the listing of the directory
1750  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1751 
1752  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1753  iFileEnd = listOfFiles.end();
1754 
1755
1756  TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
1757
1758  StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1759  uint_4 nSpectra=0;
1760  //Loop on files
1761  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1762    if (para.debuglev_>90){
1763      cout << "load file <" << *iFile << ">" << endl;
1764    }
1765    PInPersist fin(*iFile);
1766    vector<string> vec = fin.GetNameTags();
1767    list<string> listOfSpectra;
1768    //Keep only required PPF objects
1769    std::remove_copy_if(
1770                        vec.begin(), vec.end(), back_inserter(listOfSpectra),
1771                        not1(match)
1772                        );
1773   
1774    listOfSpectra.sort(stringCompare);
1775    iSpecEnd = listOfSpectra.end();
1776    //Loop of spectra matrix
1777    for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ;  ++iSpec){
1778      if (para.debuglev_>90){
1779        cout << " spactra <" << *iSpec << ">" << endl;
1780      }
1781      TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1782      fin.GetObject(aSpec,*iSpec);
1783
1784      tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
1785
1786      nSpectra++;
1787    }//eo loop on spectra in a file
1788  }//eo loop on files
1789 
1790
1791 
1792  TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1793  //Compute the median for each freq. and Channel
1794  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1795    for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
1796      TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
1797      vector<r_4> tmp;
1798      tmp0.FillTo(tmp);
1799      medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
1800    }
1801  }
1802
1803
1804  //Compute the reduced version of the mean and sigma
1805  TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1806  TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1807  reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
1808
1809
1810  sa_size_t f1320=freqToChan(1320.);
1811  sa_size_t f1345=freqToChan(1345.);
1812  sa_size_t f1355=freqToChan(1355.);
1813  sa_size_t f1380=freqToChan(1380.);
1814  //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
1815  if (para.debuglev_>9){
1816    cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
1817  }
1818  TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
1819  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1820    double meanMed1;
1821    double sigmaMed1;
1822    TVector<r_4> band1;
1823    band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
1824    MeanSigma(band1,meanMed1,sigmaMed1);
1825    double meanMed2;
1826    double sigmaMed2;
1827    TVector<r_4> band2;
1828    band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
1829    MeanSigma(band2,meanMed2,sigmaMed2);
1830    meanMed(iCh) = (meanMed1+meanMed2)*0.5;
1831  } 
1832  meanMed.Print(cout);
1833  cout << endl;
1834
1835 
1836  //Compute the sigma in the range 1320MHz-1380MHz
1837  if (para.debuglev_>9){
1838    cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
1839  }
1840  TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
1841  sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1842  sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1843
1844  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1845    double meanSigma;
1846    double sigmaSigma;
1847    TVector<r_4> band;
1848    band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
1849    MeanSigma(band,meanSigma,sigmaSigma);
1850    meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
1851    sigmaMed(iCh) = meanSigma;
1852  }
1853  sigmaMed.Print(cout);
1854  cout << endl;
1855
1856 
1857 
1858  if (para.debuglev_>9){
1859    cout << "Compute medianOfSpectraNorm" << endl;
1860  }
1861  TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
1862  for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1863    medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
1864    medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
1865  }
1866
1867 
1868
1869  {//Save the result
1870    stringstream tmp;
1871    tmp << nSpectra;
1872    string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1873    cout << "Save median based on " <<  nSpectra << " cycles " << endl;
1874    POutPersist fos(fileName);
1875
1876    string tag = "median";
1877    fos << PPFNameTag(tag) << medianOfSpectra;
1878
1879    tag = "medianNorm";
1880    fos << PPFNameTag(tag) << medianOfSpectraNorm;
1881   
1882
1883    tag = "meanmedred";
1884    fos << PPFNameTag(tag) << meanRedMtx;
1885    tag = "sigmamedred";
1886    fos << PPFNameTag(tag) << sigmaRedMtx;
1887    tag = "cycleVsfreq";
1888   
1889    TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1890    fos << PPFNameTag(tag) << tarr;
1891  }
1892}
1893
1894//-------------------------------------------------------
1895int main(int narg, char* arg[]) {
1896 
1897  int rc = 0; //return code
1898  string msg; //message used in Exceptions
1899
1900
1901
1902  //default value for initial parameters (see Para structure on top of the file)
1903  string debuglev = "0";
1904  string action = "meanDiffOnOff";
1905  string inputPath = "."; 
1906  string outputPath = "."; 
1907  string sourceName = "Abell85";
1908  string ppfFile;
1909  string nSliceInFreq = "32";
1910  string typeOfCalib="perRun";
1911  string calibFreq = "1346";
1912  string calibBandFreq="6.25";
1913  string mxcycles;
1914
1915  //  bool okarg=false;
1916  int ka=1;
1917  while (ka<narg) {
1918    if (strcmp(arg[ka],"-h")==0) {
1919      cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1920           << " -mxcycles <number> (max. number of cycles to be treated)\n"
1921           << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1922           << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1923           << " -src [Abell85]\n" 
1924           << " -inPath [.]|<top_root_dir of the ppf file>\n" 
1925           << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
1926           << " -outPath [.]|<dir of the output> \n"
1927           << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1928           << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1929           << " -debug <level> "
1930           << endl;
1931      return 0;
1932    }
1933    else if (strcmp(arg[ka],"-debug")==0) {
1934      debuglev=arg[ka+1];
1935      ka+=2;
1936    }
1937    else if (strcmp(arg[ka],"-act")==0) {
1938      action=arg[ka+1];
1939      ka+=2;
1940    }
1941    else if (strcmp(arg[ka],"-calibfreq")==0) {
1942      calibFreq=arg[ka+1];
1943      ka+=2;
1944    }   
1945    else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1946      calibBandFreq=arg[ka+1];
1947      ka+=2;
1948    }   
1949    else if (strcmp(arg[ka],"-mxcycles")==0) {
1950      mxcycles=arg[ka+1];
1951      ka+=2;
1952    }   
1953    else if (strcmp(arg[ka],"-inPath")==0) {
1954      inputPath=arg[ka+1];
1955      ka+=2;
1956    }
1957    else if (strcmp(arg[ka],"-outPath")==0) {
1958      outputPath=arg[ka+1];
1959      ka+=2;
1960    }
1961    else if (strcmp(arg[ka],"-src")==0) {
1962      sourceName=arg[ka+1];
1963      ka+=2;
1964    }
1965    else if (strcmp(arg[ka],"-ppfFile")==0) {
1966      ppfFile=arg[ka+1];
1967      ka+=2;
1968    }
1969    else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1970      nSliceInFreq=arg[ka+1];
1971      ka+=2;
1972    }
1973    else ka++;
1974  }//eo while
1975
1976  para.debuglev_   = atoi(debuglev.c_str());
1977  para.inPath_     = inputPath;
1978  para.outPath_    = outputPath;
1979  para.sourceName_ = sourceName;
1980  para.ppfFile_    = ppfFile;
1981  para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1982  para.calibFreq_   = calibFreq;
1983  para.calibBandFreq_ = calibBandFreq;
1984  para.rcalibFreq_   = atof(calibFreq.c_str());
1985  para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1986  if (mxcycles != "") {
1987    para.maxNumberCycles_ = atoi(mxcycles.c_str());
1988  } else {
1989    para.maxNumberCycles_ = std::numeric_limits<int>::max();
1990  }
1991
1992  cout << "Dump Initial parameters ............" << endl;
1993  cout << " action = " << action << "\n"
1994       << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1995       << " inputPath = " << para.inPath_  << "\n" 
1996       << " outputPath = " <<  para.outPath_ << "\n"
1997       << " sourceName = " << para.sourceName_ << "\n"
1998       << " ppfFile = " <<  para.ppfFile_ << "\n"
1999       << " nSliceInFreq = " << para.nSliceInFreq_  << "\n"
2000       << " calibFreq = " <<  para.calibFreq_ << "\n"
2001       << " calibBandFreq = " <<  para.calibBandFreq_ << "\n"
2002       << " debuglev = "  << para.debuglev_   << "\n";
2003  cout << "...................................." << endl;
2004
2005  if ( "" == ppfFile ) {
2006    cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
2007         << endl;
2008    return 999;
2009  }
2010
2011
2012  try {
2013
2014//     int b,e;
2015//     char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
2016//     printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
2017
2018    if ( action == "meanRawDiffOnOff" ) {
2019      meanRawDiffOnOffCycles();
2020    } else if (action == "medianRawDiffOnOff") {
2021      medianRawDiffOnOffCycles();
2022    } else if (action == "meanCalibBAODiffOnOff") {
2023      meanCalibBAODiffOnOffCycles();
2024    } else {
2025      msg = "Unknown action " + action;
2026      throw msg;
2027    }
2028
2029
2030  }  catch (std::exception& sex) {
2031    cerr << "mergeRawOnOff.cc std::exception :"  << (string)typeid(sex).name() 
2032         << "\n msg= " << sex.what() << endl;
2033    rc = 78;
2034  }
2035  catch ( string str ) {
2036    cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
2037  }
2038  catch (...) {
2039    cerr << "mergeRawOnOff.cc catched unknown (...) exception  " << endl; 
2040    rc = 79; 
2041  } 
2042
2043  return 0;
2044
2045}
Note: See TracBrowser for help on using the repository browser.