source: BAORadio/AmasNancay/v5/mergeAnaFiles.cc @ 675

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

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