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

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

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