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

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

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