source: Sophya/trunk/SophyaPI/PIext/servnobjm.cc@ 3435

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

Ajout fonctions tetphi2mollX/Y() et longlat2mollX/Y() ds les fichiers d'expression ntuple generes par piapp - Reza 11/12/2007

File size: 44.6 KB
RevLine 
[165]1#include <stdio.h>
2#include <stdlib.h>
3#include <ctype.h>
4
[295]5#include <typeinfo>
[2322]6#include <iostream>
[165]7#include <string>
8#include <list>
9#include <map>
10
[2615]11#include "sopnamsp.h"
[165]12#include "strutil.h"
13
14#include "nobjmgr.h"
15#include "servnobjm.h"
[330]16#include "nomgadapter.h"
[165]17#include "pistdimgapp.h"
18
[333]19#include "fct1dfit.h"
20#include "fct2dfit.h"
21
[544]22#ifdef SANS_EVOLPLANCK
[333]23#include "matrix.h"
24#include "cvector.h"
[544]25#else
26#include "tmatrix.h"
27#include "tvector.h"
[2922]28#include "matharr.h"
[584]29#include "pitvmaad.h"
[544]30#endif
31
[333]32#include "ntuple.h"
33#include "cimage.h"
34
[165]35#include "histos.h"
36#include "histos2.h"
[2792]37#include "hisprof.h"
[165]38#include "ntuple.h"
[2792]39#include "datatable.h"
[165]40
[1905]41#include "piyfxdrw.h"
[326]42#include "pisurfdr.h"
[333]43
44#include "pintuple.h"
45#include "pintup3d.h"
46
[326]47#include "pipodrw.h"
[165]48
49
[326]50
[165]51/* --Methode-- */
[2491]52Services2NObjMgr::Services2NObjMgr(string& tmpdir)
[165]53{
[1276]54SetTmpDir(tmpdir);
[333]55mImgapp = NULL;
[2491]56mOmg = NULL;
[171]57dynlink = NULL;
[165]58}
59
60/* --Methode-- */
61Services2NObjMgr::~Services2NObjMgr()
62{
63CloseDLL();
[2491]64if (mOmg) delete mOmg;
[165]65}
66
67/* --Methode-- */
[295]68void Services2NObjMgr::RegisterClass(AnyDataObj* o, NObjMgrAdapter* oa)
[165]69{
[295]70ObjAdaptList::iterator it;
71for(it = objadaplist.begin(); it != objadaplist.end(); it++)
[495]72#ifdef SANS_EVOLPLANCK
[295]73 if (typeid(*o) == typeid(*((*it).obj))) THROW(dupIdErr);
[495]74#else
75 if (typeid(*o) == typeid(*((*it).obj)))
76 throw(DuplicateIdExc("Services2NObjMgr::RegisterClass() - Duplicate class"));
77#endif
[295]78dataobj_adapter oba;
79oba.obj = o;
80oba.obja = oa;
81objadaplist.push_back(oba);
82}
[165]83
[295]84/* --Methode-- */
85NObjMgrAdapter* Services2NObjMgr::GetAdapter(AnyDataObj* o)
86{
87ObjAdaptList::iterator it;
88for(it = objadaplist.begin(); it != objadaplist.end(); it++)
89 if (typeid(*o) == typeid(*((*it).obj))) return((*it).obja->Clone(o));
90return(new NObjMgrAdapter(o));
91}
[165]92
[295]93/* --Methode-- */
[1276]94void Services2NObjMgr::SetTmpDir(string const & tmpdir)
95{
96TmpDir = tmpdir;
97PDynLinkMgr::SetTmpDir(tmpdir);
98return;
99}
100
101/* --Methode-- */
[333]102void Services2NObjMgr::PlotFunc(string const & expfunc, string & nom, double xmin, double xmax, int np, string dopt)
[295]103{
[165]104FILE *fip;
[326]105string fname = TmpDir + "func1_pia_dl.c";
106string cmd;
107int rc;
108
109if (!mImgapp) return;
110
[2755]111// Pour synchronisation d'execution simultanee
[2762]112ZSync zs(mutx_dynlink); zs.NOp();
[2755]113
[326]114cmd = "rm -f " + fname;
115rc = system(cmd.c_str());
116// printf("PlotFunc_Do> %s (Rc=%d)\n", cmd.c_str(), rc);
117
118if ((fip = fopen(fname.c_str(), "w")) == NULL) {
119 string sn = fname;
[449]120 cout << "Services2NObjMgr/PlotFunc_Error: fopen( " << sn << endl;
[326]121 return;
122 }
123
124// constitution du fichier a compiler
125fputs("#include <math.h> \n", fip);
126fputs("double func1_pia_dl_func(double x) \n{\n", fip);
127fprintf(fip,"return(%s); \n}\n", expfunc.c_str());
128fclose(fip);
129
130string func = "func1_pia_dl_func";
131DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func);
132if (!f) return;
[333]133PlotFunc(f, nom, xmin, xmax, np, dopt);
[326]134CloseDLL();
135return;
136}
137
138/* --Methode-- */
[333]139void Services2NObjMgr::PlotFunc2D(string const & expfunc, string & nom, double xmin, double xmax,
140 double ymin, double ymax, int npx, int npy, string dopt)
[326]141{
142FILE *fip;
143string fname = TmpDir + "func2_pia_dl.c";
144string cmd;
145int rc;
146
147if (!mImgapp) return;
[2755]148// Pour synchronisation d'execution simultanee
[2762]149ZSync zs(mutx_dynlink); zs.NOp();
[326]150
151cmd = "rm " + fname;
152rc = system(cmd.c_str());
153// printf("PlotFunc2D_Do> %s (Rc=%d)\n", cmd.c_str(), rc);
154
155if ((fip = fopen(fname.c_str(), "w")) == NULL) {
156 string sn = fname;
[449]157 cout << "Services2NObjMgr/PlotFunc2D_Error: fopen( " << sn << endl;
[326]158 return;
159 }
160
161// constitution du fichier a compiler
162fputs("#include <math.h> \n", fip);
163fputs("double func2_pia_dl_func(double x, double y) \n{\n", fip);
164fprintf(fip,"return(%s); \n}\n", expfunc.c_str());
165fclose(fip);
166
167string func = "func2_pia_dl_func";
168DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func);
169if (!f) return;
[333]170PlotFunc2D(f, nom, xmin, xmax, ymin, ymax, npx, npy, dopt);
[326]171CloseDLL();
172return;
173}
174
175/* --Methode-- */
[333]176void Services2NObjMgr::PlotFuncFrCFile(string const & fname, string const & func, string & nom,
177 double xmin, double xmax, int np, string dopt)
[326]178{
[2755]179// Pour synchronisation d'execution simultanee
[2762]180ZSync zs(mutx_dynlink); zs.NOp();
[326]181DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func);
182if (!f) return;
[333]183PlotFunc(f, nom, xmin, xmax, np, dopt);
[326]184CloseDLL();
185return;
186}
187
188/* --Methode-- */
[333]189void Services2NObjMgr::PlotFunc2DFrCFile(string const & fname, string const & func, string & nom,
190 double xmin, double xmax, double ymin, double ymax, int npx, int npy, string dopt)
[326]191{
[2755]192// Pour synchronisation d'execution simultanee
[2762]193ZSync zs(mutx_dynlink); zs.NOp();
[326]194DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func);
195if (!f) return;
[333]196PlotFunc2D(f, nom, xmin, xmax, ymin, ymax, npx, npy, dopt);
[326]197CloseDLL();
198return;
199}
200
201/* --Methode-- */
[333]202void Services2NObjMgr::PlotFunc(DlFunctionOfX f, string & nom, double xmin, double xmax, int np, string dopt)
[326]203{
204if (!mImgapp) return;
205
206int k;
207if (np < 1) np = 1;
208if (xmax <= xmin) xmax = xmin+1.;
209Vector* vpy = new Vector(np);
210
211double xx;
212double dx = (xmax-xmin)/np;
[500]213
214try {
215 for(k=0; k<np; k++) { xx = xmin+dx*k; (*vpy)(k) = f(xx); }
216}
[495]217#ifdef SANS_EVOLPLANCK
[500]218CATCH(merr) {
[326]219 fflush(stdout);
220 cout << endl;
221 cerr << endl;
222 string es = PeidaExc(merr);
223 cerr << "Services2NObjMgr::PlotFunc() Exception :" << merr << es;
224 delete vpy;
225 vpy = NULL;
226 } ENDTRY;
[495]227#else
228catch ( PException exc ) {
229 fflush(stdout);
230 cout << endl;
231 cerr << endl;
232 cerr << "Services2NObjMgr::PlotFunc() Exception :" << exc.Msg() << endl;
233 delete vpy;
234 vpy = NULL;
235}
236#endif
[326]237
238if (vpy) {
[344]239 string titre;
240 if (nom.length() < 1) {
241 titre = "Function f(x)";
[466]242 nom = "/autoc/f_x";
[344]243 }
244 else titre = nom;
[2924]245
246 MyObjMgr()->AddObj(vpy, nom);
247 if ( (dopt == "nodisp") || (dopt == "nodisplay") ) return;
248
[344]249 P1DArrayAdapter* vya = new POVectorAdapter(vpy, false);
[326]250 vya->DefineXCoordinate(xmin, (xmax-xmin)/np);
251 PIYfXDrawer* dr = new PIYfXDrawer(vya, NULL, true) ;
[1971]252 dopt = "thinline " + dopt;
253 int rsid = mImgapp->DispScDrawer(dr, titre, dopt);
[2924]254 MyObjMgr()->AddWRsId(nom, rsid);
255}
[326]256return;
257}
258
259/* --Methode-- */
[333]260void Services2NObjMgr::PlotFunc2D(DlFunctionOfXY f, string & nom, double xmin, double xmax, double ymin, double ymax,
[326]261 int npx, int npy, string dopt)
262{
263if (!mImgapp) return;
264
[2781]265if (npx < 2) npx = 2;
266if (npy < 2) npy = 2;
267//if (npx > 250) npx = 250;
268//if (npy > 250) npy = 250;
[326]269if (xmax <= xmin) xmax = xmin+1.;
270if (ymax <= ymin) ymax = ymin+1.;
271
272Matrix* mtx = new Matrix(npy, npx);
273
274int i,j;
275double xx, yy;
276double dx = (xmax-xmin)/npx;
277double dy = (ymax-ymin)/npy;
278// printf(" -- DBG -- %d %d , %g %g , %g %g \n", npx, npy, xmin, xmax, ymin, ymax);
[500]279try {
[326]280 for(j=0; j<npy; j++) {
281 yy = ymin+dy*j;
282 for(i=0; i<npx; i++) {
283 xx = xmin+dx*i;
284 (*mtx)(j, i) = f(xx, yy);
285 }
286 }
[500]287}
288#ifdef SANS_EVOLPLANCK
289CATCH(merr) {
[326]290 fflush(stdout);
291 cout << endl;
292 cerr << endl;
293 string es = PeidaExc(merr);
294 cerr << "Services2NObjMgr::PlotFunc2D() Exception :" << merr << es;
295 delete mtx; mtx = NULL;
296 } ENDTRY;
[495]297#else
298catch ( PException exc ) {
299 fflush(stdout);
300 cout << endl;
301 cerr << endl;
302 cerr << "Services2NObjMgr::PlotFunc2D() Exception :" << exc.Msg() << endl;
303 delete mtx; mtx = NULL;
304}
305#endif
[326]306
307if (mtx) {
[344]308 string titre;
309 if (nom.length() < 1) {
310 titre = "Function f(x,y)";
[466]311 nom = "/autoc/f2d_xy";
[344]312 }
313 else titre = nom;
[3375]314
315 MyObjMgr()->AddObj(mtx, nom);
316 if ( (dopt == "nodisp") || (dopt == "nodisplay") ) return;
317
[2781]318 int rsid = 0;
[344]319 P2DArrayAdapter* arr = new POMatrixAdapter(mtx, false);
[326]320 arr->DefineXYCoordinates(xmin, ymin, dx, dy);
[2781]321 if ( (npx <= 200) && (npy <= 200) ) {
322 PISurfaceDrawer* sdr = new PISurfaceDrawer(arr, true, true, true);
323 rsid = mImgapp->Disp3DDrawer(sdr, titre, dopt);
324 }
325 else rsid = mImgapp->DispImage(arr, titre, dopt);
[3375]326 MyObjMgr()->AddWRsId(nom, rsid);
327}
[2781]328
[326]329return;
330}
331
332/* --Methode-- */
[2180]333void Services2NObjMgr::ExpVal(string expval,string resultvarname)
334{
335 // Fill C-file to be executed
336 FILE *fip;
337 string func = "eval_pia_dl_func";
338 string fname = TmpDir + func; fname += ".c";
339 string cmd = "rm -f " + fname; system(cmd.c_str());
340 if((fip=fopen(fname.c_str(), "w"))==NULL) {
341 cout << "Services2NObjMgr/EvalExp_Error: fopen("<<fname<<")"<<endl;
342 return;
343 }
344 fprintf(fip,"#include <math.h>\n");
345 fprintf(fip,"double %s(double ___dummy_variable___) \n{\n",func.c_str());
346 // Add all variables already declared
[2491]347 DVList& varlist = MyObjMgr()->GetVarList();
[2180]348 DVList::ValList::const_iterator it;
349 for(it = varlist.Begin(); it != varlist.End(); it++) {
350#ifdef SANS_EVOLPLANCK
351 MuTyV mtv = (*it).second;
352 double value = (double)(mtv);
353#else
354 double value = (double)((*it).second.elval);
355#endif
356 string name_var = (*it).first;
357 fprintf(fip,"double %s = %.17f;\n",name_var.c_str(),value);
358 }
359 fprintf(fip,"return %s;\n}\n",expval.c_str());
360 fclose(fip);
361
362 // Dynamically link function
363 DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname,func);
364 if(!f) {
365 cout<<"Services2NObjMgr/EvalExp_Error: linking DlFunctionOfX"<<endl;
366 cout<<"...expval = "<<expval<<endl;
367 return;
368 }
369
370 // Evaluate function and close dynamic link
371 double result;
372 try {
373 result = f(0.);
374 CloseDLL();
375 } catch ( ... ) {
376 cout<<"Services2NObjMgr/EvalExp_Error: Arithmetic exception"<<endl;
377 CloseDLL();
378 return;
379 }
380
381 // Eventually store the result into variable or just print it
382 if(resultvarname.size()>0) {
[2491]383 if(MyObjMgr()->HasVar(resultvarname)) MyObjMgr()->DeleteVar(resultvarname);
[2180]384 char str[512];
385 if(result==0.) sprintf(str,"%f",result);
386 else {
387 double lr = log10(fabs(result));
388 if(lr>-17. && lr<17.) sprintf(str,"%17f",result);
389 else sprintf(str,"%.17e",result);
390 }
[2491]391 MyObjMgr()->SetVar(resultvarname,(string)str);
[2180]392 } else cout<<result<<" = "<<expval<<endl;
393
394}
395
396/* --Methode-- */
[333]397void Services2NObjMgr::DisplayPoints2D(string& nom, string& expx, string& expy,
398 string& experrx, string& experry,
[357]399 string& expcut, string dopt, string loop)
[333]400{
401NObjMgrAdapter* obja=NULL;
[2491]402obja = MyObjMgr()->GetObjAdapter(nom);
[333]403if (obja == NULL) {
[449]404 cout << "Services2NObjMgr::DisplayPoints2D() Error , No such object " << nom << endl;
[333]405 return;
406 }
407if (!mImgapp) return;
408
409// Creation NTuple
410char* ntn[4] = {"expx","expy","expex","expey",};
411NTuple* nt = NULL;
412bool haserr = false;
413
414if ( (experrx.length() > 0 ) && (experry.length() > 0 ) ) { haserr = true; nt = new NTuple(4, ntn); }
415else { haserr = false; experrx = experry = "0."; nt = new NTuple(2, ntn); }
416
[357]417ComputeExpressions(obja, expx, expy, experrx, experry, expcut, loop, nt, NULL, NULL);
[333]418
419if (nt->NEntry() < 1) {
420 cout << "Services2NObjMgr::DisplayPoints2D() Warning Zero points satisfy cut !" << endl;
421 delete nt;
422 return;
423 }
424
425// nt->Show();
426// nt->Print(0,10);
427PINTuple* pin = new PINTuple(nt, true);
428pin->SelectXY(ntn[0], ntn[1]);
429if ( haserr ) pin->SelectErrBar(ntn[2], ntn[3]);
430
[1971]431dopt = "defline " + dopt;
[333]432string titre = nom + ":" + expy + "%" + expx;
[1971]433mImgapp->DispScDrawer( (PIDrawer*)pin, titre, dopt);
[333]434return;
435}
436
437/* --Methode-- */
438void Services2NObjMgr::DisplayPoints3D(string& nom, string& expx, string& expy,
[357]439 string& expz, string& expcut, string dopt, string loop)
[333]440{
441NObjMgrAdapter* obja=NULL;
[2491]442obja = MyObjMgr()->GetObjAdapter(nom);
[333]443if (obja == NULL) {
[449]444 cout << "Services2NObjMgr::DisplayPoints3D() Error , No such object " << nom << endl;
[333]445 return;
446 }
447if (!mImgapp) return;
448
449char* ntn[3] = {"expx","expy","expz"};
450NTuple* nt = new NTuple(3,ntn); // Creation NTuple
451
452string expwt = "1.";
[357]453ComputeExpressions(obja, expx, expy, expz, expwt, expcut, loop, nt, NULL, NULL);
[333]454
455if (nt->NEntry() < 1) {
456 cout << "Services2NObjMgr::DisplayPoints3D() Warning Zero points satisfy cut !" << endl;
457 delete nt;
458 return;
459 }
[2494]460//DBG nt->Show();
461//DBG nt->Print(0,10);
[333]462PINTuple3D* pin = new PINTuple3D(nt, true);
463pin->SelectXYZ(ntn[0], ntn[1], ntn[2]);
[1971]464dopt = "defline " + dopt;
[333]465
466// Pour plot a partir de DispScDrawer
467// string nomdisp = "_NT3D_";
468// mImgapp->DispScDrawer( (PIDrawer*)pin, nomdisp, opt);
469// Pour plot a partir de Disp3DDrawer
470string titre = nom + ":" + expy + "%" + expx;
[1971]471mImgapp->Disp3DDrawer(pin, titre, dopt);
[333]472
473return;
474}
475
476/* --Methode-- */
477void Services2NObjMgr::DisplayPoints2DW(string& nom, string& expx, string& expy,
[3279]478 string& expwt, string& expcut, bool fgcolidx, string dopt, string loop)
[333]479{
480NObjMgrAdapter* obja=NULL;
[2491]481obja = MyObjMgr()->GetObjAdapter(nom);
[333]482if (obja == NULL) {
[449]483 cout << "Services2NObjMgr::DisplayPoints2DW() Error , No such object " << nom << endl;
[333]484 return;
485 }
486if (!mImgapp) return;
487
488char* ntn[3] = {"expx","expy","expw"};
489NTuple* nt = new NTuple(3,ntn); // Creation NTuple
490
491string exp = "1.";
[357]492ComputeExpressions(obja, expx, expy, expwt, exp, expcut, loop, nt, NULL, NULL);
[333]493
494if (nt->NEntry() < 1) {
495 cout << "Services2NObjMgr::DisplayPoints2DW() Warning Zero points satisfy cut !" << endl;
496 delete nt;
497 return;
498 }
499
500PINTuple* pin = new PINTuple(nt, true);
501pin->SelectXY(ntn[0], ntn[1]);
[3279]502if (fgcolidx) pin->SelectColorByIndex(ntn[2]); // Use the weight expression as a color index
503else pin->SelectWt(ntn[2]);
[333]504
505string titre = nom + ":" + expwt + "_" + expy + "%" + expx ;
[1971]506mImgapp->DispScDrawer( (PIDrawer*)pin, titre, dopt);
[333]507return;
508}
509
510/* --Methode-- */
[2999]511void Services2NObjMgr::DisplayPoints3DW(string& nom, string& expx, string& expy, string& expz,
512 string& expwt, string& expcut, string dopt, string loop)
513{
514NObjMgrAdapter* obja=NULL;
515obja = MyObjMgr()->GetObjAdapter(nom);
516if (obja == NULL) {
517 cout << "Services2NObjMgr::DisplayPoints3D() Error , No such object " << nom << endl;
518 return;
519 }
520if (!mImgapp) return;
521
522char* ntn[4] = {"expx","expy","expz","expw"};
523NTuple* nt = new NTuple(4,ntn); // Creation NTuple
524
525ComputeExpressions(obja, expx, expy, expz, expwt, expcut, loop, nt, NULL, NULL);
526
527if (nt->NEntry() < 1) {
528 cout << "Services2NObjMgr::DisplayPoints3DW() Warning Zero points satisfy cut !" << endl;
529 delete nt;
530 return;
531 }
532//DBG nt->Show();
533//DBG nt->Print(0,10);
534PINTuple3D* pin = new PINTuple3D(nt, true);
535pin->SelectXYZ(ntn[0], ntn[1], ntn[2]);
536pin->SelectWt(ntn[3]);
537pin->UseColorScale(true);
538dopt = "defline " + dopt;
539
540// Pour plot a partir de DispScDrawer
541// string nomdisp = "_NT3D_";
542// mImgapp->DispScDrawer( (PIDrawer*)pin, nomdisp, opt);
543// Pour plot a partir de Disp3DDrawer
544string titre = nom + ":" + expy + "%" + expx + " W=" + expwt;
545mImgapp->Disp3DDrawer(pin, titre, dopt);
546
547return;
548}
549
550/* --Methode-- */
[333]551void Services2NObjMgr::ProjectH1(string& nom, string& expx, string& expwt,
[357]552 string& expcut, string& nomh1, string dopt, string loop)
[333]553{
554NObjMgrAdapter* obja=NULL;
[2491]555obja = MyObjMgr()->GetObjAdapter(nom);
[333]556if (obja == NULL) {
[449]557 cout << "Services2NObjMgr::ProjectH1() Error , No such object " << nom << endl;
[333]558 return;
[2681]559}
[333]560if (!mImgapp) return;
561
562Histo* h1 = NULL;
563NTuple* nt = NULL;
564AnyDataObj* oh = NULL;
[2679]565bool h1_already_exist = false;
[2491]566if (nomh1.length() > 0) oh=MyObjMgr()->GetObj(nomh1);
[333]567else nomh1 = "/tmp/projh1d";
[2679]568if ( (oh != NULL) && (typeid(*oh) == typeid(Histo)) ) {
569 h1 = (Histo*)oh;
570 h1_already_exist = true;
571 // Pas de remise a zero ! h1->Zero();
572} else {
[333]573 char* ntn[2]= {"hxval", "hwt"};
574 nt = new NTuple(2,ntn); // Creation NTuple
[2681]575}
[333]576string expz = "0.";
[357]577ComputeExpressions(obja, expx, expwt, expz, expwt, expcut, loop, nt, h1, NULL);
[333]578
579if ((!h1) && (!nt)) return;
580if (!h1) {
581 if (nt->NEntry() < 1) {
582 cout << "Services2NObjMgr::ProjectH1() Warning Zero points satisfy cut !" << endl;
583 delete nt;
584 return;
[2681]585 }
[333]586 double xmin, xmax;
587 nt->GetMinMax(0, xmin, xmax);
588 h1 = new Histo(xmin, xmax, 100);
589 int k;
590 float* xn;
591 for(k=0; k<nt->NEntry(); k++) {
592 xn = nt->GetVec(k);
593 h1->Add(xn[0], xn[1]);
[2681]594 }
[333]595 delete nt;
[2491]596 MyObjMgr()->AddObj(h1, nomh1);
[2681]597}
[333]598
[2681]599if(!h1_already_exist || dopt.size()>0) MyObjMgr()->DisplayObj(nomh1, dopt);
[333]600return;
601}
602
603/* --Methode-- */
604void Services2NObjMgr::ProjectH2(string& nom, string& expx, string& expy, string& expwt,
[357]605 string& expcut, string& nomh2, string dopt, string loop)
[333]606{
607NObjMgrAdapter* obja=NULL;
[2491]608obja = MyObjMgr()->GetObjAdapter(nom);
[333]609if (obja == NULL) {
[449]610 cout << "Services2NObjMgr::ProjectH2() Error , No such object " << nom << endl;
[333]611 return;
[2681]612}
[333]613if (!mImgapp) return;
614
615Histo2D* h2 = NULL;
616NTuple* nt = NULL;
617AnyDataObj* oh = NULL;
[2679]618bool h2_already_exist = false;
[2491]619if (nomh2.length() > 0) oh=MyObjMgr()->GetObj(nomh2);
[333]620else nomh2 = "/tmp/projh2d";
[2679]621if ( (oh != NULL) && (typeid(*oh) == typeid(Histo2D)) ) {
622 h2 = (Histo2D*)oh;
623 h2_already_exist = true;
624 // Pas de remise a zero ! h2->Zero();
625} else {
[333]626 char* ntn[3]= {"hxval", "hyval", "hwt"};
627 nt = new NTuple(3,ntn); // Creation NTuple
[2681]628}
[333]629string expz = "0.";
[357]630ComputeExpressions(obja, expx, expy, expwt, expwt, expcut, loop, nt, NULL, h2);
[333]631
632if ((!h2) && (!nt)) return;
633if (!h2) {
634 if (nt->NEntry() < 1) {
635 cout << "Services2NObjMgr::ProjectH2() Warning Zero points satisfy cut !" << endl;
636 delete nt;
637 return;
[2681]638 }
[333]639 double xmin, xmax, ymin, ymax;
640 nt->GetMinMax(0, xmin, xmax);
[466]641 nt->GetMinMax(1, ymin, ymax);
[333]642 h2 = new Histo2D(xmin, xmax, 50, ymin, ymax, 50);
643 int k;
644 float* xn;
645 for(k=0; k<nt->NEntry(); k++) {
646 xn = nt->GetVec(k);
647 h2->Add(xn[0], xn[1], xn[2]);
[2681]648 }
[333]649 delete nt;
[2491]650 MyObjMgr()->AddObj(h2, nomh2);
[2681]651}
[333]652
[2681]653if(!h2_already_exist || dopt.size()>0) MyObjMgr()->DisplayObj(nomh2, dopt);
[333]654return;
655
656}
657
658/* --Methode-- cmv 13/10/98 */
659void Services2NObjMgr::ProjectHProf(string& nom, string& expx, string& expy, string& expwt,
[357]660 string& expcut, string& nomprof, string dopt, string loop)
661// Pour remplir un ``GeneralFitData'' a partir de divers objets:
[333]662//| nom = nom de l'objet a projeter dans un HProf.
663//| expx = expression X de definition du bin.
664//| expy = expression Y a additionner dans le bin.
665//| expwt = expression W du poids a additionner.
666//| expcut = expression du test de selection.
667//| nomprof = nom du HProf engendre (optionnel). Si l'objet n'existe pas
668//| les limites Xmin,Xmax sont calculees automatiquement.
669//| sinon ce sont celles de l'objet preexistant.
670//| opt = options generales pour le display.
671{
672NObjMgrAdapter* obja=NULL;
[2491]673obja = MyObjMgr()->GetObjAdapter(nom);
[333]674if (obja == NULL) {
[449]675 cout << "Services2NObjMgr::ProjectHProf() Error , No such object " << nom << endl;
[333]676 return;
[2681]677}
[333]678if (!mImgapp) return;
679
680HProf* hprof = NULL;
681NTuple* nt = NULL;
682AnyDataObj* oh = NULL;
[2679]683bool hp_already_exist = false;
[2491]684if (nomprof.length() > 0) oh=MyObjMgr()->GetObj(nomprof);
[333]685else nomprof = "/tmp/projprof";
[2679]686if( (oh!=NULL) && (typeid(*oh) == typeid(HProf)) ) {
687 hprof = (HProf*)oh;
688 hp_already_exist = true;
689} else {
[333]690 char* ntn[3]= {"hxval", "hyval", "hwt"};
691 nt = new NTuple(3,ntn); // Creation NTuple
692}
693string expz = "0.";
[357]694ComputeExpressions(obja, expx, expy, expwt, expwt, expcut, loop, nt, NULL, NULL, hprof);
[333]695
696if((!hprof) && (!nt)) return;
697if(!hprof) {
698 if (nt->NEntry() < 1) {
699 cout << "Services2NObjMgr::ProjectHProf() Warning Zero points satisfy cut !" << endl;
700 delete nt;
701 return;
702 }
[1091]703 r_8 xmin, xmax;
[333]704 nt->GetMinMax(0, xmin, xmax);
705 hprof = new HProf(xmin, xmax, 100);
706 int k;
707 float* xn;
708 for(k=0; k<nt->NEntry(); k++) {
709 xn = nt->GetVec(k);
710 hprof->Add(xn[0], xn[1], xn[2]);
[2681]711 }
[333]712 delete nt;
[2491]713 MyObjMgr()->AddObj(hprof, nomprof);
[2681]714}
[1090]715hprof->UpdateHisto();
[333]716
[2681]717if(!hp_already_exist || dopt.size()>0) MyObjMgr()->DisplayObj(nomprof, dopt);
[333]718return;
719}
720
[357]721
[333]722/* --Methode-- */
[357]723void Services2NObjMgr::FillVect(string& nom, string& expx, string& expv,
724 string& expcut, string& nomvec, string dopt, string loop)
[333]725{
726NObjMgrAdapter* obja=NULL;
[2491]727obja = MyObjMgr()->GetObjAdapter(nom);
[333]728if (obja == NULL) {
[357]729 cout << "Services2NObjMgr::FillVect() Error , No such object: " << nom << endl;
[333]730 return;
731 }
732if (!mImgapp) return;
733
[357]734Vector* v1 = NULL;
735AnyDataObj* ov = NULL;
[2491]736ov=MyObjMgr()->GetObj(nomvec);
[357]737if (ov != NULL) v1 = dynamic_cast<Vector *>(ov);
738if (v1 == NULL) {
739 cout << "Services2NObjMgr::FillVect() Error , No such object or not a vector: " << nomvec << endl;
740 return;
741 }
742
743char* ntn[2]= {"vi", "vv"};
744NTuple* nt = new NTuple(2,ntn); // Creation NTuple
745
746string expz = "0.";
747ComputeExpressions(obja, expx, expv, expz, expz, expcut, loop, nt);
748
749if (!nt) return;
750if (nt->NEntry() < 1) {
751 cout << "Services2NObjMgr::FillVect() Warning Zero points satisfy cut !" << endl;
752 delete nt;
753 return;
754 }
755
756 int i,k;
757 double* xn;
758 for(k=0; k<nt->NEntry(); k++) {
759 xn = nt->GetLineD(k);
[2180]760 i = int(xn[0]+0.5);
[357]761 if ( (i < 0) || i >= v1->NElts() ) continue;
762 (*v1)(i) = xn[1];
763 }
764 delete nt;
765
766
[2491]767MyObjMgr()->DisplayObj(nomvec, dopt);
[357]768return;
769}
770
771/* --Methode-- */
772void Services2NObjMgr::FillMatx(string& nom, string& expx, string& expy, string& expv,
773 string& expcut, string& nommtx, string dopt, string loop)
774{
775NObjMgrAdapter* obja=NULL;
[2491]776obja = MyObjMgr()->GetObjAdapter(nom);
[357]777if (obja == NULL) {
778 cout << "Services2NObjMgr::FillMatx() Error , No such objet " << nom << endl;
779 return;
780 }
781if (!mImgapp) return;
782
783Matrix* mtx = NULL;
784AnyDataObj* om = NULL;
[2491]785om=MyObjMgr()->GetObj(nommtx);
[357]786if (om != NULL) mtx = dynamic_cast<Matrix *>(om);
787if (mtx == NULL) {
788 cout << "Services2NObjMgr::FillMatx() Error , No such object or not a matrix " << nommtx << endl;
789 return;
790 }
791
792char* ntn[3]= {"mi", "mj", "mv"};
793NTuple* nt = new NTuple(3,ntn); // Creation NTuple
794
795string expz = "0.";
796ComputeExpressions(obja, expx, expy, expv, expz, expcut, loop, nt);
797
798if (!nt) return;
799if (nt->NEntry() < 1) {
800 cout << "Services2NObjMgr::FillMatx() Warning Zero points satisfy cut !" << endl;
801 delete nt;
802 return;
803 }
804
805 int ic, jl, k;
806 double* xn;
807 for(k=0; k<nt->NEntry(); k++) {
808 xn = nt->GetLineD(k);
[2180]809 ic = int(xn[0]+0.5);
810 jl = int(xn[1]+0.5);
[357]811 if ( (ic < 0) || ic >= mtx->NCol() ) continue;
812 if ( (jl < 0) || jl >= mtx->NRows() ) continue;
813 (*mtx)(jl, ic) = xn[2];
814 }
815 delete nt;
816
817
[2491]818MyObjMgr()->DisplayObj(nommtx, dopt);
[357]819return;
820
821}
822
823/* --Methode-- */
824void Services2NObjMgr::ExpressionToVector(string& nom, string& expx, string& expcut,
825 string& nomvec, string dopt, string loop)
826{
827NObjMgrAdapter* obja=NULL;
[2491]828obja = MyObjMgr()->GetObjAdapter(nom);
[357]829if (obja == NULL) {
[449]830 cout << "Services2NObjMgr::ExpressionToVector() Error , No such object " << nom << endl;
[357]831 return;
832 }
833if (!mImgapp) return;
834
[333]835NTuple* nt = NULL;
[2924]836// if (nomvec.length() < 1) nomvec = "/tmp/expvec";
[333]837
838char* ntn[2]= {"vecval", "vecwt"};
839nt = new NTuple(1,ntn); // Creation NTuple
840
841string expwt = "1.";
[1486]842string expz = "0.";
843string dumexpcut = expcut; if(dumexpcut.size()<=0) dumexpcut = "1.";
844ComputeExpressions(obja,expx,expz,expz,expwt,dumexpcut,loop,nt,NULL,NULL);
[333]845
846if (!nt) return;
847if (nt->NEntry() < 1) {
[357]848 cout << "Services2NObjMgr::ExpressionToVector() Warning Zero points satisfy cut !" << endl;
[333]849 delete nt;
850 return;
851 }
852
853Vector* vec = new Vector(nt->NEntry());
854int k;
855float* xn;
856for(k=0; k<nt->NEntry(); k++) {
857 xn = nt->GetVec(k);
858 (*vec)(k) = xn[0];
859 }
860delete nt;
[2922]861if (nomvec.size() > 0) {
862 MyObjMgr()->AddObj(vec, nomvec);
863 MyObjMgr()->DisplayObj(nomvec, dopt);
864}
865else {
866// On calcule et on affiche mean/sigma + min/max
867 double min, max, mean, sigma;
868 vec->MinMax(min, max);
869 MathArray<r_8> ma;
870 ma.MeanSigma(*vec, mean, sigma);
[2924]871 cout << "Size= " << vec->Size() << " Mean= " << mean << " Sigma= " << sigma
[2922]872 << " Min= " << min << " Max= " << max << endl;
873 delete vec;
874}
[333]875return;
876}
877
878/* --Methode-- */
[447]879void Services2NObjMgr::NtFromASCIIFile(string& nom,string& filename,double def_val)
[2792]880// Pour remplir un NTuple/DataTable "nom" existant a partir du fichier
[447]881// ASCII table "filename". Si il y a plus de variables dans le
882// ntuple que dans le fichier "filename",
[2792]883// les sur-numeraires numeriques sont mises a "def_val" par defaut.
[447]884{
[2491]885AnyDataObj* mobj = MyObjMgr()->GetObj(nom);
[2792]886if(mobj == NULL) {
887 cout<<"NtFromASCIIFile() Error, object "<<nom<<" not existing"<<endl;
888 return;
889}
890
891NTuple* nt = NULL;
892DataTable* dt = NULL;
893if(typeid(*mobj) == typeid(NTuple)) nt = (NTuple*) mobj;
894else if(typeid(*mobj) == typeid(DataTable)) dt = (DataTable*) mobj;
895
896if(nt==NULL && dt==NULL) {
897 cout<<"NtFromASCIIFile() Error, object "<<nom<<" not an NTuple nor a DataTable"<<endl;
898 return;
899}
[447]900if (!mImgapp) return;
901
[2792]902if(nt) nt->FillFromASCIIFile(filename, def_val);
[2833]903else if(dt) {
904 ifstream is(filename.c_str());
905 dt->FillFromASCIIFile(is);
906}
[447]907return;
908}
909
910/* --Methode-- */
[333]911void Services2NObjMgr::FillNT(string& nom, string& expx, string& expy, string& expz,
[357]912 string& expt, string& expcut, string& nomnt, string loop)
[333]913{
914NObjMgrAdapter* obja=NULL;
[2491]915obja = MyObjMgr()->GetObjAdapter(nom);
[333]916if (obja == NULL) {
[449]917 cout << "Services2NObjMgr::FillNT() Error , No such object " << nom << endl;
[333]918 return;
919 }
920if (!mImgapp) return;
921
922bool fgnnt = false;
923NTuple* nt = NULL;
924AnyDataObj* oh = NULL;
[2491]925if (nomnt.length() > 0) oh=MyObjMgr()->GetObj(nomnt);
[333]926else nomnt = "/tmp/fillnt";
927if ( (oh != NULL) && (typeid(*oh) == typeid(NTuple)) ) {
928 nt = (NTuple*)oh;
929 if (nt->NVar() > 10) {
[449]930 cout << "Services2NObjMgr::FillNT() Warning , Max 10 var in NTuple -> new NTuple" << endl;
[333]931 nt = NULL;
932 }
933 }
934if (nt == NULL) {
935 char* ntn[4]= {"x", "y","z","t"};
936 nt = new NTuple(4,ntn); // Creation NTuple
937 fgnnt = true;
938 }
939
[357]940ComputeExpressions(obja, expx, expy, expz, expt, expcut, loop, nt, NULL, NULL);
[333]941
[2491]942if (fgnnt) MyObjMgr()->AddObj(nt, nomnt);
[333]943return;
944
945}
946
947/* --Methode-- */
948void Services2NObjMgr::FillNTFrCFile(string & nom, string const & fname,
[357]949 string const & funcname, string & nomnt, string loop)
[333]950{
951if (!mImgapp) return;
952
953NObjMgrAdapter* obja=NULL;
[2491]954obja = MyObjMgr()->GetObjAdapter(nom);
[333]955if (obja == NULL) {
956 cout << "Services2NObjMgr::FillNTFrCFile( " << nom << "...) No such object" <<endl;
957 return;
958 }
[344]959bool adel = true;
960NTupleInterface* objnt = obja->GetNTupleInterface(adel);
[333]961if (objnt == NULL) {
[449]962 cout << "Services2NObjMgr::FillNTFrCFile( " << nom << "...) Not an NTupleInterface !" <<endl;
[333]963 return;
964 }
965
[2755]966// Pour synchronisation d'execution simultanee
[2762]967ZSync zs(mutx_dynlink); zs.NOp();
[2755]968
[333]969NTLoopExprFunc f = (NTLoopExprFunc)LinkFunctionFromFile(fname, funcname);
970if (!f) {
971 cerr << "Services2NObjMgr::FillNTFrCFile Error Creation NTLoopExprFunc" << endl;
[344]972 if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]973 return;
974 }
975
976bool fgnnt = false;
977NTuple* nt = NULL;
978if (nomnt.length() > 0) {
979 AnyDataObj* oh = NULL;
[2491]980 oh=MyObjMgr()->GetObj(nomnt);
[333]981 if ( (oh != NULL) && (typeid(*oh) == typeid(NTuple)) ) {
982 nt = (NTuple*)oh;
983 if (nt->NVar() > 10) {
[449]984 cout << "Services2NObjMgr::FillNTFrCFile() Warning , Max 10 var in NTuple -> new NTuple" << endl;
[333]985 nt = NULL;
986 }
987 }
988 if (nt == NULL) {
989 char* ntn[4]= {"x", "y","z","t"};
990 nt = new NTuple(4,ntn); // Creation NTuple
991 fgnnt = true;
992 }
993 }
994
995double xnt[10];
996
997int i,k;
[2668]998for(i=0; i<10; i++) xnt[i] = 0.;
[333]999
[2419]1000int_8 k1,k2,dk;
[357]1001k1 = 0; k2 = objnt->NbLines(); dk = 1;
1002DecodeLoopParameters(loop, k1, k2, dk);
1003if (k1 < 0) k1 = 0;
1004if (k2 < 0) k2 = objnt->NbLines();
[2419]1005if (k2 > (int_8)objnt->NbLines()) k2 = objnt->NbLines();
[357]1006if (dk <= 0) dk = 1;
1007
[500]1008try {
[333]1009 double* xn;
[2419]1010 int_8 kstart = k1, kend = k2;
1011 for(k=kstart; k<kend; k+=dk) {
[333]1012 xn = objnt->GetLineD(k);
[2419]1013 if (f((int_8_exprf)k, xn, xnt, xnt+1, xnt+2, xnt+3, (int_8_exprf)kstart,(int_8_exprf) kend) != 0) {
[2668]1014 if (nt) nt->Fill(xnt);
[333]1015 }
1016 }
1017 }
[500]1018#ifdef SANS_EVOLPLANCK
[333]1019CATCH(merr) {
1020 fflush(stdout);
1021 cout << endl;
1022 cerr << endl;
1023 string es = PeidaExc(merr);
1024 cerr << "Services2NObjMgr::FillNTFrCFile() Exception :" << merr << es;
1025 } ENDTRY;
[500]1026#else
1027catch ( PException exc ) {
1028 fflush(stdout);
1029 cout << endl;
1030 cerr << endl;
1031 cerr << "Services2NObjMgr::FillNTFrCFile() Exception :" << exc.Msg() << endl;
1032}
1033#endif
[344]1034
1035if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]1036CloseDLL();
1037
[2491]1038if (fgnnt) MyObjMgr()->AddObj(nt, nomnt);
[333]1039return;
1040}
1041
1042/* --Methode-- */
1043void Services2NObjMgr::PrepareNTExpressionCFile(string & nom, string const & fname,
1044 string const & funcname)
1045{
1046NObjMgrAdapter* obja=NULL;
[2491]1047obja = MyObjMgr()->GetObjAdapter(nom);
[333]1048if (obja == NULL) {
1049 cout << "Services2NObjMgr::PrepareNTExpressionCFile( " << nom << "...) No such object" <<endl;
1050 return;
1051 }
[344]1052bool adel = true;
1053NTupleInterface* objnt = obja->GetNTupleInterface(adel);
[333]1054if (objnt == NULL) {
1055 cout << "Services2NObjMgr::PrepareNTExpressionCFile( " << nom
1056 << "...) No NTupleInterface !" <<endl;
1057 return;
1058 }
1059string vardec = objnt->VarList_C("_xnti_");
1060
1061FILE *fip;
1062if ((fip = fopen(fname.c_str(), "w")) == NULL) {
1063 cout << "Services2NObjMgr::PrepareNTExpressionCFile()_Error: fopen " << fname << endl;
[344]1064 if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]1065 return;
1066 }
[344]1067
[333]1068// constitution du fichier des decalarations des variables de l'interface NTuple
1069fputs("#include <stdlib.h> \n", fip);
1070fputs("#include <stdio.h> \n", fip);
1071fputs("#include <math.h> \n\n", fip);
[344]1072
[1469]1073fputs("/* ------ Compare bits on double --------- */ \n", fip);
[2419]1074fputs("typedef long long int_8_exprf;\n", fip);
1075fputs("int_8_exprf BitCmp64(double v,int_8_exprf flg)\n", fip);
1076fputs("{return ((int_8_exprf)((v<0.) ? v-0.1 : v+0.1))&flg;}\n", fip);
[344]1077fputs("/* ------ Some random number generators --------- */ \n", fip);
[1931]1078fputs("#if defined(__ppc__) && defined(__MACH__) \n",fip);
1079fputs("#include <limits.h> \n", fip);
1080fputs("#define drand48() ((double)(random())/LONG_MAX) \n",fip);
1081fputs("#endif \n",fip);
[344]1082fputs("#define frand01() ( (float) drand48() ) \n", fip);
1083fputs("#define drand01() drand48() \n", fip);
1084fputs("#define rand01() drand48() \n", fip);
1085fputs("#define frandpm1() ( 2. * frand01() - 1.) \n", fip);
1086fputs("#define drandpm1() ( 2. * drand01() - 1.) \n", fip);
1087fputs("#define randpm1() ( 2. * drand01() - 1.) \n", fip);
1088fputs("double NorRand(void) \n", fip);
1089fputs(" { \n double x,A,B; \n LAB10: \n A = drand01(); \n", fip);
1090fputs(" if ( A == 0. ) goto LAB10; \n B = drand01(); \n", fip);
[345]1091fputs(" x = sqrt(-2.*log(A))*cos(2.*M_PI*B); \n", fip);
[344]1092fputs(" return(x); \n } \n", fip);
1093fputs("#define GauRand() NorRand() \n", fip);
1094fputs("#define gaurand() NorRand() \n\n", fip);
1095
[3435]1096fputs("/* Angle conversions */\n", fip);
1097fputs("double deg2rad(double x) { return (x*M_PI/180.); } \n", fip);
1098fputs("double rad2deg(double x) { return (x*180./M_PI); } \n", fip);
1099fputs("/* Theta(0..Pi, phi(0..2Pi) conversion to Molleweide X,Y */\n", fip);
1100fputs("double tetphi2mollX(double teta, double phi) { \n", fip);
1101fputs(" if (phi<=M_PI) return(180.-sin(acos(2*teta/M_PI-1.))*((M_PI-phi)*180./M_PI)); \n", fip);
1102fputs(" else return(180.+sin(acos(2.*teta/M_PI-1.))*((phi-M_PI)*180./M_PI)); } \n", fip);
1103fputs("double tetphi2mollY(double teta) { \n", fip);
1104fputs(" return (90.-(teta*180./M_PI)); } \n\n", fip);
1105fputs("/* Longitude(0..360), Latitude(-90..90) conversion to Molleweide X,Y */\n", fip);
1106fputs("double longlat2mollX(double lng, double lat) { \n", fip);
1107fputs(" if (lng<=180.) return(180.-sin(acos(lat/90.))*(180.-lng)); \n", fip);
1108fputs(" else return(180.+sin(acos(lat/90.))*(lng-180.)); } \n", fip);
1109fputs("double longlat2mollY(double lat) { return lat; } \n\n\n", fip);
1110
[2419]1111fputs("/* NTupleInterface Variable declaration - Generated by piapp */\n", fip);
1112fputs("/* -- Services2NObjMgr::PrepareNTExpressionCFile() -- */ \n", fip);
1113fputs("/* _nl line number or sequential index : _nstart <= _nl < _nend */ \n\n", fip);
1114fprintf(fip,"int %s(int_8_exprf _nl, double* _xnti_, double* _rx_, double* _ry_, double* _rz_, \n",
[333]1115 funcname.c_str());
[2419]1116fprintf(fip," double* _rt_, int_8_exprf _nstart, int_8_exprf _nend) \n");
[333]1117fprintf(fip, "{ \n %s \n", vardec.c_str());
1118fputs(" if (!1) { /* Cut Expression failed */ \n", fip);
1119fputs(" *_rx_ = *_ry_ = *_rz_ = *_rt_ = 0.; return(0);", fip);
1120fputs(" } \n /* Cut expression satisfied */ \n", fip);
1121fputs(" *_rx_ = 1.; \n *_ry_ = 1.; \n *_rz_ = 1.; \n *_rt_ = 1.; \n", fip);
1122fputs(" return(1); \n} \n", fip);
1123
1124fclose(fip);
[344]1125
1126if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]1127return;
1128}
1129
1130/* --Methode-- cmv 13/10/98 */
1131void Services2NObjMgr::FillGFD(string& nom, string& expx, string& expy, string& expz,
[357]1132 string& experr, string& expcut, string& nomgfd, string loop)
[333]1133// Pour remplir un ``GeneralFitData'' a partir de divers objets:
1134//| nom = nom de l'objet a transcrire selon 1D: Z=f(X) ou 2D: Z=f(X,Y) .
1135//| Vector,Matrix,Histo,HProf,Histo2D,Image<T>,StarList,NTuple,GeneralFitData
1136//| expx = expression X du GeneralFitData (1er abscisse)
1137//| expy = expression Y du GeneralFitData (2sd abscisse si non "", Z=f(X,Y))
1138//| expz = expression Z du GeneralFitData (valeur de l'ordonnee)
1139//| experr = expression de l'erreur sur l'ordonnee Z
1140//| expcut = expression du test de selection
1141//| nomgfd = nom du GeneralFitData engendre (optionnel)
1142{
1143NObjMgrAdapter* obja=NULL;
[2491]1144obja = MyObjMgr()->GetObjAdapter(nom);
[333]1145if (obja == NULL) {
[449]1146 cout << "Services2NObjMgr::FillGFD() Error , No such object "<<nom<<endl;
[333]1147 return;
1148 }
1149if(!mImgapp) return;
1150
1151// 2D ou 3D?
1152int nvar = 2;
1153if(expy.length()<=0) {nvar = 1; expy = "0.";}
1154
1155// Creation NTuple Buffer
1156char* ntn[4]= {"x","y","f","e"};
1157NTuple*nt = new NTuple(4,ntn);
1158
1159// Remplissage NTuple buffer
[357]1160ComputeExpressions(obja, expx, expy, expz, experr, expcut, loop, nt, NULL, NULL);
[333]1161if(nt->NEntry() < 1)
1162 {cout<<"Services2NObjMgr::FillGFD() Warning Zero points satisfy cut !"<<endl;
1163 delete nt; return;}
1164
1165//Remplissage de la structure GeneraFitData
1166if (nt->NEntry() <= 0) {
1167 cout<<"Services2NObjMgr::FillGFD() Warning - NData= " << nt->NEntry() << endl;
1168 delete nt;
1169 return;
1170 }
1171
1172GeneralFitData* gfd = new GeneralFitData(nvar,nt->NEntry(),0);
1173int k;
1174float* xn;
1175for(k=0; k<nt->NEntry(); k++) {
1176 xn = nt->GetVec(k);
1177 gfd->AddData(xn,xn[2],xn[3]);
1178}
1179
1180// Menage et table d'objets
1181delete nt;
[2491]1182MyObjMgr()->AddObj(gfd, nomgfd);
[333]1183return;
1184}
1185
[1067]1186/* --Methode-- cmv 12/07/00 */
1187void Services2NObjMgr::FillGFDfrVec(string nomgfd,string namx,string namy,string namz,string name)
1188// Pour remplir un ``GeneralFitData'' a partir de vecteurs
1189//| gdfrvec nomgd X Y ! !
1190//| gdfrvec nomgd X Y ! EY
1191//| gdfrvec nomgd X Y Z !
1192//| gdfrvec nomgd X Y Z EZ
1193//| - nomgfd = nom du generaldata a remplir
1194//| - namx = nom du vecteur contenant les valeurs X
1195//| - namy = nom du vecteur contenant les valeurs Y
1196//| - namz = nom du vecteur contenant les valeurs Z (ou "!")
1197//| - name = nom du vecteur contenant les valeurs des erreurs EY ou EZ
1198{
1199// Decodage des noms des vecteurs pour le remplissage du generaldata
1200if(nomgfd=="!" || nomgfd.length()<1)
1201 {cout<<"FillGFDfrVec_Error: bad GenaralData name "<<nomgfd<<endl; return;}
1202if(namx=="!" || namx.length()<1)
1203 {cout<<"FillGFDfrVec_Error: bad X vector name "<<namx<<endl; return;}
1204if(namy=="!" || namy.length()<1)
1205 {cout<<"FillGFDfrVec_Error: bad Y vector name "<<namy<<endl; return;}
1206if(namz.length()<1) namz = "!";
1207if(name.length()<1) name = "!";
1208int nvar = 0;
1209if(namz=="!") nvar = 1; else nvar = 2;
[333]1210
[1067]1211// Identify data
1212NamedObjMgr omg;
1213AnyDataObj* mobj = NULL;
1214Vector* v;
1215int nel = 0;
1216r_8 *x=NULL, *y=NULL, *z=NULL,* ez=NULL;
1217
1218if( (mobj=omg.GetObj(namx)) == NULL) {
1219 cout<<"FillGFDfrVec_Error: unknown X object "<<namx<<endl; return;
1220} else {
1221 v = (Vector*) mobj; x = v->Data(); nel=v->NElts();
1222}
1223
1224if( (mobj=omg.GetObj(namy)) == NULL) {
1225 cout<<"FillGFDfrVec_Error: unknown Y object "<<namy<<endl; return;
1226} else {
1227 v = (Vector*) mobj; y = z = v->Data(); if(v->NElts()<nel) nel=v->NElts();
1228}
1229
1230if( nvar==2 && (mobj=omg.GetObj(namz)) == NULL) {
1231 cout<<"FillGFDfrVec_Error: unknown Z object "<<namz<<endl; return;
1232} else {
1233 v = (Vector*) mobj; z = v->Data(); if(v->NElts()<nel) nel=v->NElts();
1234}
1235
1236if(name!="!") {
1237 if( (mobj=omg.GetObj(name)) == NULL) {
1238 cout<<"FillGFDfrVec_Error: unknown EZ object "<<name<<endl; return;
1239 } else {
1240 v = (Vector*) mobj; ez = v->Data(); if(v->NElts()<nel) nel=v->NElts();
1241 }
1242}
1243
1244if(nel<=0)
1245 {cout<<"FillGFDfrVec_Error: bad number of elements "<<nel<<endl; return;}
1246
1247// Create GeneralData and fill it with vectors
1248GeneralFitData* gfd = new GeneralFitData(nvar,nel+5,0);
1249if(nvar==1) gfd->SetData1(nel,x,z,ez); // On remplit Y=f(X)
1250else gfd->SetData2(nel,x,y,z,ez); // On remplit Z=f(X,y)
1251
1252// Menage et table d'objets
1253if( omg.GetObj(nomgfd) != NULL ) omg.DelObj(nomgfd);
[2491]1254MyObjMgr()->AddObj(gfd,nomgfd);
[1067]1255return;
1256}
1257
[333]1258/* --Methode-- */
1259void Services2NObjMgr::ComputeExpressions(NObjMgrAdapter* obja, string& expx,
[357]1260 string& expy, string& expz, string& expt, string& expcut, string& loop,
[333]1261 NTuple* nt, Histo* h1, Histo2D* h2, HProf* hp)
1262{
1263if (obja == NULL) return;
[344]1264bool adel = true;
1265NTupleInterface* objnt = obja->GetNTupleInterface(adel);
[333]1266if (objnt == NULL) return;
1267string vardec = objnt->VarList_C("_zz6qi_");
1268
[2755]1269// Pour synchronisation d'execution simultanee
[2762]1270ZSync zs(mutx_dynlink); zs.NOp();
[2755]1271
[333]1272PlotExprFunc f = LinkExprFunc(vardec, expx, expy, expz, expt, expcut);
1273if (!f) {
1274 cerr << "Services2NObjMgr::::ComputeExpressions() Error Creation PlotExprFunc " << endl;
[344]1275 if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]1276 return;
1277 }
1278
1279double xnt[10];
1280
[2419]1281int_8 k;
1282for(k=0; k<10; k++) xnt[k] = 0.;
1283int_8 k1,k2,dk;
[357]1284k1 = 0; k2 = objnt->NbLines(); dk = 1;
1285DecodeLoopParameters(loop, k1, k2, dk);
1286if (k1 < 0) k1 = 0;
1287if (k2 < 0) k2 = objnt->NbLines();
[2419]1288if (k2 > (int_8)objnt->NbLines()) k2 = objnt->NbLines();
[357]1289if (dk <= 0) dk = 1;
[500]1290
1291try {
[333]1292 double* xn;
[357]1293 for(k=k1; k<k2; k += dk) {
[333]1294 xn = objnt->GetLineD(k);
[2419]1295 if (f((int_8_exprf)k,xn, xnt, xnt+1, xnt+2, xnt+3) != 0) {
[2668]1296 if (nt) nt->Fill(xnt);
[333]1297 if (h1) h1->Add(xnt[0], xnt[3]);
1298 if (h2) h2->Add(xnt[0], xnt[1], xnt[3]);
1299 if (hp) hp->Add(xnt[0], xnt[1], xnt[3]);
1300 }
1301 }
1302 }
[500]1303#ifdef SANS_EVOLPLANCK
[333]1304CATCH(merr) {
1305 fflush(stdout);
1306 cout << endl;
1307 cerr << endl;
1308 string es = PeidaExc(merr);
1309 cerr << "Services2NObjMgr::ComputeExpressions() Exception :" << merr << es;
1310 } ENDTRY;
[500]1311#else
1312catch ( PException exc ) {
1313 fflush(stdout);
1314 cout << endl;
1315 cerr << endl;
1316 cerr << "Services2NObjMgr::ComputeExpressions() Exception :" << exc.Msg() << endl;
1317}
1318#endif
[333]1319
[344]1320if (adel) delete objnt; // Delete de l'objet NTupleInterface si necessaire
[333]1321// Fermeture du fichier .so
1322CloseDLL();
1323return;
1324}
1325
1326
1327/* --Methode-- */
1328PlotExprFunc Services2NObjMgr::LinkExprFunc(string& vardec, string& expx, string& expy,
1329 string& expz, string& expt, string& cut)
1330{
1331FILE *fip;
1332string fname = TmpDir + "expf_pia_dl.c";
1333string cmd;
1334int rc;
1335
1336cmd = "rm -f " + fname;
1337rc = system(cmd.c_str());
1338//DBG printf("LinkExprFunc_Do> %s (Rc=%d)\n", cmd.c_str(), rc);
1339
1340if ((fip = fopen(fname.c_str(), "w")) == NULL) {
1341 string sn = fname;
[449]1342 cout << "Services2NObjMgr/LinkExprFunc_Error: fopen( " << sn << endl;
[333]1343 return(NULL);
1344 }
1345
1346// constitution du fichier a compiler
1347fputs("#include <stdlib.h> \n", fip);
1348fputs("#include <math.h> \n", fip);
[344]1349
[1469]1350fputs("/* ------ Compare bits on double --------- */ \n", fip);
[2419]1351fputs("typedef long long int_8_exprf;\n", fip);
1352fputs("int_8_exprf BitCmp64(double v,int_8_exprf flg)\n", fip);
1353fputs("{return ((int_8_exprf)((v<0.) ? v-0.1 : v+0.1))&flg;}\n", fip);
[344]1354fputs("/* ------ Some random number generators --------- */ \n", fip);
[1931]1355fputs("#if defined(__ppc__) && defined(__MACH__) \n",fip);
1356fputs("#include <limits.h> \n", fip);
1357fputs("#define drand48() ((double)(random())/LONG_MAX) \n",fip);
1358fputs("#endif \n",fip);
[344]1359fputs("#define frand01() ( (float) drand48() ) \n", fip);
1360fputs("#define drand01() drand48() \n", fip);
1361fputs("#define rand01() drand48() \n", fip);
1362fputs("#define frandpm1() ( 2. * frand01() - 1.) \n", fip);
1363fputs("#define drandpm1() ( 2. * drand01() - 1.) \n", fip);
1364fputs("#define randpm1() ( 2. * drand01() - 1.) \n", fip);
1365fputs("double NorRand(void) \n", fip);
1366fputs(" { \n double x,A,B; \n LAB10: \n A = drand01(); \n", fip);
1367fputs(" if ( A == 0. ) goto LAB10; \n B = drand01(); \n", fip);
[345]1368fputs(" x = sqrt(-2.*log(A))*cos(2.*M_PI*B); \n", fip);
[344]1369fputs(" return(x); \n } \n", fip);
1370fputs("#define GauRand() NorRand() \n", fip);
1371fputs("#define gaurand() NorRand() \n\n", fip);
1372
[3435]1373fputs("/* Angle conversions */\n", fip);
1374fputs("double deg2rad(double x) { return (x*M_PI/180.); } \n", fip);
1375fputs("double rad2deg(double x) { return (x*180./M_PI); } \n", fip);
1376fputs("/* Theta(0..Pi, phi(0..2Pi) conversion to Molleweide X,Y */\n", fip);
1377fputs("double tetphi2mollX(double teta, double phi) { \n", fip);
1378fputs(" if (phi<=M_PI) return(180.-sin(acos(2*teta/M_PI-1.))*((M_PI-phi)*180./M_PI)); \n", fip);
1379fputs(" else return(180.+sin(acos(2.*teta/M_PI-1.))*((phi-M_PI)*180./M_PI)); } \n", fip);
1380fputs("double tetphi2mollY(double teta) { \n", fip);
1381fputs(" return (90.-(teta*180./M_PI)); } \n\n", fip);
1382fputs("/* Longitude(0..360), Latitude(-90..90) conversion to Molleweide X,Y */\n", fip);
1383fputs("double longlat2mollX(double lng, double lat) { \n", fip);
1384fputs(" if (lng<=180.) return(180.-sin(acos(lat/90.))*(180.-lng)); \n", fip);
1385fputs(" else return(180.+sin(acos(lat/90.))*(lng-180.)); } \n", fip);
1386fputs("double longlat2mollY(double lat) { return lat; } \n\n\n", fip);
1387
[2419]1388fputs("int expf_pia_dl_func(int_8_exprf _nl, double* _zz6qi_, double* _rx_6q_, double* _ry_6q_, double* _rz_6q_, double* _rt_6q_) \n{\n", fip);
[333]1389fprintf(fip,"%s \n", vardec.c_str());
1390fprintf(fip, "if (!(%s)) { *_rx_6q_ = *_ry_6q_ = *_rz_6q_ = *_rt_6q_ = 0.; return(0); } \n", cut.c_str());
1391fprintf(fip, "*_rx_6q_ = %s ; \n", expx.c_str());
1392fprintf(fip, "*_ry_6q_ = %s ; \n", expy.c_str());
1393fprintf(fip, "*_rz_6q_ = %s ; \n", expz.c_str());
1394fprintf(fip, "*_rt_6q_ = %s ; \n", expt.c_str());
1395fputs("return(1); \n} \n", fip);
1396fclose(fip);
1397string func = "expf_pia_dl_func";
1398return((PlotExprFunc)LinkFunctionFromFile(fname, func));
1399}
1400
1401
1402/* --Methode-- */
1403DlFunction Services2NObjMgr::LinkFunctionFromFile(string const & fname, string const & funcname)
1404{
1405// Le link dynamique
1406CloseDLL();
1407dynlink = PDynLinkMgr::BuildFromCFile(fname);
1408if (dynlink == NULL) {
1409 cerr << "Services2NObjMgr/LinkFunctionFromFile_Erreur: Erreur creation/Ouverture SO " << endl;
1410 return(NULL);
1411 }
1412
1413DlFunction retfunc = dynlink->GetFunction(funcname);
1414if (retfunc == NULL) {
1415 string sn = funcname;
1416 cerr << "Services2NObjMgr/LinkExprFunc_Erreur: Erreur linking " << sn << endl;
1417 CloseDLL();
1418 return(NULL);
1419 }
1420else return(retfunc);
1421}
1422
1423/* --Methode-- */
1424void Services2NObjMgr::CloseDLL()
1425{
1426if (dynlink) delete dynlink; dynlink = NULL;
1427}
1428
[1319]1429/* --Methode-- */
1430int Services2NObjMgr::ExecuteCommand(string line)
1431{
1432 if (mImgapp == NULL) return(99);
1433 return(mImgapp->CmdInterpreter()->Interpret(line));
1434}
1435
[357]1436// Fonction static
1437/* --Methode-- */
[2419]1438void Services2NObjMgr::DecodeLoopParameters(string& loop, int_8& i1, int_8& i2, int_8& di)
[357]1439{
1440// Decode des paramatres de boucle for(int i=i1; i<i2; i+=di) specifies
1441// sous forme i1[:i2[:di]]
1442// cout << "LoopParam() " << loop << " I1=" << i1 << " I2=" << i2 << " DI=" << di;
1443size_t l = loop.length();
1444if (l < 1) return;
1445size_t p = loop.find(':');
[2419]1446if (p >= l) { i1 = atol(loop.c_str()); return; }
1447i1 = atol(loop.substr(0, p).c_str());
[357]1448string aa = loop.substr(p+1);
1449p = aa.find(':');
1450if (p < aa.length() ) {
[2419]1451 i2 = atol(aa.substr(0,p).c_str());
1452 di = atol(aa.substr(p+1).c_str());
[357]1453 }
[2419]1454else i2 = atol(aa.c_str());
[357]1455// cout << "-> I1= " << i1 << " I2= " << i2 << " DI= " << di << endl;
1456return;
1457}
[333]1458
1459/* --Methode-- */
[165]1460string Services2NObjMgr::FileName2Name(string const & fn)
1461{
1462
1463char fsep[2] = {FILESEP, '\0'};
1464char tsep[2] = {'.', '\0'};
1465size_t p = fn.find_last_of(fsep);
[194]1466size_t l = fn.length();
1467if (p >= l) p = 0;
1468else p++;
1469size_t q = fn.find_first_of(tsep,p);
1470if (q < p) q = l;
[165]1471return(fn.substr(p,q-p));
1472}
1473
1474
[546]1475
[165]1476
[333]1477// SANS_EVOLPLANCK Attention !
[1105]1478#ifdef SANS_EVOLPLANCK
[333]1479#include "pclassids.h"
1480
1481/* --Methode-- */
1482char* Services2NObjMgr::PClassIdToClassName(int cid)
1483{
1484switch (cid) {
1485 case ClassId_Poly1 :
1486 return("Poly1");
1487 case ClassId_Poly2 :
1488 return("Poly2");
1489 case ClassId_Matrix :
1490 return("Matrix");
1491 case ClassId_Vector :
1492 return("Vector");
1493
1494 case ClassId_DVList :
1495 return("DVList");
1496
1497 case ClassId_Histo1D :
1498 return("Histo1D");
1499 case ClassId_Histo2D :
1500 return("Histo2D");
1501 case ClassId_HProf :
1502 return("HProf");
[2605]1503 case ClassId_HistoErr :
1504 return("HistoErr");
[333]1505 case ClassId_NTuple :
1506 return("NTuple");
[361]1507 case ClassId_XNTuple :
1508 return("XNTuple");
[333]1509 case ClassId_GeneralFitData :
1510 return("GeneralFitData");
1511
1512 case ClassId_Image :
1513 return("RzImage");
1514 case ClassId_Image + kuint_1 :
1515 return("ImageU1");
1516 case ClassId_Image + kint_1 :
1517 return("ImageI1");
1518 case ClassId_Image + kuint_2 :
1519 return("ImageU2");
1520 case ClassId_Image + kint_2 :
1521 return("ImageI2");
1522 case ClassId_Image + kuint_4 :
1523 return("ImageU4");
1524 case ClassId_Image + kint_4 :
1525 return("ImageI4");
1526 case ClassId_Image + kr_4 :
1527 return("ImageR4");
1528 case ClassId_Image + kr_8 :
1529 return("ImageR8");
1530 case ClassId_ZFidu :
1531 return("ZFidu");
1532
1533 case ClassId_StarList :
1534 return("StarList");
1535 case ClassId_Transfo :
1536 return("Transfo");
1537 case ClassId_PSF :
1538 return("PSF");
1539
1540
1541// - Ajout objet PPF
1542 default:
1543 return("AnyDataObj");
1544 }
1545}
1546
[1105]1547#else
1548char* Services2NObjMgr::PClassIdToClassName(int cid)
1549{
1550 return("AnyDataObj");
1551}
1552#endif
Note: See TracBrowser for help on using the repository browser.