source: BAORadio/AmasNancay/v5/analyse.cc@ 618

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

Comment sleep in while loop
File size: 60.3 KB
Line 
1//
2// Analyse of Amas@Nancay runs by J.E Campagne (LAL)
3// Version 0: 1/6/2011
4//-----------------------------
5
6// Utilisation de SOPHYA pour faciliter les tests ...
7#include "sopnamsp.h"
8#include "machdefs.h"
9
10#include <stdlib.h>
11#include <dirent.h>
12#include <matharr.h>
13
14// include standard c/c++
15#include <iostream>
16#include <fstream>
17#include <string>
18#include <vector>
19#include <map>
20#include <functional>
21#include <algorithm>
22#include <numeric>
23#include <list>
24#include <exception>
25
26// include sophya mesure ressource CPU/memoire ...
27#include "resusage.h"
28#include "ctimer.h"
29#include "timing.h"
30#include "timestamp.h"
31#include "strutilxx.h"
32#include "ntuple.h"
33#include "fioarr.h"
34#include "tarrinit.h"
35#include "histinit.h"
36#include "fitsioserver.h"
37#include "fiosinit.h"
38#include "ppersist.h"
39
40
41const sa_size_t NUMBER_OF_CHANNELS = 2;
42const sa_size_t NUMBER_OF_FREQ = 8192;
43const r_4 LOWER_FREQUENCY = 1250.0; //MHz
44const r_4 TOTAL_BANDWIDTH = 250.0; //MHz
45
46//decalration of non class members functions
47extern "C" {
48 int Usage(bool);
49}
50
51
52//----------------------------------------------------------
53//Utility fonctions
54class median_of_empty_list_exception:public std::exception{
55 virtual const char* what() const throw() {
56 return "Attempt to take the median of an empty list of numbers. "
57 "The median of an empty list is undefined.";
58 }
59};
60template<class RandAccessIter>
61double median(RandAccessIter begin, RandAccessIter end)
62 throw(median_of_empty_list_exception){
63 if(begin == end){ throw median_of_empty_list_exception(); }
64 std::size_t size = end - begin;
65 std::size_t middleIdx = size/2;
66 RandAccessIter target = begin + middleIdx;
67 std::nth_element(begin, target, end);
68
69 if(size % 2 != 0){ //Odd number of elements
70 return *target;
71 }else{ //Even number of elements
72 double a = *target;
73 RandAccessIter targetNeighbor= target-1;
74 std::nth_element(begin, targetNeighbor, end);
75 return (a+*targetNeighbor)/2.0;
76 }
77}
78
79//-------------
80
81// Function for deleting pointers in map.
82template<class A, class B>
83struct DeleteMapFntor
84{
85 // Overloaded () operator.
86 // This will be called by for_each() function.
87 bool operator()(pair<A,B> x) const
88 {
89 // Assuming the second item of map is to be
90 // deleted. Change as you wish.
91 delete x.second;
92 return true;
93 }
94};
95//-----
96bool compare(const pair<string,r_4>& i, const pair<string,r_4>& j) {
97 return i.second < j.second;
98}
99//-----
100sa_size_t round_sa(r_4 r) {
101 return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
102}
103//-----
104string StringToLower(string strToConvert){
105 //change each element of the string to lower case
106 for(unsigned int i=0;i<strToConvert.length();i++) {
107 strToConvert[i] = tolower(strToConvert[i]);
108 }
109 return strToConvert;//return the converted string
110}
111//-----
112//JEC 22/9/11 comparison, not case sensitive to sort File liste START
113bool stringCompare( const string &left, const string &right ){
114 if( left.size() < right.size() )
115 return true;
116 ////////////POSSIBLY A BUG return false;
117 for( string::const_iterator lit = left.begin(), rit = right.begin(); lit != left.end() && rit != right.end(); ++lit, ++rit )
118 if( tolower( *lit ) < tolower( *rit ) )
119 return true;
120 else if( tolower( *lit ) > tolower( *rit ) )
121 return false;
122 return false; ///////TO BE FIXED
123}//JEC 22/9/11 comparison, not case sensitive to sort File liste END
124//-----
125list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
126 list<string> theList;
127
128
129 DIR *dip;
130 struct dirent *dit;
131 string msg; string fileName;
132 string fullFileName;
133 size_t found;
134
135 if ((dip=opendir(dir.c_str())) == NULL ) {
136 msg = "opendir failed on directory "+dir;
137 throw msg;
138 }
139 while ( (dit = readdir(dip)) != NULL ) {
140 fileName = dit->d_name;
141 found=fileName.find(filePettern);
142 if (found != string::npos) {
143 fullFileName = dir + "/";
144 fullFileName += fileName;
145 theList.push_back(fullFileName);
146 }
147 }//eo while
148 if (closedir(dip) == -1) {
149 msg = "closedir failed on directory "+dir;
150 throw msg;
151 }
152
153 //JEC 22/9/11 START
154 theList.sort(stringCompare);
155 //JEC 22/9/11 END
156
157 return theList;
158}
159
160
161//Declaration of local classes
162//----------------------------------------------
163//Process Interface
164class IProcess {
165public:
166 IProcess() {}
167 virtual ~IProcess() {}
168 virtual int processCmd() throw(string) =0;
169};
170//------------
171//Common Process
172class ProcessBase : public IProcess {
173public:
174 ProcessBase();
175 virtual ~ProcessBase();
176 void SetInputPath(const string& inputPath) {inputPath_ = inputPath;}
177 void SetOutputPath(const string& outputPath) {outputPath_ = outputPath;}
178 void SetSourceName(const string& sourceName) {sourceName_ = sourceName;}
179 void SetDateOfRun(const string& dateOfRun) {dateOfRun_ = dateOfRun;}
180 void SetSpectraDirectory(const string& spectraDirectory) {spectraDirectory_ = spectraDirectory;}
181 void SetTypeOfFile(const string& typeOfFile) { typeOfFile_ = typeOfFile; }
182 void SetNumCycle(const string& numcycle) {numcycle_ = numcycle; }
183 void SetScaFileName(const string& scaFileName) { scaFileName_ =scaFileName; }
184
185 void SetDebugLevel(const string& debuglev) {
186 debuglev_ = atoi(debuglev.c_str());
187 }
188
189 virtual int processCmd() throw(string);
190
191protected:
192 string inputPath_;
193 string outputPath_;
194 string sourceName_;
195 string dateOfRun_;
196 string spectraDirectory_;
197 string typeOfFile_;
198
199 string numcycle_; //cycle numbers format="first,last"
200 sa_size_t ifirstCycle_;
201 sa_size_t ilastCycle_;
202
203
204 uint_4 debuglev_;
205 string scaFileName_;
206 NTuple* scaTuple_;
207 map<sa_size_t,sa_size_t> idCycleInTuple_;
208};
209ProcessBase::ProcessBase() {
210 scaTuple_ = 0;
211}
212ProcessBase::~ProcessBase() {
213 if (scaTuple_) delete scaTuple_;
214 scaTuple_ = 0;
215}
216//------------
217//Process ON/OFF data
218//------------
219class ProcessONOFFData : public ProcessBase {
220protected:
221 string freqBAOCalibration_;//string MHz
222public:
223 ProcessONOFFData(){}
224 virtual ~ProcessONOFFData(){}
225
226 void SetFreqBAOCalibration(const string& freqBAOCalibration) {
227 freqBAOCalibration_ = freqBAOCalibration;
228 }
229
230 virtual int processCmd() throw(string);
231};
232
233
234//------------
235//Process ON/OFF Raw data
236//------------
237class ProcessONOFFRawData : public ProcessBase {
238
239public:
240 ProcessONOFFRawData(){}
241 virtual ~ProcessONOFFRawData(){}
242
243 virtual int processCmd() throw(string);
244};
245
246//------------
247//Process ON/OFF data treated with the mspec (specmfib) => no filtering at all: paquets & freq.
248//------------
249class ProcessONOFFMeanData : public ProcessBase {
250
251public:
252 ProcessONOFFMeanData(){}
253 virtual ~ProcessONOFFMeanData(){}
254
255 virtual int processCmd() throw(string);
256};
257
258
259//------------
260//Process Gain
261//------------
262class ProcessGain : public ProcessBase {
263protected:
264 string mode_; //mode of data taken for gain computation On || Off
265public:
266 ProcessGain(){}
267 virtual ~ProcessGain(){}
268
269 void SetMode(const string& mode) {mode_ = mode;}
270
271 virtual int processCmd() throw(string);
272};
273//------------
274//Process Calibration
275//------------
276class ProcessCalibration : public ProcessBase {
277protected:
278 string option_; //option of calibration procedure
279 string freqBAOCalibration_;//string MHz
280 r_4 valfreqBAOCalibration_; //value MHz
281 string bandWidthBAOCalibration_;//string MHz
282 r_4 valbandWidthBAOCalibration_;//value MHz
283
284 sa_size_t lowerFreqBin_;
285 sa_size_t upperFreqBin_;
286
287public:
288 ProcessCalibration() {}
289 virtual ~ProcessCalibration(){}
290
291 void SetOption(const string& option) {option_ = option;}
292 void SetFreqBAOCalibration(const string& freqBAOCalibration) {
293 freqBAOCalibration_ = freqBAOCalibration;
294 valfreqBAOCalibration_ = atof(freqBAOCalibration_.c_str());
295 }
296 void SetBandWidthBAOCalibration(const string& bandWidthBAOCalibration) {
297 bandWidthBAOCalibration_ = bandWidthBAOCalibration;
298 valbandWidthBAOCalibration_ = atof(bandWidthBAOCalibration_.c_str());
299 }
300
301 void ComputeLowerUpperFreqBin();
302
303 virtual int processCmd() throw(string);
304};
305void ProcessCalibration::ComputeLowerUpperFreqBin() {
306 sa_size_t c0 = round_sa(NUMBER_OF_FREQ*(valfreqBAOCalibration_-LOWER_FREQUENCY)/TOTAL_BANDWIDTH);
307 sa_size_t dc = round_sa(NUMBER_OF_FREQ*valbandWidthBAOCalibration_/TOTAL_BANDWIDTH);
308 lowerFreqBin_ = c0-dc/2;
309 upperFreqBin_ = c0+dc/2;
310}
311
312
313//----------------------------------------------------
314//----------------------------------------------------
315int main(int narg, char* arg[])
316{
317
318 //Init process types
319 map<string,IProcess*> process;
320 process["meanONOFF"] = new ProcessONOFFMeanData(); //JEC 27/10/11
321 process["rawOnOff"] = new ProcessONOFFRawData();
322 process["dataOnOff"] = new ProcessONOFFData();
323 process["gain"] = new ProcessGain();
324 process["calib"] = new ProcessCalibration();
325
326 //Init Sophya related modules
327 // SophyaInit();
328 TArrayInitiator _inia; //nneded for TArray persistancy
329 FitsIOServerInit(); //needed for input file
330
331 //message used in Exceptions
332 string msg;
333
334 //Return code
335 int rc = 0;
336
337 //Arguments managements
338 if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
339 if (narg<11) return Usage(true);
340
341 string action;
342 string inputPath = ".";
343 string outputPath = ".";
344 string sourceName;
345 string scaFile;
346 string dateOfRun;
347 string spectraDirectory;
348 string freqBAOCalib = "";
349 string bandWidthBAOCalib = "";
350 string debuglev = "0";
351 string mode = "";
352 string numcycle;
353 string calibrationOpt = "";
354
355 string typeOfFile="medfiltmtx";
356
357
358 // bool okarg=false;
359 int ka=1;
360 while (ka<(narg-1)) {
361 cout << "Debug arglist: <" << arg[ka] <<">" << endl;
362 if (strcmp(arg[ka],"-debug")==0) {
363 debuglev=arg[ka+1];
364 ka+=2;
365 }
366 else if (strcmp(arg[ka],"-act")==0) {
367 action=arg[ka+1];
368 ka+=2;
369 }
370 else if (strcmp(arg[ka],"-inPath")==0) {
371 inputPath=arg[ka+1];
372 ka+=2;
373 }
374 else if (strcmp(arg[ka],"-outPath")==0) {
375 outputPath=arg[ka+1];
376 ka+=2;
377 }
378 else if (strcmp(arg[ka],"-source")==0) {
379 sourceName=arg[ka+1];
380 ka+=2;
381 }
382 else if (strcmp(arg[ka],"-sca")==0) {
383 scaFile=arg[ka+1];
384 ka+=2;
385 }
386 else if (strcmp(arg[ka],"-date")==0) {
387 dateOfRun=arg[ka+1];
388 ka+=2;
389 }
390 else if (strcmp(arg[ka],"-specdir")==0) {
391 spectraDirectory=arg[ka+1];
392 ka+=2;
393 }
394 else if (strcmp(arg[ka],"-specname")==0) {
395 typeOfFile=arg[ka+1];
396 ka+=2;
397 }
398 else if (strcmp(arg[ka],"-freqBAOCalib")==0) {
399 freqBAOCalib = arg[ka+1];
400 ka+=2;
401 }
402 else if (strcmp(arg[ka],"-bwBAOCalib")==0) {
403 bandWidthBAOCalib = arg[ka+1];
404 ka+=2;
405 }
406 else if (strcmp(arg[ka],"-mode")==0) {
407 mode =arg[ka+1];
408 ka+=2;
409 }
410 else if (strcmp(arg[ka],"-numcycle")==0) {
411 numcycle =arg[ka+1];
412 ka+=2;
413 }
414 else if (strcmp(arg[ka],"-calibopt")==0) {
415 calibrationOpt =arg[ka+1];
416 ka+=2;
417 }
418 else ka++;
419 }//eo while
420
421
422 //JEC 21/9/11 Give the input parameters START
423 cout << "Dump Iiitial parameters ............" << endl;
424 cout << " action = " << action << "\n"
425 << " inputPath = " << inputPath << "\n"
426 << " outputPath = " << outputPath << "\n"
427 << " sourceName = " << sourceName << "\n"
428 << " scaFile = " << scaFile << "\n"
429 << " dateOfRun = " << dateOfRun << "\n"
430 << " spectraDirectory = " << spectraDirectory << "\n"
431 << " spectraName = " << typeOfFile << "\n"
432 << " freqBAOCalib = " << freqBAOCalib << "\n"
433 << " bandWidthBAOCalib = " << bandWidthBAOCalib << "\n"
434 << " debuglev = " << debuglev << "\n"
435 << " mode = " << mode << "\n"
436 << " numcycle = " << numcycle << "\n"
437 << " calibrationOpt = " << calibrationOpt << endl;
438 cout << "...................................." << endl;
439 //JEC 21/9/11 Give the input parameters END
440
441
442 try {
443 //verification of action
444 if(process.find(action) == process.end()) {
445 msg = "action ";
446 msg += action + " not valid... FATAL";
447 rc = 999;
448 throw msg;
449 }
450
451
452 //
453 //Process initialisation...
454 //
455 try {
456 ProcessBase* procbase = dynamic_cast<ProcessBase*>(process[action]);
457 if (procbase == 0) {
458 msg= "action ";
459 msg += action + "Not a <process base> type...FATAL";
460 rc = 999;
461 throw msg;
462 }
463 procbase->SetInputPath(inputPath);
464 procbase->SetOutputPath(outputPath);
465
466 if ("" == sourceName) {
467 msg = "(FATAL) missingsourceName for action " + action;
468 Usage(true);
469 throw msg;
470 }
471 procbase->SetSourceName(sourceName);
472
473 if ("" == dateOfRun) {
474 msg = "(FATAL) missing dateOfRun for action " + action;
475 Usage(true);
476 throw msg;
477 }
478 procbase->SetDateOfRun(dateOfRun);
479
480
481 if ("" == spectraDirectory) {
482 msg = "(FATAL) missing spectraDirectory for action " + action;
483 Usage(true);
484 throw msg;
485 }
486 procbase->SetSpectraDirectory(spectraDirectory);
487
488 if ("" == scaFile) {
489 msg = "(FATAL) missing scaFile for action " + action;
490 Usage(true);
491 throw msg;
492 }
493 procbase->SetScaFileName(scaFile);
494
495 if ("" == numcycle) {
496 msg = "(FATAL) missing cycle number for action " + action;
497 Usage(true);
498 throw msg;
499 }
500 procbase->SetNumCycle(numcycle);
501
502
503 procbase->SetTypeOfFile(typeOfFile);
504
505 procbase->SetDebugLevel(debuglev);
506 }
507 catch(exception& e){
508 throw e.what();
509 }
510
511 //JEC 22/9/11 Make ON-OFF analysis WO any calibration START
512// try {
513// ProcessONOFFRawData* procRawdata = dynamic_cast<ProcessONOFFRawData*>(process[action]);
514// }
515// catch(exception& e){
516// throw e.what();
517// }
518 //JEC 22/9/11 Make ON-OFF analysis WO any calibration END
519
520
521 try {
522 ProcessONOFFData* procdata = dynamic_cast<ProcessONOFFData*>(process[action]);
523 if (procdata) {
524 if (freqBAOCalib == "") {
525 msg = "(FATAL) missing calibration BAO frequency for action " + action;
526 Usage(true);
527 throw msg;
528 }
529 procdata->SetFreqBAOCalibration(freqBAOCalib);
530 }
531 }
532 catch(exception& e){
533 throw e.what();
534 }
535
536
537 try {
538 ProcessGain* procgain = dynamic_cast<ProcessGain*>(process[action]);
539 if(procgain) {
540 if (mode == "") {
541 msg = "(FATAL) missing mode-type for action " + action;
542 Usage(true);
543 throw msg;
544 }
545 procgain->SetMode(mode);
546 }
547 }
548 catch(exception& e){
549 throw e.what();
550 }
551
552 try {
553 ProcessCalibration* proccalib = dynamic_cast<ProcessCalibration*>(process[action]);
554 if(proccalib) {
555 if (calibrationOpt == "") {
556 msg = "(FATAL) missing calibration option";
557 Usage(true);
558 throw msg;
559 }
560 if (freqBAOCalib == "") {
561 msg = "(FATAL) missing calibration BAO frequency for action " + action;
562 Usage(true);
563 throw msg;
564 }
565 if (bandWidthBAOCalib == "") {
566 msg = "(FATAL) missing calibration BAO frequency band width for action " + action;
567 Usage(true);
568 throw msg;
569 }
570 proccalib->SetOption(calibrationOpt);
571 proccalib->SetFreqBAOCalibration(freqBAOCalib);
572 proccalib->SetBandWidthBAOCalibration(bandWidthBAOCalib);
573 proccalib->ComputeLowerUpperFreqBin();
574 }
575 }
576 catch(exception& e){
577 throw e.what();
578 }
579
580 //
581 //execute command
582 //
583 rc = process[action]->processCmd();
584
585 }
586 catch (std::exception& sex) {
587 cerr << "\n analyse.cc std::exception :" << (string)typeid(sex).name()
588 << "\n msg= " << sex.what() << endl;
589 rc = 78;
590 }
591 catch ( string str ) {
592 cerr << "analyse.cc Exception raised: " << str << endl;
593 }
594 catch (...) {
595 cerr << " analyse.cc catched unknown (...) exception " << endl;
596 rc = 79;
597 }
598
599
600
601
602 cout << ">>>> analyse.cc ------- END ----------- RC=" << rc << endl;
603
604 //Delete processes
605 for_each(process.begin(),process.end(), DeleteMapFntor<string,IProcess*>());
606
607 return rc;
608}
609
610//---------------------------------------------------
611int Usage(bool flag) {
612 cout << "Analyse.cc usage...." << endl;
613 cout << "analyse -act <action_type>: meanONOFF,dataOnOff, rawOnOff, gain, calib\n"
614 << " -inPath <path for input files: default='.'>\n"
615 << " -outPath <path for output files: default='.'>\n"
616 << " -source <source name> \n"
617 << " -date <YYYYMMDD>\n"
618 << " -sca <file name scaXYZ.sum.trans>\n"
619 << " -specdir <generic directory name of spectra fits files>\n"
620 << " -specname <generic name of spectra fits files>\n"
621 << " -freqBAOCalib <freq in MHZ> freq. of calibration BAO\n"
622 << " valid for act=dataOnOff\n"
623 << " -bwBAOCalib <band width MHz> band width arround central freq. for calibration BAO\n"
624 << " valid for act=calib\n"
625 << " -mode <mode_type>:\n"
626 << " valid for act=gain, mode_type: On, Off\n"
627 << " -numcycle <number>,<number>:\n"
628 << " valid for all actions"
629 << " -calibopt <option>:\n"
630 << " valid for act=calib: indiv OR mean (NOT USED)"
631 << " -debug <number> [0=default]\n"
632 << " 1: normal print\n"
633 << " 2: save intermediate spectra\n"
634 << endl;
635 if (flag) {
636 cout << "use <path>/analyse -h for detailed instructions" << endl;
637 return 5;
638 }
639 return 1;
640}
641
642int ProcessBase::processCmd() throw(string) {
643 int rc =0;
644 string msg;
645 if(debuglev_>0)cout << "Process Base" << endl;
646 //------------------------
647 //Use the sca file informations
648 //------------------------
649 // string scaFullPathName = "./";
650 //TOBE FIXED scaFullPathName += sourceName_+"/"+dateOfRun_ + StringToLower(sourceName_)+"/";
651 string scaFullPathName = inputPath_ + "/"
652 + sourceName_+ "/" +dateOfRun_ + StringToLower(sourceName_)+ "/" + scaFileName_;
653 char* scaTupleColumnName[9] = {"cycle","stcalOn","spcalOn","stOn","spOn","stcalOff","spcalOff","stOff","spOff"};
654 scaTuple_ = new NTuple(9,scaTupleColumnName);
655 int n = scaTuple_->FillFromASCIIFile(scaFullPathName);
656 if(n<0){ //Error
657 msg = "(FATAL) NTuple error loading "+ scaFullPathName;
658 rc = 999;
659 throw msg;
660 }
661
662 if(debuglev_>1){
663 cout << "ProcessBase::processCmd: dump tuple in " << scaFullPathName << endl;
664 scaTuple_->Show(cout);
665 }
666
667
668 //Get the cycles (here consider consecutive cycles)
669 //The SCA file cannot be used as the DAQ can miss some cycles...
670 // r_8 firstCycle, lastCycle;
671 // scaTuple_->GetMinMax("cycle",firstCycle,lastCycle);
672 // ifirstCycle_ = (sa_size_t)firstCycle;
673 // ilastCycle_ = (sa_size_t)lastCycle;
674 //Analyse the string given by -numcycle command line
675 int ai1=0,ai2=0;
676 sscanf(numcycle_.c_str(),"%d,%d",&ai1,&ai2);
677 ifirstCycle_ = (sa_size_t)ai1;
678 ilastCycle_ = (sa_size_t)ai2;
679
680
681 //associate cycle number to index line in tuple
682 sa_size_t nLines = scaTuple_->NbLines();
683 for(sa_size_t iL=0; iL<nLines; ++iL){
684 idCycleInTuple_[(sa_size_t)scaTuple_->GetCell(iL,"cycle")]=iL;
685 }
686
687
688 return rc;
689}
690//----------------------------------------------
691//JEC 27/10/11
692//Process the files that are output of the specmfib -act = mspec (mean+sigma of signal files)
693//----------------------------------------------
694int ProcessONOFFMeanData::processCmd() throw(string) {
695 int rc = 0;
696 try {
697 rc = ProcessBase::processCmd();
698 }
699 catch (string s) {
700 throw s;
701 }
702
703 if(debuglev_>0)cout << "Process Data" << endl;
704 vector<string> modeList;
705 modeList.push_back("On");
706 modeList.push_back("Off");
707 vector<string>::const_iterator iMode;
708
709 uint_4 id;
710 string tag;
711
712 map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
713 map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
714
715 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
716 string mode = *iMode;
717 if(debuglev_>0)cout << "Process Mean-mspec Mode " << mode << endl;
718
719 //------------------------------------------
720 //Produce mean of the mean spectra per cycle
721 //------------------------------------------
722
723 string directoryName;
724 list<string> listOfSpecFiles;
725 list<string>::const_iterator iFile, iFileEnd;
726
727
728 //
729 //loop on cycles
730 //
731 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
732 directoryName = inputPath_ + "/"
733 + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
734 + mode + "/";
735 stringstream sicycle;
736 sicycle << icycle;
737 directoryName += spectraDirectory_ + sicycle.str() + "/";
738
739 //read directory
740 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
741
742
743 //compute mean of spectra created in a cycle
744 if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
745 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
746 iFileEnd = listOfSpecFiles.end();
747 r_4 nSpectres = 0;
748 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
749 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
750 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
751 aSpectrum >> spectre;
752 spectreMean += spectre;
753 nSpectres++;
754 }// end of for files
755 if(nSpectres>0) spectreMean /= nSpectres;
756
757 //save mean spectrum
758 spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
759 }//end of for cycles
760 }//end loop on mode for raw preocess
761
762 if(debuglev_>1) {//save mean spectra on file
763 cout << "Save mean spectra" << endl;
764 string fileName;
765
766 fileName = "./dataMean_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
767
768 POutPersist fos(fileName);
769 id=0;
770 iSpectreEnd = spectreCollect.end();
771 for (iSpectre = spectreCollect.begin();
772 iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
773 tag = "specMean";
774
775 tag += (iSpectre->first).first;
776 stringstream sid;
777 sid << (iSpectre->first).second;
778 tag += sid.str();
779 if(debuglev_>9) {
780 cout << "save tag<" << tag << ">" << endl;
781 }
782
783 fos << PPFNameTag(tag) << iSpectre->second;
784 }
785 }//end of save fits
786
787
788
789
790
791 cout << "OK ProcessONOFFMeanData finished" <<endl;
792 return rc;
793}
794//----------------------------------------------
795
796//JEC 22/9/11 Make ON-OFF analysis WO any calibration START
797//----------------------------------------------
798int ProcessONOFFRawData::processCmd() throw(string) {
799 int rc = 0;
800 try {
801 rc = ProcessBase::processCmd();
802 }
803 catch (string s) {
804 throw s;
805 }
806 if(debuglev_>0)cout << "Process Raw Data ON OFF" << endl;
807 vector<string> modeList;
808 modeList.push_back("On");
809 modeList.push_back("Off");
810 vector<string>::const_iterator iMode;
811
812 uint_4 id;
813 string tag;
814
815 //
816 //Process to get sucessively
817 //Raw Spectra,
818 //The processes are separated to allow intermediate save of results
819
820 map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
821 map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
822
823 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
824 string mode = *iMode;
825 if(debuglev_>0)cout << "Process RAW Mode " << mode << endl;
826
827 //------------------------------------------
828 //Produce Raw spectra per cycle
829 //------------------------------------------
830
831 string directoryName;
832 list<string> listOfSpecFiles;
833 list<string>::const_iterator iFile, iFileEnd;
834
835
836 //
837 //loop on cycles
838 //
839 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
840 directoryName = inputPath_ + "/"
841 + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
842 + mode + "/";
843 stringstream sicycle;
844 sicycle << icycle;
845 directoryName += spectraDirectory_ + sicycle.str() + "/";
846
847 //read directory
848 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
849
850
851 //compute mean of spectra created in a cycle
852 if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
853 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
854 iFileEnd = listOfSpecFiles.end();
855 r_4 nSpectres = 0;
856 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
857 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
858 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
859 aSpectrum >> spectre;
860 spectreMean += spectre;
861 nSpectres++;
862 }// end of for files
863 if(nSpectres>0) spectreMean /= nSpectres;
864
865 //save mean spectrum
866 spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
867 }//end of for cycles
868 }//end loop on mode for raw preocess
869
870 //JEC 23/9/11 DO IT
871 // if(debuglev_>1) {//save mean spectra on file
872 cout << "Save mean raw spectra" << endl;
873 string fileName;
874 fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
875
876 POutPersist fos(fileName);
877 id=0;
878 iSpectreEnd = spectreCollect.end();
879 for (iSpectre = spectreCollect.begin();
880 iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
881 tag = "specRaw";
882 tag += (iSpectre->first).first;
883 stringstream sid;
884 sid << (iSpectre->first).second;
885 tag += sid.str();
886 if(debuglev_>9) {
887 cout << "save tag<" << tag << ">" << endl;
888 }
889 fos << PPFNameTag(tag) << iSpectre->second;
890 }
891 // }//end of save fits
892
893
894 //------------------------------------------
895 // Perform ON-OFF
896 //------------------------------------------
897
898 map< sa_size_t, TMatrix<r_4> > diffCollect;
899 map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd;
900
901 TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero
902 r_4 nCycles = 0;
903 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
904 nCycles++;
905 TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory
906 TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory
907 TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff;
908 diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) ));
909 diffMeanOnOff += diffOnOff;
910 }//end loop on cycle
911 if(nCycles>0) diffMeanOnOff/=nCycles;
912
913 //extract channels and do the mean
914 TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0
915 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
916 meanOfChan += diffMeanOnOff.Row(iCh).Transpose();
917 }
918 meanOfChan /= (r_4)NUMBER_OF_CHANNELS;
919
920
921
922 {//save diff ON-OFF on Raw data
923 if(debuglev_>0)cout << "save ON-OFF RAW spectra" << endl;
924 string fileName;
925 fileName = "./diffOnOffRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
926 POutPersist fos(fileName);
927 iDiffEnd = diffCollect.end();
928 id = 0;
929
930 //JEC 22/9/11 Mean & Sigma in 32-bins size START
931 sa_size_t nSliceFreq = 32; //TODO: put as an input parameter option ?
932 sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq;
933 //JEC 22/9/11 Mean & Sigma in 32-bins size END
934
935 for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
936 tag = "specONOFFRaw";
937 stringstream sid;
938 sid << iDiff->first;
939 tag += sid.str();
940 fos << PPFNameTag(tag) << iDiff->second;
941
942 //JEC 22/9/11 Mean & Sigma in 32-bins size START
943 if (debuglev_>9) {
944 cout << "Debug slicing: slice/delta " << nSliceFreq << " " << deltaFreq << endl;
945 }
946 TMatrix<r_4> reducedMeanDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
947 TMatrix<r_4> reducedSigmaDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
948 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
949 sa_size_t freqLow= iSlice*deltaFreq;
950 sa_size_t freqHigh= freqLow + deltaFreq -1;
951 if (debuglev_>9) {
952 cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl;
953 }
954 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
955 if (debuglev_>9) {
956 cout << "Debug .......... Channel " << iCh;
957 }
958 TVector<r_4> reducedRow;
959 reducedRow = (iDiff->second).SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
960 double mean;
961 double sigma;
962 MeanSigma(reducedRow,mean,sigma);
963 if (debuglev_>9) {
964 cout << "mean/sigma " << mean << "/" << sigma << endl;
965 }
966 reducedMeanDiffOnOff(iCh,iSlice) = mean;
967 reducedSigmaDiffOnOff(iCh,iSlice) = sigma;
968 }//loop on Channel
969 }//loop on Freq. slice
970 tag = "redMeanONOFFRaw";
971 tag += sid.str();
972 fos << PPFNameTag(tag) << reducedMeanDiffOnOff;
973 tag = "redSigmaONOFFRaw";
974 tag += sid.str();
975 fos << PPFNameTag(tag) << reducedSigmaDiffOnOff;
976 //JEC 22/9/11 END
977
978 }//loop on ON-OFF spectre
979 //save the mean also
980 fos << PPFNameTag("specONOFFRawMean") << diffMeanOnOff;
981
982 //JEC 22/9/11 START
983 TMatrix<r_4> reducedMeanDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
984 TMatrix<r_4> reducedSigmaDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default
985 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){
986 sa_size_t freqLow= iSlice*deltaFreq;
987 sa_size_t freqHigh= freqLow + deltaFreq -1;
988 if (debuglev_>9) {
989 cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl;
990 }
991 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){
992 if (debuglev_>9) {
993 cout << "Debug .......... Channel " << iCh;
994 }
995 TVector<r_4> reducedRow;
996 reducedRow = diffMeanOnOff.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions();
997 double mean;
998 double sigma;
999 MeanSigma(reducedRow,mean,sigma);
1000 if (debuglev_>9) {
1001 cout << "mean/signa " << mean << "/" << sigma << endl;
1002 }
1003 reducedMeanDiffOnOffAll(iCh,iSlice) = mean;
1004 reducedSigmaDiffOnOffAll(iCh,iSlice) = sigma;
1005 }//loop on Channel
1006 }//loop on Freq. slice
1007 tag = "redMeanONOFFRawAll";
1008 fos << PPFNameTag(tag) << reducedMeanDiffOnOffAll;
1009 tag = "redSigmaONOFFRawAll";
1010 fos << PPFNameTag(tag) << reducedSigmaDiffOnOffAll;
1011 //JEC 22/9/11 END
1012
1013
1014
1015 fos << PPFNameTag("specONOFFRaw2ChanMean") << meanOfChan;
1016 }//end of save fits
1017
1018
1019 cout << "OK rawOnOff finished" <<endl;
1020 return rc;
1021} //ProcessONOFFRawData::processCmd
1022
1023//JEC 22/9/11 Make ON-OFF analysis WO any calibration END
1024//----------------------------------------------
1025int ProcessONOFFData::processCmd() throw(string) {
1026 int rc = 0;
1027 try {
1028 rc = ProcessBase::processCmd();
1029 }
1030 catch (string s) {
1031 throw s;
1032 }
1033 if(debuglev_>0)cout << "Process Data" << endl;
1034 vector<string> modeList;
1035 modeList.push_back("On");
1036 modeList.push_back("Off");
1037 vector<string>::const_iterator iMode;
1038
1039 uint_4 id;
1040 string tag;
1041
1042 //
1043 //Process to get sucessively
1044 //Raw Spectra,
1045 //BAO Calibrated Spectra
1046 //and RT Calibrated Spectra
1047 //The pocesses are separated to allow intermediate save of results
1048
1049 map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
1050 map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
1051
1052 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
1053 string mode = *iMode;
1054 if(debuglev_>0)cout << "Process RAW Mode " << mode << endl;
1055
1056 //------------------------------------------
1057 //Produce Raw spectra per cycle
1058 //------------------------------------------
1059
1060 string directoryName;
1061 list<string> listOfSpecFiles;
1062 list<string>::const_iterator iFile, iFileEnd;
1063
1064
1065 //
1066 //loop on cycles
1067 //
1068 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
1069 //TOBE FIXED directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
1070 directoryName = "./" + mode + "/";
1071 stringstream sicycle;
1072 sicycle << icycle;
1073 directoryName += spectraDirectory_ + sicycle.str() + "/";
1074
1075 //read directory
1076 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
1077
1078
1079 //compute mean of spectra created in a cycle
1080 if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
1081 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
1082 iFileEnd = listOfSpecFiles.end();
1083 r_4 nSpectres = 0;
1084 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
1085 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
1086 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1087 aSpectrum >> spectre;
1088 spectreMean += spectre;
1089 nSpectres++;
1090 }// end of for files
1091 if(nSpectres>0) spectreMean /= nSpectres;
1092
1093 //save mean spectrum
1094 spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
1095 }//end of for cycles
1096 }//end loop on mode for raw preocess
1097
1098 if(debuglev_>1) {//save mean spectra on file
1099 cout << "Save mean raw spectra" << endl;
1100 string fileName;
1101 //TOBE FIXED fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataRaw" + ".ppf";
1102 fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1103
1104 POutPersist fos(fileName);
1105 id=0;
1106 iSpectreEnd = spectreCollect.end();
1107 for (iSpectre = spectreCollect.begin();
1108 iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
1109 tag = "specRaw";
1110
1111 //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
1112// stringstream sid;
1113// sid << id;
1114// tag += sid.str();
1115 tag += (iSpectre->first).first;
1116 stringstream sid;
1117 sid << (iSpectre->first).second;
1118 tag += sid.str();
1119 if(debuglev_>9) {
1120 cout << "save tag<" << tag << ">" << endl;
1121 }
1122 //JEC 20/9/11 modif tag to take into account Mode and Cycle number END
1123
1124 fos << PPFNameTag(tag) << iSpectre->second;
1125 }
1126 }//end of save fits
1127
1128
1129
1130 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
1131 string mode = *iMode;
1132 if(debuglev_>0)cout << "Process CALIB BAO Mode " << mode << endl;
1133 //------------------------------------------
1134 // Correct Raw spectra for BAO calibration
1135 //------------------------------------------
1136 //Read BAO calibration files
1137 sa_size_t nr,nc; //values read
1138
1139 //JEC 20/9/11 use mean calibration coeff upon all cycles START
1140 string calibFileName = inputPath_+ "/"
1141 + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_)
1142 + "/calib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
1143 + mode + "_" + freqBAOCalibration_ + "MHz-All.txt";
1144
1145 if(debuglev_>0) cout << "Read file " << calibFileName << endl;
1146 ifstream ifs(calibFileName.c_str());
1147 if ( ! ifs.is_open() ) {
1148 rc = 999;
1149 throw calibFileName + " cannot be opened...";
1150 }
1151 TVector<r_4> calibBAOfactors;
1152 if(debuglev_>9) cout << "Debug 1" << endl;
1153 calibBAOfactors.ReadASCII(ifs,nr,nc);
1154 if(debuglev_>9){
1155 cout << "Debug 2: (nr,nc): "<< nr << "," << nc << endl;
1156 calibBAOfactors.Print(cout);
1157 }
1158
1159 //JEC 20/9/11 use mean calibration coeff upon all cycles END
1160
1161 //
1162 //spectra corrected by BAO calibration factor
1163 //-----make it different on Channels and Cycles (1/06/2011) OBSOLETE
1164 //use mean calibration coeff upon all cycles (20/6/11)
1165 //warning cycles are numbered from 1,...,N
1166 //
1167 if(debuglev_>0)cout << "do calibration..." << endl;
1168 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
1169 TMatrix<r_4> specmtx(spectreCollect[make_pair(mode,icycle)],true); //share the memory
1170
1171 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
1172 //JEC 20/9/11 use mean calibration coeff upon all cycles START
1173
1174 // specmtx( Range(iCh), Range::all() ) /= calibBAOfactors(iCh,icycle-1);
1175 specmtx( Range(iCh), Range::all() ) /= calibBAOfactors(iCh);
1176 //JEC 20/9/11 use mean calibration coeff upon all cycles END
1177 }
1178 }
1179 } //end loop mode for BAO calib
1180
1181 cout << "save calibrated BAO spectra" << endl;
1182 string fileName;
1183 //TO BE FIXED fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAOCalib" + ".ppf";
1184 fileName = "./dataBAOCalib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1185
1186 POutPersist fos(fileName);
1187 iSpectreEnd = spectreCollect.end();
1188 id=0;
1189 for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
1190
1191 tag = "specBAOCalib";
1192 //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
1193 // stringstream sid;
1194 // sid << id;
1195 // tag += sid.str();
1196 tag += (iSpectre->first).first;
1197 stringstream sid;
1198 sid << (iSpectre->first).second;
1199 tag += sid.str();
1200 if(debuglev_>9) {
1201 cout << "save tag<" << tag << ">" << endl;
1202 }
1203 //JEC 20/9/11 modif tag to take into account Mode and Cycle number END
1204
1205 fos << PPFNameTag(tag) << iSpectre->second;
1206 }
1207
1208 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
1209 string mode = *iMode;
1210 if(debuglev_>0)cout << "Process CALIB RT Mode " << mode << endl;
1211 //------------------------------------------
1212 // Correct BAO calib spectra for RT calibration
1213 //------------------------------------------
1214 //Very Preliminary May-June 11
1215 //coef RT @ 1346MHz Ouest - Est associees a Ch 0 et 1
1216 r_4 calibRT[NUMBER_OF_CHANNELS] = {27.724, 22.543};
1217 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
1218 TMatrix<r_4> specmtx(spectreCollect[make_pair(mode,icycle)],true); //share the memory
1219 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
1220 specmtx( Range(iCh), Range::all() ) *= calibRT[iCh];
1221 }
1222 }
1223 }//end loop on mode RT calib
1224
1225 {//save mean spectra BAO & RT calibrated on file
1226 if(debuglev_>0)cout << "save calibrated BAO & RT spectra" << endl;
1227 string fileName;
1228 //TO BE FIXED fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAORTCalib" + ".ppf";
1229 fileName = "./dataBAORTCalib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1230 POutPersist fos(fileName);
1231 iSpectreEnd = spectreCollect.end();
1232 id = 0;
1233 for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
1234 tag = "specBAORTCalib";
1235 //JEC 20/9/11 modif tag to take into account Mode and Cycle number START
1236// stringstream sid;
1237// sid << id;
1238// tag += sid.str();
1239 tag += (iSpectre->first).first;
1240 stringstream sid;
1241 sid << (iSpectre->first).second;
1242 tag += sid.str();
1243 if(debuglev_>9) {
1244 cout << "save tag<" << tag << ">" << endl;
1245 }
1246 //JEC 20/9/11 modif tag to take into account Mode and Cycle number END
1247 fos << PPFNameTag(tag) << iSpectre->second;
1248 }
1249 }//end of save fits
1250
1251 //------------------------------------------
1252 // Perform ON-OFF
1253 //------------------------------------------
1254
1255 map< sa_size_t, TMatrix<r_4> > diffCollect;
1256 map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd;
1257
1258 TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero
1259 r_4 nCycles = 0;
1260 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
1261 nCycles++;
1262 TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory
1263 TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory
1264 TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff;
1265 diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) ));
1266 diffMeanOnOff += diffOnOff;
1267 }//end loop on cycle
1268 if(nCycles>0) diffMeanOnOff/=nCycles;
1269
1270 //exctract channels and do the mean
1271 TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0
1272 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
1273 meanOfChan += diffMeanOnOff.Row(iCh).Transpose();
1274 }
1275 meanOfChan /= (r_4)NUMBER_OF_CHANNELS;
1276
1277
1278
1279 {//save diff ON-OFF BAO & RT calibrated
1280 if(debuglev_>0)cout << "save ON-OFF spectra" << endl;
1281 string fileName;
1282 //TO BE FIXED fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "diffOnOff" + ".ppf";
1283 fileName = "./diffOnOff_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1284 POutPersist fos(fileName);
1285 iDiffEnd = diffCollect.end();
1286 id = 0;
1287 for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
1288 tag = "specONOFF";
1289 stringstream sid;
1290 //JEC 20/9/11 sid << id;
1291 sid << iDiff->first;
1292 tag += sid.str();
1293 if(debuglev_>9) {
1294 cout << "save tag<" << tag << ">" << endl;
1295 }
1296 fos << PPFNameTag(tag) << iDiff->second;
1297 }
1298 //save the mean also
1299 fos << PPFNameTag("specONOFFMean") << diffMeanOnOff;
1300 fos << PPFNameTag("specONOFF2ChanMean") << meanOfChan;
1301 }//end of save fits
1302
1303 cout << "OK dataOnOff finished";
1304
1305 return rc;
1306} //ProcessONOFFData::processCmd
1307//
1308//----------------------------------------------
1309//
1310int ProcessGain::processCmd() throw(string) {
1311 int rc = 0;
1312 string msg = "";
1313
1314 try {
1315 rc = ProcessBase::processCmd();
1316 }
1317 catch (string s) {
1318 throw s;
1319 }
1320 if(debuglev_>0)cout << "Process Gain" << endl;
1321
1322 string directoryName;
1323 //TOBE FIXED directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode_ + "/";
1324 //JEC 6/09/2011 numcycle_ repalced by ifirstCycle_ according to definition of numcycle_ and the fact that for GAIN 1 cycle is involved
1325 stringstream thegaincycle;
1326 thegaincycle << ifirstCycle_;
1327 directoryName = inputPath_ + "/"
1328 + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
1329 + mode_ + "/" + spectraDirectory_ + thegaincycle.str() + "/";
1330
1331 list<string> listOfSpecFiles;
1332 list<string>::const_iterator iFile, iFileEnd;
1333 //read directory
1334
1335 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
1336
1337 //Mean power computed over the N channels to select the spectra for gain computation
1338 //The threshold is computed "on-line" as 1% of the difference (max power -min power) over the
1339 // the min power.
1340 //A possible alternative is to set an absolute value...
1341 if(debuglev_>0)cout << "compute mean poser spectra for files in " << directoryName << endl;
1342 iFileEnd = listOfSpecFiles.end();
1343 map<string, r_4> meanpowerCollect;
1344 //JEC 21/9/11 add meanpower for each Channels START
1345 map<string, r_4> meanPowerPerChanCollect;
1346 //JEC 21/9/11 add meanpower for each Channels END
1347
1348 map<string, r_4>::const_iterator iMeanPow, iMeanPowEnd;
1349
1350 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
1351 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
1352 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1353 aSpectrum >> spectre;
1354 meanpowerCollect[*iFile] = spectre.Sum()/spectre.Size();
1355 //JEC 21/9/11 add meanpower for each Channels START
1356 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
1357 TVector<r_4> specChan(NUMBER_OF_FREQ);
1358 specChan = spectre.Row(iCh).Transpose();
1359 stringstream tmp;
1360 tmp << iCh;
1361 string tag = *iFile + "_" + tmp.str();
1362 meanPowerPerChanCollect[tag] = specChan.Sum()/specChan.Size();
1363 }
1364 //JEC 21/9/11 add meanpower for each Channels END
1365 }//end of for files
1366 pair<string, r_4> minelement = *min_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
1367 pair<string, r_4> maxelement = *max_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
1368 if(debuglev_>1){
1369 cout << "meanpowerCollect has " << meanpowerCollect.size() << " spectra registered" << endl;
1370 cout << "find min mean power "<<minelement.second << " for " << minelement.first << endl;
1371 cout << "find max mean power "<<maxelement.second << " for " << maxelement.first << endl;
1372 }
1373 r_4 threshold = minelement.second + 0.01*(maxelement.second-minelement.second);
1374 if(debuglev_>1){
1375 cout << "threshold found at " << threshold <<endl;
1376 }
1377
1378 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
1379 r_4 nSpectres = 0;
1380 iMeanPowEnd = meanpowerCollect.end();
1381 for (iMeanPow = meanpowerCollect.begin(); iMeanPow != iMeanPowEnd; ++iMeanPow) {
1382 if ( iMeanPow->second <= threshold ) {
1383 //TODO avoid the reloading of the file may speed up
1384 FitsInOutFile aSpectrum(iMeanPow->first,FitsInOutFile::Fits_RO);
1385 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
1386 aSpectrum >> spectre;
1387 spectreMean += spectre;
1388 nSpectres++;
1389 }
1390 }
1391 if(debuglev_>1){
1392 cout << "Number of files selected for gain " << nSpectres <<endl;
1393 }
1394 if(nSpectres>0) {
1395 spectreMean /= nSpectres;
1396 } else {
1397 stringstream tmp;
1398 tmp << threshold;
1399 msg = "Gain: cannot find a spectra above threshold value =" + tmp.str() + " ... FATAL";
1400 throw msg;
1401 }
1402
1403 //Save gain spectra
1404 {
1405 //use ! to override the file: special features of cfitsio library
1406 string fileName;
1407 //TOBE FIXED fileName = "!./" + sourceName_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
1408 fileName = "!"+ outputPath_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
1409 if(debuglev_>1){
1410 cout << "save gains in " << fileName << endl;
1411 }
1412 FitsInOutFile fos(fileName, FitsInOutFile::Fits_Create);
1413 fos << spectreMean;
1414 }
1415 //save mean power values
1416 {
1417 vector<r_4> tmp;
1418 //JEC 21/9/11 add meanpower for each Channels START
1419 vector<r_4> tmpCh0; //for Chan 0
1420 vector<r_4> tmpCh1; //for Chan 1
1421 //JEC 21/9/11 add meanpower for each Channels END
1422 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
1423 if (debuglev_>9) {
1424 cout << "Gain: save mean power of file: " << *iFile << endl;
1425 }
1426 tmp.push_back(meanpowerCollect[*iFile]);
1427 //JEC 21/9/11 add meanpower for each Channels START
1428 stringstream tmp0;
1429 tmp0 << (sa_size_t)0;
1430 string tag0 = *iFile + "_" + tmp0.str();
1431 tmpCh0.push_back(meanPowerPerChanCollect[tag0]);
1432 if (NUMBER_OF_CHANNELS>1){
1433 stringstream tmp1;
1434 tmp1 << (sa_size_t)1;
1435 string tag1 = *iFile + "_" + tmp1.str();
1436 tmpCh1.push_back(meanPowerPerChanCollect[tag1]);
1437 }
1438 //JEC 21/9/11 add meanpower for each Channels END
1439 }
1440 string fileName;
1441 //TOBE FIXED fileName = "./" + sourceName_ + "/gain_monitor_" + dateOfRun_
1442 fileName = outputPath_ + "/gain_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1443 POutPersist fos(fileName);
1444 TVector<r_4> tvtmp(tmp);
1445 fos.PutObject(tvtmp,"gainmoni"); //Attention initialement le tag etait "monitor"...
1446 //JEC 21/9/11 add meanpower for each Channels START
1447 TVector<r_4> tvtmp0(tmpCh0);
1448 fos.PutObject(tvtmp0,"gainmoni0");
1449 if (NUMBER_OF_CHANNELS>1){
1450 TVector<r_4> tvtmp1(tmpCh1);
1451 fos.PutObject(tvtmp1,"gainmoni1");
1452 }
1453 //JEC 21/9/11 add meanpower for each Channels END
1454 }
1455
1456 cout << "OK gain finished" <<endl;
1457 return rc;
1458}//ProcessGain::processCmd
1459//
1460//----------------------------------------------
1461//
1462int ProcessCalibration::processCmd() throw(string) {
1463 int rc = 0;
1464 string msg = "";
1465
1466 try {
1467 rc = ProcessBase::processCmd();
1468 }
1469 catch (string s) {
1470 throw s;
1471 }
1472 if(debuglev_>0)cout << "Process Calibration with option "<< option_ << endl;
1473
1474 vector<string> modeList;
1475 modeList.push_back("On");
1476 modeList.push_back("Off");
1477
1478 r_8 t0absolute = -1; //TIMEWIN of the first spectra used
1479 r_8 timeTag0 = -1; //MEANTT, mean TIME TAG of the first paquet window
1480 bool first = true; // for initialisation
1481
1482 // Power Tuple
1483 // mode, cycle, time, {Ch0, ... ,ChN} mean poser in the range [f0-Bw/2, f0+Bw/2]
1484 vector<string> varPowerTupleName; //ntuple variable naming
1485 varPowerTupleName.push_back("mode");
1486 varPowerTupleName.push_back("cycle");
1487 varPowerTupleName.push_back("time");
1488 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
1489 stringstream tmp;
1490 tmp << iCh;
1491 varPowerTupleName.push_back("Ch"+tmp.str());
1492 }
1493 r_4 valPowerTuple[varPowerTupleName.size()];
1494 NTuple powerEvolution(varPowerTupleName);
1495
1496
1497 //-----------------
1498 //Start real process
1499 //------------------
1500 for (size_t iMode = 0; iMode < modeList.size(); ++iMode) {
1501 string mode = modeList[iMode];
1502 if(debuglev_>0)cout << "Process Calibration for Mode " << mode << endl;
1503
1504 //initialize the start of each calibration procedure given by the RT SCA file
1505 //see ProcessBase::processCmd for the structure of the sca file
1506 string scaStartCalibName = "stcal"+mode;
1507 sa_size_t idStartCalib = scaTuple_->ColumnIndex(scaStartCalibName);
1508
1509 string directoryName;
1510 list<string> listOfSpecFiles;
1511 list<string>::const_iterator iFile, iFileEnd;
1512 //
1513 //loop on cycles
1514 //
1515 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
1516
1517 directoryName = inputPath_ + "/"
1518 + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
1519 stringstream sicycle;
1520 sicycle << icycle;
1521 directoryName += spectraDirectory_ + sicycle.str() + "/";
1522
1523 //read directory
1524 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
1525
1526 iFileEnd = listOfSpecFiles.end();
1527 DVList header;
1528 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
1529
1530 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
1531 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
1532 aSpectrum.GetHeaderRecords(header);
1533 aSpectrum >> spectre;
1534
1535 if(first){//initialise the timer
1536 first = false;
1537 t0absolute = header.GetD("TIMEWIN")/1000.;
1538 timeTag0 = header.GetD("MEANTT");
1539 if (debuglev_>1){
1540 cout << "debug Header of " << *iFile << endl;
1541 cout << "TIMEWIN = " << setprecision(12) << t0absolute << " "
1542 << "MEANTT = " << setprecision(12) << timeTag0 << endl;
1543 }
1544 }//end init timer
1545
1546 //Set time of current file
1547 //Add to the non-precise absolute origin an precise increment
1548 r_4 timeTag = t0absolute + (header.GetD("MEANTT") - timeTag0);
1549 int index=0;
1550 valPowerTuple[index++] = iMode;
1551 valPowerTuple[index++] = icycle;
1552 valPowerTuple[index++] = timeTag-scaTuple_->GetCell(idCycleInTuple_[icycle],idStartCalib); //take the RT time start as refernce
1553
1554 //Compute the mean power of the two channel (separatly) in the range [f-bw/2, f+bw/2]
1555 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
1556 TMatrix<r_4> tmp(spectre(Range(iCh),Range(lowerFreqBin_,upperFreqBin_)),false);
1557 r_4 mean = tmp.Sum()/(r_4)tmp.Size();
1558 valPowerTuple[index++] = mean;
1559 }//end of channel loop
1560
1561 //Fill Tuple
1562 powerEvolution.Fill(valPowerTuple);
1563 }//end of files loop
1564 }//end of cycles loop
1565 }//end of mode loop
1566
1567 //store power evolution Tuple
1568 if(debuglev_>0){
1569 cout << "ProcessCalibration::processCmd: dump powerEvolution tuple" << endl;
1570 powerEvolution.Show(cout);
1571 }
1572 //
1573 //Compute Calibration Coefficients as function of mode/cycle/channels
1574 //
1575
1576 //Tsys,Incremant Intensity... Tuple
1577 //mode mode cycle [(tsys0,dint0),...,(tsysN,dintN)]
1578 vector<string> varCalibTupleName;
1579 varCalibTupleName.push_back("mode");
1580 varCalibTupleName.push_back("cycle");
1581 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
1582 stringstream tmp;
1583 tmp << iCh;
1584 varCalibTupleName.push_back("tsys"+tmp.str());
1585 varCalibTupleName.push_back("dint"+tmp.str());
1586 }
1587 r_4 valCalibTuple[varCalibTupleName.size()];
1588 NTuple calibEvolution(varCalibTupleName);
1589
1590 // select time E [twin0,twin1] U [twin2,twin3] for Tsys
1591 // time unit = second
1592 const r_4 twin0 = -3.;
1593 const r_4 twin1 = -1.;
1594 const r_4 twin2 = 6.;
1595 const r_4 twin3 = 8;
1596 const r_4 thresholdFactor = 0.20; //80% of the diff. Max-Min intensity
1597
1598 sa_size_t inModeIdx = powerEvolution.ColumnIndex("mode");
1599 sa_size_t inCycleIdx= powerEvolution.ColumnIndex("cycle");
1600 sa_size_t inTimeIdx = powerEvolution.ColumnIndex("time");
1601 sa_size_t inOffsetCh0 = powerEvolution.ColumnIndex("Ch0"); //nb Ch0 position in the powerEvolution tuple
1602 if(debuglev_>1) cout << "DEBUG: in Idx: ("
1603 << inModeIdx << ", "
1604 << inCycleIdx << ", "
1605 << inTimeIdx << ", "
1606 << inOffsetCh0 << ")"
1607 << endl;
1608
1609
1610 size_t outModeIdx = calibEvolution.ColumnIndex("mode");
1611 size_t outCycleIdx= calibEvolution.ColumnIndex("cycle");
1612 size_t outOffsetCh0 = calibEvolution.ColumnIndex("tsys0"); //nb Ch0 position in the monitor tuple
1613 size_t outNDataPerCh= 2;
1614 if(debuglev_>1) cout << "DEBUG: out Idx: ("
1615 << outModeIdx << ", "
1616 << outCycleIdx << ", "
1617 << outOffsetCh0 << ")"
1618 << endl;
1619
1620 sa_size_t nPowerEvolEntry = powerEvolution.NEntry();
1621 double minMode;
1622 double maxMode;
1623 powerEvolution.GetMinMax("mode",minMode,maxMode);
1624 double minCycleNum;
1625 double maxCycleNum;
1626 powerEvolution.GetMinMax("cycle",minCycleNum,maxCycleNum);
1627 if(debuglev_>1) cout << "DEBUG mode ("<<minMode<<","<<maxMode<<")\n"
1628 << "cycle ("<<minCycleNum<<","<<maxCycleNum<<")"
1629 << endl;
1630
1631 r_4* aPowerEvolEntry; // a ligne of the powerEvolution tuple
1632// r_4* aPowerEvolEntry = new r_4[powerEvolution.NbColumns()]; // a ligne of the powerEvolution tuple
1633
1634 r_4 minMean;
1635 r_4 maxMean;
1636
1637 for (size_t iMode = (size_t)minMode; iMode <= (size_t)maxMode; iMode++){
1638 for (size_t iCycle = (size_t)minCycleNum; iCycle <= (size_t)maxCycleNum; iCycle++ ){
1639 if(debuglev_>1) cout<<"DEBUG >>>>>>> mode/cycle: " << iMode << "/" << iCycle << endl;
1640
1641 valCalibTuple[outModeIdx]=iMode;
1642 valCalibTuple[outCycleIdx]=iCycle;
1643 //
1644 //Compute the Min && Max value of each Ch<i> data
1645 if(debuglev_>1) cout<<"DEBUG compute Min and Max for each channels" << endl;
1646 //
1647 TVector<r_4> chMean[NUMBER_OF_CHANNELS];
1648 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1649 chMean[iCh].SetSize(nPowerEvolEntry);
1650 }
1651 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1652 aPowerEvolEntry = powerEvolution.GetVec(ie);
1653 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1654 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1655 chMean[iCh](ie) = aPowerEvolEntry[iCh+inOffsetCh0];
1656 }//end cut
1657 }
1658 }//eo loop on tuple entries
1659 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1660 //
1661 //select values to compute background Tsys
1662 if(debuglev_>1) {
1663 cout<< "DEBUG >>>> Ch[" << iCh << "]" << endl;
1664 cout<< "DEBUG select values to compute background Tsys " << endl;
1665 }
1666 //
1667
1668 std::vector<r_4> lowerInt;
1669 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1670 aPowerEvolEntry = powerEvolution.GetVec(ie);
1671 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1672 r_4 time = aPowerEvolEntry[inTimeIdx];
1673 if ( (time >= twin0 && time <= twin1) ||
1674 (time >= twin2 && time <= twin3)
1675 ) lowerInt.push_back(aPowerEvolEntry[iCh+inOffsetCh0]);
1676 }//end cut
1677 } //end loop entries
1678 //
1679 //compute the Tsys
1680 if(debuglev_>1) cout <<"DEBUG compute Tsys" << endl;
1681 //
1682 std::nth_element(lowerInt.begin(),
1683 lowerInt.begin()+lowerInt.size()/2,
1684 lowerInt.end());
1685 r_4 tsys = *(lowerInt.begin()+lowerInt.size()/2);
1686 if(debuglev_>1) cout << "DEBUG tsys["<<iCh<<"] : " << tsys <<endl;
1687 //
1688 //set the threshold for DAB detection
1689 //
1690 chMean[iCh].MinMax(minMean,maxMean);
1691 minMean = (tsys > minMean) ? tsys : minMean; //pb of empty vector
1692 if(debuglev_>1) cout << "DEBUG min["<<iCh<<"] : " << minMean
1693 << " max["<<iCh<<"] : " << maxMean
1694 <<endl;
1695 r_4 deltathres = thresholdFactor * (maxMean-minMean);
1696 r_4 thres = tsys + deltathres;
1697 if(debuglev_>1) cout << "DEBUG thres["<<iCh<<"] : " << thres <<endl;
1698 //
1699 //collect upper part of the DAB
1700 if(debuglev_>1) cout <<"DEBUG collect upper part of the DAB" << endl;
1701 //
1702 std::vector<r_4> upperInt;
1703 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1704 aPowerEvolEntry = powerEvolution.GetVec(ie);
1705 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1706 r_4 mean = aPowerEvolEntry[iCh+inOffsetCh0];
1707 if (mean >= thres) upperInt.push_back(mean);
1708 }//end cut
1709 }//eo loop on channels
1710 //
1711 //compute adjacent differences to detect the 2 DAB levels
1712 //
1713 if(debuglev_>1) cout << "(DEBUG )size upper [" << iCh << "]: " << upperInt.size() <<endl;
1714 std::vector<r_4> upperIntAdjDiff(upperInt.size());
1715 std::adjacent_difference(upperInt.begin(),
1716 upperInt.end(),
1717 upperIntAdjDiff.begin());
1718 //
1719 //Search the 2 DAB states
1720 if(debuglev_>1) cout<<"DEBUG Search the 2 DAB states" << endl;
1721 //
1722 std::vector<r_4> upIntState[2];
1723 int state=-1;
1724 for (size_t i=1;i<upperInt.size();i++) {//skip first value
1725 if (fabs(upperIntAdjDiff[i])<upperInt[i]*0.010) {
1726 if(state == -1) state=0;
1727 if(state == -2) state=1;
1728 upIntState[state].push_back(upperInt[i]);
1729 } else {
1730 if (state == 0) state=-2;
1731 }
1732 }
1733 //
1734 //Take the mean of the median values of each levels
1735 if(debuglev_>1)cout << "DEBUG mean of the median values of each levels" << endl;
1736 //
1737 r_4 meanUpper = 0;
1738 r_4 nval = 0;
1739 for (int i=0;i<2;i++) {
1740 if (!upIntState[i].empty()) {
1741// std::nth_element(upIntState[i].begin(),
1742// upIntState[i].begin()+upIntState[i].size()/2,
1743// upIntState[i].end());
1744// meanUpper += *(upIntState[i].begin()+upIntState[i].size()/2);
1745 meanUpper += median(upIntState[i].begin(),upIntState[i].end());
1746 nval++;
1747 }
1748 }
1749 meanUpper /= nval;
1750 //
1751 //Finaly the increase of power due to the DAB is:
1752 //
1753 r_4 deltaInt = meanUpper - tsys;
1754 if(debuglev_>1) cout << "DEBUG deltaInt["<<iCh<<"] : " << deltaInt <<endl;
1755 //
1756 //Save for monitoring and final calibration operations
1757 //
1758 valCalibTuple[outOffsetCh0+outNDataPerCh*iCh] = tsys;
1759 valCalibTuple[outOffsetCh0+outNDataPerCh*iCh+1] = deltaInt;
1760 }//end loop on channels
1761 if(debuglev_>1) cout<<"DEBUG Fill the calibEvolution tuple" << endl;
1762 calibEvolution.Fill(valCalibTuple);
1763 }//eo loop on Cycles
1764 }//eo loop on Modes
1765 //
1766 //store calibration evolution Tuple
1767 //
1768 if(debuglev_>0){
1769 cout << "ProcessCalibration::processCmd: dump calibEvolution tuple" << endl;
1770 calibEvolution.Show(cout);
1771 }
1772 //
1773 //Compute the means per channel of the calibration coefficients (deltaInt)
1774 //and save cycle based Calibration Ctes in file named as
1775 // <source>-<date>-<mode>-<fcalib>MHz-Ch<channel>Cycles.txt
1776 // format <cycle> <coef>
1777 //the means are stored in files
1778 // <source>-<date>-<mode>-<fcalib>MHz-All.txt
1779 // format <channel> <coef>
1780 //
1781 sa_size_t nModes = (sa_size_t)(maxMode - minMode) + 1;
1782 string calibCoefFileName;
1783 ofstream calibMeanCoefFile[nModes]; //Mean over cycles
1784 ofstream calibCoefFile[nModes][NUMBER_OF_CHANNELS]; //cycle per cycle
1785 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1786 //The file for the Means Coeff.
1787 calibCoefFileName = outputPath_ + "/calib_"
1788 + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
1789 + modeList[iMode] + "_"
1790 + freqBAOCalibration_ + "MHz-All.txt";
1791 calibMeanCoefFile[iMode].open(calibCoefFileName.c_str());
1792 //The file for the cycle per cycle Coeff.
1793 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1794 stringstream chString;
1795 chString << iCh;
1796 calibCoefFileName = outputPath_ + "/calib_"
1797 + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
1798 + modeList[iMode] + "_"
1799 + freqBAOCalibration_ + "MHz-Ch" + chString.str() + "Cycles.txt";
1800 calibCoefFile[iMode][iCh].open(calibCoefFileName.c_str());
1801 }
1802 }
1803
1804 r_4* aCalibEvolEntry;
1805 sa_size_t nCalibEvolEntry = calibEvolution.NEntry();
1806
1807 TMatrix<r_4> meanCalibCoef(nModes,NUMBER_OF_CHANNELS); //by default init to 0
1808 TMatrix<r_4> nData4Mean(nModes,NUMBER_OF_CHANNELS);
1809
1810 //READ the calibration tuple, fill the array for mean and write to ascii file
1811 for (sa_size_t ie=0; ie<nCalibEvolEntry; ie++){
1812 aCalibEvolEntry = calibEvolution.GetVec(ie);
1813 if(debuglev_>1){
1814 cout << "DEBUG mode/cycle/Coef "
1815 << aCalibEvolEntry[outModeIdx] << " "
1816 << aCalibEvolEntry[outCycleIdx] << " ";
1817 }
1818 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1819 if(debuglev_>1) cout << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1] << " "; // for debug
1820 sa_size_t currentMode = (sa_size_t)(aCalibEvolEntry[outModeIdx] - minMode); //ok even if minMode <> 0
1821 nData4Mean(currentMode,iCh)++;
1822 meanCalibCoef(currentMode,iCh) += aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1];
1823
1824 calibCoefFile[currentMode][iCh] << aCalibEvolEntry[outCycleIdx] << " "
1825 << setprecision(12)
1826 << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1]
1827 << endl;
1828 }
1829 if(debuglev_>1) cout << endl; //for debug
1830 }
1831
1832 //Perform means: true means that div per 0 is treated (slower but safer)
1833 meanCalibCoef.Div(nData4Mean,true);
1834
1835 if(debuglev_>0){
1836 cout << "DEBUG nData4Mean" << endl;
1837 nData4Mean.Print(cout);
1838 cout << "DEBUG meanCalibCoef (averaged)" << endl;
1839 meanCalibCoef.Print(cout);
1840 }
1841
1842 //Save means Coef
1843 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1844 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1845 calibMeanCoefFile[iMode] << setprecision(12)
1846 << meanCalibCoef(iMode,iCh)
1847 << endl;
1848 }
1849 }
1850
1851 //Save Monitor File
1852 {
1853 string fileName = outputPath_ + "/calib_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1854 POutPersist calibMonitorFile(fileName);
1855 calibMonitorFile << PPFNameTag("powermoni") << powerEvolution;
1856 calibMonitorFile << PPFNameTag("calibmoni") << calibEvolution;
1857 }
1858
1859 //Clean & Return
1860 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1861 calibMeanCoefFile[iMode].close();
1862 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1863 calibCoefFile[iMode][iCh].close();
1864 }
1865 }
1866
1867 cout << "Ok calibration finished" << endl;
1868 return rc;
1869}//ProcessCalibration::processCmd
1870
Note: See TracBrowser for help on using the repository browser.