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

Last change on this file since 562 was 562, checked in by campagne, 14 years ago

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