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

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

move to trunk the previous HEAD (jec)
File size: 13.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 = 10000;
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//--------------------------------------------------------------
58char *regexp (const char *string, const char *patrn, int *begin, int *end) {
59 int i, w=0, len;
60 char *word = NULL;
61 regex_t rgT;
62 regmatch_t match;
63 regcomp(&rgT,patrn,REG_EXTENDED);
64 if ((regexec(&rgT,string,1,&match,0)) == 0) {
65 *begin = (int)match.rm_so;
66 *end = (int)match.rm_eo;
67 len = *end-*begin;
68 word=(char*)malloc(len+1);
69 for (i=*begin; i<*end; i++) {
70 word[w] = string[i];
71 w++; }
72 word[w]=0;
73 }
74 regfree(&rgT);
75 return word;
76}
77//-------
78sa_size_t round_sa(r_4 r) {
79 return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
80}
81//-----
82string StringToLower(string strToConvert){
83 //change each element of the string to lower case
84 for(unsigned int i=0;i<strToConvert.length();i++) {
85 strToConvert[i] = tolower(strToConvert[i]);
86 }
87 return strToConvert;//return the converted string
88}
89//-----
90bool stringCompare( const string &left, const string &right ){
91 if( left.size() < right.size() )
92 return true;
93 for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
94 if( tolower( *lit ) < tolower( *rit ) )
95 return true;
96 else if( tolower( *lit ) > tolower( *rit ) )
97 return false;
98 return false;
99}
100//-----
101list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
102 list<string> theList;
103
104
105 DIR *dip;
106 struct dirent *dit;
107 string msg; string fileName;
108 string fullFileName;
109 size_t found;
110
111 if ((dip=opendir(dir.c_str())) == NULL ) {
112 msg = "opendir failed on directory "+dir;
113 throw msg;
114 }
115 while ( (dit = readdir(dip)) != NULL ) {
116 fileName = dit->d_name;
117 found=fileName.find(filePettern);
118 if (found != string::npos) {
119 fullFileName = dir + "/";
120 fullFileName += fileName;
121 theList.push_back(fullFileName);
122 }
123 }//eo while
124 if (closedir(dip) == -1) {
125 msg = "closedir failed on directory "+dir;
126 throw msg;
127 }
128
129 theList.sort(stringCompare);
130
131 return theList;
132
133}
134//
135class StringMatch : public unary_function<string,bool> {
136public:
137 StringMatch(const string& pattern): pattern_(pattern){}
138 bool operator()(const string& aStr) const {
139
140
141 int b,e;
142 regexp(aStr.c_str(),pattern_.c_str(),&b,&e);
143
144// cout << "investigate " << aStr << " to find " << pattern_
145// << "[" <<b<<","<<e<<"]"
146// << endl;
147
148
149 if (b != 0) return false;
150 if (e != (int)aStr.size()) return false;
151 return true;
152
153 }
154private:
155 string pattern_;
156};
157//-------------------------------------------------------
158//Rebin in frequence + compute mean and sigma
159void reduceSpectra(const TMatrix<r_4>& specMtxInPut,
160 TMatrix<r_4>& meanMtx,
161 TMatrix<r_4>& sigmaMtx) {
162 sa_size_t nSliceFreq = para.nSliceInFreq_;
163 sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq;
164 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
165 sa_size_t freqLow= iSlice*deltaFreq;
166 sa_size_t freqHigh= freqLow + deltaFreq -1;
167 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
168 TVector<r_4> reducedRow;
169 reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
170 double mean;
171 double sigma;
172 MeanSigma(reducedRow,mean,sigma);
173 meanMtx(iCh,iSlice) = mean;
174 sigmaMtx(iCh,iSlice) = sigma;
175 }//eo loop on channels
176 }//eo loop on slices
177}
178//-------------------------------------------------------
179//Compute the mean of Raw spectra and also the mean/sigma of rebinned spectra
180//Used like:
181//
182void meanOnCycles() throw(string) {
183 list<string> listOfFiles;
184 string directoryName;
185 directoryName = para.inPath_ + "/" + para.sourceName_;
186
187 //Make the listing of the directory
188 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
189
190 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
191 iFileEnd = listOfFiles.end();
192
193 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
194 TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
195 uint_4 nSpectra=0;
196 //Loop on files
197 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
198 if (para.debuglev_>90){
199 cout << "load file <" << *iFile << ">" << endl;
200 }
201 PInPersist fin(*iFile);
202 vector<string> vec = fin.GetNameTags();
203 list<string> listOfSpectra;
204 //Keep only required PPF objects
205 std::remove_copy_if(
206 vec.begin(), vec.end(), back_inserter(listOfSpectra),
207 not1(match)
208 );
209
210 iSpecEnd = listOfSpectra.end();
211 listOfSpectra.sort(stringCompare);
212 //Loop of spectra matrix
213 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
214 if (para.debuglev_>90){
215 cout << " spactra <" << *iSpec << ">" << endl;
216 }
217 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
218 fin.GetObject(aSpec,*iSpec);
219 //How to see if the GetObject is ok?? Ask Reza
220 nSpectra++;
221 meanOfSpectra+=aSpec;
222 }//eo loop on spectra in a file
223 }//eo loop on files
224
225 //Normalisation
226 if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
227
228 //Compute the reduced version of the mean and sigma
229 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
230 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
231 reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
232
233 {//Save the result
234 stringstream tmp;
235 tmp << nSpectra;
236 string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
237 cout << "Save mean based on " << nSpectra << " cycles " << endl;
238 POutPersist fos(fileName);
239
240 string tag = "mean";
241 fos << PPFNameTag(tag) << meanOfSpectra;
242 tag = "meanred";
243 fos << PPFNameTag(tag) << meanRedMtx;
244 tag = "sigmared";
245 fos << PPFNameTag(tag) << sigmaRedMtx;
246 }
247}
248//-------------------------------------------------------
249//Compute the median of Raw spectra and also the mean/sigma of rebinned spectra
250//Used like:
251//
252void medianOnCycles() throw(string) {
253 list<string> listOfFiles;
254 string directoryName;
255 directoryName = para.inPath_ + "/" + para.sourceName_;
256
257 //Make the listing of the directory
258 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
259
260 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
261 iFileEnd = listOfFiles.end();
262
263
264 TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough...
265
266 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
267 uint_4 nSpectra=0;
268 //Loop on files
269 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
270 if (para.debuglev_>90){
271 cout << "load file <" << *iFile << ">" << endl;
272 }
273 PInPersist fin(*iFile);
274 vector<string> vec = fin.GetNameTags();
275 list<string> listOfSpectra;
276 //Keep only required PPF objects
277 std::remove_copy_if(
278 vec.begin(), vec.end(), back_inserter(listOfSpectra),
279 not1(match)
280 );
281
282 iSpecEnd = listOfSpectra.end();
283 listOfSpectra.sort(stringCompare);
284 //Loop of spectra matrix
285 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
286 if (para.debuglev_>90){
287 cout << " spactra <" << *iSpec << ">" << endl;
288 }
289 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
290 fin.GetObject(aSpec,*iSpec);
291
292 if ((sa_size_t)nSpectra < TOTAL_NUM_CYCLES ) {
293 //STORE THE spectra
294 // tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)).Show();
295 // aSpec.Show();
296 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
297// aSpec.Show();
298// cout << "deb 1" << endl;
299// tmp = aSpec;
300// cout << "deb 2" << endl;
301 } else {
302 stringstream tmp;
303 tmp<<nSpectra;
304 throw "FATAL: SIZE ERROR: too many cycles:"+tmp.str();
305 }
306 nSpectra++;
307 }//eo loop on spectra in a file
308 }//eo loop on files
309
310 //Resize
311
312 nSpectra--;//For further selection in Range
313 if (para.debuglev_>90){
314 cout << "Cycles < 0," << nSpectra << ">" << endl;
315 }
316
317
318 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
319 //Compute the median for each freq. and Channel
320 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
321 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
322 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra)).CompactAllDimensions());
323 vector<r_4> tmp;
324 tmp0.FillTo(tmp);
325 nth_element(tmp.begin(),tmp.begin()+tmp.size()/2,tmp.end());
326 medianOfSpectra(iCh,freq) = *(tmp.begin()+tmp.size()/2);
327 }
328 }
329
330
331 //Compute the reduced version of the mean and sigma
332 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
333 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
334 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
335
336 {//Save the result
337 stringstream tmp;
338 tmp << nSpectra;
339 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
340 cout << "Save median based on " << (nSpectra+1) << " cycles " << endl;
341 POutPersist fos(fileName);
342
343 string tag = "median";
344 fos << PPFNameTag(tag) << medianOfSpectra;
345 tag = "meanmedred";
346 fos << PPFNameTag(tag) << meanRedMtx;
347 tag = "sigmamedred";
348 fos << PPFNameTag(tag) << sigmaRedMtx;
349 }
350}
351
352//-------------------------------------------------------
353int main(int narg, char* arg[]) {
354
355 int rc = 0; //return code
356 string msg; //message used in Exceptions
357
358
359
360 //default value for initial parameters (see Para structure on top of the file)
361 string debuglev = "0";
362 string action = "meanDiffOnOff";
363 string inputPath = ".";
364 string outputPath = ".";
365 string sourceName = "Abell85";
366 string ppfFile;
367 string nSliceInFreq = "32";
368
369
370 // bool okarg=false;
371 int ka=1;
372 while (ka<(narg-1)) {
373 if (strcmp(arg[ka],"-debug")==0) {
374 debuglev=arg[ka+1];
375 ka+=2;
376 }
377 else if (strcmp(arg[ka],"-act")==0) {
378 action=arg[ka+1];
379 ka+=2;
380 }
381 else if (strcmp(arg[ka],"-inPath")==0) {
382 inputPath=arg[ka+1];
383 ka+=2;
384 }
385 else if (strcmp(arg[ka],"-outPath")==0) {
386 outputPath=arg[ka+1];
387 ka+=2;
388 }
389 else if (strcmp(arg[ka],"-src")==0) {
390 sourceName=arg[ka+1];
391 ka+=2;
392 }
393 else if (strcmp(arg[ka],"-ppfFile")==0) {
394 ppfFile=arg[ka+1];
395 ka+=2;
396 }
397 else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
398 nSliceInFreq=arg[ka+1];
399 ka+=2;
400 }
401 else ka++;
402 }//eo while
403
404 para.debuglev_ = atoi(debuglev.c_str());
405 para.inPath_ = inputPath;
406 para.outPath_ = outputPath;
407 para.sourceName_ = sourceName;
408 para.ppfFile_ = ppfFile;
409 para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
410
411 cout << "Dump Initial parameters ............" << endl;
412 cout << " action = " << action << "\n"
413 << " inputPath = " << para.inPath_ << "\n"
414 << " outputPath = " << para.outPath_ << "\n"
415 << " sourceName = " << para.sourceName_ << "\n"
416 << " ppfFile = " << para.ppfFile_ << "\n"
417 << " nSliceInFreq = " << para.nSliceInFreq_ << "\n"
418 << " debuglev = " << para.debuglev_ << "\n";
419 cout << "...................................." << endl;
420
421 if ( "" == ppfFile ) {
422 cerr << "mergeRawOnOff.cc: you have forgotten ppfFile option"
423 << endl;
424 return 999;
425 }
426
427
428 try {
429
430// int b,e;
431// char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
432// printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
433
434 if ( action == "meanDiffOnOff" ) {
435 meanOnCycles();
436 } else if (action == "medianDiffOnOff") {
437 medianOnCycles();
438 } else {
439 msg = "Unknown action " + action;
440 throw msg;
441 }
442
443
444 } catch (std::exception& sex) {
445 cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name()
446 << "\n msg= " << sex.what() << endl;
447 rc = 78;
448 }
449 catch ( string str ) {
450 cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
451 }
452 catch (...) {
453 cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl;
454 rc = 79;
455 }
456
457 return 0;
458
459}
Note: See TracBrowser for help on using the repository browser.