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

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

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