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

Last change on this file since 611 was 611, checked in by torrento, 14 years ago

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