source: Sophya/trunk/ArchTOIPipe/TestPipes/simofftst.cc@ 2650

Last change on this file since 2650 was 2014, checked in by ansari, 23 years ago

Amelioration de calcul d'offset (possibilite de coupure sur coordonnees
galactique)
Nouveau processeur nettoyeur (SimpleCleaner)
Ajout programme de traitement aksj02.cc , donnees Archeops KS1/KS3

Reza 28/5/2002
File size: 11.7 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5// $Id: simofftst.cc,v 1.3 2002-05-28 17:30:45 ansari Exp $
6
7/* Test de processeurs ds simtoipr.cc - Reza Avril 2001
8
9---------------- Exemple d'appel ---------------------
10csh> simofftst -start 104385384 -end 104399964 -range -500,500 \
11 -intoi boloMuV_26 -wtoi 8192 -wdegli 1024 \
12 inputbolo.fits filt.fits xx.ppf
13# Avec Filtre
14csh> simofftst -start 104385384 -end 104399964 -range -500,500 \
15 -intoi boloMuV_26 -wtoi 8192 -wdegli 1024 \
16 inputbolo.fits filtfft.fits spectre.ppf
17csh> $exedir/simofftst -start $sn0 -end $sn1 -intoi $cleanname \
18 -wtoi 8192 -wdegli 1024 -degli 3.,2.,3,2,0 \
19 -soff 128,9,3 -gfilt 16,2.5 -prstat \
20 -soffnt $outppf -bgalcut -20.,20. -bgal $glat \
21 -bgalfile $point $cleanfile $outfile
22
23*/
24
25
26
27#include "machdefs.h"
28#include <math.h>
29#include <unistd.h>
30
31#include "toimanager.h"
32#include "cgt.h"
33#include "fitstoirdr.h"
34#include "fitstoiwtr.h"
35#include "simtoipr.h"
36#include "simoffset.h"
37#include "nooppr.h"
38#include "timing.h"
39#include "histinit.h"
40#include "psighand.h"
41#include <stdexcept>
42
43void Usage(bool fgerr)
44{
45 cout << endl;
46 if (fgerr) {
47 cout << " simofftst : Argument Error ! simofftst -h for usage " << endl;
48 exit(1);
49 }
50 else {
51 cout << "\n Usage : simofftst [-intoi toiname] [-start snb] [-end sne] \n"
52 << " [-dbg] [-wtoi sz] [-wdegli sz] [-range min,max] \n"
53 << " [-degli ns,ns2,mxpt,mnpt,wrec] [-soff wsz,nptfit,degpol] \n"
54 << " [-soffnt PPFName] [-gfilt wsz,sigma] \n"
55 << " [-bgalcut bmin,bmax] [-bgal toiname] [-bgalfile pointFitsName] \n"
56 << " [-nooutflg] [-nowrtms] [-nowrticd] [-prstat] [-useseqbuff] \n"
57 << " inFitsName outFitsName \n"
58 << " -dbg : sets TOISeqBuffered debug level to 1 \n"
59 << " -start snb : sets the start sample num \n"
60 << " -end sne : sets the end sample num \n"
61 << " -range min,max : sets the acceptable range for intoi \n"
62 << " default= -16000,16000\n"
63 << " -intoi toiName : select input TOI name (def bolo)\n"
64 << " -wtoi sz : sets TOISeqBuff buffer size (def= 8192)\n"
65 << " -wdegli sz : sets deglitcher window size (def= 512) \n"
66 << " -degli ns,ns2,maxnpt,minnpt,wrec : sets deglitcher parameters\n"
67 << " -soff wsz,nptfit,degpol : set SimpleOffsetEstimator parameters\n"
68 << " -soffnt PPFName : Writes out SimpleOffsetEstimator NTuple \n"
69 << " -gfilt wsz,sigma : Activate gaussian filter (par= window size, sigma) \n"
70 // << " -wfft sz : Activate Fourier filter and sets its width \n"
71 // << " -keepfft n : Keeps n spectra (if Fourier filter activated) \n"
72 // << " -killfreq bf,nharm,sigf : kills nharm frequency, basefreq=bf \n"
73 // << " with a gaussian with width=sigf \n"
74 << " -nooutflg : Don't write flags to output FITS \n"
75 << " -wrtms : Write mean/sigma TOI's from Deglitcher \n"
76 << " -wrticd : Write incopie,degli TOI's from deglitcher \n"
77 << " -prstat : PrintStat with ProcSampleCounter \n"
78 << " -useseqbuff : Use TOISeqBuffered TOI's (default: TOISegmented) \n"
79 << endl;
80 }
81 if (fgerr) exit(1);
82}
83
84int main(int narg, char** arg) {
85
86 if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
87
88 cout << "simofftst starting - Decoding arguments " << " narg=" << narg << endl;
89
90 bool fgdbg = false;
91 bool fgsetstart = false;
92
93 bool fgprstat = false;
94 bool fgsegmented = true;
95
96 int dbglev = 0;
97
98 int wtoi = 8192;
99 int wdegli = 512;
100
101 int wfft = 4096;
102 int keepfft = 0;
103 int nmax = 10;
104 int istart = 0;
105 int iend = 0;
106
107 double range_min = -16000;
108 double range_max = 16000.;
109
110 // Deglitcher parameters
111 double degli_ns1 = 3.;
112 double degli_ns2 = 1.5;
113 int degli_minnpt = 2;
114 int degli_maxnpt = 4;
115 int degli_wrec = 10;
116
117 // SimpleOffsetEstimator parameters
118 int soff_mwsz = 256;
119 int soff_nptfit = 5;
120 int soff_degpol = 2;
121 bool soff_nt = false;
122 string soff_nt_ppfname = "";
123
124 // Gaussian filter parameters
125 bool gfilt_fg = false;
126 int gfilt_wsz = 16;
127 double gfilt_sigma = 2.;
128
129 // Coordinate file for galactic cut
130 string pointfile = "pointgal.fits";
131 double bmin = -20;
132 double bmax = 20.;
133 string bgaltoi = "bgal";
134 bool bgalcut = false;
135
136 // File names
137 string infile;
138 string outfile;
139 string outppfname;
140 string intoi = "bolo";
141
142 bool fg_wrtflag = true;
143 bool fg_wrtms = false;
144 bool fg_wrticd = false;
145
146
147 bool fg_f_filt = false;
148 bool fg_killfreq = false;
149 int bf_killfreq = 1;
150 int nharm_killfreq = 1;
151 double sigf_killfreq = 1.;
152
153 if (narg < 4) Usage(true);
154 int ko=1;
155 // decoding arguments
156 for(int ia=1; ia<narg; ia++) {
157 if (strcmp(arg[ia],"-start") == 0) {
158 if (ia == narg-1) Usage(true); // -start est suivi d'un argument
159 istart = atoi(arg[ia+1]); ia++;
160 fgsetstart = true;
161 }
162 else if (strcmp(arg[ia],"-end") == 0) {
163 if (ia == narg-1) Usage(true);
164 iend = atoi(arg[ia+1]); ia++;
165 }
166 else if (strcmp(arg[ia],"-wtoi") == 0) {
167 if (ia == narg-1) Usage(true);
168 wtoi = atoi(arg[ia+1]); ia++;
169 }
170 else if (strcmp(arg[ia],"-wdegli") == 0) {
171 if (ia == narg-1) Usage(true);
172 wdegli = atoi(arg[ia+1]); ia++;
173 }
174 else if (strcmp(arg[ia],"-degli") == 0) {
175 if (ia == narg-1) Usage(true);
176 sscanf(arg[ia+1],"%lf,%lf,%d,%d,%d", &degli_ns1, &degli_ns2,
177 &degli_maxnpt, &degli_minnpt, &degli_wrec);
178 ia++;
179 }
180 else if (strcmp(arg[ia],"-soff") == 0) {
181 if (ia == narg-1) Usage(true);
182 sscanf(arg[ia+1],"%d,%d,%d", &soff_mwsz, &soff_nptfit, &soff_degpol);
183 ia++;
184 }
185 else if (strcmp(arg[ia],"-soffnt") == 0) {
186 soff_nt = true;
187 soff_nt_ppfname = arg[ia+1]; ia++;
188 }
189 else if (strcmp(arg[ia],"-gfilt") == 0) {
190 if (ia == narg-1) Usage(true);
191 sscanf(arg[ia+1],"%d,%lf", &gfilt_wsz, &gfilt_sigma);
192 gfilt_fg = true;
193 ia++;
194 }
195 else if (strcmp(arg[ia],"-bgalcut") == 0) {
196 if (ia == narg-1) Usage(true);
197 sscanf(arg[ia+1],"%lf,%lf", &bmin, &bmax);
198 bgalcut = true;
199 ia++;
200 }
201 else if (strcmp(arg[ia],"-bgalfile") == 0) {
202 if (ia == narg-1) Usage(true);
203 pointfile = arg[ia+1];
204 ia++;
205 }
206 else if (strcmp(arg[ia],"-bgal") == 0) {
207 if (ia == narg-1) Usage(true);
208 bgaltoi = arg[ia+1];
209 ia++;
210 }
211 else if (strcmp(arg[ia],"-wfft") == 0) {
212 if (ia == narg-1) Usage(true);
213 wfft = atoi(arg[ia+1]); ia++;
214 fg_f_filt = true;
215 }
216 else if (strcmp(arg[ia],"-keepfft") == 0) {
217 if (ia == narg-1) Usage(true);
218 keepfft = atoi(arg[ia+1]); ia++;
219 }
220 else if (strcmp(arg[ia],"-killfreq") == 0) {
221 if (ia == narg-1) Usage(true);
222 sscanf(arg[ia+1],"%d,%d,%lf",&bf_killfreq, &nharm_killfreq, &sigf_killfreq);
223 fg_killfreq = true;
224 ia++;
225 }
226 else if (strcmp(arg[ia],"-range") == 0) {
227 if (ia == narg-1) Usage(true);
228 sscanf(arg[ia+1],"%lf,%lf",&range_min, &range_max);
229 ia++;
230 }
231 else if (strcmp(arg[ia],"-intoi") == 0) {
232 if (ia == narg-1) Usage(true);
233 intoi = arg[ia+1]; ia++;
234 }
235 else if (strcmp(arg[ia],"-dbg") == 0) fgdbg = true;
236 else if (strcmp(arg[ia],"-nooutflg") == 0) fg_wrtflag = false;
237 else if (strcmp(arg[ia],"-wrtms") == 0) fg_wrtms = true;
238 else if (strcmp(arg[ia],"-wrticd") == 0) fg_wrticd = true;
239
240 else if (strcmp(arg[ia],"-prstat") == 0) fgprstat = true;
241 else if (strcmp(arg[ia],"-useseqbuff") == 0) fgsegmented = false;
242
243 else { ko = ia; break; } // Debut des noms
244 }
245
246 if (iend < istart) iend = istart+wtoi*(nmax+5);
247 if ((narg-ko) < 2) Usage(true);
248 infile = arg[ko];
249 outfile = arg[ko+1];
250 // outppfname = arg[ko+2];
251
252 cout << " Initializing SOPHYA ... " << endl;
253 SophyaInit();
254 SophyaConfigureSignalhandling(true);
255
256 InitTim();
257
258 cout << ">>>> simofftst: Infile= " << infile << " outFile="
259 << outfile << endl;
260 cout << ">>>> Window Size WTOI= " << wtoi << " WDEGLI= " << wdegli
261 << " iStart= " << istart << " iEnd= " << iend << endl;
262 cout << ">>>> InTOIName= " << intoi << endl;
263
264 try {
265 TOIManager* mgr = TOIManager::getManager();
266
267 // mgr->setRequestedSample(11680920,11710584);
268 // mgr->setRequestedSample(104121000, 104946120);
269 if (fgsetstart)
270 mgr->setRequestedSample(istart, iend);
271
272 cout << "> Creating FITSTOIReader object - InFile=" << infile << endl;
273 FITSTOIReader r(infile);
274
275 cout << "> Creating SimpleDeglitcher() " << endl;
276 SimpleDeglitcher degl(wdegli);
277 degl.SetDetectionParam(degli_ns1, degli_ns2, degli_maxnpt,
278 degli_minnpt, degli_wrec);
279 cout << " Setting Range for deglitcher: " << range_min
280 << " - " << range_max << endl;
281 degl.SetRange(range_min, range_max);
282
283 cout << "> Creating SimpleDeglitcher(" << soff_mwsz << "," << soff_nptfit
284 << "," << soff_degpol << ")" << endl;
285 SimpleOffsetEstimator offes(soff_mwsz, soff_nptfit, soff_degpol);
286 offes.SavePolyNTuple(soff_nt, soff_nt_ppfname);
287
288 if (gfilt_fg)
289 cout << "> Creating a GaussianFilter SimpleFilter Object " << endl;
290 double G_sigma = gfilt_sigma;
291 double G_a = 1./(G_sigma*sqrt(M_PI*2.));
292 SimpleFilter filt(gfilt_wsz, SimpleFilter::GaussFilter, G_a, G_sigma);
293
294
295 FITSTOIReader* rbgal = NULL;
296 if (bgalcut) { // if Galactic cut
297 cout << "> Creating bgal FITSTOIReader object - InFile=" << pointfile
298 << " bgaltoiname= " << bgaltoi << endl;
299 cout << " offes.SetBGalCut( " << bmin << "," << bmax << ")" << endl;
300 rbgal = new FITSTOIReader(pointfile);
301 offes.SetBGalCut(bmin, bmax);
302 }
303
304 cout << "> Creating FITSTOIWriter OutFitsName= " << outfile << endl;
305 FITSTOIWriter w(outfile);
306
307 CGT plombier(fgsegmented, wtoi);
308 plombier.SetDebugLevel(dbglev);
309
310 cout << "> Connecting Processors through plombier ... " << endl;
311 string inname = "in";
312 plombier.Connect(r, intoi, degl, inname);
313 plombier.Connect(degl, "out", offes, "in");
314 if (gfilt_fg) {
315 plombier.Connect(offes, "out", filt, "in");
316 plombier.Connect(filt, "out", w, "degoffil", "", 0, fg_wrtflag);
317 }
318 else
319 plombier.Connect(offes, "out", w, "degoff", "", 0, fg_wrtflag);
320
321 if (bgalcut) {
322 inname = "bgal";
323 plombier.Connect(*rbgal, bgaltoi, offes, inname);
324 }
325
326 if (fg_wrticd) {
327 plombier.Connect(degl, "incopie", w, "incopie");
328 plombier.Connect(offes, "incopie", w, "indegl");
329 }
330 if (fg_wrtms) {
331 plombier.Connect(degl, "mean", w, "meandegl");
332 plombier.Connect(degl, "sigma", w, "sigdegl");
333 }
334
335 cout << "> Plombier status before start" << endl;
336 cout << plombier ;
337
338 cout << offes;
339
340 PrtTim("starting processors");
341 plombier.Start();
342
343
344 // ------------------- Impression continu de stat ------------------------
345 if (fgprstat) {
346 ProcSampleCounter<SimpleOffsetEstimator> stats(offes);
347 stats.InfoMessage() = "simofftst/Info";
348 stats.PrintStats();
349 }
350 // -----------------------------------------------------------------------
351
352 cout << degl;
353 cout << offes;
354 if (gfilt_fg) cout << filt;
355 cout << w;
356
357 mgr->joinAll();
358 PrtTim("End threads");
359 if (bgalcut) delete rbgal;
360
361 }
362 catch (PThrowable & exc) {
363 cerr << "\n simofftst: Catched Exception \n" << (string)typeid(exc).name()
364 << " - Msg= " << exc.Msg() << endl;
365 }
366 catch (const std::exception & sex) {
367 cerr << "\n simofftst: Catched std::exception \n"
368 << (string)typeid(sex).name() << endl;
369 }
370 catch (...) {
371 cerr << "\n simofftst: some other exception was caught ! " << endl;
372 }
373
374 return(0);
375}
Note: See TracBrowser for help on using the repository browser.