source: BAORadio/AmasNancay/analyse.cc@ 526

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

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