source: BAORadio/AmasNancay/analyse.cc@ 513

Last change on this file since 513 was 507, checked in by campagne, 14 years ago

New import
File size: 42.4 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>: data, 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=data\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 fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAOCalib" + ".ppf";
735 POutPersist fos(fileName);
736 iSpectreEnd = spectreCollect.end();
737 id=0;
738 for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
739 tag = "specBAOCalib";
740 stringstream sid;
741 sid << id;
742 tag += sid.str();
743 fos << PPFNameTag(tag) << (*iSpectre).second;
744 }
745 }//end of save fits
746
747
748 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
749 string mode = *iMode;
750 if(debuglev_>0)cout << "Process CALIB RT Mode " << mode << endl;
751 //------------------------------------------
752 // Correct BAO calib spectra for RT calibration
753 //------------------------------------------
754 //Very Preliminary May-June 11
755 //coef RT @ 1346MHz Ouest - Est associees a Ch 0 et 1
756 r_4 calibRT[NUMBER_OF_CHANNELS] = {27.724, 22.543};
757 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
758 TMatrix<r_4> specmtx(spectreCollect[make_pair(mode,icycle)],true); //share the memory
759 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){
760 specmtx( Range(iCh), Range::all() ) *= calibRT[iCh];
761 }
762 }
763 }//end loop on mode RT calib
764
765 {//save mean spectra BAO & RT calibrated on file
766 if(debuglev_>0)cout << "save calibrated BAO & RT spectra" << endl;
767 string fileName;
768 fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataBAORTCalib" + ".ppf";
769 POutPersist fos(fileName);
770 iSpectreEnd = spectreCollect.end();
771 id = 0;
772 for (iSpectre = spectreCollect.begin();iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
773 tag = "specBAORTCalib";
774 stringstream sid;
775 sid << id;
776 tag += sid.str();
777 fos << PPFNameTag(tag) << (*iSpectre).second;
778 }
779 }//end of save fits
780
781 //------------------------------------------
782 // Perform ON-OFF
783 //------------------------------------------
784
785 map< sa_size_t, TMatrix<r_4> > diffCollect;
786 map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd;
787
788 TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero
789 r_4 nCycles = 0;
790 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
791 nCycles++;
792 TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory
793 TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory
794 TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff;
795 diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) ));
796 diffMeanOnOff += diffOnOff;
797 }//end loop on cycle
798 if(nCycles>0) diffMeanOnOff/=nCycles;
799
800 //exctract channels and do the mean
801 TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0
802 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) {
803 meanOfChan += diffMeanOnOff.Row(iCh).Transpose();
804 }
805 meanOfChan /= (r_4)NUMBER_OF_CHANNELS;
806
807
808
809 {//save diff ON-OFF BAO & RT calibrated
810 if(debuglev_>0)cout << "save ON-OFF spectra" << endl;
811 string fileName;
812 fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "diffOnOff" + ".ppf";
813 POutPersist fos(fileName);
814 iDiffEnd = diffCollect.end();
815 id = 0;
816 for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) {
817 tag = "specONOFF";
818 stringstream sid;
819 sid << id;
820 tag += sid.str();
821 fos << PPFNameTag(tag) << (*iDiff).second;
822 }
823 //save the mean also
824 fos << PPFNameTag("specONOFFMean") << diffMeanOnOff;
825 fos << PPFNameTag("specONOFF2ChanMean") << meanOfChan;
826 }//end of save fits
827
828 return rc;
829} //ProcessONOFFData::processCmd
830//
831//----------------------------------------------
832//
833int ProcessGain::processCmd() throw(string) {
834 int rc = 0;
835 string msg = "";
836
837 try {
838 rc = ProcessBase::processCmd();
839 }
840 catch (string s) {
841 throw s;
842 }
843 if(debuglev_>0)cout << "Process Gain" << endl;
844
845 string directoryName;
846 //TOBE FIXED directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode_ + "/";
847 //JEC 6/09/2011 numcycle_ repalced by ifirstCycle_ according to definition of numcycle_ and the fact that for GAIN 1 cycle is involved
848 stringstream thegaincycle;
849 thegaincycle << ifirstCycle_;
850 directoryName = inputPath_ + "/"
851 + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
852 + mode_ + "/" + spectraDirectory_ + thegaincycle.str() + "/";
853
854 list<string> listOfSpecFiles;
855 list<string>::const_iterator iFile, iFileEnd;
856 //read directory
857
858 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
859
860 //Mean power computed over the N channels to select the spectra for gain computation
861 //The threshold is computed "on-line" as 1% of the difference (max power -min power) over the
862 // the min power.
863 //A possible alternative is to set an absolute value...
864 if(debuglev_>0)cout << "compute mean poser spectra for files in " << directoryName << endl;
865 iFileEnd = listOfSpecFiles.end();
866 map<string, r_4> meanpowerCollect;
867 map<string, r_4>::const_iterator iMeanPow, iMeanPowEnd;
868
869 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
870 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
871 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
872 aSpectrum >> spectre;
873 meanpowerCollect[*iFile] = spectre.Sum()/spectre.Size();
874 }//end of for files
875 pair<string, r_4> minelement = *min_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
876 pair<string, r_4> maxelement = *max_element(meanpowerCollect.begin(),meanpowerCollect.end(),compare);
877 if(debuglev_>1){
878 cout << "meanpowerCollect has " << meanpowerCollect.size() << " spectra registered" << endl;
879 cout << "find min mean power "<<minelement.second << " for " << minelement.first << endl;
880 cout << "find max mean power "<<maxelement.second << " for " << maxelement.first << endl;
881 }
882 r_4 threshold = minelement.second + 0.01*(maxelement.second-minelement.second);
883 if(debuglev_>1){
884 cout << "threshold found at " << threshold <<endl;
885 }
886
887 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
888 r_4 nSpectres = 0;
889 iMeanPowEnd = meanpowerCollect.end();
890 for (iMeanPow = meanpowerCollect.begin(); iMeanPow != iMeanPowEnd; ++iMeanPow) {
891 if ( iMeanPow->second <= threshold ) {
892 //TODO avoid the reloading of the file may speed up
893 FitsInOutFile aSpectrum(iMeanPow->first,FitsInOutFile::Fits_RO);
894 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
895 aSpectrum >> spectre;
896 spectreMean += spectre;
897 nSpectres++;
898 }
899 }
900 if(debuglev_>1){
901 cout << "Number of files selected for gain " << nSpectres <<endl;
902 }
903 if(nSpectres>0) {
904 spectreMean /= nSpectres;
905 } else {
906 stringstream tmp;
907 tmp << threshold;
908 msg = "Gain: cannot find a spectra above threshold value =" + tmp.str() + " ... FATAL";
909 throw msg;
910 }
911
912 //Save gain spectra
913 {
914 //use ! to override the file: special features of cfitsio library
915 string fileName;
916 //TOBE FIXED fileName = "!./" + sourceName_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
917 fileName = "!"+ outputPath_ + "/gain_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".fits";
918 if(debuglev_>1){
919 cout << "save gains in " << fileName << endl;
920 }
921 FitsInOutFile fos(fileName, FitsInOutFile::Fits_Create);
922 fos << spectreMean;
923 }
924 //save mean power values
925 {
926 vector<r_4> tmp;
927 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
928 tmp.push_back(meanpowerCollect[*iFile]);
929 }
930 string fileName;
931 //TOBE FIXED fileName = "./" + sourceName_ + "/gain_monitor_" + dateOfRun_
932 fileName = outputPath_ + "/gain_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
933 POutPersist fos(fileName);
934 TVector<r_4> tvtmp(tmp);
935 fos.PutObject(tvtmp,"gainmoni");
936 }
937
938
939 cout << "OK gain finished" <<endl;
940 return rc;
941}//ProcessGain::processCmd
942//
943//----------------------------------------------
944//
945int ProcessCalibration::processCmd() throw(string) {
946 int rc = 0;
947 string msg = "";
948
949 try {
950 rc = ProcessBase::processCmd();
951 }
952 catch (string s) {
953 throw s;
954 }
955 if(debuglev_>0)cout << "Process Calibration with option "<< option_ << endl;
956
957 vector<string> modeList;
958 modeList.push_back("On");
959 modeList.push_back("Off");
960
961 r_8 t0absolute = -1; //TIMEWIN of the first spectra used
962 r_8 timeTag0 = -1; //MEANTT, mean TIME TAG of the first paquet window
963 bool first = true; // for initialisation
964
965 // Power Tuple
966 // mode, cycle, time, {Ch0, ... ,ChN} mean poser in the range [f0-Bw/2, f0+Bw/2]
967 vector<string> varPowerTupleName; //ntuple variable naming
968 varPowerTupleName.push_back("mode");
969 varPowerTupleName.push_back("cycle");
970 varPowerTupleName.push_back("time");
971 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
972 stringstream tmp;
973 tmp << iCh;
974 varPowerTupleName.push_back("Ch"+tmp.str());
975 }
976 r_4 valPowerTuple[varPowerTupleName.size()];
977 NTuple powerEvolution(varPowerTupleName);
978
979
980 //-----------------
981 //Start real process
982 //------------------
983 for (size_t iMode = 0; iMode < modeList.size(); ++iMode) {
984 string mode = modeList[iMode];
985 if(debuglev_>0)cout << "Process Calibration for Mode " << mode << endl;
986
987 //initialize the start of each calibration procedure given by the RT SCA file
988 //see ProcessBase::processCmd for the structure of the sca file
989 string scaStartCalibName = "stcal"+mode;
990 sa_size_t idStartCalib = scaTuple_->ColumnIndex(scaStartCalibName);
991
992 string directoryName;
993 list<string> listOfSpecFiles;
994 list<string>::const_iterator iFile, iFileEnd;
995 //
996 //loop on cycles
997 //
998 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
999
1000 directoryName = inputPath_ + "/"
1001 + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";
1002 stringstream sicycle;
1003 sicycle << icycle;
1004 directoryName += spectraDirectory_ + sicycle.str() + "/";
1005
1006 //read directory
1007 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_);
1008
1009 iFileEnd = listOfSpecFiles.end();
1010 DVList header;
1011 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
1012
1013 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) {
1014 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
1015 aSpectrum.GetHeaderRecords(header);
1016 aSpectrum >> spectre;
1017
1018 if(first){//initialise the timer
1019 first = false;
1020 t0absolute = header.GetD("TIMEWIN")/1000.;
1021 timeTag0 = header.GetD("MEANTT");
1022 if (debuglev_>1){
1023 cout << "debug Header of " << *iFile << endl;
1024 cout << "TIMEWIN = " << setprecision(12) << t0absolute << " "
1025 << "MEANTT = " << setprecision(12) << timeTag0 << endl;
1026 }
1027 }//end init timer
1028
1029 //Set time of current file
1030 //Add to the non-precise absolute origin an precise increment
1031 r_4 timeTag = t0absolute + (header.GetD("MEANTT") - timeTag0);
1032 int index=0;
1033 valPowerTuple[index++] = iMode;
1034 valPowerTuple[index++] = icycle;
1035 valPowerTuple[index++] = timeTag-scaTuple_->GetCell(idCycleInTuple_[icycle],idStartCalib); //take the RT time start as refernce
1036
1037 //Compute the mean power of the two channel (separatly) in the range [f-bw/2, f+bw/2]
1038 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
1039 TMatrix<r_4> tmp(spectre(Range(iCh),Range(lowerFreqBin_,upperFreqBin_)),false);
1040 r_4 mean = tmp.Sum()/(r_4)tmp.Size();
1041 valPowerTuple[index++] = mean;
1042 }//end of channel loop
1043
1044 //Fill Tuple
1045 powerEvolution.Fill(valPowerTuple);
1046 }//end of files loop
1047 }//end of cycles loop
1048 }//end of mode loop
1049
1050 //store power evolution Tuple
1051 if(debuglev_>0){
1052 cout << "ProcessCalibration::processCmd: dump powerEvolution tuple" << endl;
1053 powerEvolution.Show(cout);
1054 }
1055 //
1056 //Compute Calibration Coefficients as function of mode/cycle/channels
1057 //
1058
1059 //Tsys,Incremant Intensity... Tuple
1060 //mode mode cycle [(tsys0,dint0),...,(tsysN,dintN)]
1061 vector<string> varCalibTupleName;
1062 varCalibTupleName.push_back("mode");
1063 varCalibTupleName.push_back("cycle");
1064 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS;++iCh){
1065 stringstream tmp;
1066 tmp << iCh;
1067 varCalibTupleName.push_back("tsys"+tmp.str());
1068 varCalibTupleName.push_back("dint"+tmp.str());
1069 }
1070 r_4 valCalibTuple[varCalibTupleName.size()];
1071 NTuple calibEvolution(varCalibTupleName);
1072
1073 // select time E [twin0,twin1] U [twin2,twin3] for Tsys
1074 // time unit = second
1075 const r_4 twin0 = -3.;
1076 const r_4 twin1 = -1.;
1077 const r_4 twin2 = 6.;
1078 const r_4 twin3 = 8;
1079 const r_4 thresholdFactor = 0.80; //80% of the diff. Max-Min intensity
1080
1081 sa_size_t inModeIdx = powerEvolution.ColumnIndex("mode");
1082 sa_size_t inCycleIdx= powerEvolution.ColumnIndex("cycle");
1083 sa_size_t inTimeIdx = powerEvolution.ColumnIndex("time");
1084 sa_size_t inOffsetCh0 = powerEvolution.ColumnIndex("Ch0"); //nb Ch0 position in the powerEvolution tuple
1085 if(debuglev_>1) cout << "DEBUG: in Idx: ("
1086 << inModeIdx << ", "
1087 << inCycleIdx << ", "
1088 << inTimeIdx << ", "
1089 << inOffsetCh0 << ")"
1090 << endl;
1091
1092
1093 size_t outModeIdx = calibEvolution.ColumnIndex("mode");
1094 size_t outCycleIdx= calibEvolution.ColumnIndex("cycle");
1095 size_t outOffsetCh0 = calibEvolution.ColumnIndex("tsys0"); //nb Ch0 position in the monitor tuple
1096 size_t outNDataPerCh= 2;
1097 if(debuglev_>1) cout << "DEBUG: out Idx: ("
1098 << outModeIdx << ", "
1099 << outCycleIdx << ", "
1100 << outOffsetCh0 << ")"
1101 << endl;
1102
1103 sa_size_t nPowerEvolEntry = powerEvolution.NEntry();
1104 double minMode;
1105 double maxMode;
1106 powerEvolution.GetMinMax("mode",minMode,maxMode);
1107 double minCycleNum;
1108 double maxCycleNum;
1109 powerEvolution.GetMinMax("cycle",minCycleNum,maxCycleNum);
1110 if(debuglev_>1) cout << "DEBUG mode ("<<minMode<<","<<maxMode<<")\n"
1111 << "cycle ("<<minCycleNum<<","<<maxCycleNum<<")"
1112 << endl;
1113
1114 r_4* aPowerEvolEntry; // a ligne of the powerEvolution tuple
1115// r_4* aPowerEvolEntry = new r_4[powerEvolution.NbColumns()]; // a ligne of the powerEvolution tuple
1116
1117 r_4 minMean;
1118 r_4 maxMean;
1119
1120 for (size_t iMode = (size_t)minMode; iMode <= (size_t)maxMode; iMode++){
1121 for (size_t iCycle = (size_t)minCycleNum; iCycle <= (size_t)maxCycleNum; iCycle++ ){
1122 if(debuglev_>1) cout<<"DEBUG >>>>>>> mode/cycle: " << iMode << "/" << iCycle << endl;
1123
1124 valCalibTuple[outModeIdx]=iMode;
1125 valCalibTuple[outCycleIdx]=iCycle;
1126 //
1127 //Compute the Min && Max value of each Ch<i> data
1128 if(debuglev_>1) cout<<"DEBUG compute Min and Max for each channels" << endl;
1129 //
1130 TVector<r_4> chMean[NUMBER_OF_CHANNELS];
1131 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1132 chMean[iCh].SetSize(nPowerEvolEntry);
1133 }
1134 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1135 aPowerEvolEntry = powerEvolution.GetVec(ie);
1136 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1137 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1138 chMean[iCh](ie) = aPowerEvolEntry[iCh+inOffsetCh0];
1139 }//end cut
1140 }
1141 }//eo loop on tuple entries
1142 for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;iCh++){
1143 //
1144 //select values to compute background Tsys
1145 if(debuglev_>1) {
1146 cout<< "DEBUG >>>> Ch[" << iCh << "]" << endl;
1147 cout<< "DEBUG select values to compute background Tsys " << endl;
1148 }
1149 //
1150
1151 std::vector<r_4> lowerInt;
1152 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1153 aPowerEvolEntry = powerEvolution.GetVec(ie);
1154 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1155 r_4 time = aPowerEvolEntry[inTimeIdx];
1156 if ( (time >= twin0 && time <= twin1) ||
1157 (time >= twin2 && time <= twin3)
1158 ) lowerInt.push_back(aPowerEvolEntry[iCh+inOffsetCh0]);
1159 }//end cut
1160 } //end loop entries
1161 //
1162 //compute the Tsys
1163 if(debuglev_>1) cout <<"DEBUG compute Tsys" << endl;
1164 //
1165 std::nth_element(lowerInt.begin(),
1166 lowerInt.begin()+lowerInt.size()/2,
1167 lowerInt.end());
1168 r_4 tsys = *(lowerInt.begin()+lowerInt.size()/2);
1169 if(debuglev_>1) cout << "DEBUG tsys["<<iCh<<"] : " << tsys <<endl;
1170 //
1171 //set the threshold for DAB detection
1172 //
1173 chMean[iCh].MinMax(minMean,maxMean);
1174 minMean = (tsys > minMean) ? tsys : minMean; //pb of empty vector
1175 if(debuglev_>1) cout << "DEBUG min["<<iCh<<"] : " << minMean
1176 << " max["<<iCh<<"] : " << maxMean
1177 <<endl;
1178 r_4 deltathres = thresholdFactor * (maxMean-minMean);
1179 r_4 thres = tsys + deltathres;
1180 if(debuglev_>1) cout << "DEBUG thres["<<iCh<<"] : " << thres <<endl;
1181 //
1182 //collect upper part of the DAB
1183 if(debuglev_>1) cout <<"DEBUG collect upper part of the DAB" << endl;
1184 //
1185 std::vector<r_4> upperInt;
1186 for (sa_size_t ie=0; ie<nPowerEvolEntry; ie++){
1187 aPowerEvolEntry = powerEvolution.GetVec(ie);
1188 if ( (size_t)aPowerEvolEntry[inModeIdx] == iMode && (size_t)aPowerEvolEntry[inCycleIdx] == iCycle ){
1189 r_4 mean = aPowerEvolEntry[iCh+inOffsetCh0];
1190 if (mean >= thres) upperInt.push_back(mean);
1191 }//end cut
1192 }//eo loop on channels
1193 //
1194 //compute adjacent differences to detect the 2 DAB levels
1195 //
1196 if(debuglev_>1) cout << "(DEBUG )size upper [" << iCh << "]: " << upperInt.size() <<endl;
1197 std::vector<r_4> upperIntAdjDiff(upperInt.size());
1198 std::adjacent_difference(upperInt.begin(),
1199 upperInt.end(),
1200 upperIntAdjDiff.begin());
1201 //
1202 //Search the 2 DAB states
1203 if(debuglev_>1) cout<<"DEBUG Search the 2 DAB states" << endl;
1204 //
1205 std::vector<r_4> upIntState[2];
1206 int state=-1;
1207 for (size_t i=1;i<upperInt.size();i++) {//skip first value
1208 if (fabs(upperIntAdjDiff[i])<upperInt[i]*0.010) {
1209 if(state == -1) state=0;
1210 if(state == -2) state=1;
1211 upIntState[state].push_back(upperInt[i]);
1212 } else {
1213 if (state == 0) state=-2;
1214 }
1215 }
1216 //
1217 //Take the mean of the median values of each levels
1218 if(debuglev_>1)cout << "DEBUG mean of the median values of each levels" << endl;
1219 //
1220 r_4 meanUpper = 0;
1221 r_4 nval = 0;
1222 for (int i=0;i<2;i++) {
1223 if (!upIntState[i].empty()) {
1224 std::nth_element(upIntState[i].begin(),
1225 upIntState[i].begin()+upIntState[i].size()/2,
1226 upIntState[i].end());
1227 meanUpper += *(upIntState[i].begin()+upIntState[i].size()/2);
1228 nval++;
1229 }
1230 }
1231 meanUpper /= nval;
1232 //
1233 //Finaly the increase of power due to the DAB is:
1234 //
1235 r_4 deltaInt = meanUpper - tsys;
1236 if(debuglev_>1) cout << "DEBUG deltaInt["<<iCh<<"] : " << deltaInt <<endl;
1237 //
1238 //Save for monitoring and final calibration operations
1239 //
1240 valCalibTuple[outOffsetCh0+outNDataPerCh*iCh] = tsys;
1241 valCalibTuple[outOffsetCh0+outNDataPerCh*iCh+1] = deltaInt;
1242 }//end loop on channels
1243 if(debuglev_>1) cout<<"DEBUG Fill the calibEvolution tuple" << endl;
1244 calibEvolution.Fill(valCalibTuple);
1245 }//eo loop on Cycles
1246 }//eo loop on Modes
1247 //
1248 //store calibration evolution Tuple
1249 //
1250 if(debuglev_>0){
1251 cout << "ProcessCalibration::processCmd: dump calibEvolution tuple" << endl;
1252 calibEvolution.Show(cout);
1253 }
1254 //
1255 //Compute the means per channel of the calibration coefficients (deltaInt)
1256 //and save cycle based Calibration Ctes in file named as
1257 // <source>-<date>-<mode>-<fcalib>MHz-Ch<channel>Cycles.txt
1258 // format <cycle> <coef>
1259 //the means are stored in files
1260 // <source>-<date>-<mode>-<fcalib>MHz-All.txt
1261 // format <channel> <coef>
1262 //
1263 sa_size_t nModes = (sa_size_t)(maxMode - minMode) + 1;
1264 string calibCoefFileName;
1265 ofstream calibMeanCoefFile[nModes]; //Mean over cycles
1266 ofstream calibCoefFile[nModes][NUMBER_OF_CHANNELS]; //cycle per cycle
1267 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1268 //The file for the Means Coeff.
1269 calibCoefFileName = outputPath_ + "/calib_"
1270 + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
1271 + modeList[iMode] + "_"
1272 + freqBAOCalibration_ + "MHz-All.txt";
1273 calibMeanCoefFile[iMode].open(calibCoefFileName.c_str());
1274 //The file for the cycle per cycle Coeff.
1275 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1276 stringstream chString;
1277 chString << iCh;
1278 calibCoefFileName = outputPath_ + "/calib_"
1279 + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"
1280 + modeList[iMode] + "_"
1281 + freqBAOCalibration_ + "MHz-Ch" + chString.str() + "Cycles.txt";
1282 calibCoefFile[iMode][iCh].open(calibCoefFileName.c_str());
1283 }
1284 }
1285
1286 r_4* aCalibEvolEntry;
1287 sa_size_t nCalibEvolEntry = calibEvolution.NEntry();
1288
1289 TMatrix<r_4> meanCalibCoef(nModes,NUMBER_OF_CHANNELS); //by default init to 0
1290 TMatrix<r_4> nData4Mean(nModes,NUMBER_OF_CHANNELS);
1291
1292 //READ the calibration tuple, fill the array for mean and write to ascii file
1293 for (sa_size_t ie=0; ie<nCalibEvolEntry; ie++){
1294 aCalibEvolEntry = calibEvolution.GetVec(ie);
1295 if(debuglev_>1){
1296 cout << "DEBUG mode/cycle/Coef "
1297 << aCalibEvolEntry[outModeIdx] << " "
1298 << aCalibEvolEntry[outCycleIdx] << " ";
1299 }
1300 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1301 if(debuglev_>1) cout << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1] << " "; // for debug
1302 sa_size_t currentMode = (sa_size_t)(aCalibEvolEntry[outModeIdx] - minMode); //ok even if minMode <> 0
1303 nData4Mean(currentMode,iCh)++;
1304 meanCalibCoef(currentMode,iCh) += aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1];
1305
1306 calibCoefFile[currentMode][iCh] << aCalibEvolEntry[outCycleIdx] << " "
1307 << setprecision(12)
1308 << aCalibEvolEntry[outOffsetCh0+outNDataPerCh*iCh+1]
1309 << endl;
1310 }
1311 if(debuglev_>1) cout << endl; //for debug
1312 }
1313
1314 //Perform means: true means that div per 0 is treated (slower but safer)
1315 meanCalibCoef.Div(nData4Mean,true);
1316
1317 if(debuglev_>0){
1318 cout << "DEBUG nData4Mean" << endl;
1319 nData4Mean.Print(cout);
1320 cout << "DEBUG meanCalibCoef (averaged)" << endl;
1321 meanCalibCoef.Print(cout);
1322 }
1323
1324 //Save means Coef
1325 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1326 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1327 calibMeanCoefFile[iMode] << setprecision(12)
1328 << meanCalibCoef(iMode,iCh)
1329 << endl;
1330 }
1331 }
1332
1333 //Save Monitor File
1334 {
1335 string fileName = outputPath_ + "/calib_monitor_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
1336 POutPersist calibMonitorFile(fileName);
1337 calibMonitorFile << PPFNameTag("powermoni") << powerEvolution;
1338 calibMonitorFile << PPFNameTag("calibmoni") << calibEvolution;
1339 }
1340
1341 //Clean & Return
1342 for (sa_size_t iMode=0; iMode<nModes; iMode++){
1343 calibMeanCoefFile[iMode].close();
1344 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){
1345 calibCoefFile[iMode][iCh].close();
1346 }
1347 }
1348
1349 cout << "Ok calibration finished" << endl;
1350 return rc;
1351}//ProcessCalibration::processCmd
1352
Note: See TracBrowser for help on using the repository browser.