source: BAORadio/AmasNancay/trunk/mergeRawOnOff.cc@ 553

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

include image (jec)
File size: 16.4 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} para;
55//--------------------------------------------------------------
56//Utility functions
57
58sa_size_t freqTochan(r_4 f){
59 return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ);
60}
61//---------
62 class median_of_empty_list_exception:public std::exception{
63 virtual const char* what() const throw() {
64 return "Attempt to take the median of an empty list of numbers. "
65 "The median of an empty list is undefined.";
66 }
67 };
68 template<class RandAccessIter>
69 double median(RandAccessIter begin, RandAccessIter end)
70 throw(median_of_empty_list_exception){
71 if(begin == end){ throw median_of_empty_list_exception(); }
72 std::size_t size = end - begin;
73 std::size_t middleIdx = size/2;
74 RandAccessIter target = begin + middleIdx;
75 std::nth_element(begin, target, end);
76
77 if(size % 2 != 0){ //Odd number of elements
78 return *target;
79 }else{ //Even number of elements
80 double a = *target;
81 RandAccessIter targetNeighbor= target-1;
82 std::nth_element(begin, targetNeighbor, end);
83 return (a+*targetNeighbor)/2.0;
84 }
85 }
86
87//--------------------------------------------------------------
88char *regexp (const char *string, const char *patrn, int *begin, int *end) {
89 int i, w=0, len;
90 char *word = NULL;
91 regex_t rgT;
92 regmatch_t match;
93 regcomp(&rgT,patrn,REG_EXTENDED);
94 if ((regexec(&rgT,string,1,&match,0)) == 0) {
95 *begin = (int)match.rm_so;
96 *end = (int)match.rm_eo;
97 len = *end-*begin;
98 word=(char*)malloc(len+1);
99 for (i=*begin; i<*end; i++) {
100 word[w] = string[i];
101 w++; }
102 word[w]=0;
103 }
104 regfree(&rgT);
105 return word;
106}
107//-------
108sa_size_t round_sa(r_4 r) {
109 return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
110}
111//-----
112string StringToLower(string strToConvert){
113 //change each element of the string to lower case
114 for(unsigned int i=0;i<strToConvert.length();i++) {
115 strToConvert[i] = tolower(strToConvert[i]);
116 }
117 return strToConvert;//return the converted string
118}
119//-----
120bool stringCompare( const string &left, const string &right ){
121 if( left.size() < right.size() )
122 return true;
123 for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
124 if( tolower( *lit ) < tolower( *rit ) )
125 return true;
126 else if( tolower( *lit ) > tolower( *rit ) )
127 return false;
128 return false;
129}
130//-----
131list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
132 list<string> theList;
133
134
135 DIR *dip;
136 struct dirent *dit;
137 string msg; string fileName;
138 string fullFileName;
139 size_t found;
140
141 if ((dip=opendir(dir.c_str())) == NULL ) {
142 msg = "opendir failed on directory "+dir;
143 throw msg;
144 }
145 while ( (dit = readdir(dip)) != NULL ) {
146 fileName = dit->d_name;
147 found=fileName.find(filePettern);
148 if (found != string::npos) {
149 fullFileName = dir + "/";
150 fullFileName += fileName;
151 theList.push_back(fullFileName);
152 }
153 }//eo while
154 if (closedir(dip) == -1) {
155 msg = "closedir failed on directory "+dir;
156 throw msg;
157 }
158
159 theList.sort(stringCompare);
160
161 return theList;
162
163}
164//
165class StringMatch : public unary_function<string,bool> {
166public:
167 StringMatch(const string& pattern): pattern_(pattern){}
168 bool operator()(const string& aStr) const {
169
170
171 int b,e;
172 regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
173
174// cout << "investigate " << aStr << " to find " << pattern_
175// << "[" <<b<<","<<e<<"]"
176// << endl;
177
178
179 if (b != 0) return false;
180 if (e != (int)aStr.size()) return false;
181 return true;
182
183 }
184private:
185 string pattern_;
186};
187//-------------------------------------------------------
188//Rebin in frequence + compute mean and sigma
189void reduceSpectra(const TMatrix<r_4>& specMtxInPut,
190 TMatrix<r_4>& meanMtx,
191 TMatrix<r_4>& sigmaMtx) {
192 sa_size_t nSliceFreq = para.nSliceInFreq_;
193 sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq;
194 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
195 sa_size_t freqLow= iSlice*deltaFreq;
196 sa_size_t freqHigh= freqLow + deltaFreq -1;
197 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
198 TVector<r_4> reducedRow;
199 reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
200 double mean;
201 double sigma;
202 MeanSigma(reducedRow,mean,sigma);
203 meanMtx(iCh,iSlice) = mean;
204 sigmaMtx(iCh,iSlice) = sigma;
205 }//eo loop on channels
206 }//eo loop on slices
207}
208//-------------------------------------------------------
209//Compute the mean of Raw spectra and also the mean/sigma of rebinned spectra
210//Used like:
211//
212void meanOnCycles() throw(string) {
213 list<string> listOfFiles;
214 string directoryName;
215 directoryName = para.inPath_ + "/" + para.sourceName_;
216
217 //Make the listing of the directory
218 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
219
220 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
221 iFileEnd = listOfFiles.end();
222
223 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
224 TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
225 uint_4 nSpectra=0;
226 //Loop on files
227 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
228 if (para.debuglev_>90){
229 cout << "load file <" << *iFile << ">" << endl;
230 }
231 PInPersist fin(*iFile);
232 vector<string> vec = fin.GetNameTags();
233 list<string> listOfSpectra;
234 //Keep only required PPF objects
235 std::remove_copy_if(
236 vec.begin(), vec.end(), back_inserter(listOfSpectra),
237 not1(match)
238 );
239
240 iSpecEnd = listOfSpectra.end();
241 listOfSpectra.sort(stringCompare);
242 //Loop of spectra matrix
243 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
244 if (para.debuglev_>90){
245 cout << " spactra <" << *iSpec << ">" << endl;
246 }
247 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
248 fin.GetObject(aSpec,*iSpec);
249 //How to see if the GetObject is ok?? Ask Reza
250 nSpectra++;
251 meanOfSpectra+=aSpec;
252 }//eo loop on spectra in a file
253 }//eo loop on files
254
255 //Normalisation
256 if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
257
258 //Compute the reduced version of the mean and sigma
259 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
260 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
261 reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
262
263 {//Save the result
264 stringstream tmp;
265 tmp << nSpectra;
266 string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
267 cout << "Save mean based on " << nSpectra << " cycles " << endl;
268 POutPersist fos(fileName);
269
270 string tag = "mean";
271 fos << PPFNameTag(tag) << meanOfSpectra;
272 tag = "meanred";
273 fos << PPFNameTag(tag) << meanRedMtx;
274 tag = "sigmared";
275 fos << PPFNameTag(tag) << sigmaRedMtx;
276 }
277}
278//-------------------------------------------------------
279//Compute the median of Raw spectra and also the mean/sigma of rebinned spectra
280//Used like:
281//
282void medianOnCycles() throw(string) {
283 list<string> listOfFiles;
284 string directoryName;
285 directoryName = para.inPath_ + "/" + para.sourceName_;
286
287 //Make the listing of the directory
288 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
289
290 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
291 iFileEnd = listOfFiles.end();
292
293
294 TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough...
295
296 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
297 uint_4 nSpectra=0;
298 //Loop on files
299 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
300 if (para.debuglev_>90){
301 cout << "load file <" << *iFile << ">" << endl;
302 }
303 PInPersist fin(*iFile);
304 vector<string> vec = fin.GetNameTags();
305 list<string> listOfSpectra;
306 //Keep only required PPF objects
307 std::remove_copy_if(
308 vec.begin(), vec.end(), back_inserter(listOfSpectra),
309 not1(match)
310 );
311
312 iSpecEnd = listOfSpectra.end();
313 listOfSpectra.sort(stringCompare);
314 //Loop of spectra matrix
315 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ; ++iSpec){
316 if (para.debuglev_>90){
317 cout << " spactra <" << *iSpec << ">" << endl;
318 }
319 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
320 fin.GetObject(aSpec,*iSpec);
321
322 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
323
324 nSpectra++;
325 }//eo loop on spectra in a file
326 }//eo loop on files
327
328
329
330 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
331 //Compute the median for each freq. and Channel
332 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
333 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
334 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
335 vector<r_4> tmp;
336 tmp0.FillTo(tmp);
337 medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
338 }
339 }
340
341
342 //Compute the reduced version of the mean and sigma
343 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
344 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
345 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
346
347
348 sa_size_t f1320=freqTochan(1320.);
349 sa_size_t f1345=freqTochan(1345.);
350 sa_size_t f1355=freqTochan(1355.);
351 sa_size_t f1380=freqTochan(1380.);
352 //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
353 if (para.debuglev_>9){
354 cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
355 }
356 TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
357 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
358 double meanMed1;
359 double sigmaMed1;
360 TVector<r_4> band1;
361 band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
362 MeanSigma(band1,meanMed1,sigmaMed1);
363 double meanMed2;
364 double sigmaMed2;
365 TVector<r_4> band2;
366 band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
367 MeanSigma(band2,meanMed2,sigmaMed2);
368 meanMed(iCh) = (meanMed1+meanMed2)*0.5;
369 }
370 meanMed.Print(cout);
371
372 //Compute the sigma in the range 1320MHz-1380MHz
373 if (para.debuglev_>9){
374 cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
375 }
376 TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
377 sa_size_t redf1320=(1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
378 sa_size_t redf1380=(1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_;
379
380 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
381 double meanSigma;
382 double sigmaSigma;
383 TVector<r_4> band;
384 band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
385 MeanSigma(band,meanSigma,sigmaSigma);
386 meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
387 sigmaMed(iCh) = meanSigma;
388 }
389 sigmaMed.Print(cout);
390
391
392 if (para.debuglev_>9){
393 cout << "Compute medianOfSpectraNorm" << endl;
394 }
395 TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
396 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
397 medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
398 medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
399 }
400
401
402
403 {//Save the result
404 stringstream tmp;
405 tmp << nSpectra;
406 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
407 cout << "Save median based on " << nSpectra << " cycles " << endl;
408 POutPersist fos(fileName);
409
410 string tag = "median";
411 fos << PPFNameTag(tag) << medianOfSpectra;
412
413 tag = "medianNorm";
414 fos << PPFNameTag(tag) << medianOfSpectraNorm;
415
416
417 tag = "meanmedred";
418 fos << PPFNameTag(tag) << meanRedMtx;
419 tag = "sigmamedred";
420 fos << PPFNameTag(tag) << sigmaRedMtx;
421 tag = "cycleVsfreq";
422
423 TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
424 fos << PPFNameTag(tag) << tarr;
425 }
426}
427
428//-------------------------------------------------------
429int main(int narg, char* arg[]) {
430
431 int rc = 0; //return code
432 string msg; //message used in Exceptions
433
434
435
436 //default value for initial parameters (see Para structure on top of the file)
437 string debuglev = "0";
438 string action = "meanDiffOnOff";
439 string inputPath = ".";
440 string outputPath = ".";
441 string sourceName = "Abell85";
442 string ppfFile;
443 string nSliceInFreq = "32";
444
445
446 // bool okarg=false;
447 int ka=1;
448 while (ka<(narg-1)) {
449 if (strcmp(arg[ka],"-h")==0) {
450 cout << "Usage: -act [meanDiffOnOff]|medianDiffOnOff\n"
451 << " -src [Abell85]\n -inPath [.]|<top_root_dir of the ppf file>\n"
452 << " (ex. /sps/baoradio/AmasNancay/JEC/\n "
453 << " -outPath [.]|<dir of the output> \n"
454 << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
455 << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)"
456 << " -debug <level> "
457 << endl;
458 return 0;
459 }
460 else if (strcmp(arg[ka],"-debug")==0) {
461 debuglev=arg[ka+1];
462 ka+=2;
463 }
464 else if (strcmp(arg[ka],"-act")==0) {
465 action=arg[ka+1];
466 ka+=2;
467 }
468 else if (strcmp(arg[ka],"-inPath")==0) {
469 inputPath=arg[ka+1];
470 ka+=2;
471 }
472 else if (strcmp(arg[ka],"-outPath")==0) {
473 outputPath=arg[ka+1];
474 ka+=2;
475 }
476 else if (strcmp(arg[ka],"-src")==0) {
477 sourceName=arg[ka+1];
478 ka+=2;
479 }
480 else if (strcmp(arg[ka],"-ppfFile")==0) {
481 ppfFile=arg[ka+1];
482 ka+=2;
483 }
484 else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
485 nSliceInFreq=arg[ka+1];
486 ka+=2;
487 }
488 else ka++;
489 }//eo while
490
491 para.debuglev_ = atoi(debuglev.c_str());
492 para.inPath_ = inputPath;
493 para.outPath_ = outputPath;
494 para.sourceName_ = sourceName;
495 para.ppfFile_ = ppfFile;
496 para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
497
498 cout << "Dump Initial parameters ............" << endl;
499 cout << " action = " << action << "\n"
500 << " inputPath = " << para.inPath_ << "\n"
501 << " outputPath = " << para.outPath_ << "\n"
502 << " sourceName = " << para.sourceName_ << "\n"
503 << " ppfFile = " << para.ppfFile_ << "\n"
504 << " nSliceInFreq = " << para.nSliceInFreq_ << "\n"
505 << " debuglev = " << para.debuglev_ << "\n";
506 cout << "...................................." << endl;
507
508 if ( "" == ppfFile ) {
509 cerr << "mergeRawOnOff.cc: you have forgotten ppfFile option"
510 << endl;
511 return 999;
512 }
513
514
515 try {
516
517// int b,e;
518// char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
519// printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
520
521 if ( action == "meanDiffOnOff" ) {
522 meanOnCycles();
523 } else if (action == "medianDiffOnOff") {
524 medianOnCycles();
525 } else {
526 msg = "Unknown action " + action;
527 throw msg;
528 }
529
530
531 } catch (std::exception& sex) {
532 cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name()
533 << "\n msg= " << sex.what() << endl;
534 rc = 78;
535 }
536 catch ( string str ) {
537 cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
538 }
539 catch (...) {
540 cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl;
541 rc = 79;
542 }
543
544 return 0;
545
546}
Note: See TracBrowser for help on using the repository browser.