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

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

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