source: BAORadio/AmasNancay/v3/mergeAnaFiles.cc@ 654

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

improve mean ON-OFF with Calib BAO (jec)
File size: 35.7 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//
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 meanCalibBAODiffOnOffCycles() 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 char* onffTupleName[7]={"cycle",
270 "onoffRaw0","onoffRaw1",
271 "onoffRun0","onoffRun1",
272 "onoffCycle0","onoffCycle1"};
273 NTuple onoffevolution(7,onffTupleName);
274 r_4 xnt[7];
275
276 //Lower and Higher freq. bin to perform mean follow up
277 sa_size_t chLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
278 sa_size_t chHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
279 if (para.debuglev_>90){
280 cout << "freq. band for follow up [" << chLow << ", "<< chHigh << "]" << endl;
281 }
282
283 //Loop on files/run
284
285 int totalNumberCycles=0; //total number of cycles for normalisation
286 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
287 if (para.debuglev_>90){
288 cout << "load file <" << *iFile << ">" << endl;
289 }
290
291 vector<string> tokens;
292 split(*iFile,"_",tokens);
293 string dateOfRun = tokens[1];
294 if (para.debuglev_>90){
295 cout << "date <" << dateOfRun << ">" << endl;
296 }
297 vector<string> tokens2;
298 split(tokens[2],".",tokens2);
299 string srcLower = tokens2[0];
300
301
302
303 PInPersist fin(*iFile);
304 vector<string> vec = fin.GetNameTags();
305
306 vector<string> modeList;
307 modeList.push_back("On");
308 modeList.push_back("Off");
309 vector<string>::const_iterator iMode;
310
311 map<string, list<int> > cycleModeCollect;
312
313 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
314 list<string> listOfSpectra;
315 //Keep only required PPF objects
316 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
317 std::remove_copy_if(
318 vec.begin(), vec.end(), back_inserter(listOfSpectra),
319 not1(StringMatch(matchstr))
320 );
321
322 listOfSpectra.sort(stringCompare);
323 iSpecEnd = listOfSpectra.end();
324
325 matchstr = "[0-9]+";
326 //Loop of spectra matrix
327 list<int> listOfCycles;
328 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
329 int b,e;
330 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
331 if (para.debuglev_>90){
332 cout << " spactra <" << *iSpec << ">" << endl;
333 cout << " cycle " << iSpec->substr(b) << endl;
334 }
335 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
336 }//end loop spectra
337 cycleModeCollect[*iMode] = listOfCycles;
338 }//end of mode
339
340 //Take the Intersection of the list Of cycles in mode Off and On
341 list<int> commonCycles;
342 set_intersection(cycleModeCollect["On"].begin(),
343 cycleModeCollect["On"].end(),
344 cycleModeCollect["Off"].begin(),
345 cycleModeCollect["Off"].end(),
346 back_inserter(commonCycles)
347 );
348
349 if (para.debuglev_>90){
350 cout << "Liste of cycles common to On & Off: <";
351 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
352 cout << *i << " ";
353 }
354 cout << ">" << endl;
355 }
356
357 //
358 //Load BAO Calibration factors "per Cycle and Channels"
359 //Compute the mean per Cycle to
360 // fill factors "per Run and Channels" with the same cycle list
361 //
362 //
363 //TODO improve the code....
364
365 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
366 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
367 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
368 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
369
370 string calibFileName;
371 ifstream ifs;
372 sa_size_t nr,nc; //values read
373
374 //OFF Cycle per Channel
375 calibFileName = directoryName + "/"
376 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
377 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
378 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
379 ifs.open(calibFileName.c_str());
380 if ( ! ifs.is_open() ) {
381
382 throw calibFileName + " cannot be opened...";
383 }
384 calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
385 if(para.debuglev_>9){
386 cout << "(nr,nc): "<< nr << "," << nc << endl;
387 calibBAOfactors_Off_Cycle_Ch0.Print(cout);
388 cout << endl;
389 }
390
391 TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
392 calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
393 {//Compute the mean
394 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
395 double mean,sigma;
396 MeanSigma(coef,mean,sigma);
397 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
398 }
399 if(para.debuglev_>9){
400 cout << "Fill calib. with mean value " << endl;
401 calibBAOfactors_Off_Run_Ch0.Print(cout);
402 cout << endl;
403 }
404 ifs.close();
405
406 //
407 calibFileName = directoryName + "/"
408 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
409 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
410 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
411 ifs.open(calibFileName.c_str());
412 if ( ! ifs.is_open() ) {
413
414 throw calibFileName + " cannot be opened...";
415 }
416 calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
417 if(para.debuglev_>9){
418 cout << "(nr,nc): "<< nr << "," << nc << endl;
419 calibBAOfactors_Off_Cycle_Ch1.Print(cout);
420 cout << endl;
421 }
422 TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
423 calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
424 {//Compute the mean
425 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
426 double mean,sigma;
427 MeanSigma(coef,mean,sigma);
428 // cout << "Mean: " << mean << " sigma " << sigma << endl;
429 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
430 }
431 if(para.debuglev_>9){
432 cout << "Fill calib. with mean value " << endl;
433 calibBAOfactors_Off_Run_Ch1.Print(cout);
434 cout << endl;
435 }
436 ifs.close();
437
438 //ON Cycle per Channel
439 calibFileName = directoryName + "/"
440 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
441 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
442 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
443 ifs.open(calibFileName.c_str());
444 if ( ! ifs.is_open() ) {
445
446 throw calibFileName + " cannot be opened...";
447 }
448 calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
449 if(para.debuglev_>9){
450 cout << "(nr,nc): "<< nr << "," << nc << endl;
451 calibBAOfactors_On_Cycle_Ch0.Print(cout);
452 cout << endl;
453 }
454
455 TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
456 calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
457 {//Compute the mean
458 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
459 double mean,sigma;
460 MeanSigma(coef,mean,sigma);
461 // cout << "Mean: " << mean << " sigma " << sigma << endl;
462 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
463 }
464 if(para.debuglev_>9){
465 cout << "Fill calib. with mean value " << endl;
466 calibBAOfactors_On_Run_Ch0.Print(cout);
467 cout << endl;
468 }
469 ifs.close();
470
471
472 calibFileName = directoryName + "/"
473 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
474 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
475 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
476 ifs.open(calibFileName.c_str());
477 if ( ! ifs.is_open() ) {
478 throw calibFileName + " cannot be opened...";
479 }
480 calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
481 if(para.debuglev_>9){
482 cout << "(nr,nc): "<< nr << "," << nc << endl;
483 calibBAOfactors_On_Cycle_Ch1.Print(cout);
484 cout << endl;
485 }
486 TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
487 calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
488 {//Compute the mean
489 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
490 double mean,sigma;
491 MeanSigma(coef,mean,sigma);
492 // cout << "Mean: " << mean << " sigma " << sigma << endl;
493 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
494 }
495 if(para.debuglev_>9){
496 cout << "Fill calib. with mean value " << endl;
497 calibBAOfactors_On_Run_Ch1.Print(cout);
498 cout << endl;
499 }
500
501 ifs.close();
502
503 //link <cycle> - <calibration coefficient>
504 //We cannot rely on identical cycle list of the OFF and ON calibration
505 map<int,r_4> calibBAO_Off_Run_Ch0;
506 map<int,r_4> calibBAO_Off_Run_Ch1;
507 map<int,r_4> calibBAO_On_Run_Ch0;
508 map<int,r_4> calibBAO_On_Run_Ch1;
509
510 map<int,r_4> calibBAO_Off_Cycle_Ch0;
511 map<int,r_4> calibBAO_Off_Cycle_Ch1;
512 map<int,r_4> calibBAO_On_Cycle_Ch0;
513 map<int,r_4> calibBAO_On_Cycle_Ch1;
514
515 //per Run based BAO coefficients
516 nr = calibBAOfactors_Off_Run_Ch0.NRows();
517 for (sa_size_t ir=0; ir<nr; ++ir){
518// cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
519// << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
520
521 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
522 = calibBAOfactors_Off_Run_Ch0(ir,1);
523 calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
524 = calibBAOfactors_Off_Cycle_Ch0(ir,1);
525 calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
526 = calibBAOfactors_Off_Run_Ch1(ir,1);
527 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
528 = calibBAOfactors_Off_Cycle_Ch1(ir,1);
529 }//eo loop on coef Off
530
531 nr = calibBAOfactors_On_Run_Ch0.NRows();
532 for (sa_size_t ir=0; ir<nr; ++ir){
533 calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
534 = calibBAOfactors_On_Run_Ch0(ir,1);
535 calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
536 = calibBAOfactors_On_Cycle_Ch0(ir,1);
537 calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
538 = calibBAOfactors_On_Run_Ch1(ir,1);
539 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
540 = calibBAOfactors_On_Cycle_Ch1(ir,1);
541 }//eo loop on coef On
542
543 //Loop on cyles
544 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
545
546 //Look if the cycle has been calibrated...
547 bool isCycleCalibrated =
548 ( calibBAO_On_Run_Ch0.count(*ic) *
549 calibBAO_On_Run_Ch1.count(*ic) *
550 calibBAO_Off_Run_Ch0.count(*ic) *
551 calibBAO_Off_Run_Ch1.count(*ic) *
552 calibBAO_On_Cycle_Ch0.count(*ic) *
553 calibBAO_On_Cycle_Ch1.count(*ic) *
554 calibBAO_Off_Cycle_Ch0.count(*ic) *
555 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
556
557 if(para.debuglev_>9) {
558 cout << "Calibration coefficients for cycle "<<*ic << endl;
559 if (isCycleCalibrated) {
560 cout << "Cycle calibrated " << endl;
561 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
562 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
563 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
564 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
565 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
566 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
567 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
568 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
569 } else {
570 cout << "Cycle NOT calibrated " << endl;
571 }
572 }//debug
573
574
575 if ( ! isCycleCalibrated ) continue;
576
577 string ppftag;
578 //load ON phase
579 stringstream cycle;
580 cycle << *ic;
581
582 ppftag = "specRawOn"+cycle.str();
583 TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
584 fin.GetObject(aSpecOn,ppftag);
585
586 TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
587 aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
588 aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
589
590 TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
591 aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
592 aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
593
594 ppftag = "specRawOff"+cycle.str();
595 TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
596 fin.GetObject(aSpecOff,ppftag);
597
598 TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
599 aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
600 aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
601
602 TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
603 aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
604 aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
605
606
607 //Perform the difference ON-OFF with the different calibration options
608 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
609 meanDiffONOFF_noCalib += diffOnOff_noCalib;
610
611 TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
612 meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
613
614 TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
615 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
616
617 //
618 totalNumberCycles++;
619 //Fill NTuple
620 xnt[0] = totalNumberCycles;
621
622 TVector<r_4> meanInRange_noCalib(NUMBER_OF_CHANNELS);
623 meanInRange(diffOnOff_noCalib,chLow,chHigh,meanInRange_noCalib);
624 xnt[1] = meanInRange_noCalib(0);
625 xnt[2] = meanInRange_noCalib(1);
626
627 TVector<r_4> meanInRange_perRunCalib(NUMBER_OF_CHANNELS);
628 meanInRange(diffOnOff_perRunCalib,chLow,chHigh,meanInRange_perRunCalib);
629 xnt[3] = meanInRange_perRunCalib(0);
630 xnt[4] = meanInRange_perRunCalib(1);
631
632 TVector<r_4> meanInRange_perCycleCalib(NUMBER_OF_CHANNELS);
633 meanInRange(diffOnOff_perCycleCalib,chLow,chHigh,meanInRange_perCycleCalib);
634 xnt[5] = meanInRange_perCycleCalib(0);
635 xnt[6] = meanInRange_perCycleCalib(1);
636
637 onoffevolution.Fill(xnt);
638
639 //Quit if enough
640 if (totalNumberCycles >= para.maxNumberCycles_) break;
641
642 }//eo loop on cycles
643 if (totalNumberCycles >= para.maxNumberCycles_) break;
644
645 }//eo loop on spectra in a file
646 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
647 //Normalisation
648 if(totalNumberCycles > 0){
649 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles;
650 meanDiffONOFF_perRunCalib /= (r_4)totalNumberCycles;
651 meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
652 }
653
654 //Compute the reduced version of the mean and sigma
655 TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
656 TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
657 reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
658
659 TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
660 TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
661 reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
662
663 TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
664 TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
665 reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
666
667 {//save the results
668 stringstream tmp;
669 tmp << totalNumberCycles;
670 string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
671
672 POutPersist fos(fileName);
673 cout << "Save output in " << fileName << endl;
674
675 string tag = "meanNoCalib";
676 fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
677 tag = "meanPerRunCalib";
678 fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
679 tag = "meanPerCycleCalib";
680 fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
681
682 tag = "redmeanNoCalib";
683 fos << PPFNameTag(tag) << meanRedMtx_noCalib;
684 tag = "redsigmaNoCalib";
685 fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
686
687 tag = "redmeanPerRunCalib";
688 fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
689 tag = "redsigmaPerRunCalib";
690 fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
691
692 tag = "redmeanPerCycleCalib";
693 fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
694 tag = "redsigmaPerCycleCalib";
695 fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
696
697 tag = "onoffevol";
698 fos << PPFNameTag(tag) << onoffevolution;
699 }//end of save
700}
701//-------------------------------------------------------
702//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
703//Used like:
704//
705void meanRawDiffOnOffCycles() throw(string) {
706 list<string> listOfFiles;
707 string directoryName;
708 directoryName = para.inPath_ + "/" + para.sourceName_;
709
710 //Make the listing of the directory
711 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
712
713 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
714 iFileEnd = listOfFiles.end();
715
716 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
717 TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
718 uint_4 nSpectra=0;
719 //Loop on files
720 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
721 if (para.debuglev_>90){
722 cout << "load file <" << *iFile << ">" << endl;
723 }
724 PInPersist fin(*iFile);
725 vector<string> vec = fin.GetNameTags();
726 list<string> listOfSpectra;
727 //Keep only required PPF objects
728 std::remove_copy_if(
729 vec.begin(), vec.end(), back_inserter(listOfSpectra),
730 not1(match)
731 );
732
733 listOfSpectra.sort(stringCompare);
734 iSpecEnd = listOfSpectra.end();
735 //Loop of spectra matrix
736 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
737 if (para.debuglev_>90){
738 cout << " spactra <" << *iSpec << ">" << endl;
739 }
740 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
741 fin.GetObject(aSpec,*iSpec);
742 //How to see if the GetObject is ok?? Ask Reza
743 nSpectra++;
744 meanOfSpectra+=aSpec;
745 }//eo loop on spectra in a file
746 }//eo loop on files
747
748 //Normalisation
749 if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
750
751 //Compute the reduced version of the mean and sigma
752 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
753 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
754 reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
755
756 {//Save the result
757 stringstream tmp;
758 tmp << nSpectra;
759 string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
760 cout << "Save mean based on " << nSpectra << " cycles " << endl;
761 POutPersist fos(fileName);
762
763 string tag = "mean";
764 fos << PPFNameTag(tag) << meanOfSpectra;
765 tag = "meanred";
766 fos << PPFNameTag(tag) << meanRedMtx;
767 tag = "sigmared";
768 fos << PPFNameTag(tag) << sigmaRedMtx;
769 }
770}
771//-------------------------------------------------------
772//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
773//Used like:
774//
775void medianRawDiffOnOffCycles() 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
787 TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
788
789 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
790 uint_4 nSpectra=0;
791 //Loop on files
792 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
793 if (para.debuglev_>90){
794 cout << "load file <" << *iFile << ">" << endl;
795 }
796 PInPersist fin(*iFile);
797 vector<string> vec = fin.GetNameTags();
798 list<string> listOfSpectra;
799 //Keep only required PPF objects
800 std::remove_copy_if(
801 vec.begin(), vec.end(), back_inserter(listOfSpectra),
802 not1(match)
803 );
804
805 listOfSpectra.sort(stringCompare);
806 iSpecEnd = listOfSpectra.end();
807 //Loop of spectra matrix
808 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ; ++iSpec){
809 if (para.debuglev_>90){
810 cout << " spactra <" << *iSpec << ">" << endl;
811 }
812 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
813 fin.GetObject(aSpec,*iSpec);
814
815 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
816
817 nSpectra++;
818 }//eo loop on spectra in a file
819 }//eo loop on files
820
821
822
823 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
824 //Compute the median for each freq. and Channel
825 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
826 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
827 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
828 vector<r_4> tmp;
829 tmp0.FillTo(tmp);
830 medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
831 }
832 }
833
834
835 //Compute the reduced version of the mean and sigma
836 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
837 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
838 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
839
840
841 sa_size_t f1320=freqToChan(1320.);
842 sa_size_t f1345=freqToChan(1345.);
843 sa_size_t f1355=freqToChan(1355.);
844 sa_size_t f1380=freqToChan(1380.);
845 //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
846 if (para.debuglev_>9){
847 cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
848 }
849 TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
850 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
851 double meanMed1;
852 double sigmaMed1;
853 TVector<r_4> band1;
854 band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
855 MeanSigma(band1,meanMed1,sigmaMed1);
856 double meanMed2;
857 double sigmaMed2;
858 TVector<r_4> band2;
859 band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
860 MeanSigma(band2,meanMed2,sigmaMed2);
861 meanMed(iCh) = (meanMed1+meanMed2)*0.5;
862 }
863 meanMed.Print(cout);
864 cout << endl;
865
866
867 //Compute the sigma in the range 1320MHz-1380MHz
868 if (para.debuglev_>9){
869 cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
870 }
871 TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
872 sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
873 sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
874
875 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
876 double meanSigma;
877 double sigmaSigma;
878 TVector<r_4> band;
879 band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
880 MeanSigma(band,meanSigma,sigmaSigma);
881 meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
882 sigmaMed(iCh) = meanSigma;
883 }
884 sigmaMed.Print(cout);
885 cout << endl;
886
887
888
889 if (para.debuglev_>9){
890 cout << "Compute medianOfSpectraNorm" << endl;
891 }
892 TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
893 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
894 medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
895 medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
896 }
897
898
899
900 {//Save the result
901 stringstream tmp;
902 tmp << nSpectra;
903 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
904 cout << "Save median based on " << nSpectra << " cycles " << endl;
905 POutPersist fos(fileName);
906
907 string tag = "median";
908 fos << PPFNameTag(tag) << medianOfSpectra;
909
910 tag = "medianNorm";
911 fos << PPFNameTag(tag) << medianOfSpectraNorm;
912
913
914 tag = "meanmedred";
915 fos << PPFNameTag(tag) << meanRedMtx;
916 tag = "sigmamedred";
917 fos << PPFNameTag(tag) << sigmaRedMtx;
918 tag = "cycleVsfreq";
919
920 TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
921 fos << PPFNameTag(tag) << tarr;
922 }
923}
924
925//-------------------------------------------------------
926int main(int narg, char* arg[]) {
927
928 int rc = 0; //return code
929 string msg; //message used in Exceptions
930
931
932
933 //default value for initial parameters (see Para structure on top of the file)
934 string debuglev = "0";
935 string action = "meanDiffOnOff";
936 string inputPath = ".";
937 string outputPath = ".";
938 string sourceName = "Abell85";
939 string ppfFile;
940 string nSliceInFreq = "32";
941 string typeOfCalib="perRun";
942 string calibFreq = "1346";
943 string calibBandFreq="6.25";
944 string mxcycles;
945
946 // bool okarg=false;
947 int ka=1;
948 while (ka<narg) {
949 if (strcmp(arg[ka],"-h")==0) {
950 cout << "Usage: -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
951 << " -mxcycles <number> (max. number of cycles to be treated)\n"
952 << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
953 << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
954 << " -src [Abell85]\n -inPath [.]|<top_root_dir of the ppf file>\n"
955 << " (ex. /sps/baoradio/AmasNancay/JEC/\n "
956 << " -outPath [.]|<dir of the output> \n"
957 << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
958 << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
959 << " -debug <level> "
960 << endl;
961 return 0;
962 }
963 else if (strcmp(arg[ka],"-debug")==0) {
964 debuglev=arg[ka+1];
965 ka+=2;
966 }
967 else if (strcmp(arg[ka],"-act")==0) {
968 action=arg[ka+1];
969 ka+=2;
970 }
971 else if (strcmp(arg[ka],"-calibfreq")==0) {
972 calibFreq=arg[ka+1];
973 ka+=2;
974 }
975 else if (strcmp(arg[ka],"-calibbandfreq")==0) {
976 calibBandFreq=arg[ka+1];
977 ka+=2;
978 }
979 else if (strcmp(arg[ka],"-mxcycles")==0) {
980 mxcycles=arg[ka+1];
981 ka+=2;
982 }
983 else if (strcmp(arg[ka],"-inPath")==0) {
984 inputPath=arg[ka+1];
985 ka+=2;
986 }
987 else if (strcmp(arg[ka],"-outPath")==0) {
988 outputPath=arg[ka+1];
989 ka+=2;
990 }
991 else if (strcmp(arg[ka],"-src")==0) {
992 sourceName=arg[ka+1];
993 ka+=2;
994 }
995 else if (strcmp(arg[ka],"-ppfFile")==0) {
996 ppfFile=arg[ka+1];
997 ka+=2;
998 }
999 else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1000 nSliceInFreq=arg[ka+1];
1001 ka+=2;
1002 }
1003 else ka++;
1004 }//eo while
1005
1006 para.debuglev_ = atoi(debuglev.c_str());
1007 para.inPath_ = inputPath;
1008 para.outPath_ = outputPath;
1009 para.sourceName_ = sourceName;
1010 para.ppfFile_ = ppfFile;
1011 para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1012 para.calibFreq_ = calibFreq;
1013 para.calibBandFreq_ = calibBandFreq;
1014 para.rcalibFreq_ = atof(calibFreq.c_str());
1015 para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1016 if (mxcycles != "") {
1017 para.maxNumberCycles_ = atoi(mxcycles.c_str());
1018 } else {
1019 para.maxNumberCycles_ = std::numeric_limits<int>::max();
1020 }
1021
1022 cout << "Dump Initial parameters ............" << endl;
1023 cout << " action = " << action << "\n"
1024 << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1025 << " inputPath = " << para.inPath_ << "\n"
1026 << " outputPath = " << para.outPath_ << "\n"
1027 << " sourceName = " << para.sourceName_ << "\n"
1028 << " ppfFile = " << para.ppfFile_ << "\n"
1029 << " nSliceInFreq = " << para.nSliceInFreq_ << "\n"
1030 << " calibFreq = " << para.calibFreq_ << "\n"
1031 << " calibBandFreq = " << para.calibBandFreq_ << "\n"
1032 << " debuglev = " << para.debuglev_ << "\n";
1033 cout << "...................................." << endl;
1034
1035 if ( "" == ppfFile ) {
1036 cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1037 << endl;
1038 return 999;
1039 }
1040
1041
1042 try {
1043
1044// int b,e;
1045// char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1046// printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1047
1048 if ( action == "meanRawDiffOnOff" ) {
1049 meanRawDiffOnOffCycles();
1050 } else if (action == "medianRawDiffOnOff") {
1051 medianRawDiffOnOffCycles();
1052 } else if (action == "meanCalibBAODiffOnOff") {
1053 meanCalibBAODiffOnOffCycles();
1054 } else {
1055 msg = "Unknown action " + action;
1056 throw msg;
1057 }
1058
1059
1060 } catch (std::exception& sex) {
1061 cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name()
1062 << "\n msg= " << sex.what() << endl;
1063 rc = 78;
1064 }
1065 catch ( string str ) {
1066 cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1067 }
1068 catch (...) {
1069 cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl;
1070 rc = 79;
1071 }
1072
1073 return 0;
1074
1075}
Note: See TracBrowser for help on using the repository browser.