source: Sophya/trunk/AddOn/TAcq/branap.cc@ 3993

Last change on this file since 3993 was 3993, checked in by ansari, 14 years ago

Ajout possibilite de faire un DataTable (NTuple) avec la puissance spectrale dans plusieurs bandes en fonction du temps, Reza 13/05/2011

File size: 12.3 KB
Line 
1#include <stdlib.h>
2#include <string.h>
3
4#include "branap.h"
5#include "minifits.h"
6#include "strutilxx.h"
7#include "sopnamsp.h"
8
9//--------------------------------------------------------------
10// Projet BAORadio - (C) LAL/IRFU 2008-2010
11// Classe de gestion des parametres programmes d'analyse
12//--------------------------------------------------------------
13
14/* --Methode-- */
15BRAnaParam::BRAnaParam(uint_4 nmean, uint_4 nzon, uint_4 npaqz)
16{
17 outpath_="./";
18 fgfitsout_=false;
19 nmean_=nmean;
20 nbloc_=1;
21 imin_=imax_=0;
22 istep_=1;
23 rdsamefc_=true; // read paquets with same frame counter
24 freqmin_=0; freqmax_=-1;
25 nbinfreq_=1;
26 paqsize_=16424;
27 nzones_=nzon;
28 npaqinzone_=npaqz;
29 fgdatafft_=false; fgforceraworfft_=false;
30 fgsinglechannel_=false; fgforcesingleortwochan_=false;
31 prtlevel_=0;
32 prtmodulo_=50000;
33 nbcalgrp_=1;
34 nthreads_=1;
35
36 spec_win_sz_=0;
37 spw_ext_sz_=0;
38 nbmax_specwfiles_=0;
39
40 vmin_=0.; vmax_=9e99;
41 nbands_=0; bandfirst_ = bandlast_ = 0;
42 fgdtpaq_ = fgdtms_ = false;
43
44 fgfreqfilter_=false; //JEC 1/2/11
45 medhalfwidth_=50; //JEC 6/4/11
46
47 gainfile_="";
48
49 proctimeduration_=9.e9;
50 fgtimeselect_=false;
51}
52
53/* --Methode-- */
54int BRAnaParam::DecodeArgs(int narg, char* arg[])
55{
56 if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
57 if (narg<5) return Usage(true);
58
59 bool okarg=false;
60 int ka=1;
61 while (ka<(narg-1)) {
62 if (strcmp(arg[ka],"-act")==0) {
63 action_=arg[ka+1];
64 ka+=2;
65 }
66 else if (strcmp(arg[ka],"-out")==0) {
67 outpath_=arg[ka+1];
68 size_t lenp=outpath_.size();
69 if ((lenp>0)&&(outpath_[lenp-1]!='/')) outpath_+='/';
70 ka+=2;
71 }
72 else if (strcmp(arg[ka],"-fitsout")==0) {
73 fgfitsout_=true;
74 ka++;
75 }
76 else if (strcmp(arg[ka],"-nmean")==0) {
77 nmean_=atoi(arg[ka+1]);
78 ka+=2;
79 }
80 else if (strcmp(arg[ka],"-nbloc")==0) {
81 nbloc_=atoi(arg[ka+1]);
82 ka+=2;
83 }
84 else if (strcmp(arg[ka],"-freq")==0) {
85 sscanf(arg[ka+1],"%d,%d,%d",&freqmin_,&freqmax_,&nbinfreq_);
86 ka+=2;
87 }
88 else if (strcmp(arg[ka],"-zones")==0) {
89 int nzon=4;
90 int npaqz=128;
91 sscanf(arg[ka+1],"%d,%d",&nzon,&npaqz);
92 nzones_=nzon; npaqinzone_=npaqz;
93 ka+=2;
94 }
95 else if (strcmp(arg[ka],"-prt")==0) {
96 sscanf(arg[ka+1],"%d,%ld",&prtlevel_,&prtmodulo_);
97 ka+=2;
98 }
99 else if (strcmp(arg[ka],"-nthr")==0) {
100 nthreads_=atoi(arg[ka+1]);
101 ka+=2;
102 }
103 else if (strcmp(arg[ka],"-nvcal")==0) {
104 nbcalgrp_=atoi(arg[ka+1]);
105 ka+=2;
106 }
107 else if (strcmp(arg[ka],"-nosfc")==0) {
108 rdsamefc_=false;
109 ka++;
110 }
111 else if (strcmp(arg[ka],"-singlechan")==0) {
112 fgsinglechannel_=true; fgforcesingleortwochan_ = true;
113 ka++;
114 }
115 else if (strcmp(arg[ka],"-twochan")==0) {
116 fgsinglechannel_=false; fgforcesingleortwochan_ = true;
117 ka++;
118 }
119 else if (strcmp(arg[ka],"-rawdata")==0) {
120 fgforceraworfft_=true; fgdatafft_ = false;
121 ka++;
122 }
123 else if (strcmp(arg[ka],"-fftdata")==0) {
124 fgforceraworfft_=true; fgdatafft_ = true;
125 ka++;
126 }
127 else if (strcmp(arg[ka],"-varcut")==0) {
128 sscanf(arg[ka+1],"%lg,%lg",&vmin_,&vmax_);
129 ka+=2;
130 }
131 else if (strcmp(arg[ka],"-nband")==0) {
132 sscanf(arg[ka+1],"%d,%d,%d",&nbands_,&bandfirst_,&bandlast_);
133 ka+=2;
134 }
135 else if (strcmp(arg[ka],"-fdtpaq")==0) {
136 fgdtpaq_=true;
137 ka++;
138 }
139 else if (strcmp(arg[ka],"-fdtms")==0) {
140 fgdtms_=true;
141 ka++;
142 }
143 else if (strcmp(arg[ka],"-freqfilter")==0) {
144 fgfreqfilter_=true; ka++;
145 if (strcmp(arg[ka],"-")!=0) medhalfwidth_=atof(arg[ka]);
146 ka++;
147 }
148 //-tmproc hh:mm:ss,nseconds
149 else if (strcmp(arg[ka],"-tmproc")==0) {
150 int ah=0,am=0;
151 double as=0, als=0;
152 sscanf(arg[ka+1],"%d:%d:%lg,%lg",&ah,&am,&as,&als);
153 fgtimeselect_=true; proctimeduration_=als;
154 proctimestart_.SetHour(ah,am,as);
155 ka+=2;
156 }
157 else if (strcmp(arg[ka],"-gain")==0) {
158 gainfile_=arg[ka+1];
159 ka+=2;
160 }
161 else if (strcmp(arg[ka],"-tspwin")==0) {
162 int ai1=0,ai2=0,ai3=0;
163 sscanf(arg[ka+1],"%d,%d,%d",&ai1,&ai2,&ai3);
164 spec_win_sz_=ai1; spw_ext_sz_=ai2; nbmax_specwfiles_=ai3;
165 ka+=2;
166 }
167 else if (strcmp(arg[ka],"-in")==0) {
168 if ((narg-ka)<4) {
169 cout << " BRAnaParam::DecodeArgs() / Argument error " << endl;
170 return Usage(true);
171 }
172 sscanf(arg[ka+1],"%d,%d,%d",&imin_,&imax_,&istep_); ka+=2;
173 while(ka<(narg-1)) {
174 string inpath = arg[ka];
175 size_t lenp=inpath.size();
176 if (lenp<1) inpath="./";
177 if ((lenp>0)&&(inpath[lenp-1]!='/')) inpath+='/';
178 vector<string> fiblist;
179 string sa1 = arg[ka+1];
180 FillVStringFrString(sa1, fiblist, ',');
181 char dbuff[32];
182 for(size_t i=0; i<fiblist.size(); i++) {
183 sprintf(dbuff,"Fiber%d/",(int)atoi(fiblist[i].c_str()));
184 dirlist_.push_back(inpath+dbuff);
185 }
186 ka += 2;
187 }
188 okarg=true;
189 }
190 else ka++;
191 }
192
193 if (!okarg) {
194 cout << " BRAnaParam::DecodeArgs() / Argument error " << endl;
195 return Usage(true);
196 }
197 return 0;
198}
199
200/* --Methode-- */
201int BRAnaParam::Usage(bool fgshort)
202{
203 cout << " --- BRAnaParam : Reading/Processing BAORadio FITS files parameters " << endl;
204 cout << " Usage: prgname [-act ACT] [-out OutPath] [-fitsout] \n"
205 << " [-nmean NMean] [-zones NZones,nPaqinZone] \n"
206 << " [-nbloc NBloc] [-freq NumFreqMin,NumFreqMax,NBinFreq] \n"
207 << " [-prt lev,modulo] [-nvcal n] [-nthr n] [-nosfc]\n"
208 << " [-singlechan] [-twochan] [-fftdata] [-rawdata] \n"
209 << " [-freqfilter medhw] [-gain filename] [-varcut min,max] [-nband nband,first,last] \n"
210 << " [-freqfilter] [-gain filename] [-varcut min,max] [-nband nband,first,last] \n"
211 << " [-tmproc hh:mm:ss,nseconds] [-fdtpaq] [-fdtms] [-tspwin wsz,extsz,nfiles] \n"
212 << " -in Imin,Imax,Istep InPath FiberList [InPath2 FiberList2 InPath3 FiberList3 ...] \n" << endl;
213 if (fgshort) {
214 cout << " prgname -h for detailed instructions" << endl;
215 return 5;
216 }
217 cout << " -act Action: cube3d or vis or viscktt or mspec \n"
218 << " cube3d: create 3D fits cubes \n"
219 << " vis: compute visibilites (vismfib program) \n"
220 << " viscktt: compute visibilities and check TimeTag/FrameCounter (vismfib program)\n"
221 << " mspec: compute and save mean spectra for each channel \n"
222 << " bproc: run BRBaseProcessor for debug/printing (use -prt) \n"
223 << " -out OutPath: Output directory name \n"
224 << " -fitsout : Force FITS format for output files \n"
225 << " -nmean NMean: Number of packet used for spectra/visibility computation \n"
226 << " -nbloc NBloc: Number of MemMgr blocs in output file\n"
227 << " -zones NZones,NbPaqinZone : Number of Zones and number of paquets in one zone (RAcqMemZoneMgr) \n"
228 << " -freq NumFreqMin,NumFreqMax,NBinFreq : Frequency zone and number of bins \n"
229 << " -prt lev,modulo : Print level (0,1,2...) and print counter modulo (10000, 50000 ...) \n"
230 << " -nvcal n : number of BRVisibilityCalculator objects running in parallel in BRVisCalcGroup (default=1) \n"
231 << " -nthr n : number of threads for parallel execution in BRVisibilityCalculator (default=1) \n"
232 << " -nosfc : Don't force reading with SAME FrameCounter \n"
233 << " -singlechan : Force one channel per fiber \n"
234 << " -twochan : Force two channels per fiber \n"
235 << " -rawdata : Force raw data mode (firmware raw) \n"
236 << " -fftdata : Force FFT data mode (firmware fft) \n"
237 << " -varcut min,max : min-max cut on variance \n"
238 << " -nband nband,first,last: numbers of freq. bands and first and last bands used for cuts\n"
239 << " -tmproc hh:mm:ss,nseconds : processing time window definition \n"
240 << " -fdtpaq : force per paquet data table filling (specmfib) \n"
241 << " -fdtms : force time averaged/filtered data table filling; use -freq for defining bands (specmfib)\n"
242 << " -freqfilter medhw: force median filtering on the frequencies \n"
243 << " with half window width medhw (use - for default=" << medhalfwidth_ << ") \n"
244 << " -gain filename : spectral response fits file name \n"
245 << " -tspwin wsz,extsz,nfiles : spectra time (paquet no) window (ex: -tspwin 120,4,5) \n"
246 << " wsz,extsz= window,extension size; nfiles= maximum number of windows saved \n"
247 << " -in : input files/directory definition : \n"
248 << " Imin,Imax,Istep: fits files signalII.fits Imin<=II<=Imax Istep=increment \n"
249 << " InPath: Input directerory fits files in InPath/FiberJJ directory \n"
250 << " FiberList: List of fiber numbers (JJ - Ex: 2,1,4 ) \n" << endl;
251 return 1;
252}
253
254/* --Fonction-- */
255int BRAnaParam::PaqSizeFromFits()
256{
257 uint_4 paqsz, npaqf;
258 char flnm[1024];
259 sprintf(flnm,"%s/signal%d.fits",dirlist_[0].c_str(),imin_);
260 bool fgdatafft_in_fits=false;
261 bool fgsinglechan_in_fits=false;
262 SOPHYA::TimeStamp tmstart;
263 int rc = DecodeMiniFitsHeader(flnm,paqsize_, npaqinfile_,fgdatafft_in_fits, fgsinglechan_in_fits,tmstart);
264 if (!fgforcesingleortwochan_) fgsinglechannel_=fgsinglechan_in_fits;
265 if (!fgforceraworfft_) fgdatafft_=fgdatafft_in_fits;
266 if (fgtimeselect_) {
267 int year,month,day;
268 tmstart.GetDate(year,month,day);
269 proctimestart_.SetDate(year,month,day);
270 proctimeend_.Set(proctimestart_.ToDays()+proctimeduration_/86400.);
271 }
272 return rc;
273}
274
275/* --Fonction-- */
276ostream& BRAnaParam::Print(ostream& os)
277{
278 os << " BRAnaParam::Print() dirlist_.size()=" << dirlist_.size() << " Input directories: " << endl;
279 for(size_t k=0; k< dirlist_.size(); k++)
280 cout << k+1 << " : " << dirlist_[k] << endl;
281 cout << " IMin= " << imin_ << " IMax= " << imax_ << " IStep= " << istep_
282 << ((rdsamefc_)?" SameFrameCounter read mode":" AllOKPaquets read mode ") << endl;
283 if (fgtimeselect_) {
284 cout << " Processing time window, StartTime=" << proctimestart_ << " duration= " << proctimeduration_
285 << " EndTime=" << proctimeend_ << endl;
286 }
287 cout << " OutPath= " << outpath_ << (fgfitsout_?" force FITS output":" PPF output") << endl;
288 cout << " Action=" << action_ << " NMean=" << nmean_ << " NBloc=" << nbloc_ << endl;
289 cout << " FreqMin= " << freqmin_ << " FreqMax= " << freqmax_ << " NBinFreq= " << nbinfreq_ << endl;
290 if (fgdtpaq_) cout<< " Fill Per Paquet DadaTable ";
291 if (fgdtms_) cout<< " Fill Time Averaged/Filtered binned power DataTable ";
292 if (!fgdtpaq_&&!fgdtms_) cout << " NO DataTable " << endl;
293 else cout << endl;
294 cout << " GainFileName=" << gainfile_ << " Variance: Min= " << vmin_ << " Max= " << vmax_
295 << " Bands: N=" << nbands_ << " First=" << bandfirst_ << " Last=" << bandlast_
296 << endl;
297 if (fgfreqfilter_) {
298 cout << " force frequence median filtering (action gain), half width =" << medhalfwidth_
299 << endl;
300 } else {
301 cout<<" NO freq. filtering" << endl;
302 }
303 cout << " Spectra TimeWindow (Nb.Paquets) Size=" << spec_win_sz_ << " ExtensionSize=" << spw_ext_sz_
304 << " MaxNbFile=" << nbmax_specwfiles_ << endl;
305 cout << " PaqSize=" << paqsize_ << " - NZones=" << nzones_ << " NPaqZone=" << npaqinzone_
306 << " PrtLevel=" << prtlevel_ << " PrtCntModulo=" << prtmodulo_ << endl;
307 cout << " AcqDataMode: " << ((fgdatafft_)?" Data_FFT " : " Data_Raw " )
308 << ((fgsinglechannel_)?" SingleChannel " : " TwoChannels " ) << endl;
309 cout << " NbVisibCalculator in Group: " << nbcalgrp_ << " with N//threads: " << nthreads_ << endl;
310
311
312 return os;
313}
314
315/* --Fonction-- */
316int BRAnaParam::DecodeMiniFitsHeader(const char* filename, uint_4& paqsz, uint_4& npaq,
317 bool& fgdatafft, bool& fgsinglechannel, SOPHYA::TimeStamp& tmstart)
318{
319 cout << " DecodeMiniFitsHeader - Opening file: " << filename << endl;
320 MiniFITSFile mff(filename, MF_Read);
321 string acqmode=mff.GetKeyValue("ACQMODE");
322 cout << "DecodeMiniFitsHeader()... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
323 << " NAxis2=" << mff.NAxis2() << " AcqMode=" << acqmode << endl;
324 paqsz = mff.NAxis1();
325 npaq = mff.NAxis2();
326 if (acqmode.substr(0,3)=="fft") fgdatafft=true;
327 if (acqmode.find("1c") < acqmode.length()) fgsinglechannel=true;
328 string fkvs=mff.GetKeyValue("DATEOBS");
329 if (fkvs.length()>0) tmstart.Set(fkvs);
330 fkvs=mff.GetKeyValue("TMSTART");
331 if (fkvs.length()>0) tmstart.Set(fkvs);
332 return 0;
333}
334
Note: See TracBrowser for help on using the repository browser.