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

Last change on this file since 615 was 615, checked in by torrento, 14 years ago

Added calibcoeffNtp function
File size: 71.0 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//AST 18/11/11
276//-------------------------------------------------------
277//Make n-tuple with calibration coefficients
278//
279void calibCoeffNtp() throw(string) {
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 static const int NINFO=19;
291 char* calibTupleName[NINFO]={"date","run","cycle"
292 ,"meancoeffoff0","meancoeffoff1"
293 ,"meancoeffon0","meancoeffon1"
294 ,"coeffoff0","coeffoff1"
295 ,"coeffon0","coeffon1"
296 ,"invmeancoeffoff0","invmeancoeffoff1"
297 ,"invmeancoeffon0","invmeancoeffon1"
298 ,"invcoeffoff0","invcoeffoff1"
299 ,"invcoeffon0","invcoeffon1"
300 };
301 NTuple calibcoeffs(NINFO,calibTupleName);
302 r_4 xnt[NINFO];
303
304 int totalNumberRuns=0; //total number of runs
305 int totalNumberCycles=0; //total number of cycles for normalisation
306 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
307 if (para.debuglev_>90){
308 cout << "load file <" << *iFile << ">" << endl;
309 }
310
311 vector<string> tokens;
312 split(*iFile,"_",tokens);
313 string dateOfRun = tokens[1];
314 r_4 rdateOfRun = atoi(dateOfRun.c_str()) - 2011*1.e4;
315
316 if (para.debuglev_>90){
317 cout << "date <" << dateOfRun << ">" << endl;
318 cout << "rdate <" << rdateOfRun << ">" << endl;
319 }
320 vector<string> tokens2;
321 split(tokens[2],".",tokens2);
322 string srcLower = tokens2[0];
323
324 PInPersist fin(*iFile);
325 vector<string> vec = fin.GetNameTags();
326
327 vector<string> modeList;
328 modeList.push_back("On");
329 modeList.push_back("Off");
330 vector<string>::const_iterator iMode;
331
332 map<string, list<int> > cycleModeCollect;
333
334 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
335 list<string> listOfSpectra;
336 //Keep only required PPF objects
337 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
338 std::remove_copy_if(
339 vec.begin(), vec.end(), back_inserter(listOfSpectra),
340 not1(StringMatch(matchstr))
341 );
342
343 listOfSpectra.sort(stringCompare);
344 iSpecEnd = listOfSpectra.end();
345
346 matchstr = "[0-9]+";
347 //Loop of spectra matrix
348 list<int> listOfCycles;
349 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
350 int b,e;
351 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
352 if (para.debuglev_>90){
353 cout << " spectra <" << *iSpec << ">" << endl;
354 cout << " cycle " << iSpec->substr(b) << endl;
355 }
356 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
357 }//end loop spectra
358 cycleModeCollect[*iMode] = listOfCycles;
359 }//end of mode
360
361 //Take the Intersection of the list Of cycles in mode Off and On
362 list<int> commonCycles;
363 set_intersection(cycleModeCollect["On"].begin(),
364 cycleModeCollect["On"].end(),
365 cycleModeCollect["Off"].begin(),
366 cycleModeCollect["Off"].end(),
367 back_inserter(commonCycles)
368 );
369
370 if (para.debuglev_>90){
371 cout << "Liste of cycles common to On & Off: <";
372 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
373 cout << *i << " ";
374 }
375 cout << ">" << endl;
376 }
377
378 //
379 //Load BAO Calibration factors "per Cycle and Channels"
380 //Compute the mean per Cycle to
381 // fill factors "per Run and Channels" with the same cycle list
382 //
383 //
384 //TODO improve the code....
385
386 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
387 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
388 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
389 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
390
391 string calibFileName;
392 ifstream ifs;
393 sa_size_t nr,nc; //values read
394
395 //OFF Cycle per Channel
396 calibFileName = directoryName + "/"
397 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
398 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
399 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
400 ifs.open(calibFileName.c_str());
401 if ( ! ifs.is_open() ) {
402 throw calibFileName + " cannot be opened...";
403 }
404 calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
405 if(para.debuglev_>9){
406 cout << "(nr,nc): "<< nr << "," << nc << endl;
407 calibBAOfactors_Off_Cycle_Ch0.Print(cout);
408 cout << endl;
409 }
410
411 TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
412 calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
413 {//Compute the mean
414 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
415 double mean,sigma;
416 MeanSigma(coef,mean,sigma);
417 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
418 }
419 if(para.debuglev_>9){
420 cout << "Fill calib. with mean value " << endl;
421 calibBAOfactors_Off_Run_Ch0.Print(cout);
422 cout << endl;
423 }
424
425 TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch0(nr,nc);
426 invcalibBAOfactors_Off_Cycle_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
427 {//Compute the inverse value
428 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
429 invcalibBAOfactors_Off_Cycle_Ch0.Column(1).SetCst(1);
430 invcalibBAOfactors_Off_Cycle_Ch0.Column(1).Div(coef);
431 if(para.debuglev_>9){
432 cout << "Fill calib. with inverse value " << endl;
433 invcalibBAOfactors_Off_Cycle_Ch0.Print(cout);
434 cout << endl;
435 }
436 }
437 TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch0(nr,nc);
438 invcalibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
439 {//Compute the inverse mean
440 TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
441 double mean,sigma;
442 MeanSigma(coef,mean,sigma);
443 invcalibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
444 }
445 if(para.debuglev_>9){
446 cout << "Fill calib. with inverse mean value " << endl;
447 invcalibBAOfactors_Off_Run_Ch0.Print(cout);
448 cout << endl;
449 }
450 ifs.close();
451
452 //
453 calibFileName = directoryName + "/"
454 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
455 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
456 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
457 ifs.open(calibFileName.c_str());
458 if ( ! ifs.is_open() ) {
459
460 throw calibFileName + " cannot be opened...";
461 }
462 calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
463 if(para.debuglev_>9){
464 cout << "(nr,nc): "<< nr << "," << nc << endl;
465 calibBAOfactors_Off_Cycle_Ch1.Print(cout);
466 cout << endl;
467 }
468 TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
469 calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
470 {//Compute the mean
471 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
472 double mean,sigma;
473 MeanSigma(coef,mean,sigma);
474 // cout << "Mean: " << mean << " sigma " << sigma << endl;
475 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
476 }
477 if(para.debuglev_>9){
478 cout << "Fill calib. with mean value " << endl;
479 calibBAOfactors_Off_Run_Ch1.Print(cout);
480 cout << endl;
481 }
482 TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch1(nr,nc);
483 invcalibBAOfactors_Off_Cycle_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
484 {//Compute the inverse value
485 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
486 invcalibBAOfactors_Off_Cycle_Ch1.Column(1).SetCst(1);
487 invcalibBAOfactors_Off_Cycle_Ch1.Column(1).Div(coef);
488 if(para.debuglev_>9){
489 cout << "Fill calib. with inverse value " << endl;
490 invcalibBAOfactors_Off_Cycle_Ch1.Print(cout);
491 cout << endl;
492 }
493 }
494 TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch1(nr,nc);
495 invcalibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
496 {//Compute the inverse mean
497 TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
498 double mean,sigma;
499 MeanSigma(coef,mean,sigma);
500 invcalibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
501 }
502 if(para.debuglev_>9){
503 cout << "Fill calib. with inverse mean value " << endl;
504 invcalibBAOfactors_Off_Run_Ch1.Print(cout);
505 cout << endl;
506 }
507 ifs.close();
508
509 //ON Cycle per Channel
510 calibFileName = directoryName + "/"
511 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
512 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
513 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
514 ifs.open(calibFileName.c_str());
515 if ( ! ifs.is_open() ) {
516
517 throw calibFileName + " cannot be opened...";
518 }
519 calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
520 if(para.debuglev_>9){
521 cout << "(nr,nc): "<< nr << "," << nc << endl;
522 calibBAOfactors_On_Cycle_Ch0.Print(cout);
523 cout << endl;
524 }
525
526 TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
527 calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
528 {//Compute the mean
529 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
530 double mean,sigma;
531 MeanSigma(coef,mean,sigma);
532 // cout << "Mean: " << mean << " sigma " << sigma << endl;
533 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
534 }
535 if(para.debuglev_>9){
536 cout << "Fill calib. with mean value " << endl;
537 calibBAOfactors_On_Run_Ch0.Print(cout);
538 cout << endl;
539 }
540
541 TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch0(nr,nc);
542 invcalibBAOfactors_On_Cycle_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
543 {//Compute the inverse value
544 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
545 invcalibBAOfactors_On_Cycle_Ch0.Column(1).SetCst(1);
546 invcalibBAOfactors_On_Cycle_Ch0.Column(1).Div(coef);
547 if(para.debuglev_>9){
548 cout << "Fill calib. with inverse value " << endl;
549 invcalibBAOfactors_On_Cycle_Ch0.Print(cout);
550 cout << endl;
551 }
552 }
553 TMatrix<r_4> invcalibBAOfactors_On_Run_Ch0(nr,nc);
554 invcalibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
555 {//Compute the inverse mean
556 TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
557 double mean,sigma;
558 MeanSigma(coef,mean,sigma);
559 invcalibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
560 }
561 if(para.debuglev_>9){
562 cout << "Fill calib. with inverse mean value " << endl;
563 invcalibBAOfactors_On_Run_Ch0.Print(cout);
564 cout << endl;
565 }
566 ifs.close();
567
568 calibFileName = directoryName + "/"
569 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
570 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
571 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
572 ifs.open(calibFileName.c_str());
573 if ( ! ifs.is_open() ) {
574 throw calibFileName + " cannot be opened...";
575 }
576 calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
577 if(para.debuglev_>9){
578 cout << "(nr,nc): "<< nr << "," << nc << endl;
579 calibBAOfactors_On_Cycle_Ch1.Print(cout);
580 cout << endl;
581 }
582 TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
583 calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
584 {//Compute the mean
585 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
586 double mean,sigma;
587 MeanSigma(coef,mean,sigma);
588 // cout << "Mean: " << mean << " sigma " << sigma << endl;
589 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
590 }
591 if(para.debuglev_>9){
592 cout << "Fill calib. with mean value " << endl;
593 calibBAOfactors_On_Run_Ch1.Print(cout);
594 cout << endl;
595 }
596
597 TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch1(nr,nc);
598 invcalibBAOfactors_On_Cycle_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
599 {//Compute the inverse value
600 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
601 invcalibBAOfactors_On_Cycle_Ch1.Column(1).SetCst(1);
602 invcalibBAOfactors_On_Cycle_Ch1.Column(1).Div(coef);
603 if(para.debuglev_>9){
604 cout << "Fill calib. with inverse value " << endl;
605 invcalibBAOfactors_On_Cycle_Ch1.Print(cout);
606 cout << endl;
607 }
608 }
609 TMatrix<r_4> invcalibBAOfactors_On_Run_Ch1(nr,nc);
610 invcalibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
611 {//Compute the inverse mean
612 TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
613 double mean,sigma;
614 MeanSigma(coef,mean,sigma);
615 invcalibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
616 }
617 if(para.debuglev_>9){
618 cout << "Fill calib. with inverse mean value " << endl;
619 invcalibBAOfactors_On_Run_Ch1.Print(cout);
620 cout << endl;
621 }
622 ifs.close();
623
624 //link <cycle> - <calibration coefficient>
625 //We cannot rely on identical cycle list of the OFF and ON calibration
626 map<int,r_4> calibBAO_Off_Run_Ch0;
627 map<int,r_4> calibBAO_Off_Run_Ch1;
628 map<int,r_4> calibBAO_On_Run_Ch0;
629 map<int,r_4> calibBAO_On_Run_Ch1;
630
631 map<int,r_4> calibBAO_Off_Cycle_Ch0;
632 map<int,r_4> calibBAO_Off_Cycle_Ch1;
633 map<int,r_4> calibBAO_On_Cycle_Ch0;
634 map<int,r_4> calibBAO_On_Cycle_Ch1;
635
636 map<int,r_4> invcalibBAO_Off_Run_Ch0;
637 map<int,r_4> invcalibBAO_Off_Run_Ch1;
638 map<int,r_4> invcalibBAO_On_Run_Ch0;
639 map<int,r_4> invcalibBAO_On_Run_Ch1;
640
641 map<int,r_4> invcalibBAO_Off_Cycle_Ch0;
642 map<int,r_4> invcalibBAO_Off_Cycle_Ch1;
643 map<int,r_4> invcalibBAO_On_Cycle_Ch0;
644 map<int,r_4> invcalibBAO_On_Cycle_Ch1;
645
646 //per Run based BAO coefficients
647 nr = calibBAOfactors_Off_Run_Ch0.NRows();
648 for (sa_size_t ir=0; ir<nr; ++ir){
649// cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
650// << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
651
652 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
653 = calibBAOfactors_Off_Run_Ch0(ir,1);
654 calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
655 = calibBAOfactors_Off_Cycle_Ch0(ir,1);
656 calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
657 = calibBAOfactors_Off_Run_Ch1(ir,1);
658 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
659 = calibBAOfactors_Off_Cycle_Ch1(ir,1);
660
661 invcalibBAO_Off_Run_Ch0[(int)invcalibBAOfactors_Off_Run_Ch0(ir,0)]
662 = invcalibBAOfactors_Off_Run_Ch0(ir,1);
663 invcalibBAO_Off_Cycle_Ch0[(int)invcalibBAOfactors_Off_Cycle_Ch0(ir,0)]
664 = invcalibBAOfactors_Off_Cycle_Ch0(ir,1);
665 invcalibBAO_Off_Run_Ch1[(int)invcalibBAOfactors_Off_Run_Ch1(ir,0)]
666 = invcalibBAOfactors_Off_Run_Ch1(ir,1);
667 invcalibBAO_Off_Cycle_Ch1[(int)invcalibBAOfactors_Off_Cycle_Ch1(ir,0)]
668 = invcalibBAOfactors_Off_Cycle_Ch1(ir,1);
669 }//eo loop on coef Off
670
671 nr = calibBAOfactors_On_Run_Ch0.NRows();
672 for (sa_size_t ir=0; ir<nr; ++ir){
673 calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
674 = calibBAOfactors_On_Run_Ch0(ir,1);
675 calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
676 = calibBAOfactors_On_Cycle_Ch0(ir,1);
677 calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
678 = calibBAOfactors_On_Run_Ch1(ir,1);
679 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
680 = calibBAOfactors_On_Cycle_Ch1(ir,1);
681
682 invcalibBAO_On_Run_Ch0[(int)invcalibBAOfactors_On_Run_Ch0(ir,0)]
683 = invcalibBAOfactors_On_Run_Ch0(ir,1);
684 invcalibBAO_On_Cycle_Ch0[(int)invcalibBAOfactors_On_Cycle_Ch0(ir,0)]
685 = invcalibBAOfactors_On_Cycle_Ch0(ir,1);
686 invcalibBAO_On_Run_Ch1[(int)invcalibBAOfactors_On_Run_Ch1(ir,0)]
687 = invcalibBAOfactors_On_Run_Ch1(ir,1);
688 invcalibBAO_On_Cycle_Ch1[(int)invcalibBAOfactors_On_Cycle_Ch1(ir,0)]
689 = invcalibBAOfactors_On_Cycle_Ch1(ir,1);
690 }//eo loop on coef On
691
692 //Loop on cyles
693 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
694
695 //Look if the cycle has been calibrated...
696 bool isCycleCalibrated =
697 ( calibBAO_On_Run_Ch0.count(*ic) *
698 calibBAO_On_Run_Ch1.count(*ic) *
699 calibBAO_Off_Run_Ch0.count(*ic) *
700 calibBAO_Off_Run_Ch1.count(*ic) *
701 calibBAO_On_Cycle_Ch0.count(*ic) *
702 calibBAO_On_Cycle_Ch1.count(*ic) *
703 calibBAO_Off_Cycle_Ch0.count(*ic) *
704 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
705
706 if(para.debuglev_>9) {
707 cout << "Calibration coefficients for cycle "<<*ic << endl;
708 if (isCycleCalibrated) {
709 cout << "Cycle calibrated " << endl;
710 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
711 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
712 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
713 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
714 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
715 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
716 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
717 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
718 } else {
719 cout << "Cycle NOT calibrated " << endl;
720 }
721 }//debug
722
723 if ( ! isCycleCalibrated ) continue;
724
725 totalNumberCycles++;
726 //Fill NTuple
727 xnt[0] = rdateOfRun;
728 xnt[1] = totalNumberRuns;
729 xnt[2] = totalNumberCycles;
730 xnt[3] = calibBAO_Off_Run_Ch0[*ic];
731 xnt[4] = calibBAO_Off_Run_Ch1[*ic];
732 xnt[5] = calibBAO_On_Run_Ch0[*ic];
733 xnt[6] = calibBAO_On_Run_Ch1[*ic];
734 xnt[7] = calibBAO_Off_Cycle_Ch0[*ic];
735 xnt[8] = calibBAO_Off_Cycle_Ch1[*ic];
736 xnt[9] = calibBAO_On_Cycle_Ch0[*ic];
737 xnt[10] = calibBAO_On_Cycle_Ch1[*ic];
738 xnt[11] = invcalibBAO_Off_Run_Ch0[*ic];
739 xnt[12] = invcalibBAO_Off_Run_Ch1[*ic];
740 xnt[13] = invcalibBAO_On_Run_Ch0[*ic];
741 xnt[14] = invcalibBAO_On_Run_Ch1[*ic];
742 xnt[15] = invcalibBAO_Off_Cycle_Ch0[*ic];
743 xnt[16] = invcalibBAO_Off_Cycle_Ch1[*ic];
744 xnt[17] = invcalibBAO_On_Cycle_Ch0[*ic];
745 xnt[18] = invcalibBAO_On_Cycle_Ch1[*ic];
746
747 calibcoeffs.Fill(xnt);
748
749 //Quit if enough
750 if (totalNumberCycles >= para.maxNumberCycles_) break;
751
752 }//eo loop on cycles
753 if (totalNumberCycles >= para.maxNumberCycles_) break;
754 totalNumberRuns++;
755 }//eo loop on spectra in a file
756 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
757
758 {//save the results
759 stringstream tmp;
760 tmp << totalNumberCycles;
761 string fileName = para.outPath_+"/calibcoeffTuple_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
762
763 POutPersist fos(fileName);
764 cout << "Save output in " << fileName << endl;
765
766 string tag = "calib";
767 fos << PPFNameTag(tag) << calibcoeffs;
768
769 }//end of save
770
771}
772//-------------------------------------------------------
773//Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra
774//
775void meanCalibBAODiffOnOffCycles() throw(string) {
776
777 list<string> listOfFiles;
778 string directoryName;
779 directoryName = para.inPath_ + "/" + para.sourceName_;
780
781 //Make the listing of the directory
782 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
783
784 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
785 iFileEnd = listOfFiles.end();
786
787 //mean ON-OFF over the list of cycles
788 TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //set to 0
789 TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //set to 0
790 TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //set to 0
791 TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //set to 0
792 static const int NINFO=25;
793 char* onffTupleName[NINFO]={"cycle"
794 ,"onoffRaw0","onoffRaw1"
795 ,"onoffRun0","onoffRun1"
796 ,"onoffCycle0","onoffCycle1"
797 ,"onoffRaw01420","onoffRaw11420"
798 ,"onoffRun01420","onoffRun11420"
799 ,"onoffCycle01420","onoffCycle11420"
800 ,"onoffRaw01420side","onoffRaw11420side"
801 ,"onoffRun01420side","onoffRun11420side"
802 ,"onoffCycle01420side","onoffCycle11420side"
803 ,"onoffRaw0f14101415","onoffRaw1f14101415"
804 ,"offRaw0f14101415","offRaw1f14101415"
805 ,"onRaw0f14101415","onRaw1f14101415"
806 };
807 NTuple onoffevolution(NINFO,onffTupleName);
808 r_4 xnt[NINFO];
809
810 //Lower and Higher freq. arround the Calibration Freq. bin to perform mean follow up
811 sa_size_t chCalibLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));
812 sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
813 //Lower and Higher freq. just arround 1420.4MHz Freq. bin to perform mean follow up
814 sa_size_t ch1420Low;
815 sa_size_t ch1420High;
816 if (para.sourceName_ == "Abell85") {
817 ch1420Low = freqToChan(1420.4-0.2); // Abell85
818 ch1420High = freqToChan(1420.4+0.2);
819 } else if (para.sourceName_ == "Abell1205") {
820 ch1420Low = freqToChan(1420.4-0.3); // Abell1205
821 ch1420High = freqToChan(1420.4+0.2);
822 } else if (para.sourceName_ == "Abell2440") {
823 ch1420Low = freqToChan(1420.4); // Abell2440
824 ch1420High = freqToChan(1420.4+0.3);
825 } else {
826 ch1420Low = freqToChan(1420.4-0.2); // Abell85
827 ch1420High = freqToChan(1420.4+0.2);
828 }
829 //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
830 sa_size_t ch1420aLow = freqToChan(1418);
831 sa_size_t ch1420aHigh = freqToChan(1419);
832 sa_size_t ch1420bLow = freqToChan(1422);
833 sa_size_t ch1420bHigh = freqToChan(1423);
834
835 //1400-1420Mhz
836 sa_size_t ch1400 = freqToChan(1400);
837 // sa_size_t ch1405 = freqToChan(1400);
838 sa_size_t ch1410 = freqToChan(1410);
839 sa_size_t ch1415 = freqToChan(1415);
840 sa_size_t ch1420 = freqToChan(1420);
841
842 if (para.debuglev_>0){
843 cout << "freq. band for follow up [" << chCalibLow << ", "<< chCalibHigh << "]" << endl;
844 cout << "freq. band for follow up [" << ch1420Low << ", "<< ch1420High << "]" << endl;
845 cout << "freq. band for follow up [" << ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
846 cout << "freq. band for follow up [" << ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
847 }
848
849 //Loop on files/run
850
851 int totalNumberCycles=0; //total number of cycles for normalisation
852 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
853 if (para.debuglev_>90){
854 cout << "load file <" << *iFile << ">" << endl;
855 }
856
857 vector<string> tokens;
858 split(*iFile,"_",tokens);
859 string dateOfRun = tokens[1];
860 if (para.debuglev_>90){
861 cout << "date <" << dateOfRun << ">" << endl;
862 }
863 vector<string> tokens2;
864 split(tokens[2],".",tokens2);
865 string srcLower = tokens2[0];
866
867
868
869 PInPersist fin(*iFile);
870 vector<string> vec = fin.GetNameTags();
871
872 vector<string> modeList;
873 modeList.push_back("On");
874 modeList.push_back("Off");
875 vector<string>::const_iterator iMode;
876
877 map<string, list<int> > cycleModeCollect;
878
879 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
880 list<string> listOfSpectra;
881 //Keep only required PPF objects
882 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
883 std::remove_copy_if(
884 vec.begin(), vec.end(), back_inserter(listOfSpectra),
885 not1(StringMatch(matchstr))
886 );
887
888 listOfSpectra.sort(stringCompare);
889 iSpecEnd = listOfSpectra.end();
890
891 matchstr = "[0-9]+";
892 //Loop of spectra matrix
893 list<int> listOfCycles;
894 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
895 int b,e;
896 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
897 if (para.debuglev_>90){
898 cout << " spactra <" << *iSpec << ">" << endl;
899 cout << " cycle " << iSpec->substr(b) << endl;
900 }
901 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
902 }//end loop spectra
903 cycleModeCollect[*iMode] = listOfCycles;
904 }//end of mode
905
906 //Take the Intersection of the list Of cycles in mode Off and On
907 list<int> commonCycles;
908 set_intersection(cycleModeCollect["On"].begin(),
909 cycleModeCollect["On"].end(),
910 cycleModeCollect["Off"].begin(),
911 cycleModeCollect["Off"].end(),
912 back_inserter(commonCycles)
913 );
914
915 if (para.debuglev_>90){
916 cout << "Liste of cycles common to On & Off: <";
917 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
918 cout << *i << " ";
919 }
920 cout << ">" << endl;
921 }
922
923 //
924 //Load BAO Calibration factors "per Cycle and Channels"
925 //Compute the mean per Cycle to
926 // fill factors "per Run and Channels" with the same cycle list
927 //
928 //
929 //TODO improve the code....
930
931 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
932 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
933 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
934 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
935
936 string calibFileName;
937 ifstream ifs;
938 sa_size_t nr,nc; //values read
939
940 //OFF Cycle per Channel
941 calibFileName = directoryName + "/"
942 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
943 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
944 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
945 ifs.open(calibFileName.c_str());
946 if ( ! ifs.is_open() ) {
947
948 throw calibFileName + " cannot be opened...";
949 }
950 calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
951 if(para.debuglev_>9){
952 cout << "(nr,nc): "<< nr << "," << nc << endl;
953 calibBAOfactors_Off_Cycle_Ch0.Print(cout);
954 cout << endl;
955 }
956
957 TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
958 calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
959 {//Compute the mean
960 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
961 double mean,sigma;
962 MeanSigma(coef,mean,sigma);
963 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
964 }
965 if(para.debuglev_>9){
966 cout << "Fill calib. with mean value " << endl;
967 calibBAOfactors_Off_Run_Ch0.Print(cout);
968 cout << endl;
969 }
970 ifs.close();
971
972 //
973 calibFileName = directoryName + "/"
974 + "calib_" + dateOfRun + "_" + srcLower + "_Off_"
975 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
976 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
977 ifs.open(calibFileName.c_str());
978 if ( ! ifs.is_open() ) {
979
980 throw calibFileName + " cannot be opened...";
981 }
982 calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
983 if(para.debuglev_>9){
984 cout << "(nr,nc): "<< nr << "," << nc << endl;
985 calibBAOfactors_Off_Cycle_Ch1.Print(cout);
986 cout << endl;
987 }
988 TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
989 calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
990 {//Compute the mean
991 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
992 double mean,sigma;
993 MeanSigma(coef,mean,sigma);
994 // cout << "Mean: " << mean << " sigma " << sigma << endl;
995 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
996 }
997 if(para.debuglev_>9){
998 cout << "Fill calib. with mean value " << endl;
999 calibBAOfactors_Off_Run_Ch1.Print(cout);
1000 cout << endl;
1001 }
1002 ifs.close();
1003
1004 //ON Cycle per Channel
1005 calibFileName = directoryName + "/"
1006 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
1007 + para.calibFreq_ +"MHz-Ch0Cycles.txt";
1008 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
1009 ifs.open(calibFileName.c_str());
1010 if ( ! ifs.is_open() ) {
1011
1012 throw calibFileName + " cannot be opened...";
1013 }
1014 calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
1015 if(para.debuglev_>9){
1016 cout << "(nr,nc): "<< nr << "," << nc << endl;
1017 calibBAOfactors_On_Cycle_Ch0.Print(cout);
1018 cout << endl;
1019 }
1020
1021 TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
1022 calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
1023 {//Compute the mean
1024 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
1025 double mean,sigma;
1026 MeanSigma(coef,mean,sigma);
1027 // cout << "Mean: " << mean << " sigma " << sigma << endl;
1028 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
1029 }
1030 if(para.debuglev_>9){
1031 cout << "Fill calib. with mean value " << endl;
1032 calibBAOfactors_On_Run_Ch0.Print(cout);
1033 cout << endl;
1034 }
1035 ifs.close();
1036
1037
1038 calibFileName = directoryName + "/"
1039 + "calib_" + dateOfRun + "_" + srcLower + "_On_"
1040 + para.calibFreq_ +"MHz-Ch1Cycles.txt";
1041 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
1042 ifs.open(calibFileName.c_str());
1043 if ( ! ifs.is_open() ) {
1044 throw calibFileName + " cannot be opened...";
1045 }
1046 calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
1047 if(para.debuglev_>9){
1048 cout << "(nr,nc): "<< nr << "," << nc << endl;
1049 calibBAOfactors_On_Cycle_Ch1.Print(cout);
1050 cout << endl;
1051 }
1052 TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
1053 calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
1054 {//Compute the mean
1055 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
1056 double mean,sigma;
1057 MeanSigma(coef,mean,sigma);
1058 // cout << "Mean: " << mean << " sigma " << sigma << endl;
1059 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
1060 }
1061 if(para.debuglev_>9){
1062 cout << "Fill calib. with mean value " << endl;
1063 calibBAOfactors_On_Run_Ch1.Print(cout);
1064 cout << endl;
1065 }
1066
1067 ifs.close();
1068
1069 //link <cycle> - <calibration coefficient>
1070 //We cannot rely on identical cycle list of the OFF and ON calibration
1071 map<int,r_4> calibBAO_Off_Run_Ch0;
1072 map<int,r_4> calibBAO_Off_Run_Ch1;
1073 map<int,r_4> calibBAO_On_Run_Ch0;
1074 map<int,r_4> calibBAO_On_Run_Ch1;
1075
1076 map<int,r_4> calibBAO_Off_Cycle_Ch0;
1077 map<int,r_4> calibBAO_Off_Cycle_Ch1;
1078 map<int,r_4> calibBAO_On_Cycle_Ch0;
1079 map<int,r_4> calibBAO_On_Cycle_Ch1;
1080
1081 //per Run based BAO coefficients
1082 nr = calibBAOfactors_Off_Run_Ch0.NRows();
1083 for (sa_size_t ir=0; ir<nr; ++ir){
1084// cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
1085// << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
1086
1087 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
1088 = calibBAOfactors_Off_Run_Ch0(ir,1);
1089 calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
1090 = calibBAOfactors_Off_Cycle_Ch0(ir,1);
1091 calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
1092 = calibBAOfactors_Off_Run_Ch1(ir,1);
1093 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
1094 = calibBAOfactors_Off_Cycle_Ch1(ir,1);
1095 }//eo loop on coef Off
1096
1097 nr = calibBAOfactors_On_Run_Ch0.NRows();
1098 for (sa_size_t ir=0; ir<nr; ++ir){
1099 calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
1100 = calibBAOfactors_On_Run_Ch0(ir,1);
1101 calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
1102 = calibBAOfactors_On_Cycle_Ch0(ir,1);
1103 calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
1104 = calibBAOfactors_On_Run_Ch1(ir,1);
1105 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
1106 = calibBAOfactors_On_Cycle_Ch1(ir,1);
1107 }//eo loop on coef On
1108
1109 //Loop on cyles
1110 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1111
1112 //Look if the cycle has been calibrated...
1113 bool isCycleCalibrated =
1114 ( calibBAO_On_Run_Ch0.count(*ic) *
1115 calibBAO_On_Run_Ch1.count(*ic) *
1116 calibBAO_Off_Run_Ch0.count(*ic) *
1117 calibBAO_Off_Run_Ch1.count(*ic) *
1118 calibBAO_On_Cycle_Ch0.count(*ic) *
1119 calibBAO_On_Cycle_Ch1.count(*ic) *
1120 calibBAO_Off_Cycle_Ch0.count(*ic) *
1121 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
1122
1123 if(para.debuglev_>9) {
1124 cout << "Calibration coefficients for cycle "<<*ic << endl;
1125 if (isCycleCalibrated) {
1126 cout << "Cycle calibrated " << endl;
1127 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
1128 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
1129 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
1130 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
1131 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
1132 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
1133 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
1134 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
1135 } else {
1136 cout << "Cycle NOT calibrated " << endl;
1137 }
1138 }//debug
1139
1140
1141 if ( ! isCycleCalibrated ) continue;
1142
1143 string ppftag;
1144 //load ON phase
1145 stringstream cycle;
1146 cycle << *ic;
1147
1148 ppftag = "specRawOn"+cycle.str();
1149 TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1150 fin.GetObject(aSpecOn,ppftag);
1151
1152 TMatrix<r_4> aSpecOn_BAOCalibRun(aSpecOn,false);
1153 aSpecOn_BAOCalibRun(Range(0),Range::all()) /= calibBAO_On_Run_Ch0[*ic];
1154 aSpecOn_BAOCalibRun(Range(1),Range::all()) /= calibBAO_On_Run_Ch1[*ic];
1155
1156 TMatrix<r_4> aSpecOn_BAOCalibCycle(aSpecOn,false);
1157 aSpecOn_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_On_Cycle_Ch0[*ic];
1158 aSpecOn_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_On_Cycle_Ch1[*ic];
1159
1160 //Load OFF phase
1161 ppftag = "specRawOff"+cycle.str();
1162 TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1163 fin.GetObject(aSpecOff,ppftag);
1164
1165 TMatrix<r_4> aSpecOff_BAOCalibRun(aSpecOff,false);
1166 aSpecOff_BAOCalibRun(Range(0),Range::all()) /= calibBAO_Off_Run_Ch0[*ic];
1167 aSpecOff_BAOCalibRun(Range(1),Range::all()) /= calibBAO_Off_Run_Ch1[*ic];
1168
1169 TMatrix<r_4> aSpecOff_BAOCalibCycle(aSpecOff,false);
1170 aSpecOff_BAOCalibCycle(Range(0),Range::all()) /= calibBAO_Off_Cycle_Ch0[*ic];
1171 aSpecOff_BAOCalibCycle(Range(1),Range::all()) /= calibBAO_Off_Cycle_Ch1[*ic];
1172
1173
1174 //Perform the difference ON-OFF with the different calibration options
1175 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1176 meanDiffONOFF_noCalib += diffOnOff_noCalib;
1177
1178 //JEC 29/10/11 add ON-OFF/OFF
1179 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1180 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1181 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1182 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
1183
1184 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
1185 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1186
1187 TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
1188 meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
1189
1190 TMatrix<r_4> diffOnOff_perCycleCalib = aSpecOn_BAOCalibCycle - aSpecOff_BAOCalibCycle;
1191 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
1192
1193 //
1194 totalNumberCycles++;
1195 //Fill NTuple
1196 xnt[0] = totalNumberCycles;
1197
1198 //Follow up arround the Calibration Freq.
1199 TVector<r_4> meanInRange_CalibFreq_noCalib(NUMBER_OF_CHANNELS);
1200 meanInRange(diffOnOff_noCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_noCalib);
1201 xnt[1] = meanInRange_CalibFreq_noCalib(0);
1202 xnt[2] = meanInRange_CalibFreq_noCalib(1);
1203
1204 TVector<r_4> meanInRange_CalibFreq_perRunCalib(NUMBER_OF_CHANNELS);
1205 meanInRange(diffOnOff_perRunCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perRunCalib);
1206 xnt[3] = meanInRange_CalibFreq_perRunCalib(0);
1207 xnt[4] = meanInRange_CalibFreq_perRunCalib(1);
1208
1209 TVector<r_4> meanInRange_CalibFreq_perCycleCalib(NUMBER_OF_CHANNELS);
1210 meanInRange(diffOnOff_perCycleCalib,chCalibLow,chCalibHigh,meanInRange_CalibFreq_perCycleCalib);
1211 xnt[5] = meanInRange_CalibFreq_perCycleCalib(0);
1212 xnt[6] = meanInRange_CalibFreq_perCycleCalib(1);
1213
1214
1215 //Follow up arround the 1420.4MHz Freq.
1216 TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1217 meanInRange(diffOnOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1218 xnt[7] = meanInRange_1420Freq_noCalib(0);
1219 xnt[8] = meanInRange_1420Freq_noCalib(1);
1220
1221 TVector<r_4> meanInRange_1420Freq_perRunCalib(NUMBER_OF_CHANNELS);
1222 meanInRange(diffOnOff_perRunCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perRunCalib);
1223 xnt[9] = meanInRange_1420Freq_perRunCalib(0);
1224 xnt[10] = meanInRange_1420Freq_perRunCalib(1);
1225
1226 TVector<r_4> meanInRange_1420Freq_perCycleCalib(NUMBER_OF_CHANNELS);
1227 meanInRange(diffOnOff_perCycleCalib,ch1420Low,ch1420High,meanInRange_1420Freq_perCycleCalib);
1228 xnt[11] = meanInRange_1420Freq_perCycleCalib(0);
1229 xnt[12] = meanInRange_1420Freq_perCycleCalib(1);
1230
1231
1232 //Follow up below the 1420.4MHz Freq.
1233 TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1234 meanInRange(diffOnOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1235 TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1236 meanInRange(diffOnOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1237
1238 xnt[13] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1239 xnt[14] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1240
1241
1242 TVector<r_4> meanInRange_1420aFreq_perRun(NUMBER_OF_CHANNELS);
1243 meanInRange(diffOnOff_perRunCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perRun);
1244 TVector<r_4> meanInRange_1420bFreq_perRun(NUMBER_OF_CHANNELS);
1245 meanInRange(diffOnOff_perRunCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perRun);
1246
1247 xnt[15] = (meanInRange_1420aFreq_perRun(0) + meanInRange_1420bFreq_perRun(0))/2.;
1248 xnt[16] = (meanInRange_1420aFreq_perRun(1) + meanInRange_1420bFreq_perRun(1))/2.;
1249
1250
1251 TVector<r_4> meanInRange_1420aFreq_perCycle(NUMBER_OF_CHANNELS);
1252 meanInRange(diffOnOff_perCycleCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_perCycle);
1253 TVector<r_4> meanInRange_1420bFreq_perCycle(NUMBER_OF_CHANNELS);
1254 meanInRange(diffOnOff_perCycleCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_perCycle);
1255
1256 xnt[17] = (meanInRange_1420aFreq_perCycle(0) + meanInRange_1420bFreq_perCycle(0))/2.;
1257 xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
1258
1259
1260 //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1261 TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
1262 meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
1263 xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
1264 xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
1265
1266 //JEC 18/11/11 follow up the 1400-1420MHz OFF only
1267 TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1268 aSpecOffovOff.Div(aSpecOffFitltered);
1269
1270 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1271 meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1272
1273 xnt[21] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1274 xnt[22] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1275
1276 TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1277 aSpecOnovOff.Div(aSpecOffFitltered);
1278
1279 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1280 meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1281
1282 xnt[23] = meanInRange_1410a1415Freq_ON_noCalib(0);
1283 xnt[24] = meanInRange_1410a1415Freq_ON_noCalib(1);
1284
1285 //store infos to Ntuple
1286 onoffevolution.Fill(xnt);
1287
1288 //Quit if enough
1289 if (totalNumberCycles >= para.maxNumberCycles_) break;
1290
1291 }//eo loop on cycles
1292 if (totalNumberCycles >= para.maxNumberCycles_) break;
1293
1294 }//eo loop on spectra in a file
1295 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1296 //Normalisation
1297 if(totalNumberCycles > 0){
1298 //JEC 29/10 add ON-OFF/OFF
1299 meanDiffONOFFovOFF_noCalib /= (r_4)totalNumberCycles;
1300 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles;
1301 meanDiffONOFF_perRunCalib /= (r_4)totalNumberCycles;
1302 meanDiffONOFF_perCycleCalib /= (r_4)totalNumberCycles;
1303 }
1304
1305 //Compute the reduced version of the mean and sigma
1306 TMatrix<r_4> meanRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1307 TMatrix<r_4> sigmaRedMtx_noCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1308 reduceSpectra(meanDiffONOFF_noCalib,meanRedMtx_noCalib,sigmaRedMtx_noCalib);
1309
1310 TMatrix<r_4> meanRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1311 TMatrix<r_4> sigmaRedMtx_perRunCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1312 reduceSpectra(meanDiffONOFF_perRunCalib,meanRedMtx_perRunCalib,sigmaRedMtx_perRunCalib);
1313
1314 TMatrix<r_4> meanRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1315 TMatrix<r_4> sigmaRedMtx_perCycleCalib(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1316 reduceSpectra(meanDiffONOFF_perCycleCalib,meanRedMtx_perCycleCalib,sigmaRedMtx_perCycleCalib);
1317
1318 {//save the results
1319 stringstream tmp;
1320 tmp << totalNumberCycles;
1321 string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1322
1323 POutPersist fos(fileName);
1324 cout << "Save output in " << fileName << endl;
1325
1326 string tag = "meanNoCalib";
1327 fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1328
1329 //JEC 29/10/11
1330 tag = "meanOvOffNoCalib";
1331 fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1332
1333 tag = "meanPerRunCalib";
1334 fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
1335 tag = "meanPerCycleCalib";
1336 fos << PPFNameTag(tag) << meanDiffONOFF_perCycleCalib;
1337
1338 tag = "redmeanNoCalib";
1339 fos << PPFNameTag(tag) << meanRedMtx_noCalib;
1340 tag = "redsigmaNoCalib";
1341 fos << PPFNameTag(tag) << sigmaRedMtx_noCalib;
1342
1343 tag = "redmeanPerRunCalib";
1344 fos << PPFNameTag(tag) << meanRedMtx_perRunCalib;
1345 tag = "redsigmaPerRunCalib";
1346 fos << PPFNameTag(tag) << sigmaRedMtx_perRunCalib;
1347
1348 tag = "redmeanPerCycleCalib";
1349 fos << PPFNameTag(tag) << meanRedMtx_perCycleCalib;
1350 tag = "redsigmaPerCycleCalib";
1351 fos << PPFNameTag(tag) << sigmaRedMtx_perCycleCalib;
1352
1353 tag = "onoffevol";
1354 fos << PPFNameTag(tag) << onoffevolution;
1355 }//end of save
1356}
1357//JEC 14/11/11 New meanRawDiffOnOffCycles START
1358//-------------------------------------------------------
1359//Compute the mean of Diff ON-OFF/OFF from Raw spectra
1360//Used like:
1361//
1362void meanRawDiffOnOffCycles() throw(string) {
1363 list<string> listOfFiles;
1364 string directoryName;
1365 directoryName = para.inPath_ + "/" + para.sourceName_;
1366
1367 //Make the listing of the directory
1368 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1369
1370 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1371 iFileEnd = listOfFiles.end();
1372
1373 //mean ON-OFF over the list of cycles
1374 TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/GAIN
1375 TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
1376 TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // ON/Filtered_OFF
1377 TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
1378
1379 //Tuple only for RAW things to follow
1380 static const int NINFO=11;
1381 char* onffTupleName[NINFO]={"cycle"
1382 ,"onoffRaw01420","onoffRaw11420"
1383 ,"onoffRaw01420side","onoffRaw11420side"
1384 ,"onoffRaw0f14101415","onoffRaw1f14101415"
1385 ,"offRaw0f14101415","offRaw1f14101415"
1386 ,"onRaw0f14101415","onRaw1f14101415"
1387 };
1388 NTuple onoffevolution(NINFO,onffTupleName);
1389 r_4 xnt[NINFO];
1390
1391 //Lower and Higher freq. just arround 1420.4MHz Freq. bin to perform mean follow up
1392 sa_size_t ch1420Low = freqToChan(1420.4-0.2);
1393 sa_size_t ch1420High = freqToChan(1420.4+0.2);
1394
1395 //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
1396 sa_size_t ch1420aLow = freqToChan(1418);
1397 sa_size_t ch1420aHigh = freqToChan(1419);
1398 sa_size_t ch1420bLow = freqToChan(1422);
1399 sa_size_t ch1420bHigh = freqToChan(1423);
1400
1401 //1400-1420Mhz
1402 sa_size_t ch1400 = freqToChan(1400);
1403 // sa_size_t ch1405 = freqToChan(1400);
1404 sa_size_t ch1410 = freqToChan(1410);
1405 sa_size_t ch1415 = freqToChan(1415);
1406 sa_size_t ch1420 = freqToChan(1420);
1407
1408 if (para.debuglev_>0){
1409 cout << "freq. band for follow up [" << ch1420Low << ", "<< ch1420High << "]" << endl;
1410 cout << "freq. band for follow up [" << ch1420aLow << ", "<< ch1420aHigh << "]" << endl;
1411 cout << "freq. band for follow up [" << ch1420bLow << ", "<< ch1420bHigh << "]" << endl;
1412 }
1413
1414 int totalNumberCycles=0; //total number of cycles
1415
1416 //Loop on files
1417 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1418 if (para.debuglev_>90){
1419 }
1420 PInPersist fin(*iFile);
1421 vector<string> vec = fin.GetNameTags();
1422
1423 vector<string> modeList;
1424 modeList.push_back("On");
1425 modeList.push_back("Off");
1426 vector<string>::const_iterator iMode;
1427
1428 map<string, list<int> > cycleModeCollect;
1429
1430 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
1431 list<string> listOfSpectra;
1432 //Keep only required PPF objects
1433 string matchstr = "specRaw"+(*iMode)+"[0-9]+";
1434 std::remove_copy_if(
1435 vec.begin(), vec.end(), back_inserter(listOfSpectra),
1436 not1(StringMatch(matchstr))
1437 );
1438
1439 listOfSpectra.sort(stringCompare);
1440 iSpecEnd = listOfSpectra.end();
1441
1442 matchstr = "[0-9]+";
1443 //Loop of spectra matrix
1444 list<int> listOfCycles;
1445 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){
1446 int b,e;
1447 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
1448 if (para.debuglev_>90){
1449 cout << " spactra <" << *iSpec << ">" << endl;
1450 cout << " cycle " << iSpec->substr(b) << endl;
1451 }
1452 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
1453 }//end loop spectra
1454 cycleModeCollect[*iMode] = listOfCycles;
1455 }//end of mode
1456
1457 //Take the Intersection of the list Of cycles in mode Off and On
1458 list<int> commonCycles;
1459 set_intersection(cycleModeCollect["On"].begin(),
1460 cycleModeCollect["On"].end(),
1461 cycleModeCollect["Off"].begin(),
1462 cycleModeCollect["Off"].end(),
1463 back_inserter(commonCycles)
1464 );
1465
1466 if (para.debuglev_>90){
1467 cout << "Liste of cycles common to On & Off: <";
1468 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
1469 cout << *i << " ";
1470 }
1471 cout << ">" << endl;
1472 }
1473
1474 //Loop on cyles
1475 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
1476
1477 string ppftag;
1478 //load ON phase
1479 stringstream cycle;
1480 cycle << *ic;
1481 ppftag = "specRawOn"+cycle.str();
1482 TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1483 fin.GetObject(aSpecOn,ppftag);
1484
1485 //Load OFF phase
1486 ppftag = "specRawOff"+cycle.str();
1487 TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1488 fin.GetObject(aSpecOff,ppftag);
1489
1490 //Perform the difference ON-OFF
1491 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
1492 meanDiffONOFF_noCalib += diffOnOff_noCalib;
1493
1494 //JEC 29/10/11 add ON-OFF/OFF
1495 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
1496 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1497 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
1498 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
1499
1500 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
1501 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
1502
1503
1504 //JEC 15/11/11 add ON/OFF and OFF/OFF
1505 TMatrix<r_4> onOvOff(aSpecOn,false);
1506 onOvOff.Div(aSpecOffFitltered);
1507 meanONovOFF_noCalib += onOvOff;
1508
1509 TMatrix<r_4> offOvOff(aSpecOff,false);
1510 offOvOff.Div(aSpecOffFitltered);
1511 meanOFFovOFF_noCalib += offOvOff;
1512
1513 totalNumberCycles++;
1514
1515
1516 //Fill NTuple
1517 xnt[0] = totalNumberCycles;
1518
1519
1520 //Follow up arround the 1420.4MHz Freq.
1521 TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
1522 meanInRange(diffOnOffOvOff_noCalib,ch1420Low,ch1420High,meanInRange_1420Freq_noCalib);
1523 xnt[1] = meanInRange_1420Freq_noCalib(0);
1524 xnt[2] = meanInRange_1420Freq_noCalib(1);
1525
1526
1527 //Follow up below the 1420.4MHz Freq.
1528 TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
1529 meanInRange(diffOnOffOvOff_noCalib,ch1420aLow,ch1420aHigh,meanInRange_1420aFreq_noCalib);
1530 TVector<r_4> meanInRange_1420bFreq_noCalib(NUMBER_OF_CHANNELS);
1531 meanInRange(diffOnOffOvOff_noCalib,ch1420bLow,ch1420bHigh,meanInRange_1420bFreq_noCalib);
1532
1533 xnt[3] = (meanInRange_1420aFreq_noCalib(0) + meanInRange_1420bFreq_noCalib(0))/2.;
1534 xnt[4] = (meanInRange_1420aFreq_noCalib(1) + meanInRange_1420bFreq_noCalib(1))/2.;
1535
1536
1537 //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
1538 TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
1539 meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
1540 xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
1541 xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
1542
1543 //JEC 18/11/11 follow up the 1400-1420MHz OFF only
1544 TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
1545 aSpecOffovOff.Div(aSpecOffFitltered);
1546
1547 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
1548 meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
1549
1550 xnt[7] = meanInRange_1410a1415Freq_OFF_noCalib(0);
1551 xnt[8] = meanInRange_1410a1415Freq_OFF_noCalib(1);
1552
1553 TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
1554 aSpecOnovOff.Div(aSpecOffFitltered);
1555
1556 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
1557 meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
1558
1559 xnt[9] = meanInRange_1410a1415Freq_ON_noCalib(0);
1560 xnt[10] = meanInRange_1410a1415Freq_ON_noCalib(1);
1561
1562 //store infos to Ntuple
1563 onoffevolution.Fill(xnt);
1564
1565 //Quit if enough
1566 if (totalNumberCycles >= para.maxNumberCycles_) break;
1567 }//end of cycles
1568
1569 if (totalNumberCycles >= para.maxNumberCycles_) break;
1570
1571 }//end files
1572 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
1573 //Normalisation
1574 if(totalNumberCycles > 0){
1575 meanDiffONOFFovOFF_noCalib /= (r_4)totalNumberCycles;
1576 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles;
1577 meanONovOFF_noCalib /= (r_4)totalNumberCycles;
1578 meanOFFovOFF_noCalib /= (r_4)totalNumberCycles;
1579 }
1580
1581 {//save results
1582 stringstream tmp;
1583 tmp << totalNumberCycles;
1584 string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1585
1586 POutPersist fos(fileName);
1587 cout << "Save output in " << fileName << endl;
1588
1589 string tag = "meanNoCalib";
1590 fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
1591
1592 tag = "meanOvOffNoCalib";
1593 fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
1594
1595 tag = "meanOnovOffNoCalib";
1596 fos << PPFNameTag(tag) << meanONovOFF_noCalib;
1597
1598 tag = "meanOffovOffNoCalib";
1599 fos << PPFNameTag(tag) << meanOFFovOFF_noCalib;
1600
1601 tag = "onoffevol";
1602 fos << PPFNameTag(tag) << onoffevolution;
1603
1604 }//end save
1605}
1606//JEC 14/11/11 New meanRawDiffOnOffCycles END
1607//-------------------------------------------------------
1608//JEC 14/11/11 Obsolete START
1609//-------------------------------------------------------
1610//Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1611//Used like:
1612//
1613// void meanRawDiffOnOffCycles() throw(string) {
1614// list<string> listOfFiles;
1615// string directoryName;
1616// directoryName = para.inPath_ + "/" + para.sourceName_;
1617
1618// //Make the listing of the directory
1619// listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1620
1621// list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1622// iFileEnd = listOfFiles.end();
1623
1624// StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1625// TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1626// uint_4 nSpectra=0;
1627// //Loop on files
1628// for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1629// if (para.debuglev_>90){
1630// cout << "load file <" << *iFile << ">" << endl;
1631// }
1632// PInPersist fin(*iFile);
1633// vector<string> vec = fin.GetNameTags();
1634// list<string> listOfSpectra;
1635// //Keep only required PPF objects
1636// std::remove_copy_if(
1637// vec.begin(), vec.end(), back_inserter(listOfSpectra),
1638// not1(match)
1639// );
1640
1641// listOfSpectra.sort(stringCompare);
1642// iSpecEnd = listOfSpectra.end();
1643// //Loop of spectra matrix
1644// for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){
1645// if (para.debuglev_>90){
1646// cout << " spactra <" << *iSpec << ">" << endl;
1647// }
1648// TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1649// fin.GetObject(aSpec,*iSpec);
1650// //How to see if the GetObject is ok?? Ask Reza
1651// nSpectra++;
1652// meanOfSpectra+=aSpec;
1653// }//eo loop on spectra in a file
1654// }//eo loop on files
1655
1656// //Normalisation
1657// if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra);
1658
1659// //Compute the reduced version of the mean and sigma
1660// TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1661// TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1662// reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx);
1663
1664// {//Save the result
1665// stringstream tmp;
1666// tmp << nSpectra;
1667// string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1668// cout << "Save mean based on " << nSpectra << " cycles " << endl;
1669// POutPersist fos(fileName);
1670
1671// string tag = "mean";
1672// fos << PPFNameTag(tag) << meanOfSpectra;
1673// tag = "meanred";
1674// fos << PPFNameTag(tag) << meanRedMtx;
1675// tag = "sigmared";
1676// fos << PPFNameTag(tag) << sigmaRedMtx;
1677// }
1678// }
1679//JEC 14/11/11 Obsolete END
1680//-------------------------------------------------------
1681//Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra
1682//Used like:
1683//
1684void medianRawDiffOnOffCycles() throw(string) {
1685 list<string> listOfFiles;
1686 string directoryName;
1687 directoryName = para.inPath_ + "/" + para.sourceName_;
1688
1689 //Make the listing of the directory
1690 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
1691
1692 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
1693 iFileEnd = listOfFiles.end();
1694
1695
1696 TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...
1697
1698 StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects
1699 uint_4 nSpectra=0;
1700 //Loop on files
1701 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
1702 if (para.debuglev_>90){
1703 cout << "load file <" << *iFile << ">" << endl;
1704 }
1705 PInPersist fin(*iFile);
1706 vector<string> vec = fin.GetNameTags();
1707 list<string> listOfSpectra;
1708 //Keep only required PPF objects
1709 std::remove_copy_if(
1710 vec.begin(), vec.end(), back_inserter(listOfSpectra),
1711 not1(match)
1712 );
1713
1714 listOfSpectra.sort(stringCompare);
1715 iSpecEnd = listOfSpectra.end();
1716 //Loop of spectra matrix
1717 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < para.maxNumberCycles_ ; ++iSpec){
1718 if (para.debuglev_>90){
1719 cout << " spactra <" << *iSpec << ">" << endl;
1720 }
1721 TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1722 fin.GetObject(aSpec,*iSpec);
1723
1724 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec;
1725
1726 nSpectra++;
1727 }//eo loop on spectra in a file
1728 }//eo loop on files
1729
1730
1731
1732 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1733 //Compute the median for each freq. and Channel
1734 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1735 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){
1736 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions());
1737 vector<r_4> tmp;
1738 tmp0.FillTo(tmp);
1739 medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end());
1740 }
1741 }
1742
1743
1744 //Compute the reduced version of the mean and sigma
1745 TMatrix<r_4> meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1746 TMatrix<r_4> sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_);
1747 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx);
1748
1749
1750 sa_size_t f1320=freqToChan(1320.);
1751 sa_size_t f1345=freqToChan(1345.);
1752 sa_size_t f1355=freqToChan(1355.);
1753 sa_size_t f1380=freqToChan(1380.);
1754 //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]
1755 if (para.debuglev_>9){
1756 cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl;
1757 }
1758 TVector<r_4>meanMed(NUMBER_OF_CHANNELS);
1759 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1760 double meanMed1;
1761 double sigmaMed1;
1762 TVector<r_4> band1;
1763 band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions();
1764 MeanSigma(band1,meanMed1,sigmaMed1);
1765 double meanMed2;
1766 double sigmaMed2;
1767 TVector<r_4> band2;
1768 band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions();
1769 MeanSigma(band2,meanMed2,sigmaMed2);
1770 meanMed(iCh) = (meanMed1+meanMed2)*0.5;
1771 }
1772 meanMed.Print(cout);
1773 cout << endl;
1774
1775
1776 //Compute the sigma in the range 1320MHz-1380MHz
1777 if (para.debuglev_>9){
1778 cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl;
1779 }
1780 TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS);
1781 sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1782 sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_);
1783
1784 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1785 double meanSigma;
1786 double sigmaSigma;
1787 TVector<r_4> band;
1788 band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions();
1789 MeanSigma(band,meanSigma,sigmaSigma);
1790 meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra
1791 sigmaMed(iCh) = meanSigma;
1792 }
1793 sigmaMed.Print(cout);
1794 cout << endl;
1795
1796
1797
1798 if (para.debuglev_>9){
1799 cout << "Compute medianOfSpectraNorm" << endl;
1800 }
1801 TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data...
1802 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
1803 medianOfSpectraNorm.Row(iCh) -= meanMed(iCh);
1804 medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh);
1805 }
1806
1807
1808
1809 {//Save the result
1810 stringstream tmp;
1811 tmp << nSpectra;
1812 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
1813 cout << "Save median based on " << nSpectra << " cycles " << endl;
1814 POutPersist fos(fileName);
1815
1816 string tag = "median";
1817 fos << PPFNameTag(tag) << medianOfSpectra;
1818
1819 tag = "medianNorm";
1820 fos << PPFNameTag(tag) << medianOfSpectraNorm;
1821
1822
1823 tag = "meanmedred";
1824 fos << PPFNameTag(tag) << meanRedMtx;
1825 tag = "sigmamedred";
1826 fos << PPFNameTag(tag) << sigmaRedMtx;
1827 tag = "cycleVsfreq";
1828
1829 TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1)));
1830 fos << PPFNameTag(tag) << tarr;
1831 }
1832}
1833
1834//-------------------------------------------------------
1835int main(int narg, char* arg[]) {
1836
1837 int rc = 0; //return code
1838 string msg; //message used in Exceptions
1839
1840
1841
1842 //default value for initial parameters (see Para structure on top of the file)
1843 string debuglev = "0";
1844 string action = "meanDiffOnOff";
1845 string inputPath = ".";
1846 string outputPath = ".";
1847 string sourceName = "Abell85";
1848 string ppfFile;
1849 string nSliceInFreq = "32";
1850 string typeOfCalib="perRun";
1851 string calibFreq = "1346";
1852 string calibBandFreq="6.25";
1853 string mxcycles;
1854
1855 // bool okarg=false;
1856 int ka=1;
1857 while (ka<narg) {
1858 if (strcmp(arg[ka],"-h")==0) {
1859 cout << "Usage: -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
1860 << " -mxcycles <number> (max. number of cycles to be treated)\n"
1861 << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
1862 << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
1863 << " -src [Abell85]\n"
1864 << " -inPath [.]|<top_root_dir of the ppf file>\n"
1865 << " (ex. /sps/baoradio/AmasNancay/JEC/\n "
1866 << " -outPath [.]|<dir of the output> \n"
1867 << " -nSliceInFreq [32]|<number of bin in freq. to cumulate>\n"
1868 << " -ppfFile <generic name of the input ppf files> (ex. diffOnOffRaw)\n"
1869 << " -debug <level> "
1870 << endl;
1871 return 0;
1872 }
1873 else if (strcmp(arg[ka],"-debug")==0) {
1874 debuglev=arg[ka+1];
1875 ka+=2;
1876 }
1877 else if (strcmp(arg[ka],"-act")==0) {
1878 action=arg[ka+1];
1879 ka+=2;
1880 }
1881 else if (strcmp(arg[ka],"-calibfreq")==0) {
1882 calibFreq=arg[ka+1];
1883 ka+=2;
1884 }
1885 else if (strcmp(arg[ka],"-calibbandfreq")==0) {
1886 calibBandFreq=arg[ka+1];
1887 ka+=2;
1888 }
1889 else if (strcmp(arg[ka],"-mxcycles")==0) {
1890 mxcycles=arg[ka+1];
1891 ka+=2;
1892 }
1893 else if (strcmp(arg[ka],"-inPath")==0) {
1894 inputPath=arg[ka+1];
1895 ka+=2;
1896 }
1897 else if (strcmp(arg[ka],"-outPath")==0) {
1898 outputPath=arg[ka+1];
1899 ka+=2;
1900 }
1901 else if (strcmp(arg[ka],"-src")==0) {
1902 sourceName=arg[ka+1];
1903 ka+=2;
1904 }
1905 else if (strcmp(arg[ka],"-ppfFile")==0) {
1906 ppfFile=arg[ka+1];
1907 ka+=2;
1908 }
1909 else if (strcmp(arg[ka],"-nSliceInFreq")==0) {
1910 nSliceInFreq=arg[ka+1];
1911 ka+=2;
1912 }
1913 else ka++;
1914 }//eo while
1915
1916 para.debuglev_ = atoi(debuglev.c_str());
1917 para.inPath_ = inputPath;
1918 para.outPath_ = outputPath;
1919 para.sourceName_ = sourceName;
1920 para.ppfFile_ = ppfFile;
1921 para.nSliceInFreq_ = atoi(nSliceInFreq.c_str());
1922 para.calibFreq_ = calibFreq;
1923 para.calibBandFreq_ = calibBandFreq;
1924 para.rcalibFreq_ = atof(calibFreq.c_str());
1925 para.rcalibBandFreq_ = atof(calibBandFreq.c_str());
1926 if (mxcycles != "") {
1927 para.maxNumberCycles_ = atoi(mxcycles.c_str());
1928 } else {
1929 para.maxNumberCycles_ = std::numeric_limits<int>::max();
1930 }
1931
1932 cout << "Dump Initial parameters ............" << endl;
1933 cout << " action = " << action << "\n"
1934 << " maxNumberCycles = " << para.maxNumberCycles_ << "\n"
1935 << " inputPath = " << para.inPath_ << "\n"
1936 << " outputPath = " << para.outPath_ << "\n"
1937 << " sourceName = " << para.sourceName_ << "\n"
1938 << " ppfFile = " << para.ppfFile_ << "\n"
1939 << " nSliceInFreq = " << para.nSliceInFreq_ << "\n"
1940 << " calibFreq = " << para.calibFreq_ << "\n"
1941 << " calibBandFreq = " << para.calibBandFreq_ << "\n"
1942 << " debuglev = " << para.debuglev_ << "\n";
1943 cout << "...................................." << endl;
1944
1945 if ( "" == ppfFile ) {
1946 cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option"
1947 << endl;
1948 return 999;
1949 }
1950
1951
1952 try {
1953
1954// int b,e;
1955// char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e);
1956// printf("->%s<-\n(b=%d e=%d)\n",match,b,e);
1957
1958 if ( action == "meanRawDiffOnOff" ) {
1959 meanRawDiffOnOffCycles();
1960 } else if (action == "medianRawDiffOnOff") {
1961 medianRawDiffOnOffCycles();
1962 } else if (action == "meanCalibBAODiffOnOff") {
1963 meanCalibBAODiffOnOffCycles();
1964 } else if (action == "calibCoeffNtp") {
1965 calibCoeffNtp();
1966 } else {
1967 msg = "Unknown action " + action;
1968 throw msg;
1969 }
1970
1971
1972 } catch (std::exception& sex) {
1973 cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name()
1974 << "\n msg= " << sex.what() << endl;
1975 rc = 78;
1976 }
1977 catch ( string str ) {
1978 cerr << "mergeRawOnOff.cc Exception raised: " << str << endl;
1979 }
1980 catch (...) {
1981 cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl;
1982 rc = 79;
1983 }
1984
1985 return 0;
1986
1987}
Note: See TracBrowser for help on using the repository browser.