source: BAORadio/AmasNancay/trunk/analyse.cc@ 555

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

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