source: BAORadio/AmasNancay/v1/analyse.cc@ 650

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

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