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

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

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