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

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

suivi du ON et OFF separement (jec)
File size: 48.1 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 ,"offRaw0f14001415","offRaw1f14001415"
308 ,"onRaw0f14001415","onRaw1f14001415"
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 ch1415 = freqToChan(1415);
329 sa_size_t ch1420 = freqToChan(1420);
330
331 if (para.debuglev_>0){
332 cout << "freq. band for follow up [" << chCalibLow << ", "<< chCalibHigh << "]" << endl;
333 cout << "freq. band for follow up [" << ch1420Low << ", "<< ch1420High << "]" << endl;
334 cout << "freq. band for follow up [" << ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
335 cout << "freq. band for follow up [" << ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
336 }
337
338 //Loop on files/run
339
340 int totalNumberCycles=0; //total number of cycles for normalisation
341 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
342 if (para.debuglev_>90){
343 cout << "load file <" << *iFile << ">" << endl;
344 }
345
346 vector<string> tokens;
347 split(*iFile,"_",tokens);
348 string dateOfRun = tokens[1];
349 if (para.debuglev_>90){
350 cout << "date <" << dateOfRun << ">" << endl;
351 }
352 vector<string> tokens2;
353 split(tokens[2],".",tokens2);
354 string srcLower = tokens2[0];
355
356
357
358 PInPersist fin(*iFile);
359 vector<string> vec = fin.GetNameTags();
360
361 vector<string> modeList;
362 modeList.push_back("On");
363 modeList.push_back("Off");
364 vector<string>::const_iterator iMode;
365
366 map<string, list<int> > cycleModeCollect;
367
368 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
369 list<string> listOfSpectra;
370 //Keep only required PPF objects
371 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
372 std::remove_copy_if(
373 vec.begin(), vec.end(), back_inserter(listOfSpectra),
374 not1(StringMatch(matchstr))
375 );
376
377 listOfSpectra.sort(stringCompare);
378 iSpecEnd = listOfSpectra.end();
379
380 matchstr = "[0-9]+";
381 //Loop of spectra matrix
382 list<int> listOfCycles;
383 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
384 int b,e;
385 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
386 if (para.debuglev_>90){
387 cout << " spactra <" << *iSpec << ">" << endl;
388 cout << " cycle " << iSpec->substr(b) << endl;
389 }
390 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
391 }//end loop spectra
392 cycleModeCollect[*iMode] = listOfCycles;
393 }//end of mode
394
395 //Take the Intersection of the list Of cycles in mode Off and On
396 list<int> commonCycles;
397 set_intersection(cycleModeCollect["On"].begin(),
398 cycleModeCollect["On"].end(),
399 cycleModeCollect["Off"].begin(),
400 cycleModeCollect["Off"].end(),
401 back_inserter(commonCycles)
402 );
403
404 if (para.debuglev_>90){
405 cout << "Liste of cycles common to On & Off: <";
406 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
407 cout << *i << " ";
408 }
409 cout << ">" << endl;
410 }
411
412 //
413 //Load BAO Calibration factors "per Cycle and Channels"
414 //Compute the mean per Cycle to
415 // fill factors "per Run and Channels" with the same cycle list
416 //
417 //
418 //TODO improve the code....
419
420 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
421 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
422 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
423 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
424
425 string calibFileName;
426 ifstream ifs;
427 sa_size_t nr,nc; //values read
428
429 //OFF Cycle per Channel
430 calibFileName = directoryName + "/"
431 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
432 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
433 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
434 ifs.open(calibFileName.c_str());
435 if ( ! ifs.is_open() ) {
436
437 throw calibFileName + " cannot be opened...";
438 }
439 calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
440 if(para.debuglev_>9){
441 cout << "(nr,nc): "<< nr << "," << nc << endl;
442 calibBAOfactors_Off_Cycle_Ch0.Print(cout);
443 cout << endl;
444 }
445
446 TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
447 calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
448 {//Compute the mean
449 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
450 double mean,sigma;
451 MeanSigma(coef,mean,sigma);
452 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
453 }
454 if(para.debuglev_>9){
455 cout << "Fill calib. with mean value " << endl;
456 calibBAOfactors_Off_Run_Ch0.Print(cout);
457 cout << endl;
458 }
459 ifs.close();
460
461 //
462 calibFileName = directoryName + "/"
463 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
464 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
465 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
466 ifs.open(calibFileName.c_str());
467 if ( ! ifs.is_open() ) {
468
469 throw calibFileName + " cannot be opened...";
470 }
471 calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
472 if(para.debuglev_>9){
473 cout << "(nr,nc): "<< nr << "," << nc << endl;
474 calibBAOfactors_Off_Cycle_Ch1.Print(cout);
475 cout << endl;
476 }
477 TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
478 calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
479 {//Compute the mean
480 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
481 double mean,sigma;
482 MeanSigma(coef,mean,sigma);
483 // cout << "Mean: " << mean << " sigma " << sigma << endl;
484 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
485 }
486 if(para.debuglev_>9){
487 cout << "Fill calib. with mean value " << endl;
488 calibBAOfactors_Off_Run_Ch1.Print(cout);
489 cout << endl;
490 }
491 ifs.close();
492
493 //ON Cycle per Channel
494 calibFileName = directoryName + "/"
495 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
496 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
497 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
498 ifs.open(calibFileName.c_str());
499 if ( ! ifs.is_open() ) {
500
501 throw calibFileName + " cannot be opened...";
502 }
503 calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
504 if(para.debuglev_>9){
505 cout << "(nr,nc): "<< nr << "," << nc << endl;
506 calibBAOfactors_On_Cycle_Ch0.Print(cout);
507 cout << endl;
508 }
509
510 TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
511 calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
512 {//Compute the mean
513 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
514 double mean,sigma;
515 MeanSigma(coef,mean,sigma);
516 // cout << "Mean: " << mean << " sigma " << sigma << endl;
517 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
518 }
519 if(para.debuglev_>9){
520 cout << "Fill calib. with mean value " << endl;
521 calibBAOfactors_On_Run_Ch0.Print(cout);
522 cout << endl;
523 }
524 ifs.close();
525
526
527 calibFileName = directoryName + "/"
528 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
529 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
530 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
531 ifs.open(calibFileName.c_str());
532 if ( ! ifs.is_open() ) {
533 throw calibFileName + " cannot be opened...";
534 }
535 calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
536 if(para.debuglev_>9){
537 cout << "(nr,nc): "<< nr << "," << nc << endl;
538 calibBAOfactors_On_Cycle_Ch1.Print(cout);
539 cout << endl;
540 }
541 TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
542 calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
543 {//Compute the mean
544 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
545 double mean,sigma;
546 MeanSigma(coef,mean,sigma);
547 // cout << "Mean: " << mean << " sigma " << sigma << endl;
548 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
549 }
550 if(para.debuglev_>9){
551 cout << "Fill calib. with mean value " << endl;
552 calibBAOfactors_On_Run_Ch1.Print(cout);
553 cout << endl;
554 }
555
556 ifs.close();
557
558 //link <cycle> - <calibration coefficient>
559 //We cannot rely on identical cycle list of the OFF and ON calibration
560 map<int,r_4> calibBAO_Off_Run_Ch0;
561 map<int,r_4> calibBAO_Off_Run_Ch1;
562 map<int,r_4> calibBAO_On_Run_Ch0;
563 map<int,r_4> calibBAO_On_Run_Ch1;
564
565 map<int,r_4> calibBAO_Off_Cycle_Ch0;
566 map<int,r_4> calibBAO_Off_Cycle_Ch1;
567 map<int,r_4> calibBAO_On_Cycle_Ch0;
568 map<int,r_4> calibBAO_On_Cycle_Ch1;
569
570 //per Run based BAO coefficients
571 nr = calibBAOfactors_Off_Run_Ch0.NRows();
572 for (sa_size_t ir=0; ir<nr; ++ir){
573// cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
574// << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
575
576 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
577 = calibBAOfactors_Off_Run_Ch0(ir,1);
578 calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
579 = calibBAOfactors_Off_Cycle_Ch0(ir,1);
580 calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
581 = calibBAOfactors_Off_Run_Ch1(ir,1);
582 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
583 = calibBAOfactors_Off_Cycle_Ch1(ir,1);
584 }//eo loop on coef Off
585
586 nr = calibBAOfactors_On_Run_Ch0.NRows();
587 for (sa_size_t ir=0; ir<nr; ++ir){
588 calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
589 = calibBAOfactors_On_Run_Ch0(ir,1);
590 calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
591 = calibBAOfactors_On_Cycle_Ch0(ir,1);
592 calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
593 = calibBAOfactors_On_Run_Ch1(ir,1);
594 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
595 = calibBAOfactors_On_Cycle_Ch1(ir,1);
596 }//eo loop on coef On
597
598 //Loop on cyles
599 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
600
601 //Look if the cycle has been calibrated...
602 bool isCycleCalibrated =
603 ( calibBAO_On_Run_Ch0.count(*ic) *
604 calibBAO_On_Run_Ch1.count(*ic) *
605 calibBAO_Off_Run_Ch0.count(*ic) *
606 calibBAO_Off_Run_Ch1.count(*ic) *
607 calibBAO_On_Cycle_Ch0.count(*ic) *
608 calibBAO_On_Cycle_Ch1.count(*ic) *
609 calibBAO_Off_Cycle_Ch0.count(*ic) *
610 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
611
612 if(para.debuglev_>9) {
613 cout << "Calibration coefficients for cycle "<<*ic << endl;
614 if (isCycleCalibrated) {
615 cout << "Cycle calibrated " << endl;
616 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
617 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
618 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
619 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
620 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
621 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
622 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
623 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
624 } else {
625 cout << "Cycle NOT calibrated " << endl;
626 }
627 }//debug
628
629
630 if ( ! isCycleCalibrated ) continue;
631
632 string ppftag;
633 //load ON phase
634 stringstream cycle;
635 cycle << *ic;
636
637 ppftag = "specRawOn"+cycle.str();
638 TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
639 fin.GetObject(aSpecOn,ppftag);
640
641 TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
642 aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
643 aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
644
645 TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
646 aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
647 aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
648
649 //Load OFF phase
650 ppftag = "specRawOff"+cycle.str();
651 TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
652 fin.GetObject(aSpecOff,ppftag);
653
654 TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
655 aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
656 aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
657
658 TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
659 aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
660 aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
661
662
663 //Perform the difference ON-OFF with the different calibration options
664 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
665 meanDiffONOFF_noCalib += diffOnOff_noCalib;
666
667 //JEC 29/10/11 add ON-OFF/OFF
668 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
669 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
670 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
671 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
672
673 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
674 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
675
676 TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
677 meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
678
679 TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
680 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
681
682 //
683 totalNumberCycles++;
684 //Fill NTuple
685 xnt[0] = totalNumberCycles;
686
687 //Follow up arround the Calibration Freq.
688 TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
689 meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
690 xnt[1] = meanInRange_CalibFreq_noCalib(0);
691 xnt[2] = meanInRange_CalibFreq_noCalib(1);
692
693 TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
694 meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
695 xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
696 xnt[4] = meanInRange_CalibFreq_perRunCalib(1);
697
698 TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
699 meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
700 xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
701 xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);
702
703
704 //Follow up arround the 1420.4MHz Freq.
705 TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
706 meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
707 xnt[7] = meanInRange_1420Freq_noCalib(0);
708 xnt[8] = meanInRange_1420Freq_noCalib(1);
709
710 TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
711 meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
712 xnt[9] = meanInRange_1420Freq_perRunCalib(0);
713 xnt[10] = meanInRange_1420Freq_perRunCalib(1);
714
715 TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
716 meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
717 xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
718 xnt[12] = meanInRange_1420Freq_perCycleCalib(1);
719
720
721 //Follow up below the 1420.4MHz Freq.
722 TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
723 meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
724 TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
725 meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
726
727 xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
728 xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
729
730
731 TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
732 meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
733 TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
734 meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);
735
736 xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
737 xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;
738
739
740 TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
741 meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
742 TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
743 meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);
744
745 xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
746 xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
747
748
749 //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
750 TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
751 meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
752 xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
753 xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
754
755 //JEC 18/11/11 follow up the 1400-1420MHz OFF only
756 TVector<r_4> meanInRange_1400a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
757 meanInRange(aSpecOff,ch1400,ch1415,meanInRange_1400a1415Freq_OFF_noCalib);
758
759 xnt[21] = meanInRange_1400a1415Freq_OFF_noCalib(0);
760 xnt[22] = meanInRange_1400a1415Freq_OFF_noCalib(1);
761
762 TVector<r_4> meanInRange_1400a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
763 meanInRange(aSpecOn,ch1400,ch1415,meanInRange_1400a1415Freq_ON_noCalib);
764
765 xnt[23] = meanInRange_1400a1415Freq_ON_noCalib(0);
766 xnt[24] = meanInRange_1400a1415Freq_ON_noCalib(1);
767
768 //store infos to Ntuple
769 onoffevolution.Fill(xnt);
770
771 //Quit if enough
772 if (totalNumberCycles >= para.maxNumberCycles_) break;
773
774 }//eo loop on cycles
775 if (totalNumberCycles >= para.maxNumberCycles_) break;
776
777 }//eo loop on spectra in a file
778 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
779 //Normalisation
780 if(totalNumberCycles > 0){
781 //JEC 29/10 add ON-OFF/OFF
782 meanDiffONOFFovOFF_noCalib /= (r_4)totalNumberCycles;
783 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles;
784 meanDiffONOFF_perRunCalib /= (r_4)totalNumberCycles;
785 meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
786 }
787
788 //Compute the reduced version of the mean and sigma
789 TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
790 TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
791 reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
792
793 TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
794 TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
795 reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
796
797 TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
798 TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
799 reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
800
801 {//save the results
802 stringstream tmp;
803 tmp << totalNumberCycles;
804 string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
805
806 POutPersist fos(fileName);
807 cout << "Save output in " << fileName << endl;
808
809 string tag = "meanNoCalib";
810 fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
811
812 //JEC 29/10/11
813 tag = "meanOvOffNoCalib";
814 fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
815
816 tag = "meanPerRunCalib";
817 fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
818 tag = "meanPerCycleCalib";
819 fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
820
821 tag = "redmeanNoCalib";
822 fos << PPFNameTag(tag) << meanRedMtx_noCalib;
823 tag = "redsigmaNoCalib";
824 fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
825
826 tag = "redmeanPerRunCalib";
827 fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
828 tag = "redsigmaPerRunCalib";
829 fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
830
831 tag = "redmeanPerCycleCalib";
832 fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
833 tag = "redsigmaPerCycleCalib";
834 fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
835
836 tag = "onoffevol";
837 fos << PPFNameTag(tag) << onoffevolution;
838 }//end of save
839}
840//JEC 14/11/11 New meanRawDiffOnOffCycles START
841//-------------------------------------------------------
842//Compute the mean of Diff ON-OFF/OFF from Raw spectra
843//Used like:
844//
845void meanRawDiffOnOffCycles() throw(string) {
846 list<string> listOfFiles;
847 string directoryName;
848 directoryName = para.inPath_ + "/" + para.sourceName_;
849
850 //Make the listing of the directory
851 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
852
853 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
854 iFileEnd = listOfFiles.end();
855
856 //mean ON-OFF over the list of cycles
857 TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/GAIN
858 TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
859 TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // ON/Filtered_OFF
860 TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
861
862 int totalNumberCycles=0; //total number of cycles
863
864 //Loop on files
865 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
866 if (para.debuglev_>90){
867 }
868 PInPersist fin(*iFile);
869 vector<string> vec = fin.GetNameTags();
870
871 vector<string> modeList;
872 modeList.push_back("On");
873 modeList.push_back("Off");
874 vector<string>::const_iterator iMode;
875
876 map<string, list<int> > cycleModeCollect;
877
878 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
879 list<string> listOfSpectra;
880 //Keep only required PPF objects
881 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
882 std::remove_copy_if(
883 vec.begin(), vec.end(), back_inserter(listOfSpectra),
884 not1(StringMatch(matchstr))
885 );
886
887 listOfSpectra.sort(stringCompare);
888 iSpecEnd = listOfSpectra.end();
889
890 matchstr = "[0-9]+";
891 //Loop of spectra matrix
892 list<int> listOfCycles;
893 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
894 int b,e;
895 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
896 if (para.debuglev_>90){
897 cout << " spactra <" << *iSpec << ">" << endl;
898 cout << " cycle " << iSpec->substr(b) << endl;
899 }
900 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
901 }//end loop spectra
902 cycleModeCollect[*iMode] = listOfCycles;
903 }//end of mode
904
905 //Take the Intersection of the list Of cycles in mode Off and On
906 list<int> commonCycles;
907 set_intersection(cycleModeCollect["On"].begin(),
908 cycleModeCollect["On"].end(),
909 cycleModeCollect["Off"].begin(),
910 cycleModeCollect["Off"].end(),
911 back_inserter(commonCycles)
912 );
913
914 if (para.debuglev_>90){
915 cout << "Liste of cycles common to On & Off: <";
916 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
917 cout << *i << " ";
918 }
919 cout << ">" << endl;
920 }
921
922 //Loop on cyles
923 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
924
925 string ppftag;
926 //load ON phase
927 stringstream cycle;
928 cycle << *ic;
929 ppftag = "specRawOn"+cycle.str();
930 TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
931 fin.GetObject(aSpecOn,ppftag);
932
933 //Load OFF phase
934 ppftag = "specRawOff"+cycle.str();
935 TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
936 fin.GetObject(aSpecOff,ppftag);
937
938 //Perform the difference ON-OFF
939 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
940 meanDiffONOFF_noCalib += diffOnOff_noCalib;
941
942 //JEC 29/10/11 add ON-OFF/OFF
943 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
944 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
945 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
946 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
947
948 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
949 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
950
951
952 //JEC 15/11/11 add ON/OFF and OFF/OFF
953 TMatrix<r_4> onOvOff(aSpecOn,false);
954 onOvOff.Div(aSpecOffFitltered);
955 meanONovOFF_noCalib += onOvOff;
956
957 TMatrix<r_4> offOvOff(aSpecOff,false);
958 offOvOff.Div(aSpecOffFitltered);
959 meanOFFovOFF_noCalib += offOvOff;
960
961 totalNumberCycles++;
962
963 //Quit if enough
964 if (totalNumberCycles >= para.maxNumberCycles_) break;
965 }//end of cycles
966
967 if (totalNumberCycles >= para.maxNumberCycles_) break;
968
969 }//end files
970 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
971 //Normalisation
972 if(totalNumberCycles > 0){
973 meanDiffONOFFovOFF_noCalib /= (r_4)totalNumberCycles;
974 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles;
975 meanONovOFF_noCalib /= (r_4)totalNumberCycles;
976 meanOFFovOFF_noCalib /= (r_4)totalNumberCycles;
977 }
978
979 {//save results
980 stringstream tmp;
981 tmp << totalNumberCycles;
982 string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
983
984 POutPersist fos(fileName);
985 cout << "Save output in " << fileName << endl;
986
987 string tag = "meanNoCalib";
988 fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
989
990 tag = "meanOvOffNoCalib";
991 fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
992
993 tag = "meanOnovOffNoCalib";
994 fos << PPFNameTag(tag) << meanONovOFF_noCalib;
995
996 tag = "meanOffovOffNoCalib";
997 fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;
998
999 }//end save
1000}
1001//JEC 14/11/11 New meanRawDiffOnOffCycles END
1002//-------------------------------------------------------
1003//JEC 14/11/11 Obsolete START
1004//-------------------------------------------------------
1005//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1006//Used like:
1007//
1008// void meanRawDiffOnOffCycles() throw(string) {
1009// list<string> listOfFiles;
1010// string directoryName;
1011// directoryName = para.inPath_ + "/" + para.sourceName_;
1012
1013// //Make the listing of the directory
1014// listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1015
1016// list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1017// iFileEnd = listOfFiles.end();
1018
1019// StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1020// TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1021// uint_4 nSpectra=0;
1022// //Loop on files
1023// for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1024// if (para.debuglev_>90){
1025// cout << "load file <" << *iFile << ">" << endl;
1026// }
1027// PInPersist fin(*iFile);
1028// vector<string> vec = fin.GetNameTags();
1029// list<string> listOfSpectra;
1030// //Keep only required PPF objects
1031// std::remove_copy_if(
1032// vec.begin(), vec.end(), back_inserter(listOfSpectra),
1033// not1(match)
1034// );
1035
1036// listOfSpectra.sort(stringCompare);
1037// iSpecEnd = listOfSpectra.end();
1038// //Loop of spectra matrix
1039// for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
1040// if (para.debuglev_>90){
1041// cout << " spactra <" << *iSpec << ">" << endl;
1042// }
1043// TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1044// fin.GetObject(aSpec,*iSpec);
1045// //How to see if the GetObject is ok?? Ask Reza
1046// nSpectra++;
1047// meanOfSpectra+=aSpec;
1048// }//eo loop on spectra in a file
1049// }//eo loop on files
1050
1051// //Normalisation
1052// if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
1053
1054// //Compute the reduced version of the mean and sigma
1055// TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1056// TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1057// reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
1058
1059// {//Save the result
1060// stringstream tmp;
1061// tmp << nSpectra;
1062// string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1063// cout << "Save mean based on " << nSpectra << " cycles " << endl;
1064// POutPersist fos(fileName);
1065
1066// string tag = "mean";
1067// fos << PPFNameTag(tag) << meanOfSpectra;
1068// tag = "meanred";
1069// fos << PPFNameTag(tag) << meanRedMtx;
1070// tag = "sigmared";
1071// fos << PPFNameTag(tag) << sigmaRedMtx;
1072// }
1073// }
1074//JEC 14/11/11 Obsolete END
1075//-------------------------------------------------------
1076//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1077//Used like:
1078//
1079void medianRawDiffOnOffCycles() throw(string) {
1080 list<string> listOfFiles;
1081 string directoryName;
1082 directoryName = para.inPath_ + "/" + para.sourceName_;
1083
1084 //Make the listing of the directory
1085 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1086
1087 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1088 iFileEnd = listOfFiles.end();
1089
1090
1091 TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
1092
1093 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1094 uint_4 nSpectra=0;
1095 //Loop on files
1096 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1097 if (para.debuglev_>90){
1098 cout << "load file <" << *iFile << ">" << endl;
1099 }
1100 PInPersist fin(*iFile);
1101 vector<string> vec = fin.GetNameTags();
1102 list<string> listOfSpectra;
1103 //Keep only required PPF objects
1104 std::remove_copy_if(
1105 vec.begin(), vec.end(), back_inserter(listOfSpectra),
1106 not1(match)
1107 );
1108
1109 listOfSpectra.sort(stringCompare);
1110 iSpecEnd = listOfSpectra.end();
1111 //Loop of spectra matrix
1112 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ; ++iSpec){
1113 if (para.debuglev_>90){
1114 cout << " spactra <" << *iSpec << ">" << endl;
1115 }
1116 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1117 fin.GetObject(aSpec,*iSpec);
1118
1119 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
1120
1121 nSpectra++;
1122 }//eo loop on spectra in a file
1123 }//eo loop on files
1124
1125
1126
1127 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1128 //Compute the median for each freq. and Channel
1129 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1130 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
1131 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
1132 vector<r_4> tmp;
1133 tmp0.FillTo(tmp);
1134 medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
1135 }
1136 }
1137
1138
1139 //Compute the reduced version of the mean and sigma
1140 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1141 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1142 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
1143
1144
1145 sa_size_t f1320=freqToChan(1320.);
1146 sa_size_t f1345=freqToChan(1345.);
1147 sa_size_t f1355=freqToChan(1355.);
1148 sa_size_t f1380=freqToChan(1380.);
1149 //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
1150 if (para.debuglev_>9){
1151 cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
1152 }
1153 TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
1154 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1155 double meanMed1;
1156 double sigmaMed1;
1157 TVector<r_4> band1;
1158 band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
1159 MeanSigma(band1,meanMed1,sigmaMed1);
1160 double meanMed2;
1161 double sigmaMed2;
1162 TVector<r_4> band2;
1163 band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
1164 MeanSigma(band2,meanMed2,sigmaMed2);
1165 meanMed(iCh) = (meanMed1+meanMed2)*0.5;
1166 }
1167 meanMed.Print(cout);
1168 cout << endl;
1169
1170
1171 //Compute the sigma in the range 1320MHz-1380MHz
1172 if (para.debuglev_>9){
1173 cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
1174 }
1175 TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
1176 sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1177 sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1178
1179 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1180 double meanSigma;
1181 double sigmaSigma;
1182 TVector<r_4> band;
1183 band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
1184 MeanSigma(band,meanSigma,sigmaSigma);
1185 meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
1186 sigmaMed(iCh) = meanSigma;
1187 }
1188 sigmaMed.Print(cout);
1189 cout << endl;
1190
1191
1192
1193 if (para.debuglev_>9){
1194 cout << "Compute medianOfSpectraNorm" << endl;
1195 }
1196 TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
1197 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1198 medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
1199 medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
1200 }
1201
1202
1203
1204 {//Save the result
1205 stringstream tmp;
1206 tmp << nSpectra;
1207 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1208 cout << "Save median based on " << nSpectra << " cycles " << endl;
1209 POutPersist fos(fileName);
1210
1211 string tag = "median";
1212 fos << PPFNameTag(tag) << medianOfSpectra;
1213
1214 tag = "medianNorm";
1215 fos << PPFNameTag(tag) << medianOfSpectraNorm;
1216
1217
1218 tag = "meanmedred";
1219 fos << PPFNameTag(tag) << meanRedMtx;
1220 tag = "sigmamedred";
1221 fos << PPFNameTag(tag) << sigmaRedMtx;
1222 tag = "cycleVsfreq";
1223
1224 TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1225 fos << PPFNameTag(tag) << tarr;
1226 }
1227}
1228
1229//-------------------------------------------------------
1230int main(int narg, char* arg[]) {
1231
1232 int rc = 0; //return code
1233 string msg; //message used in Exceptions
1234
1235
1236
1237 //default value for initial parameters (see Para structure on top of the file)
1238 string debuglev = "0";
1239 string action = "meanDiffOnOff";
1240 string inputPath = ".";
1241 string outputPath = ".";
1242 string sourceName = "Abell85";
1243 string ppfFile;
1244 string nSliceInFreq = "32";
1245 string typeOfCalib="perRun";
1246 string calibFreq = "1346";
1247 string calibBandFreq="6.25";
1248 string mxcycles;
1249
1250 // bool okarg=false;
1251 int ka=1;
1252 while (ka<narg) {
1253 if (strcmp(arg[ka],"-h")==0) {
1254 cout << "Usage: -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1255 << " -mxcycles <number> (max. number of cycles to be treated)\n"
1256 << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1257 << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1258 << " -src [Abell85]\n"
1259 << " -inPath [.]|<top_root_dir of the ppf file>\n"
1260 << " (ex. /sps/baoradio/AmasNancay/JEC/\n "
1261 << " -outPath [.]|<dir of the output> \n"
1262 << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1263 << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1264 << " -debug <level> "
1265 << endl;
1266 return 0;
1267 }
1268 else if (strcmp(arg[ka],"-debug")==0) {
1269 debuglev=arg[ka+1];
1270 ka+=2;
1271 }
1272 else if (strcmp(arg[ka],"-act")==0) {
1273 action=arg[ka+1];
1274 ka+=2;
1275 }
1276 else if (strcmp(arg[ka],"-calibfreq")==0) {
1277 calibFreq=arg[ka+1];
1278 ka+=2;
1279 }
1280 else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1281 calibBandFreq=arg[ka+1];
1282 ka+=2;
1283 }
1284 else if (strcmp(arg[ka],"-mxcycles")==0) {
1285 mxcycles=arg[ka+1];
1286 ka+=2;
1287 }
1288 else if (strcmp(arg[ka],"-inPath")==0) {
1289 inputPath=arg[ka+1];
1290 ka+=2;
1291 }
1292 else if (strcmp(arg[ka],"-outPath")==0) {
1293 outputPath=arg[ka+1];
1294 ka+=2;
1295 }
1296 else if (strcmp(arg[ka],"-src")==0) {
1297 sourceName=arg[ka+1];
1298 ka+=2;
1299 }
1300 else if (strcmp(arg[ka],"-ppfFile")==0) {
1301 ppfFile=arg[ka+1];
1302 ka+=2;
1303 }
1304 else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1305 nSliceInFreq=arg[ka+1];
1306 ka+=2;
1307 }
1308 else ka++;
1309 }//eo while
1310
1311 para.debuglev_ = atoi(debuglev.c_str());
1312 para.inPath_ = inputPath;
1313 para.outPath_ = outputPath;
1314 para.sourceName_ = sourceName;
1315 para.ppfFile_ = ppfFile;
1316 para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1317 para.calibFreq_ = calibFreq;
1318 para.calibBandFreq_ = calibBandFreq;
1319 para.rcalibFreq_ = atof(calibFreq.c_str());
1320 para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1321 if (mxcycles != "") {
1322 para.maxNumberCycles_ = atoi(mxcycles.c_str());
1323 } else {
1324 para.maxNumberCycles_ = std::numeric_limits<int>::max();
1325 }
1326
1327 cout << "Dump Initial parameters ............" << endl;
1328 cout << " action = " << action << "\n"
1329 << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1330 << " inputPath = " << para.inPath_ << "\n"
1331 << " outputPath = " << para.outPath_ << "\n"
1332 << " sourceName = " << para.sourceName_ << "\n"
1333 << " ppfFile = " << para.ppfFile_ << "\n"
1334 << " nSliceInFreq = " << para.nSliceInFreq_ << "\n"
1335 << " calibFreq = " << para.calibFreq_ << "\n"
1336 << " calibBandFreq = " << para.calibBandFreq_ << "\n"
1337 << " debuglev = " << para.debuglev_ << "\n";
1338 cout << "...................................." << endl;
1339
1340 if ( "" == ppfFile ) {
1341 cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1342 << endl;
1343 return 999;
1344 }
1345
1346
1347 try {
1348
1349// int b,e;
1350// char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1351// printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1352
1353 if ( action == "meanRawDiffOnOff" ) {
1354 meanRawDiffOnOffCycles();
1355 } else if (action == "medianRawDiffOnOff") {
1356 medianRawDiffOnOffCycles();
1357 } else if (action == "meanCalibBAODiffOnOff") {
1358 meanCalibBAODiffOnOffCycles();
1359 } else {
1360 msg = "Unknown action " + action;
1361 throw msg;
1362 }
1363
1364
1365 } catch (std::exception& sex) {
1366 cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name()
1367 << "\n msg= " << sex.what() << endl;
1368 rc = 78;
1369 }
1370 catch ( string str ) {
1371 cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1372 }
1373 catch (...) {
1374 cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl;
1375 rc = 79;
1376 }
1377
1378 return 0;
1379
1380}
Note: See TracBrowser for help on using the repository browser.