source: Sophya/trunk/SophyaLib/HiStats/ntuple.cc@ 975

Last change on this file since 975 was 846, checked in by ansari, 26 years ago

ajout SetVal pour I/O fichiers fits
File size: 12.6 KB
Line 
1#include <stdio.h>
2#include <string.h>
3
4#include "strutil.h"
5#include "perrors.h"
6#include "ntuple.h"
7
8
9#define BADVAL -1.e19
10#define LENNAME 8
11#define LENNAME1 (LENNAME+1)
12
13//++
14// Class NTuple
15// Lib Outils++
16// include ntuple.h
17//
18// Classe de ntuples
19//--
20//++
21// Links Parents
22// PPersist
23// NTupleInterface
24//--
25
26/* --Methode-- */
27//++
28NTuple::NTuple()
29//
30// Createur par defaut
31//--
32{
33mNVar = mNEnt = mBlk = mNBlk = 0;
34mVar = NULL;
35mVarD = NULL;
36mNames = NULL;
37mInfo = NULL;
38}
39
40
41//++
42NTuple::NTuple(int nvar, char** noms, int blk)
43//
44// Createur d'un ntuple de `nvar' variables dont les
45// noms sont dans le tableau de cahines de caracteres `noms'
46// avec `blk' d'evenements par blocks.
47//--
48{
49mNVar = mNEnt = mBlk = mNBlk = 0;
50mVar = NULL;
51mVarD = NULL;
52mNames = NULL;
53mInfo = NULL;
54if (nvar <= 0) THROW(sizeMismatchErr);
55mNVar = nvar;
56mVar = new r_4[nvar];
57mVarD = new r_8[nvar];
58if (blk < 10) blk = 10;
59mBlk = blk;
60// On prend des noms de LENNAME char pour le moment
61mNames = new char[nvar*LENNAME1];
62r_4* pt = new r_4[nvar*blk];
63mNBlk = 1;
64mPtr.push_back(pt);
65int i;
66for(i=0; i<nvar; i++)
67 { strncpy(mNames+i*LENNAME1, noms[i], LENNAME);
68 mNames[i*LENNAME1+LENNAME] = '\0'; }
69return;
70}
71
72// cmv 8/10/99
73//++
74NTuple::NTuple(const NTuple& NT)
75//
76// Createur par copie (clone).
77//--
78: mNVar(0), mNEnt(0), mBlk(0), mNBlk(0)
79, mVar(NULL), mVarD(NULL), mNames(NULL), mInfo(NULL)
80{
81if(NT.mNVar<=0) return; // cas ou NT est cree par defaut
82mNVar = NT.mNVar;
83mBlk = NT.mBlk;
84mVar = new r_4[NT.mNVar];
85mVarD = new r_8[NT.mNVar];
86mNames = new char[NT.mNVar*LENNAME1];
87
88int i;
89r_4* pt = new r_4[mNVar*mBlk];
90mNBlk = 1; mPtr.push_back(pt);
91
92for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i));
93
94if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}
95
96if(NT.mNEnt<=0) return;
97for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}
98
99return;
100}
101
102/* --Methode-- */
103//++
104NTuple::NTuple(char *flnm)
105//
106// Createur lecture fichier ppersist.
107//--
108{
109mNVar = mNEnt = mBlk = mNBlk = 0;
110mVar = NULL;
111mVarD = NULL;
112mNames = NULL;
113mInfo = NULL;
114PInPersist s(flnm);
115ObjFileIO<NTuple> fiont(this);
116fiont.Read(s);
117}
118
119/* --Methode-- */
120NTuple::~NTuple()
121{
122Clean();
123}
124
125/* --Methode-- */
126void NTuple::Clean()
127{
128if (mVar) delete[] mVar;
129if (mVarD) delete[] mVarD;
130if (mNames) delete[] mNames;
131if (mInfo) delete mInfo;
132int i;
133if(mNBlk>0) for(i=0; i<mNBlk; i++) delete[] mPtr[i];
134mPtr.erase(mPtr.begin(), mPtr.end());
135mNVar = mNEnt = mBlk = mNBlk = 0;
136mVar = NULL;
137mVarD = NULL;
138mNames = NULL;
139mInfo = NULL;
140return;
141}
142
143/* --Methode-- cmv 08/10/99 */
144//++
145NTuple& NTuple::operator = (const NTuple& NT)
146//
147// Operateur egal (clone).
148//--
149{
150if(this == &NT) return *this;
151Clean();
152if(NT.mNVar<=0) return *this; // cas ou NT est cree par defaut
153mNVar = NT.mNVar;
154mBlk = NT.mBlk;
155mVar = new r_4[NT.mNVar];
156mVarD = new r_8[NT.mNVar];
157mNames = new char[NT.mNVar*LENNAME1];
158
159int i;
160r_4* pt = new r_4[mNVar*mBlk];
161mNBlk = 1; mPtr.push_back(pt);
162
163for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i));
164
165if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}
166
167if(NT.mNEnt<=0) return *this;
168for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}
169
170// En fait il faudrait un createur par copie qui partage les donnees
171// quand l'objet est temporaire... trop complique A FAIRE ! cmv.
172return *this;
173}
174
175/* --Methode-- */
176//++
177void NTuple::Fill(r_4* x)
178//
179// Remplit le ntuple avec le tableau cd reels `x'.
180//--
181{
182int numb = mNEnt/mBlk;
183if (numb >= mNBlk) {
184 r_4* pt = new r_4[mNVar*mBlk];
185 mNBlk++;
186 mPtr.push_back(pt);
187}
188int offb = mNEnt-numb*mBlk;
189memcpy((mPtr[numb]+offb*mNVar), x, mNVar*sizeof(r_4));
190mNEnt++;
191return;
192}
193
194
195/* --Methode-- */
196//++
197float NTuple::GetVal(int n, int k) const
198//
199// Retourne la valeur de la variable `k' de l'evenement `n'.
200//--
201{
202if (n >= mNEnt) return(BADVAL);
203if ( (k < 0) || (k >= mNVar) ) return(BADVAL);
204int numb = n/mBlk;
205int offb = n-numb*mBlk;
206return(*(mPtr[numb]+offb*mNVar+k));
207}
208
209/* --Methode-- */
210//++
211void NTuple::SetVal(int n, int k, float value)
212//
213// initialise la valeur de la variable `k' de l'evenement `n'.
214//--
215{
216if (n >= mNEnt) return;
217if ( (k < 0) || (k >= mNVar) ) return;
218int numb = n/mBlk;
219int offb = n-numb*mBlk;
220*(mPtr[numb]+offb*mNVar+k) = value;
221}
222
223/* --Methode-- */
224//++
225int NTuple::IndexNom(const char* nom) const
226//
227// Retourne le numero de la variable de nom `nom'.
228//--
229{
230int i;
231for(i=0; i<mNVar; i++)
232 if ( strcmp(nom, mNames+i*LENNAME1) == 0) return(i);
233return(-1);
234}
235
236
237static char nomretour[2*LENNAME1];
238/* --Methode-- */
239//++
240char* NTuple::NomIndex(int k) const
241//
242// Retourne le nom de la variable numero 'k'
243//--
244{
245nomretour[0] = '\0';
246if ((k >= 0) && (k < mNVar)) strcpy(nomretour, mNames+k*LENNAME1);
247return(nomretour);
248}
249
250
251/* --Methode-- */
252//++
253r_4* NTuple::GetVec(int n, r_4* ret) const
254//
255// Retourne l'evenement `n' dans le vecteur `ret'.
256//--
257{
258int i;
259if (ret == NULL) ret = mVar;
260if (n >= mNEnt) {
261 for(i=0; i<mNVar; i++) ret[i] = BADVAL;
262 return(ret);
263}
264
265int numb = n/mBlk;
266int offb = n-numb*mBlk;
267memcpy(ret, (mPtr[numb]+offb*mNVar), mNVar*sizeof(r_4));
268return(ret);
269}
270
271/* --Methode-- */
272//++
273r_8* NTuple::GetVecD(int n, r_8* ret) const
274//
275// Retourne l'evenement `n' dans le vecteur `ret'.
276//--
277{
278int i;
279if (ret == NULL) ret = mVarD;
280float * fr = GetVec(n);
281for(i=0; i<mNVar; i++) ret[i] = fr[i];
282return(ret);
283}
284
285
286
287/* --Methode-- */
288//++
289DVList& NTuple::Info()
290//
291// Renvoie une référence sur l'objet DVList Associé
292//--
293{
294if (mInfo == NULL) mInfo = new DVList;
295return(*mInfo);
296}
297
298/* --Methode-- */
299//++
300void NTuple::Print(int num, int nmax) const
301//
302// Imprime `nmax' evenements a partir du numero `num'.
303//--
304{
305int i,j;
306
307printf("Num ");
308for(i=0; i<mNVar; i++) printf("%8s ", mNames+i*LENNAME1);
309putchar('\n');
310
311if (nmax <= 0) nmax = 1;
312if (num < 0) num = 0;
313nmax += num;
314if (nmax > mNEnt) nmax = mNEnt;
315for(i=num; i<nmax; i++) {
316 GetVec(i, NULL);
317 printf("%6d ", i);
318 for(j=0; j<mNVar; j++) printf("%8g ", (float)mVar[j]);
319 putchar('\n');
320}
321return;
322}
323
324/* --Methode-- */
325//++
326void NTuple::Show(ostream& os) const
327//
328// Imprime l'information generale sur le ntuple.
329//--
330{
331os << "NTuple: NVar= " << mNVar << " NEnt=" << mNEnt
332 << " (Blk Sz,Nb= " << mBlk << " ," << mNBlk << ")\n";
333os << " Variables Min Max \n";
334int i;
335double min, max;
336char buff[128];
337for(i=0; i<mNVar; i++) {
338 GetMinMax(i, min, max);
339 sprintf(buff, "%3d %16s %10lg %10lg \n", i, mNames+i*LENNAME1, min, max);
340 os << (string)buff ;
341 }
342os << endl;
343}
344
345
346/* --Methode-- */
347//++
348int NTuple::FillFromASCIIFile(string const& fn, float defval)
349//
350// Remplit le ntuple a partir d'un fichier ASCII.
351// Renvoie le nombre de lignes ajoutees.
352//--
353{
354if (NbColumns() < 1) {
355 cout << "NTuple::FillFromASCIIFile() Ntuple has " << NbColumns() << " columns" << endl;
356 return(-1);
357 }
358FILE * fip = fopen(fn.c_str(), "r");
359if (fip == NULL) {
360 cout << "NTuple::FillFromASCIIFile() Error opening file " << fn << endl;
361 return(-2);
362 }
363
364char lineb[4096];
365char *line;
366char* ccp;
367int j,kk;
368int postab, posb;
369float* xv = new float[NbColumns()];
370
371int nlread = 0;
372int nvar = NbColumns();
373// On boucle sur toutes les lignes
374while (fgets(lineb,4096,fip) != NULL) {
375 lineb[4095] = '\0';
376 j = 0; line = lineb;
377// On enleve les espaces et tab de debut
378 while ( (line[j] != '\0') && ((line[j] == ' ') || (line[j] == '\t')) ) j++;
379 line = lineb+j;
380// Il faut que le premier caractere non-espace soit un digit, ou + ou - ou .
381 if (!( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') )) continue;
382 ccp = line;
383 for(kk=0; kk<nvar; kk++) xv[kk] = defval;
384 for(kk=0; kk<nvar; kk++) {
385// Les mots sont separes par des espaces ou des tab
386 postab = posc(ccp, '\t' );
387 posb = posc(ccp, ' ' );
388 if (postab >= 0) {
389 if (posb < 0) posb = postab;
390 else if (postab < posb) posb = postab;
391 }
392 if (posb >= 0) ccp[posb] = '\0';
393 if ( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') )
394 xv[kk] = atof(ccp);
395 if (posb < 0) break;
396 ccp += posb+1; j = 0;
397 while ( (ccp[j] != '\0') && ((ccp[j] == ' ') || (ccp[j] == '\t')) ) j++;
398 ccp += j;
399 }
400 Fill(xv);
401 nlread++;
402 }
403
404delete[] xv;
405cout << "NTuple::FillFromASCIIFile( " << fn << ") " << nlread << " fill from file " << endl;
406return(nlread);
407}
408
409
410// ------- Implementation de l interface NTuple ---------
411
412/* --Methode-- */
413uint_4 NTuple::NbLines() const
414{
415return(NEntry());
416}
417/* --Methode-- */
418uint_4 NTuple::NbColumns() const
419{
420return(NVar());
421}
422
423/* --Methode-- */
424r_8 * NTuple::GetLineD(int n) const
425{
426return(GetVecD(n));
427}
428
429/* --Methode-- */
430r_8 NTuple::GetCell(int n, int k) const
431{
432return(GetVal(n, k));
433}
434
435/* --Methode-- */
436r_8 NTuple::GetCell(int n, string const & nom) const
437{
438return(GetVal(n, nom.c_str()));
439}
440
441/* --Methode-- */
442//++
443void NTuple::GetMinMax(int k, double& min, double& max) const
444//
445// Retourne le minimum et le maximum de la variable `k'.
446//--
447{
448min = 9.e19; max = -9.e19;
449if ( (k < 0) || (k >= mNVar) ) return;
450int jb,ib,i;
451double x;
452i=0;
453for(jb=0; jb< mNBlk; jb++)
454 for(ib=0; ib< mBlk; ib++) {
455 if (i >= mNEnt) break;
456 i++;
457 x = *(mPtr[jb]+ib*mNVar+k);
458 if(i==1) {min = x; max = x;}
459 if (x < min) min = x;
460 if (x > max) max = x;
461 }
462return;
463}
464
465/* --Methode-- */
466void NTuple::GetMinMax(string const & nom, double& min, double& max) const
467{
468GetMinMax(IndexNom(nom.c_str()), min, max);
469}
470
471/* --Methode-- */
472int NTuple::ColumnIndex(string const & nom) const
473{
474return(IndexNom(nom.c_str()));
475}
476
477/* --Methode-- */
478string NTuple::ColumnName(int k) const
479{
480return(NomIndex(k));
481}
482
483/* --Methode-- */
484//++
485string NTuple::VarList_C(const char* nomx) const
486//
487// Retourne une chaine de caracteres avec la declaration des noms de
488// variables. si "nomx!=NULL" , des instructions d'affectation
489// a partir d'un tableau "nomx[i]" sont ajoutees.
490//--
491{
492string rets="";
493int i;
494for(i=0; i<mNVar; i++) {
495 if ( (i%5 == 0) && (i > 0) ) rets += ";";
496 if (i%5 == 0) rets += "\ndouble ";
497 else rets += ",";
498 rets += mNames+i*LENNAME1;
499 }
500rets += "; \n";
501if (nomx) {
502 char buff[256];
503 for(i=0; i<mNVar; i++) {
504 sprintf(buff,"%s=%s[%d]; ", mNames+i*LENNAME1, nomx, i);
505 rets += buff;
506 if ( (i%3 == 0) && (i > 0) ) rets += "\n";
507 }
508 }
509
510return(rets);
511}
512
513
514/* --Methode-- */
515//++
516string NTuple::LineHeaderToString() const
517// Retourne une chaine de caracteres avec la liste des noms de
518// variables, utilisables pour une impression
519//--
520{
521char buff[32];
522string rets=" Num ";
523for(int i=0; i<mNVar; i++) {
524 sprintf(buff, "%8s ", mNames+i*LENNAME1);
525 rets += buff;
526 }
527rets += '\n';
528return(rets);
529}
530
531/* --Methode-- */
532//++
533string NTuple::LineToString(int n) const
534// Retourne une chaine de caracteres avec le contenu de la ligne "n"
535// utilisable pour une impression
536//--
537{
538char buff[32];
539double* val;
540val = GetLineD(n);
541sprintf(buff,"%6d: ",n);
542string rets=buff;
543int i;
544for(i=0; i<mNVar; i++) {
545 sprintf(buff, "%8.3g ", val[i]);
546 rets += buff;
547 }
548rets += '\n';
549return(rets);
550}
551
552
553/* --Methode-- */
554//++
555void ObjFileIO<NTuple>::WriteSelf(POutPersist& s) const
556//
557// Ecriture ppersist du ntuple.
558//--
559{
560char strg[256];
561if (dobj->mInfo) sprintf(strg, "NVar=%6d NEnt=%9d BlkSz=%6d NBlk=%6d HasInfo",
562 (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
563else sprintf(strg, "NVar=%6d NEnt=%9d BlkSz=%6d NBlk=%6d ",
564 (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
565s.PutLine(strg);
566s.PutI4(dobj->mNVar);
567s.PutBytes(dobj->mNames, dobj->mNVar*LENNAME1);
568s.PutI4(dobj->mNEnt);
569s.PutI4(dobj->mBlk);
570s.PutI4(dobj->mNBlk);
571if (dobj->mInfo) s << (*(dobj->mInfo));
572int jb;
573for(jb=0; jb<dobj->mNBlk; jb++)
574 s.PutR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk);
575return;
576}
577
578/* --Methode-- */
579//++
580void ObjFileIO<NTuple>::ReadSelf(PInPersist& s)
581//
582// Lecture ppersist du ntuple.
583//--
584{
585
586dobj->Clean();
587
588char strg[256];
589s.GetLine(strg, 255);
590// Pour savoir s'il y avait un DVList Info associe
591bool hadinfo = false;
592if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo = true;
593
594s.GetI4(dobj->mNVar);
595dobj->mNames = new char[dobj->mNVar*LENNAME1];
596dobj->mVar = new r_4[dobj->mNVar];
597dobj->mVarD = new r_8[dobj->mNVar];
598s.GetBytes(dobj->mNames, dobj->mNVar*LENNAME1);
599s.GetI4(dobj->mNEnt);
600s.GetI4(dobj->mBlk);
601s.GetI4(dobj->mNBlk);
602
603if (hadinfo) { // Lecture eventuelle du DVList Info
604 if (dobj->mInfo == NULL) dobj->mInfo = new DVList;
605 s >> (*(dobj->mInfo));
606 }
607
608int jb;
609for(jb=0; jb<dobj->mNBlk; jb++) {
610 r_4* pt = new r_4[dobj->mNVar*dobj->mBlk];
611 dobj->mPtr.push_back(pt);
612 s.GetR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk);
613}
614
615}
616
617
618#ifdef __CXX_PRAGMA_TEMPLATES__
619#pragma define_template ObjFileIO<NTuple>
620#endif
621
622#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
623template class ObjFileIO<NTuple>;
624#endif
Note: See TracBrowser for help on using the repository browser.