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

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

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