source: Sophya/trunk/FrEROS/AnaLC/anafsv.cc@ 3308

Last change on this file since 3308 was 3308, checked in by ansari, 18 years ago

Creation du module AnaLC (lecture suivi EROS avec SOPHYA) dans la base
SOPHYA - cmv+reza 22/08/2007

  • Property svn:executable set to *
File size: 20.1 KB
Line 
1
2#include <stdlib.h>
3#include <stdio.h>
4#include <ctype.h>
5#include <math.h>
6#include <iostream>
7#include <typeinfo>
8
9#include "machdefs.h"
10
11#include "fmath.h"
12
13#include "fsvcache.h"
14#include "timing.h"
15
16#include "fsvst.h"
17
18#include "pexceptions.h" // SOPHYA Exceptions
19
20#include "histats.h" // SOPHYA HiStats module include files
21#include "sopnamsp.h" // using SOPHYA namespace
22
23#include "machdefs.h"
24
25#include "tvector.h"
26#include "srandgen.h"
27#include "fioarr.h"
28#include "sopemtx.h"
29#include "pexceptions.h"
30#include "matharr.h"
31
32#include "sambainit.h"
33
34#include "lcproc.h" // declaration of LightCurveProc class
35
36
37DataTable dtclum(int numet, int nmes, TIMEINFO *tminf, MESUREU *mesu);
38DataTable dttimeinfo(int nmes, TIMEINFO *tminf);
39DataTable createdtmes(int nmes);
40void dtmes(DataTable& dtmes, int nmes, STARINFO *sti, MESUREU *mesu);
41
42
43struct resolution
44 {
45 int nbin;
46 float xmin, xmax, dxbin;
47 int *nst, *nstr;
48 float *fmoy, *fmoyr;
49 float *dflx, *sflx, *dflxe, *sflxe;
50 };
51
52typedef struct resolution SRESOL;
53
54
55SRESOL *resolinit(int nmes, int nbin, float xmin, float xmax);
56int resolst(SRESOL *sres, int nmes, STARINFO *sti, MESUREU *mesu);
57int resolend(SRESOL *sres, int nmes);
58
59
60
61
62/* ------------------------------------------------------------------- */
63
64
65
66char *Pgnm = "hanafsv";
67
68int main(int narg, char *arg[] )
69{
70char strg[128],*svf;
71SUIVIFIP *sfip;
72int numes1,numes2;
73int numet1, numet2, incn;
74float xmin, xmax;
75float ymin, ymax;
76
77int i,k,jt,rc,dnp,fgq;
78int nbstar,nbmes;
79int nstar, nmes;
80
81int prtl,prfg,clfg,ntfg;
82
83int nbin;
84float minlog,maxlog;
85
86
87GLOBINFO glinf;
88STARINFO star;
89TIMEINFO *tminf, *tmi;
90
91TIMEINFOU *tminu;
92MESURE *mes;
93MESUREU *mesu;
94
95double *TStart;
96
97SRESOL *sres;
98
99/* ................................................... */
100
101InitTim();
102
103if (narg > 1)
104 if (strcmp(arg[1],"-h") == 0)
105 {
106 printf("Usage: %s [-Flag] NomSuivi\n",Pgnm);
107 puts (" -inq : Marque une pause apres lecture TimeInfo ");
108 puts (" -act prfg,clfg,ntfg : Flags de controle ");
109 puts (" -prt prtl : Niveau d'impression (0->3)");
110 puts (" -his nbin,MinLog,MaxLog : Binning histo resolution");
111 puts (" -mes m1-m2 : Traitement des mesures m1 a m2 ");
112 puts (" -numet n1-n2,incn : Traitement des etoiles n1 a n2 (C1)");
113 puts (" -xpos min-max : Coupure min <= XPosEtoile <= max (C2)");
114 puts (" -ypos min-max : Coupure min <= YPosEtoile <= max (C3)");
115 puts (" On fait C1 & C2 & C3, Valeur par defaut = Tout");
116 puts (" prfg (PrintFlag) : 0 , 1= Glob+TimeInfo");
117 puts (" 2= +StarInfo, 3= +MesuresU (=default), 4= +Mesures");
118 puts (" clfg (Flag Courbe de lumiere) : 0(=defaut),1");
119 puts (" ntfg (Bits Flag Ntuple): B0(=defaut) TimeInfo, B1 Resol, B2 Mesures\n");
120 return(0);
121 }
122
123/* Decodage des arguments */
124numes1 = 1; numes2 = 9999999;
125numet1 = 1; numet2 = 9999999; incn = 1;
126xmin = -9.e19; xmax = 9.e19;
127ymin = -9.e19; ymax = 9.e19;
128prtl = 0; prfg = 3;
129clfg = 0; ntfg = 1;
130nbin = 30;
131minlog = 2.0; maxlog = 5.0;
132fgq = 0; /* On continue jusqu'au bout par defaut */
133
134if (narg < 2) goto PbArg;
135k = 1;
136while (*(arg[k]) == '-')
137 {
138 if ( strcmp(arg[k],"-inq") == 0) fgq = 1;
139 else if ( strcmp(arg[k],"-mes") == 0)
140 { if (k == narg-1) goto PbArg;
141 sscanf(arg[k+1],"%d-%d",&numes1,&numes2); k++; }
142 else if ( strcmp(arg[k],"-numet") == 0)
143 { if (k == narg-1) goto PbArg;
144 sscanf(arg[k+1],"%d-%d,%d",&numet1,&numet2,&incn); k++; }
145 else if ( strcmp(arg[k],"-xpos") == 0)
146 { if (k == narg-1) goto PbArg;
147 sscanf(arg[k+1],"%g-%g",&xmin,&xmax); k++; }
148 else if ( strcmp(arg[k],"-ypos") == 0)
149 { if (k == narg-1) goto PbArg;
150 sscanf(arg[k+1],"%g-%g",&ymin,&ymax); k++; }
151 else if ( strcmp(arg[k],"-prt") == 0)
152 { if (k == narg-1) goto PbArg;
153 sscanf(arg[k+1],"%d",&prtl); k++; }
154 else if ( strcmp(arg[k],"-act") == 0)
155 { if (k == narg-1) goto PbArg;
156 sscanf(arg[k+1],"%d,%d,%d",&prfg,&clfg,&ntfg); k++; }
157 else if ( strcmp(arg[k],"-his") == 0)
158 { if (k == narg-1) goto PbArg;
159 sscanf(arg[k+1],"%d,%g,%g",&nbin,&minlog,&maxlog); k++; }
160
161 if (++k >= narg) break;
162 }
163
164if (k < narg) goto Suite;
165PbArg:
166printf("%s_Erreur : Pb. arguments ( %s -h pour aide) \n", Pgnm, Pgnm);
167return(1);
168
169Suite:
170svf = arg[k];
171
172cout << " hanafsv: Args decoded, starting suivi processing ..." << endl;
173SophyaInit();
174sfip = NULL; rc = 0;
175try {
176
177 POutPersist po("LCInfo.ppf");
178 POutPersist po2("TimInfo.ppf");
179 POutPersist po3("Mesu.ppf");
180
181 if ( (sfip = SuiviOpen(svf,SUOF_RO_MEM2)) == NULL ) return(5);
182
183 nbstar = SuiviGetNbStars(sfip);
184 nbmes = SuiviGetNbMesures(sfip);
185
186 if (numes1 < 1) numes1 = 1;
187 if (numes2 < 1) numes2 = 1;
188 if (numes2 > nbmes) numes2 = nbmes;
189 if (numes2 < numes1) numes1 = numes2;
190 nmes = numes2-numes1+1;
191 if (numet1 < 1) numet1 = 1;
192 if (numet2 < 1) numet2 = 1;
193 if (numet2 > nbstar) numet2 = nbstar;
194 if (numet2 < numet1) numet1 = numet2;
195 nstar = numet2 - numet1 +1;
196
197 //--> new tminf = (TIMEINFO *)malloc(nmes*sizeof(TIMEINFO));
198 tminf = new TIMEINFO[nmes];
199 //--> new tminu = (TIMEINFOU *)malloc(nmes*sizeof(TIMEINFOU));
200 tminu = new TIMEINFOU[nmes];
201 // mes = (MESURE *)malloc(nmes*sizeof(MESURE));
202 mes = new MESURE[nmes];
203 //--> new mesu = (MESUREU *)malloc(nmes*sizeof(MESUREU));
204 mesu = new MESUREU[nmes];
205 //-> new TStart = (double *) malloc(nmes*sizeof(double));
206 TStart = new double[nmes];
207 if ( (tminf == NULL) || (tminu == NULL) || (mes == NULL) ||
208 (mesu == NULL) || (TStart == NULL))
209 { puts("hanafsv_Erreur : Pb allocation"); return(7); }
210
211 DataTable dtm = createdtmes(nmes);
212
213 LightCurveProc mylcp("lcproc.ppf"); // mylcp("lcproc.ppf", 512);
214
215 printf("\n ====> hanafsv : Lecture du fichier de suivi %s \n", svf);
216 printf(" %d <= NumMes <= %d %d <= NumEt <= %d (mod %d)\n",
217 numes1, numes2, numet1, numet2, incn);
218 printf(" %g <= XPosRef <= %g ET %g <= YPosRef <= %g \n\n",
219 xmin,xmax,ymin,ymax);
220 printf("PrtFg= %d (Level=%d) ClFg= %d NtFg= %d Nbin,Min,MaxLog= %d %g %g \n",
221 prfg, prtl, clfg, ntfg, nbin, minlog, maxlog);
222
223
224
225 sres = NULL;
226 sres = resolinit(nmes, nbin, minlog, maxlog);
227 if (sres == NULL) { puts("hanafsv_Erreur : Pb resolinit !"); goto Fin; }
228
229
230 /* Lecture GlobInfo */
231 SuiviReadGlobInfo(sfip, (char *)(&glinf));
232 if (prfg > 0) PrtGlobInfo(&glinf, prtl);
233
234 /* Boucle lecture mesure */
235 for(jt=0; jt<nmes; jt++)
236 {
237 tmi = tminf+jt;
238 rc = SuiviReadTimeInfo(sfip,jt+numes1,jt+numes1,(char *)tmi);
239 if (rc != 0)
240 { printf("hanafsv_Erreur: Pb SuiviReadTimeInfo( , %d, ...) \n", jt+numes1);
241 return(10); }
242 DecodeTim(tminf+jt, tminu+jt);
243 TStart[jt] = (double)tmi->TStart / (double)86400;
244 if (prfg > 0) PrtTimeInfo(tmi, jt+numes1, prtl);
245 if ((tminu+jt)->FgCalib < 1)
246 {
247 (mesu+jt)->Xi2 = -1.;
248 (mesu+jt)->Flux = (mesu+jt)->FluxB = 0.;
249 (mesu+jt)->ErrFlux = (mesu+jt)->ErrFluxB = 99999.0;
250 }
251 } // fin de boucle for(jt) - boucle de lecture des TimeInfo
252
253 mylcp.TimeInfo(nmes, tminf);
254
255 if (ntfg&1) {
256 DataTable TimInf = dttimeinfo(nmes, tminf);
257 cout << " Writing TimInf to po2 (TimInfo.ppf) ... " << endl;
258 po2 << PPFNameTag("dttimeinfo") << TimInf ;
259 }
260
261
262 if (fgq) {
263 puts("\n\n =====> <CR> pour continuer, <q><CR> pour stop ...");
264 gets(strg);
265 if (toupper(strg[0]) == 'Q') return(0);
266 }
267
268 dnp = numet2-numet1;
269 if (dnp > 20) dnp /= 10;
270 else dnp = 1;
271
272 for (i=numet1; i<=numet2; i+=incn ) /* Loop over stars in suivi file */
273 {
274
275 if (((i-numet1)%dnp) == 0)
276 { sprintf(strg,"LectureStar[%d] ",i); PrtTim(strg); }
277
278 rc = SuiviReadStarInfo(sfip,i,(char *)(&star));
279
280 if (rc != 0) return(10);
281
282 if ( (star.XPos < xmin) || (star.XPos > xmax) ||
283 (star.YPos < ymin) || (star.YPos > ymax) ) continue;
284
285 if (prfg > 1)
286 {
287 printf(" --------- StarInfo/Mesures Etoile [%d] ----------\n",i);
288 PrtStarInfo(&star, i, prtl);
289 }
290 if (star.NumEt < 0) continue;
291 if (star.FgRef == 0) continue;
292
293
294 for (jt=0; jt<nmes; jt++)
295 {
296 rc = SuiviReadMesures(sfip,i,jt+numes1,jt+numes1,
297 (char *)(mes+jt));
298 if (rc != 0)
299 {
300 printf("Erreur: Star[%d] Pb lecture Mesure[%d] (Rc=%d) \n",
301 i+1,jt+1,rc);
302 return(15);
303 }
304 if ((tminu+jt)->FgCalib > 0)
305 {
306 DecodeMes(mes+jt, mesu+jt);
307 Calibre_F_E(mesu+jt, tminu+jt);
308 if (prfg > 3) PrtMesure(mes+jt, jt+numes1, prtl);
309 if (prfg > 2) PrtMesureU(mesu+jt, jt+numes1, prtl);
310 }
311 else if (prfg > 2)
312 PrtMesure(mes+jt, jt+numes1, prtl);
313
314 } /* Fin de boucle sur les mesures */
315
316
317 // a star with all its measured fluxes is ready here to be processed
318 // We process our light-curve
319 mylcp.ProcessLC(i, nmes, &star, tminf, mesu);
320
321 if (clfg > 0){
322 DataTable LiCu = dtclum(i, nmes, tminf, mesu);
323 char tagbuf[32];
324 sprintf(tagbuf,"cl%d",i);
325 po << PPFNameTag(tagbuf) << LiCu ;
326 }
327
328 if (ntfg&2)
329 resolst(sres, nmes, &star, mesu);
330
331 if (ntfg&4) {
332 dtmes(dtm, nmes, &star, mesu);
333 }
334
335
336 } /* Fin de boucle sur les etoiles */
337
338Fin:
339
340 cout << " Writing dtmes to file po3 Mesu.ppf , dtm: \n" << dtm << endl;
341 po3 << PPFNameTag("dtmes") << dtm ;
342 SuiviClose(sfip);
343
344 cout << "---- End Of Loop over stars, LightCurveProc.getDTS() : \n" << mylcp.getDTS() << endl;
345
346 if ((ntfg&2) && (sres != NULL))
347 resolend(sres,nmes);
348 rc = 0;
349}
350catch (PThrowable & pex) {
351 cerr << " hanafsv/Error - Exception catched " << (string)typeid(pex).name()
352 << " - Msg= " << pex.Msg() << endl;
353 rc = 98;
354}
355catch ( ... )
356 {
357 cerr << " hanafsv/Error - Exception ... catched " << endl;
358 rc = 99;
359}
360
361
362
363return(rc);
364}
365
366
367
368
369///////////////////////// DataTabale dtclum ///////////////////////////
370DataTable dtclum(int numet, int nmes, TIMEINFO *tminf, MESUREU *mesu)
371{
372int idc,jt;
373
374 DataTable dt(100);
375 dt.AddFloatColumn("TStart");
376 dt.AddFloatColumn("Flux");
377 dt.AddFloatColumn("Xi2");
378 dt.AddFloatColumn("Fond");
379 dt.AddFloatColumn("ErrFlux");
380 dt.AddFloatColumn("SigX");
381 dt.AddFloatColumn("SugY");
382 dt.AddFloatColumn("Rho");
383 dt.AddFloatColumn("DelX");
384 dt.AddFloatColumn("DelY");
385 dt.AddFloatColumn("X");
386 dt.AddFloatColumn("Y");
387
388 double xnt[15];
389
390idc = numet;
391
392for(jt=0; jt<nmes; jt++)
393 {
394 xnt[0] = tminf[jt].TStart/86400.;
395 xnt[1] = (double) mesu[jt].Flux;
396 xnt[2] = (double) mesu[jt].Xi2;
397 xnt[3] = (double) mesu[jt].Fond;
398 xnt[4] = (double) mesu[jt].ErrFlux;
399 xnt[5] = (double) tminf[jt].SigX;
400 xnt[6] = (double) tminf[jt].SigY;
401 xnt[7] = (double) tminf[jt].Rho;
402 xnt[8] = (double) tminf[jt].DelX;
403 xnt[9] = (double) tminf[jt].DelY;
404 xnt[10] = (double) mesu[jt].X ;
405 xnt[11] = (double) mesu[jt].Y ;
406 dt.AddRow(xnt);
407 }
408
409return(dt);
410}
411
412
413
414//////////////////// DataTable dttimeinfo //////////////////////////
415DataTable dttimeinfo(int nmes, TIMEINFO *tminf)
416{
417int idn,jt,i;
418
419 DataTable dt(10000);
420 dt.AddFloatColumn("NumPhoto");
421 dt.AddFloatColumn("TStart");
422 dt.AddFloatColumn("Expose");
423 dt.AddFloatColumn("Nb0kGF");
424 dt.AddFloatColumn("Nb0kPS");
425 dt.AddFloatColumn("Fond");
426 dt.AddFloatColumn("SigFond");
427 dt.AddFloatColumn("SigX");
428 dt.AddFloatColumn("SigY");
429 dt.AddFloatColumn("Rho");
430 dt.AddFloatColumn("SgX");
431 dt.AddFloatColumn("SgY");
432 dt.AddFloatColumn("DelX");
433 dt.AddFloatColumn("DelY");
434 dt.AddFloatColumn("Absorption");
435 dt.AddFloatColumn("Airmass");
436 dt.AddFloatColumn("FgCalib");
437
438 for(i=0; i<6; i++){
439 char calname[32];
440 sprintf(calname,"Calib%d",i);
441 dt.AddFloatColumn(calname);
442 }
443
444 for(i=0; i<6; i++)
445 for(int j =0; j<3;j++){
446 char name[32];
447 sprintf(name,"Resol%d%d",j,i);
448 dt.AddFloatColumn(name);
449 }
450
451
452 double xnt[48];
453
454 idn = 300;
455 for(i=0; i<36; i++) xnt[i] = 0.;
456
457for(jt=0; jt<nmes; jt++)
458 {
459 xnt[0] = tminf[jt].NumPhoto;
460 xnt[1] = tminf[jt].TStart/86400.;
461 xnt[2] = tminf[jt].Expose;
462 xnt[3] = tminf[jt].NbOkGF;
463 xnt[4] = tminf[jt].NbOkPS;
464 xnt[5] = tminf[jt].Fond;
465 xnt[6] = tminf[jt].SigFond;
466
467 xnt[7] = tminf[jt].SigX;
468 xnt[8] = tminf[jt].SigY;
469 xnt[9] = tminf[jt].Rho;
470 xnt[10] = tminf[jt].SgX;
471 xnt[11] = tminf[jt].SgY;
472 xnt[12] = tminf[jt].DelX;
473 xnt[13] = tminf[jt].DelY;
474 xnt[14] = tminf[jt].Absorption;
475 xnt[15] = tminf[jt].AirMass;
476
477 xnt[16] = tminf[jt].FgCalib;
478 for(i=0; i<6; i++)
479 {
480 xnt[17+i] = tminf[jt].Calib[i];
481 xnt[23+i*3] = tminf[jt].Resol[0][i];
482 xnt[24+i*3] = tminf[jt].Resol[1][i];
483 xnt[25+i*3] = tminf[jt].Resol[2][i];
484 }
485 dt.AddRow(xnt);
486 }
487
488return(dt);
489}
490
491
492
493////////////////// DataTable creatdtmes ////////////////////////
494DataTable createdtmes(int nmes)
495
496{
497
498 DataTable dt(100);
499 dt.AddFloatColumn("NumEt");
500 dt.AddFloatColumn("FluxRef");
501 dt.AddFloatColumn("XPos");
502 dt.AddFloatColumn("YPos");
503 dt.AddFloatColumn("FgRef");
504 dt.AddFloatColumn("NbVois");
505 dt.AddFloatColumn("DisMin");
506 dt.AddFloatColumn("DisM2");
507
508 if (nmes > 10) nmes = 10;
509
510 char nflux[32]; char nerrflux[32]; char nfluxb[32]; char nerrfluxb[32];
511 char nfond[32]; char nxi2[32]; char nx[32]; char ny[32]; char nsigx[32];
512 char nsigy[32]; char nflx[32]; char nfnd[32]; char nxi[32];
513
514 for (int k=0; k<nmes; k++) {
515 sprintf(nflux,"Flux%d",k); sprintf(nerrflux,"ErrFlux%d",k);
516 sprintf(nfluxb,"FluxB%d",k); sprintf(nerrfluxb,"ErrFluxB%d",k);
517 sprintf(nfond,"Fond%d",k); sprintf(nxi2,"Xi2%d",k); sprintf(nx,"X%d",k);
518 sprintf(ny,"Y%d",k); sprintf(nsigx,"SigX%d",k); sprintf(nsigy,"SigY%d",k);
519 sprintf(nflx,"Flx%d",k); sprintf(nfnd,"Fnd%d",k); sprintf(nxi,"Xi%d",k);
520
521 dt.AddFloatColumn(nflux);
522 dt.AddFloatColumn(nerrflux);
523 dt.AddFloatColumn(nfluxb);
524 dt.AddFloatColumn(nerrfluxb);
525 dt.AddFloatColumn(nfond);
526 dt.AddFloatColumn(nxi2);
527 dt.AddFloatColumn(nx);
528 dt.AddFloatColumn(ny);
529 dt.AddFloatColumn(nsigx);
530 dt.AddFloatColumn(nsigy);
531 dt.AddFloatColumn(nflx);
532 dt.AddFloatColumn(nfnd);
533 dt.AddFloatColumn(nxi);
534 }
535 return dt;
536}
537
538void dtmes(DataTable& dt, int nmes, STARINFO *sti, MESUREU *mesu)
539{
540int idn,jt,off,sz;
541
542
543double xnt[200];
544
545
546//idn = 400;
547if (nmes > 10) nmes=10;
548
549xnt[0] = sti->NumEt;
550xnt[1] = sti->FluxRef;
551xnt[2] = sti->XPos;
552xnt[3] = sti->YPos;
553xnt[4] = sti->FgRef;
554xnt[5] = sti->NbVois;
555xnt[6] = sti->DisMin;
556xnt[7] = sti->DisM2;
557
558off = 7; sz = 13;
559if (nmes > 10) nmes = 10;
560for(jt=0; jt<nmes; jt++)
561 {
562 xnt[off+sz*jt+1] = mesu[jt].Flux;
563 xnt[off+sz*jt+2] = mesu[jt].ErrFlux;
564 xnt[off+sz*jt+3] = mesu[jt].FluxB;
565 xnt[off+sz*jt+4] = mesu[jt].ErrFluxB;
566 xnt[off+sz*jt+5] = mesu[jt].Fond;
567 xnt[off+sz*jt+6] = mesu[jt].Xi2;
568 xnt[off+sz*jt+7] = mesu[jt].X;
569 xnt[off+sz*jt+8] = mesu[jt].Y;
570 xnt[off+sz*jt+9] = mesu[jt].SigX;
571 xnt[off+sz*jt+10] = mesu[jt].SigY;
572 xnt[off+sz*jt+11] = mesu[jt].Flx0;
573 xnt[off+sz*jt+12] = mesu[jt].Fnd0;
574 xnt[off+sz*jt+13] = mesu[jt].Xi20;
575 }
576 dt.AddRow(xnt);
577
578return;
579}
580
581
582
583/* Nouvelle-Fonction */
584
585SRESOL *resolinit(int nmes, int nbin, float xmin, float xmax)
586
587{
588int i,j;
589SRESOL *sres;
590
591if (nbin < 1) return(NULL);
592if (xmax <= xmin) return(NULL);
593
594//-> new sres = malloc(nmes*sizeof(SRESOL));
595sres = new SRESOL[nmes];
596if (sres == NULL) return(NULL);
597for(i=0; i<nmes; i++)
598 {
599 sres[i].nbin = nbin;
600 sres[i].xmin = xmin;
601 sres[i].xmax = xmax;
602 sres[i].dxbin = (xmax-xmin)/(float)nbin;
603 //-> new sres[i].nst = malloc(nbin*sizeof(int));
604 sres[i].nst = new int[nbin];
605 if (sres[i].nst == NULL) break;
606 //-> new sres[i].nstr = malloc((nbin+2)*sizeof(int));
607 sres[i].nstr = new int[nbin+2];
608 if (sres[i].nstr == NULL)
609 { delete[] sres[i].nst; break; }
610 //->new sres[i].fmoy = malloc(nbin*sizeof(float));
611 sres[i].fmoy = new float[nbin];
612 if (sres[i].fmoy == NULL)
613 { delete[] sres[i].nst; delete[] sres[i].nstr; break; }
614 //->new sres[i].fmoyr = malloc(nbin*sizeof(float));
615 sres[i].fmoyr = new float[nbin] ;
616 if (sres[i].fmoyr == NULL)
617 { delete[] sres[i].nst; delete[] sres[i].nstr;
618 delete[] sres[i].fmoy; break; }
619 //-> new sres[i].dflx = malloc(nbin*sizeof(float));
620 sres[i].dflx = new float[nbin];
621 if (sres[i].dflx == NULL)
622 { delete[] sres[i].nst; delete[] sres[i].nstr;
623 delete[] sres[i].fmoy; delete[] sres[i].fmoyr; break; }
624 //->new sres[i].sflx = malloc(nbin*sizeof(float));
625 sres[i].sflx = new float[nbin];
626 if (sres[i].sflx == NULL)
627 { delete[] sres[i].nst; delete[] sres[i].nstr;
628 delete[] sres[i].fmoy; delete[] sres[i].fmoyr;
629 delete[] sres[i].dflx; break; }
630 //->new sres[i].dflxe = malloc(nbin*sizeof(float));
631 sres[i].dflxe = new float[nbin];
632 if (sres[i].dflxe == NULL)
633 { delete[] sres[i].nst; delete[] sres[i].nstr;
634 delete[] sres[i].fmoy; delete[] sres[i].fmoyr;
635 delete[] sres[i].dflx; delete[] sres[i].sflx; break; }
636 //->new sres[i].sflxe = malloc(nbin*sizeof(float));
637 sres[i].sflxe = new float[nbin];
638 if (sres[i].sflxe == NULL)
639 { delete[] sres[i].nst; delete[] sres[i].nstr;
640 delete[] sres[i].fmoy; delete[] sres[i].fmoyr;
641 delete[] sres[i].dflx; delete[] sres[i].sflx;
642 delete[] sres[i].dflxe; break; }
643 }
644
645if (i != nmes)
646 {
647/* BUG REZA ??? puts("resolinit_Erreur : Pb malloc/new %d %d \n", i, nmes); */
648 printf("resolinit_Erreur : Pb malloc %d %d \n", i, nmes);
649 for(j=0; j<i; j++)
650 {
651 delete[] sres[j].nst;
652 delete[] sres[j].nstr;
653 delete[] sres[j].fmoy;
654 delete[] sres[j].fmoyr;
655 delete[] sres[j].dflx;
656 delete[] sres[j].sflx;
657 delete[] sres[j].dflxe;
658 delete[] sres[j].sflxe;
659 }
660 delete[] sres;
661 return(NULL);
662 }
663for(i=0; i<nmes; i++)
664 {
665 for(j=0; j<nbin; j++)
666 {
667 sres[i].nst[j] = sres[i].nstr[j] = 0 ;
668 sres[i].fmoy[j] = sres[i].fmoyr[j] = 0. ;
669 sres[i].dflx[j] = 0. ;
670 sres[i].sflx[j] = 0. ;
671 sres[i].dflxe[j] = 0. ;
672 sres[i].sflxe[j] = 0. ;
673 }
674 }
675
676return(sres);
677}
678
679
680/* Nouvelle-Fonction */
681
682int resolst(SRESOL *sres, int nmes, STARINFO *sti, MESUREU *mesu)
683{
684int i,j,ib;
685float err, flx, flxr;
686float xlg, dx;
687
688flxr = sti->FluxRef;
689if (flxr < 1.) return(-1);
690xlg = (float)log10((double)flxr);
691
692for(j=0; j<nmes; j++)
693 {
694 ib = (int)((xlg-sres[j].xmin)/sres[j].dxbin);
695 if (ib < 0) sres[j].nstr[sres[j].nbin]++;
696 if (ib >= sres[j].nbin) sres[j].nstr[sres[j].nbin+1]++;
697 if ( (ib < 0) || (ib >= sres[j].nbin) ) continue;
698 sres[j].nstr[ib]++;
699 sres[j].fmoyr[ib] += flxr;
700 if (mesu[j].Xi2 < 0.) continue;
701 err = mesu[j].ErrFlux;
702 flx = mesu[j].Flux;
703 if (err/flxr < 0.01) err = 0.01*flxr;
704 /* if (fabs((double)(err/flx)) > (1./xlg)) continue; */
705 if (fabs((double)((flx-flxr)/err)) > 8.) continue;
706 sres[j].nst[ib]++;
707 sres[j].fmoy[ib] += flx;
708 dx = flxr-flx;
709 sres[j].dflx[ib] += dx;
710 sres[j].sflx[ib] += (dx*dx);
711 dx /= err;
712 sres[j].dflxe[ib] += dx;
713 sres[j].sflxe[ib] += (dx*dx);
714 }
715
716return(ib);
717}
718
719/* Nouvelle-Fonction */
720int resolend(SRESOL *sres, int nmes)
721
722{
723int i,j,ib,idn;
724float ns;
725float xnt[10];
726
727printf("Resolend_Info: Under/OverFlow %d %d \n",
728 sres[0].nstr[sres[0].nbin], sres[0].nstr[sres[0].nbin+1]);
729for(j=0; j<nmes; j++)
730 {
731 idn = 501+j;
732 for(ib=0; ib<sres[j].nbin; ib++)
733 {
734 xnt[0] = ib;
735 xnt[1] = sres[j].xmin+ib*sres[j].dxbin;
736 for(i=2; i<10; i++) xnt[i] = 0.;
737 if (sres[j].nstr[ib] < 1) goto Fill;
738 sres[j].fmoyr[ib] /= (float)sres[j].nstr[ib];
739 xnt[8] = sres[j].nstr[ib];
740 xnt[9] = sres[j].fmoyr[ib];
741 if (sres[j].nst[ib] < 1) goto Fill;
742 xnt[2] = ns = sres[j].nst[ib];
743 sres[j].fmoy[ib] /= ns;
744 xnt[3] = sres[j].fmoy[ib];
745 sres[j].dflx[ib] /= ns;
746 sres[j].sflx[ib] /= ns;
747 sres[j].sflx[ib] -= sres[j].dflx[ib]*sres[j].dflx[ib];
748 if (sres[j].sflx[ib] >= 0.)
749 sres[j].sflx[ib] = (float)sqrt((double)sres[j].sflx[ib]);
750 xnt[4] = sres[j].dflx[ib];
751 xnt[5] = sres[j].sflx[ib];
752 sres[j].dflxe[ib] /= ns;
753 sres[j].sflxe[ib] /= ns;
754 sres[j].sflxe[ib] -= sres[j].dflxe[ib]*sres[j].dflxe[ib];
755 if (sres[j].sflxe[ib] >= 0.)
756 sres[j].sflxe[ib] = (float)sqrt((double)sres[j].sflxe[ib]);
757 xnt[6] = sres[j].dflxe[ib];
758 xnt[7] = sres[j].sflxe[ib];
759 Fill:
760 // --- SOPHYA NTuple fill hfn_(&idn, xnt);
761 printf("Resol[%d,%d] ns,nsr %d %d fm,df,sf= %g (%g) %g %g \n",j,ib,
762 sres[j].nst[ib], sres[j].nstr[ib],
763 sres[j].fmoy[ib], sres[j].fmoyr[ib],
764 sres[j].dflx[ib], sres[j].sflx[ib]);
765 }
766 }
767
768for(j=0; j<nmes; j++)
769 {
770 delete[] sres[j].nst;
771 delete[] sres[j].nstr;
772 delete[] sres[j].fmoy;
773 delete[] sres[j].fmoyr;
774 delete[] sres[j].dflx;
775 delete[] sres[j].sflx;
776 delete[] sres[j].dflxe;
777 delete[] sres[j].sflxe;
778 }
779delete[] sres;
780
781return (0);
782
783}
Note: See TracBrowser for help on using the repository browser.