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

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

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