Changeset 333 in Sophya for trunk/SophyaPI/PIext/servnobjm.cc
- Timestamp:
- Jul 12, 1999, 1:12:29 PM (26 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaPI/PIext/servnobjm.cc
r331 r333 21 21 #include "pistdimgapp.h" 22 22 23 #include "fct1dfit.h" 24 #include "fct2dfit.h" 25 26 #include "matrix.h" 27 #include "cvector.h" 28 #include "ntuple.h" 29 #include "cimage.h" 30 23 31 #include "histos.h" 24 32 #include "histos2.h" … … 28 36 #include "piscdrawwdg.h" 29 37 #include "pisurfdr.h" 38 39 #include "pintuple.h" 40 #include "pintup3d.h" 41 30 42 #include "pipodrw.h" 31 43 … … 33 45 34 46 /* --Methode-- */ 35 Services2NObjMgr::Services2NObjMgr( PIStdImgApp* app, string& tmpdir)47 Services2NObjMgr::Services2NObjMgr(NamedObjMgr* omg, string& tmpdir) 36 48 { 37 49 TmpDir = tmpdir; 38 50 PDynLinkMgr::SetTmpDir(tmpdir); 39 mImgapp = app; 51 mImgapp = NULL; 52 mOmg = omg; 40 53 dynlink = NULL; 41 54 InitGrAttNames(); … … 71 84 72 85 /* --Methode-- */ 73 void Services2NObjMgr::Nobj_ComputeExpressions(NObjMgrAdapter* obja, string& expx, string& expy, string& expz, 74 string& expwt, string& expcut, NTuple* nt, Histo* h1, Histo2D* h2, HProf* hp) 86 void Services2NObjMgr::PlotFunc(string const & expfunc, string & nom, double xmin, double xmax, int np, string dopt) 87 { 88 FILE *fip; 89 string fname = TmpDir + "func1_pia_dl.c"; 90 string cmd; 91 int rc; 92 93 if (!mImgapp) return; 94 95 cmd = "rm -f " + fname; 96 rc = system(cmd.c_str()); 97 // printf("PlotFunc_Do> %s (Rc=%d)\n", cmd.c_str(), rc); 98 99 if ((fip = fopen(fname.c_str(), "w")) == NULL) { 100 string sn = fname; 101 cout << "Services2NObjMgr/PlotFunc_Erreur: Pb. Ouverture " << sn << endl; 102 return; 103 } 104 105 // constitution du fichier a compiler 106 fputs("#include <math.h> \n", fip); 107 fputs("double func1_pia_dl_func(double x) \n{\n", fip); 108 fprintf(fip,"return(%s); \n}\n", expfunc.c_str()); 109 fclose(fip); 110 111 string func = "func1_pia_dl_func"; 112 DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func); 113 if (!f) return; 114 PlotFunc(f, nom, xmin, xmax, np, dopt); 115 CloseDLL(); 116 return; 117 } 118 119 /* --Methode-- */ 120 void Services2NObjMgr::PlotFunc2D(string const & expfunc, string & nom, double xmin, double xmax, 121 double ymin, double ymax, int npx, int npy, string dopt) 122 { 123 FILE *fip; 124 string fname = TmpDir + "func2_pia_dl.c"; 125 string cmd; 126 int rc; 127 128 if (!mImgapp) return; 129 130 cmd = "rm " + fname; 131 rc = system(cmd.c_str()); 132 // printf("PlotFunc2D_Do> %s (Rc=%d)\n", cmd.c_str(), rc); 133 134 if ((fip = fopen(fname.c_str(), "w")) == NULL) { 135 string sn = fname; 136 cout << "Services2NObjMgr/PlotFunc2D_Erreur: Pb. Ouverture " << sn << endl; 137 return; 138 } 139 140 // constitution du fichier a compiler 141 fputs("#include <math.h> \n", fip); 142 fputs("double func2_pia_dl_func(double x, double y) \n{\n", fip); 143 fprintf(fip,"return(%s); \n}\n", expfunc.c_str()); 144 fclose(fip); 145 146 string func = "func2_pia_dl_func"; 147 DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func); 148 if (!f) return; 149 PlotFunc2D(f, nom, xmin, xmax, ymin, ymax, npx, npy, dopt); 150 CloseDLL(); 151 return; 152 } 153 154 /* --Methode-- */ 155 void Services2NObjMgr::PlotFuncFrCFile(string const & fname, string const & func, string & nom, 156 double xmin, double xmax, int np, string dopt) 157 { 158 DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func); 159 if (!f) return; 160 PlotFunc(f, nom, xmin, xmax, np, dopt); 161 CloseDLL(); 162 return; 163 } 164 165 /* --Methode-- */ 166 void Services2NObjMgr::PlotFunc2DFrCFile(string const & fname, string const & func, string & nom, 167 double xmin, double xmax, double ymin, double ymax, int npx, int npy, string dopt) 168 { 169 DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func); 170 if (!f) return; 171 PlotFunc2D(f, nom, xmin, xmax, ymin, ymax, npx, npy, dopt); 172 CloseDLL(); 173 return; 174 } 175 176 /* --Methode-- */ 177 void Services2NObjMgr::PlotFunc(DlFunctionOfX f, string & nom, double xmin, double xmax, int np, string dopt) 178 { 179 if (!mImgapp) return; 180 181 int k; 182 if (np < 1) np = 1; 183 if (xmax <= xmin) xmax = xmin+1.; 184 Vector* vpy = new Vector(np); 185 186 double xx; 187 double dx = (xmax-xmin)/np; 188 TRY { 189 for(k=0; k<np; k++) { xx = xmin+dx*k; (*vpy)(k) = f(xx); } 190 } CATCH(merr) { 191 fflush(stdout); 192 cout << endl; 193 cerr << endl; 194 string es = PeidaExc(merr); 195 cerr << "Services2NObjMgr::PlotFunc() Exception :" << merr << es; 196 delete vpy; 197 vpy = NULL; 198 } ENDTRY; 199 200 201 if (vpy) { 202 string titre = (nom.length() < 1) ? "Function f(x)" : nom; 203 P1DArrayAdapter* vya = new POVectorAdapter(vpy, true); 204 vya->DefineXCoordinate(xmin, (xmax-xmin)/np); 205 PIYfXDrawer* dr = new PIYfXDrawer(vya, NULL, true) ; 206 bool fgsr = true; 207 dopt = "thinline," + dopt; 208 int opt = DecodeDispOption(dopt, fgsr); 209 int rsid = mImgapp->DispScDrawer(dr, titre, opt); 210 if (fgsr) mImgapp->RestoreGraphicAtt(); 211 if (nom.length() > 0) { 212 mOmg->AddObj(vpy, nom); 213 mOmg->AddWRsId(nom, rsid); 214 } 215 } 216 217 return; 218 } 219 220 /* --Methode-- */ 221 void Services2NObjMgr::PlotFunc2D(DlFunctionOfXY f, string & nom, double xmin, double xmax, double ymin, double ymax, 222 int npx, int npy, string dopt) 223 { 224 if (!mImgapp) return; 225 226 if (npx < 3) npx = 3; 227 if (npy < 3) npy = 3; 228 if (npx > 250) npx = 250; 229 if (npy > 250) npy = 250; 230 if (xmax <= xmin) xmax = xmin+1.; 231 if (ymax <= ymin) ymax = ymin+1.; 232 233 Matrix* mtx = new Matrix(npy, npx); 234 235 int i,j; 236 double xx, yy; 237 double dx = (xmax-xmin)/npx; 238 double dy = (ymax-ymin)/npy; 239 // printf(" -- DBG -- %d %d , %g %g , %g %g \n", npx, npy, xmin, xmax, ymin, ymax); 240 TRY { 241 for(j=0; j<npy; j++) { 242 yy = ymin+dy*j; 243 for(i=0; i<npx; i++) { 244 xx = xmin+dx*i; 245 (*mtx)(j, i) = f(xx, yy); 246 } 247 } 248 } CATCH(merr) { 249 fflush(stdout); 250 cout << endl; 251 cerr << endl; 252 string es = PeidaExc(merr); 253 cerr << "Services2NObjMgr::PlotFunc2D() Exception :" << merr << es; 254 delete mtx; mtx = NULL; 255 } ENDTRY; 256 257 if (mtx) { 258 string titre = (nom.length() < 1) ? "Function f(x,y)" : nom; 259 bool fgsr = true; 260 int opt = DecodeDispOption(dopt, fgsr); 261 P2DArrayAdapter* arr = new POMatrixAdapter(mtx, true); 262 arr->DefineXYCoordinates(xmin, ymin, dx, dy); 263 PISurfaceDrawer* sdr = new PISurfaceDrawer(arr, true, true, true); 264 int rsid = mImgapp->Disp3DDrawer(sdr, titre, opt); 265 if (fgsr) mImgapp->RestoreGraphicAtt(); 266 if (nom.length() > 0) { 267 mOmg->AddObj(mtx, nom); 268 mOmg->AddWRsId(nom, rsid); 269 } 270 } 271 272 return; 273 } 274 275 /* --Methode-- */ 276 void Services2NObjMgr::DisplayPoints2D(string& nom, string& expx, string& expy, 277 string& experrx, string& experry, 278 string& expcut, string dopt) 279 { 280 NObjMgrAdapter* obja=NULL; 281 obja = mOmg->GetObjAdapter(nom); 282 if (obja == NULL) { 283 cout << "Services2NObjMgr::DisplayPoints2D() Error , Pas d'objet de nom " << nom << endl; 284 return; 285 } 286 if (!mImgapp) return; 287 288 // Creation NTuple 289 char* ntn[4] = {"expx","expy","expex","expey",}; 290 NTuple* nt = NULL; 291 bool haserr = false; 292 293 if ( (experrx.length() > 0 ) && (experry.length() > 0 ) ) { haserr = true; nt = new NTuple(4, ntn); } 294 else { haserr = false; experrx = experry = "0."; nt = new NTuple(2, ntn); } 295 296 ComputeExpressions(obja, expx, expy, experrx, experry, expcut, nt, NULL, NULL); 297 298 if (nt->NEntry() < 1) { 299 cout << "Services2NObjMgr::DisplayPoints2D() Warning Zero points satisfy cut !" << endl; 300 delete nt; 301 return; 302 } 303 304 // nt->Show(); 305 // nt->Print(0,10); 306 PINTuple* pin = new PINTuple(nt, true); 307 pin->SelectXY(ntn[0], ntn[1]); 308 if ( haserr ) pin->SelectErrBar(ntn[2], ntn[3]); 309 310 bool fgsr = true; 311 dopt = "defline," + dopt; 312 int opt = DecodeDispOption(dopt, fgsr); 313 string titre = nom + ":" + expy + "%" + expx; 314 mImgapp->DispScDrawer( (PIDrawer*)pin, titre, opt); 315 if (fgsr) mImgapp->RestoreGraphicAtt(); 316 return; 317 } 318 319 /* --Methode-- */ 320 void Services2NObjMgr::DisplayPoints3D(string& nom, string& expx, string& expy, 321 string& expz, string& expcut, string dopt) 322 { 323 NObjMgrAdapter* obja=NULL; 324 obja = mOmg->GetObjAdapter(nom); 325 if (obja == NULL) { 326 cout << "Services2NObjMgr::DisplayPoints3D() Error , Pas d'objet de nom " << nom << endl; 327 return; 328 } 329 if (!mImgapp) return; 330 331 char* ntn[3] = {"expx","expy","expz"}; 332 NTuple* nt = new NTuple(3,ntn); // Creation NTuple 333 334 string expwt = "1."; 335 ComputeExpressions(obja, expx, expy, expz, expwt, expcut, nt, NULL, NULL); 336 337 if (nt->NEntry() < 1) { 338 cout << "Services2NObjMgr::DisplayPoints3D() Warning Zero points satisfy cut !" << endl; 339 delete nt; 340 return; 341 } 342 nt->Show(); 343 nt->Print(0,10); 344 PINTuple3D* pin = new PINTuple3D(nt, true); 345 pin->SelectXYZ(ntn[0], ntn[1], ntn[2]); 346 bool fgsr = true; 347 dopt = "defline," + dopt; 348 int opt = DecodeDispOption(dopt, fgsr); 349 350 // Pour plot a partir de DispScDrawer 351 // string nomdisp = "_NT3D_"; 352 // mImgapp->DispScDrawer( (PIDrawer*)pin, nomdisp, opt); 353 // Pour plot a partir de Disp3DDrawer 354 string titre = nom + ":" + expy + "%" + expx; 355 mImgapp->Disp3DDrawer(pin, titre, opt); 356 357 if (fgsr) mImgapp->RestoreGraphicAtt(); 358 return; 359 } 360 361 /* --Methode-- */ 362 void Services2NObjMgr::DisplayPoints2DW(string& nom, string& expx, string& expy, 363 string& expwt, string& expcut, string dopt) 364 { 365 NObjMgrAdapter* obja=NULL; 366 obja = mOmg->GetObjAdapter(nom); 367 if (obja == NULL) { 368 cout << "Services2NObjMgr::DisplayPoints2DW() Error , Pas d'objet de nom " << nom << endl; 369 return; 370 } 371 if (!mImgapp) return; 372 373 char* ntn[3] = {"expx","expy","expw"}; 374 NTuple* nt = new NTuple(3,ntn); // Creation NTuple 375 376 string exp = "1."; 377 ComputeExpressions(obja, expx, expy, expwt, exp, expcut, nt, NULL, NULL); 378 379 if (nt->NEntry() < 1) { 380 cout << "Services2NObjMgr::DisplayPoints2DW() Warning Zero points satisfy cut !" << endl; 381 delete nt; 382 return; 383 } 384 385 PINTuple* pin = new PINTuple(nt, true); 386 pin->SelectXY(ntn[0], ntn[1]); 387 pin->SelectWt(ntn[2]); 388 389 bool fgsr = true; 390 int opt = DecodeDispOption(dopt, fgsr); 391 string titre = nom + ":" + expwt + "_" + expy + "%" + expx ; 392 mImgapp->DispScDrawer( (PIDrawer*)pin, titre, opt); 393 if (fgsr) mImgapp->RestoreGraphicAtt(); 394 return; 395 } 396 397 /* --Methode-- */ 398 void Services2NObjMgr::ProjectH1(string& nom, string& expx, string& expwt, 399 string& expcut, string& nomh1, string dopt) 400 { 401 NObjMgrAdapter* obja=NULL; 402 obja = mOmg->GetObjAdapter(nom); 403 if (obja == NULL) { 404 cout << "Services2NObjMgr::ProjectH1() Error , Pas d'objet de nom " << nom << endl; 405 return; 406 } 407 if (!mImgapp) return; 408 409 Histo* h1 = NULL; 410 NTuple* nt = NULL; 411 AnyDataObj* oh = NULL; 412 if (nomh1.length() > 0) oh=mOmg->GetObj(nomh1); 413 else nomh1 = "/tmp/projh1d"; 414 if ( (oh != NULL) && (typeid(*oh) == typeid(Histo)) ) h1 = (Histo*)oh; // Pas de remise a zero ! h1->Zero(); 415 else { 416 char* ntn[2]= {"hxval", "hwt"}; 417 nt = new NTuple(2,ntn); // Creation NTuple 418 } 419 string expz = "0."; 420 ComputeExpressions(obja, expx, expwt, expz, expwt, expcut, nt, h1, NULL); 421 422 if ((!h1) && (!nt)) return; 423 if (!h1) { 424 if (nt->NEntry() < 1) { 425 cout << "Services2NObjMgr::ProjectH1() Warning Zero points satisfy cut !" << endl; 426 delete nt; 427 return; 428 } 429 double xmin, xmax; 430 nt->GetMinMax(0, xmin, xmax); 431 h1 = new Histo(xmin, xmax, 100); 432 int k; 433 float* xn; 434 for(k=0; k<nt->NEntry(); k++) { 435 xn = nt->GetVec(k); 436 h1->Add(xn[0], xn[1]); 437 } 438 delete nt; 439 mOmg->AddObj(h1, nomh1); 440 } 441 442 mOmg->DisplayObj(nomh1, dopt); 443 return; 444 } 445 446 /* --Methode-- */ 447 void Services2NObjMgr::ProjectH2(string& nom, string& expx, string& expy, string& expwt, 448 string& expcut, string& nomh2, string dopt) 449 { 450 NObjMgrAdapter* obja=NULL; 451 obja = mOmg->GetObjAdapter(nom); 452 if (obja == NULL) { 453 cout << "Services2NObjMgr::ProjectH2() Error , Pas d'objet de nom " << nom << endl; 454 return; 455 } 456 if (!mImgapp) return; 457 458 Histo2D* h2 = NULL; 459 NTuple* nt = NULL; 460 AnyDataObj* oh = NULL; 461 if (nomh2.length() > 0) oh=mOmg->GetObj(nomh2); 462 else nomh2 = "/tmp/projh2d"; 463 if ( (oh != NULL) && (typeid(*oh) == typeid(Histo2D)) ) h2 = (Histo2D*)oh; // Pas de remise a zero ! h2->Zero(); 464 else { 465 char* ntn[3]= {"hxval", "hyval", "hwt"}; 466 nt = new NTuple(3,ntn); // Creation NTuple 467 } 468 string expz = "0."; 469 ComputeExpressions(obja, expx, expy, expwt, expwt, expcut, nt, NULL, h2); 470 471 if ((!h2) && (!nt)) return; 472 if (!h2) { 473 if (nt->NEntry() < 1) { 474 cout << "Services2NObjMgr::ProjectH2() Warning Zero points satisfy cut !" << endl; 475 delete nt; 476 return; 477 } 478 double xmin, xmax, ymin, ymax; 479 nt->GetMinMax(0, xmin, xmax); 480 nt->GetMinMax(0, ymin, ymax); 481 h2 = new Histo2D(xmin, xmax, 50, ymin, ymax, 50); 482 int k; 483 float* xn; 484 for(k=0; k<nt->NEntry(); k++) { 485 xn = nt->GetVec(k); 486 h2->Add(xn[0], xn[1], xn[2]); 487 } 488 delete nt; 489 mOmg->AddObj(h2, nomh2); 490 } 491 492 mOmg->DisplayObj(nomh2, dopt); 493 return; 494 495 } 496 497 /* --Methode-- cmv 13/10/98 */ 498 void Services2NObjMgr::ProjectHProf(string& nom, string& expx, string& expy, string& expwt, 499 string& expcut, string& nomprof, string dopt) 500 // Pour remplir un ``GeneralFitData'' a partir de divers objets: 501 //| nom = nom de l'objet a projeter dans un HProf. 502 //| expx = expression X de definition du bin. 503 //| expy = expression Y a additionner dans le bin. 504 //| expwt = expression W du poids a additionner. 505 //| expcut = expression du test de selection. 506 //| nomprof = nom du HProf engendre (optionnel). Si l'objet n'existe pas 507 //| les limites Xmin,Xmax sont calculees automatiquement. 508 //| sinon ce sont celles de l'objet preexistant. 509 //| opt = options generales pour le display. 510 { 511 NObjMgrAdapter* obja=NULL; 512 obja = mOmg->GetObjAdapter(nom); 513 if (obja == NULL) { 514 cout << "Services2NObjMgr::ProjectHProf() Error , Pas d'objet de nom " << nom << endl; 515 return; 516 } 517 if (!mImgapp) return; 518 519 HProf* hprof = NULL; 520 NTuple* nt = NULL; 521 AnyDataObj* oh = NULL; 522 if (nomprof.length() > 0) oh=mOmg->GetObj(nomprof); 523 else nomprof = "/tmp/projprof"; 524 if( (oh!=NULL) && (typeid(*oh) == typeid(HProf)) ) hprof = (HProf*)oh; 525 else { 526 char* ntn[3]= {"hxval", "hyval", "hwt"}; 527 nt = new NTuple(3,ntn); // Creation NTuple 528 } 529 string expz = "0."; 530 ComputeExpressions(obja, expx, expy, expwt, expwt, expcut, nt, NULL, NULL, hprof); 531 532 if((!hprof) && (!nt)) return; 533 if(!hprof) { 534 if (nt->NEntry() < 1) { 535 cout << "Services2NObjMgr::ProjectHProf() Warning Zero points satisfy cut !" << endl; 536 delete nt; 537 return; 538 } 539 double xmin, xmax; 540 nt->GetMinMax(0, xmin, xmax); 541 hprof = new HProf(xmin, xmax, 100); 542 int k; 543 float* xn; 544 for(k=0; k<nt->NEntry(); k++) { 545 xn = nt->GetVec(k); 546 hprof->Add(xn[0], xn[1], xn[2]); 547 } 548 delete nt; 549 mOmg->AddObj(hprof, nomprof); 550 } 551 hprof->UpdateHisto(); 552 553 mOmg->DisplayObj(nomprof, dopt); 554 return; 555 } 556 557 /* --Methode-- */ 558 void Services2NObjMgr::FillVect(string& nom, string& expx, string& expcut, 559 string& nomvec, string dopt) 560 { 561 NObjMgrAdapter* obja=NULL; 562 obja = mOmg->GetObjAdapter(nom); 563 if (obja == NULL) { 564 cout << "Services2NObjMgr::FillVect() Error , Pas d'objet de nom " << nom << endl; 565 return; 566 } 567 if (!mImgapp) return; 568 569 NTuple* nt = NULL; 570 if (nomvec.length() < 1) nomvec = "/tmp/fillvec"; 571 572 char* ntn[2]= {"vecval", "vecwt"}; 573 nt = new NTuple(1,ntn); // Creation NTuple 574 575 string expwt = "1."; 576 string expz = "0."; 577 ComputeExpressions(obja, expx, expz, expz, expwt, expcut, nt, NULL, NULL); 578 579 if (!nt) return; 580 if (nt->NEntry() < 1) { 581 cout << "Services2NObjMgr::FillVect() Warning Zero points satisfy cut !" << endl; 582 delete nt; 583 return; 584 } 585 586 Vector* vec = new Vector(nt->NEntry()); 587 int k; 588 float* xn; 589 for(k=0; k<nt->NEntry(); k++) { 590 xn = nt->GetVec(k); 591 (*vec)(k) = xn[0]; 592 } 593 delete nt; 594 mOmg->AddObj(vec, nomvec); 595 mOmg->DisplayObj(nomvec, dopt); 596 return; 597 } 598 599 /* --Methode-- */ 600 void Services2NObjMgr::FillNT(string& nom, string& expx, string& expy, string& expz, 601 string& expt, string& expcut, string& nomnt) 602 { 603 NObjMgrAdapter* obja=NULL; 604 obja = mOmg->GetObjAdapter(nom); 605 if (obja == NULL) { 606 cout << "Services2NObjMgr::FillNT() Error , Pas d'objet de nom " << nom << endl; 607 return; 608 } 609 if (!mImgapp) return; 610 611 bool fgnnt = false; 612 NTuple* nt = NULL; 613 AnyDataObj* oh = NULL; 614 if (nomnt.length() > 0) oh=mOmg->GetObj(nomnt); 615 else nomnt = "/tmp/fillnt"; 616 if ( (oh != NULL) && (typeid(*oh) == typeid(NTuple)) ) { 617 nt = (NTuple*)oh; 618 if (nt->NVar() > 10) { 619 cout << "Services2NObjMgr::FillNT() Warning , Max 10 var ds NTuple -> new NTuple" << endl; 620 nt = NULL; 621 } 622 } 623 if (nt == NULL) { 624 char* ntn[4]= {"x", "y","z","t"}; 625 nt = new NTuple(4,ntn); // Creation NTuple 626 fgnnt = true; 627 } 628 629 ComputeExpressions(obja, expx, expy, expz, expt, expcut, nt, NULL, NULL); 630 631 if (fgnnt) mOmg->AddObj(nt, nomnt); 632 return; 633 634 } 635 636 /* --Methode-- */ 637 void Services2NObjMgr::FillNTFrCFile(string & nom, string const & fname, 638 string const & funcname, string & nomnt) 639 { 640 if (!mImgapp) return; 641 642 NObjMgrAdapter* obja=NULL; 643 obja = mOmg->GetObjAdapter(nom); 644 if (obja == NULL) { 645 cout << "Services2NObjMgr::FillNTFrCFile( " << nom << "...) No such object" <<endl; 646 return; 647 } 648 NTupleInterface* objnt = obja->GetNTupleInterface(); 649 if (objnt == NULL) { 650 cout << "Services2NObjMgr::FillNTFrCFile( " << nom << "...) No NTupleInterface !" <<endl; 651 return; 652 } 653 654 NTLoopExprFunc f = (NTLoopExprFunc)LinkFunctionFromFile(fname, funcname); 655 if (!f) { 656 cerr << "Services2NObjMgr::FillNTFrCFile Error Creation NTLoopExprFunc" << endl; 657 return; 658 } 659 660 bool fgnnt = false; 661 NTuple* nt = NULL; 662 if (nomnt.length() > 0) { 663 AnyDataObj* oh = NULL; 664 oh=mOmg->GetObj(nomnt); 665 if ( (oh != NULL) && (typeid(*oh) == typeid(NTuple)) ) { 666 nt = (NTuple*)oh; 667 if (nt->NVar() > 10) { 668 cout << "Services2NObjMgr::FillNTFrCFile() Warning , Max 10 var ds NTuple -> new NTuple" << endl; 669 nt = NULL; 670 } 671 } 672 if (nt == NULL) { 673 char* ntn[4]= {"x", "y","z","t"}; 674 nt = new NTuple(4,ntn); // Creation NTuple 675 fgnnt = true; 676 } 677 } 678 679 double xnt[10]; 680 float fxnt[10]; 681 682 int i,k; 683 for(i=0; i<10; i++) fxnt[i] = xnt[i] = 0.; 684 685 686 // $CHECK$ A virer des que possible - Pb blocage application quand trop d'impression 687 // redirige - On redirige la sortie sur le terminal 688 bool red = mImgapp->HasRedirectedStdOutErr(); 689 mImgapp->RedirectStdOutErr(false); 690 691 TRY { 692 double* xn; 693 int kmax = objnt->NbLines(); 694 for(k=0; k<kmax; k++) { 695 xn = objnt->GetLineD(k); 696 if (f(xn, xnt, xnt+1, xnt+2, xnt+3, k, kmax) != 0) { 697 if (nt) { 698 for(i=0; i<4; i++) fxnt[i] = xnt[i]; 699 nt->Fill(fxnt); 700 } 701 } 702 } 703 } 704 CATCH(merr) { 705 fflush(stdout); 706 cout << endl; 707 cerr << endl; 708 string es = PeidaExc(merr); 709 cerr << "Services2NObjMgr::FillNTFrCFile() Exception :" << merr << es; 710 } ENDTRY; 711 712 CloseDLL(); 713 714 // $CHECK$ A virer des que possible On redirige la sortie sur la fenetre PIConsole 715 mImgapp->RedirectStdOutErr(red); 716 717 if (fgnnt) mOmg->AddObj(nt, nomnt); 718 return; 719 } 720 721 /* --Methode-- */ 722 void Services2NObjMgr::PrepareNTExpressionCFile(string & nom, string const & fname, 723 string const & funcname) 724 { 725 NObjMgrAdapter* obja=NULL; 726 obja = mOmg->GetObjAdapter(nom); 727 if (obja == NULL) { 728 cout << "Services2NObjMgr::PrepareNTExpressionCFile( " << nom << "...) No such object" <<endl; 729 return; 730 } 731 NTupleInterface* objnt = obja->GetNTupleInterface(); 732 if (objnt == NULL) { 733 cout << "Services2NObjMgr::PrepareNTExpressionCFile( " << nom 734 << "...) No NTupleInterface !" <<endl; 735 return; 736 } 737 string vardec = objnt->VarList_C("_xnti_"); 738 739 FILE *fip; 740 if ((fip = fopen(fname.c_str(), "w")) == NULL) { 741 cout << "Services2NObjMgr::PrepareNTExpressionCFile()_Error: fopen " << fname << endl; 742 return; 743 } 744 // constitution du fichier des decalarations des variables de l'interface NTuple 745 fputs("#include <stdlib.h> \n", fip); 746 fputs("#include <stdio.h> \n", fip); 747 fputs("#include <math.h> \n\n", fip); 748 fputs("/* NTupleInterface Variable declaration - Generated by piapp \n", fip); 749 fputs(" -- Services2NObjMgr::PrepareNTExpressionCFile() -- */ \n\n", fip); 750 fprintf(fip,"int %s(double* _xnti_, double* _rx_, double* _ry_, double* _rz_, \n", 751 funcname.c_str()); 752 fprintf(fip," double* _rt_, int _n_, int _nmax_) \n"); 753 fprintf(fip, "{ \n %s \n", vardec.c_str()); 754 fputs(" if (!1) { /* Cut Expression failed */ \n", fip); 755 fputs(" *_rx_ = *_ry_ = *_rz_ = *_rt_ = 0.; return(0);", fip); 756 fputs(" } \n /* Cut expression satisfied */ \n", fip); 757 fputs(" *_rx_ = 1.; \n *_ry_ = 1.; \n *_rz_ = 1.; \n *_rt_ = 1.; \n", fip); 758 fputs(" return(1); \n} \n", fip); 759 760 fclose(fip); 761 return; 762 } 763 764 765 /* --Methode-- cmv 13/10/98 */ 766 void Services2NObjMgr::FillGFD(string& nom, string& expx, string& expy, string& expz, 767 string& experr, string& expcut, string& nomgfd) 768 // Pour remplir un ``GeneralFitData'' a partir de divers objets: 769 //| nom = nom de l'objet a transcrire selon 1D: Z=f(X) ou 2D: Z=f(X,Y) . 770 //| Vector,Matrix,Histo,HProf,Histo2D,Image<T>,StarList,NTuple,GeneralFitData 771 //| expx = expression X du GeneralFitData (1er abscisse) 772 //| expy = expression Y du GeneralFitData (2sd abscisse si non "", Z=f(X,Y)) 773 //| expz = expression Z du GeneralFitData (valeur de l'ordonnee) 774 //| experr = expression de l'erreur sur l'ordonnee Z 775 //| expcut = expression du test de selection 776 //| nomgfd = nom du GeneralFitData engendre (optionnel) 777 { 778 NObjMgrAdapter* obja=NULL; 779 obja = mOmg->GetObjAdapter(nom); 780 if (obja == NULL) { 781 cout << "Services2NObjMgr::FillGFD() Error , Pas d'objet de nom "<<nom<<endl; 782 return; 783 } 784 if(!mImgapp) return; 785 786 // 2D ou 3D? 787 int nvar = 2; 788 if(expy.length()<=0) {nvar = 1; expy = "0.";} 789 790 // $CHECK$ - cmv calculait le nombre d'entree ndata 791 // en fonction de chaque objet - Je l'ai vire Reza 11/05/99 792 793 // Creation NTuple Buffer 794 char* ntn[4]= {"x","y","f","e"}; 795 NTuple*nt = new NTuple(4,ntn); 796 797 // Remplissage NTuple buffer 798 ComputeExpressions(obja, expx, expy, expz, experr, expcut, nt, NULL, NULL); 799 if(nt->NEntry() < 1) 800 {cout<<"Services2NObjMgr::FillGFD() Warning Zero points satisfy cut !"<<endl; 801 delete nt; return;} 802 803 //Remplissage de la structure GeneraFitData 804 if (nt->NEntry() <= 0) { 805 cout<<"Services2NObjMgr::FillGFD() Warning - NData= " << nt->NEntry() << endl; 806 delete nt; 807 return; 808 } 809 810 GeneralFitData* gfd = new GeneralFitData(nvar,nt->NEntry(),0); 811 int k; 812 float* xn; 813 for(k=0; k<nt->NEntry(); k++) { 814 xn = nt->GetVec(k); 815 gfd->AddData(xn,xn[2],xn[3]); 816 } 817 818 // Menage et table d'objets 819 delete nt; 820 mOmg->AddObj(gfd, nomgfd); 821 return; 822 } 823 824 825 ///////////////////// Fit 1D et 2D ////////////////////////// 826 /* --Function static propre aux routines de fit 1D et 2D-- cmv 13/10/98 */ 827 struct DFOptions { 828 int okres, okfun; 829 int polcx,polcy; double xc,yc; 830 double err_e, err_E; 831 double stc2; 832 int lp,lpg; 833 int i1,i2,j1,j2; 834 }; 835 typedef struct DFOptions DFOPTIONS; 836 static void DecodeFitsOptions(string par,string step,string min,string max,string opt 837 ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O); 838 void DecodeFitsOptions(string par,string step,string min,string max,string opt 839 ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O) 840 //| Pour decoder les "string" et remplir les vecteurs du fit (cf commentaires dans Fit1D) 841 { 842 // set des vecteurs et decodage des string correspondantes 843 int NParMax = 100; 844 Par.Realloc(NParMax); Step.Realloc(NParMax); 845 Min.Realloc(NParMax); Max.Realloc(NParMax); 846 { 847 Vector* v=NULL; string* s=NULL; 848 {for(int i=0;i<NParMax;i++) {Par(i)=0.; Step(i)=1.; Min(i)=1.; Max(i)=-1.;}} 849 for(int j=0;j<4;j++) { 850 if(j==0) {v=&Par; s=∥} 851 else if(j==1) {v=&Step; s=&step;} 852 else if(j==2) {v=&Min; s=&min;} 853 else if(j==3) {v=&Max; s=&max;} 854 if(s->length()>0) *s += ","; 855 for(int i=0;i<NParMax;i++) { 856 if(s->length()<=0) break; 857 sscanf(s->c_str(),"%lf",&(*v)(i)); 858 size_t p = s->find_first_of(',') + 1; 859 if(p>=s->length()) *s = ""; else *s = s->substr(p); 860 } 861 } 862 } 863 864 // Decodage de options de opt 865 O.okres = O.okfun = 0; 866 O.polcx = O.polcy = 0; 867 O.xc = O.yc = 0.; 868 O.stc2 = 1.e-3; 869 O.err_e = O.err_E = -1.; 870 O.lp = 1; O.lpg = 0; 871 O.i1 = O.j1 = O.i2 = O.j2 = -1; 872 873 if(opt.length()<=0) return; 874 opt = "," + opt + ","; 875 876 if(strstr(opt.c_str(),",r,")) O.okres = 1; // residus 877 if(strstr(opt.c_str(),",f,")) O.okfun = 1; // fonction fittee 878 if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc) 879 O.polcx = 2; // Le centrage est calcule automatiquement 880 size_t p = opt.find(",x"); 881 size_t q = opt.find_first_of(',',p+1); 882 string dum = opt.substr(p,q-p); 883 if(dum.length()>2) { 884 sscanf(dum.c_str(),",x%lf",&O.xc); 885 O.polcx = 1; // Le centrage est fixe par la valeur lue 886 } 887 } 888 if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc) 889 O.polcy = 2; // Le centrage est calcule automatiquement 890 size_t p = opt.find(",y"); 891 size_t q = opt.find_first_of(',',p+1); 892 string dum = opt.substr(p,q-p); 893 if(dum.length()>2) { 894 sscanf(dum.c_str(),",y%lf",&O.yc); 895 O.polcy = 1; // Le centrage est fixe par la valeur lue 896 } 897 } 898 if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)" 899 size_t p = opt.find(",E"); 900 size_t q = opt.find_first_of(',',p+1); 901 string dum = opt.substr(p,q-p); 902 if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&O.err_E); 903 if(O.err_E<=0.) O.err_E = 1.; 904 O.err_e=-1.; 905 } 906 if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b" 907 size_t p = opt.find(",e"); 908 size_t q = opt.find_first_of(',',p+1); 909 string dum = opt.substr(p,q-p); 910 if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&O.err_e); 911 if(O.err_e<=0.) O.err_e = 1.; 912 O.err_E=-1.; 913 } 914 if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2 915 size_t p = opt.find(",X"); 916 size_t q = opt.find_first_of(',',p+1); 917 string dum = opt.substr(p,q-p); 918 if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&O.stc2); 919 if(O.stc2<=0.) O.stc2 = 1.e-3; 920 } 921 if(strstr(opt.c_str(),",l")) { // niveau de print 922 size_t p = opt.find(",l"); 923 size_t q = opt.find_first_of(',',p+1); 924 string dum = opt.substr(p,q-p); 925 float ab; 926 if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab); 927 if(ab<0) ab = 0.; 928 O.lp = (int) ab; O.lpg = (int) 10.*(ab-O.lp); 929 } 930 if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X 931 size_t p = opt.find(",I"); 932 size_t q = opt.find_first_of(',',p+1); 933 string dum = opt.substr(p,q-p); 934 if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&O.i1,&O.i2); 935 } 936 if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y 937 size_t p = opt.find(",J"); 938 size_t q = opt.find_first_of(',',p+1); 939 string dum = opt.substr(p,q-p); 940 if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&O.j1,&O.j2); 941 } 942 return; 943 } 944 945 /* --Methode-- cmv 13/10/98 */ 946 void Services2NObjMgr::Fit12D(string& nom, string& func, 947 string par,string step,string min,string max, 948 string opt) 949 //| --------------- Fit d'objets a 1 et 2 dimensions --------------- 950 //| nom : nom de l'objet qui peut etre: 951 //| fit-1D: Vector,Histo1D,HProf ou GeneraFitData(1D) 952 //| fit-2D: Matrix,Histo2D,Imagexx ou GeneraFitData(2D) 953 //| func : pnn = fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D 954 //| : Pnn = fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D 955 //| : gnn = fit gaussienne (hauteur) + polynome de degre nn 1D 956 //| : g = fit gaussienne (hauteur) 1D 957 //| : enn = fit exponentielle + polynome de degre nn 1D 958 //| : e = fit exponentielle 1D 959 //| : Gnn = fit gaussienne (volume) + polynome de degre nn 1D 960 //| : G = fit gaussienne (volume) 1D 961 //| : = fit gaussienne+fond (volume) 2D 962 //| : Gi = fit gaussienne+fond integree (volume) 2D 963 //| : d = fit DL de gaussienne+fond (volume) 2D 964 //| : di = fit DL de gaussienne+fond integree (volume) 2D 965 //| : D = fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D 966 //| : Di = fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D 967 //| : M = fit Moffat+fond (expos=p6) (volume) 2D 968 //| : Mi = fit Moffat+fond integree (expos=p6) (volume) 2D 969 //| par : p1,...,pn valeur d'initialisation des parametres (def=0) 970 //| step : s1,...,sn valeur des steps de depart (def=1) 971 //| min : m1,...,mn valeur des minima (def=1) 972 //| max : M1,...,Mn valeur des maxima (def=-1) (max<=min : pas de limite) 973 //| opt : options "Eaa.b,eaa.b,f,r,caa.b,Xaa.b" 974 //| f = generation d'un Objet identique contenant la fonction fittee 975 //| r = generation d'un Objet identique contenant les residus 976 //| Xaa.b = aa.b valeur du DXi2 d'arret (def=1.e-3) 977 //| la.b = niveau "a.b" de print: a=niveau de print Fit1/2D 978 //| b=niveau de debug GeneralFit 979 //| Ii1/i2 numeros des bins X de l'histos utilises pour le fit [i1,i2] 980 //|2D Jj1/j2 numeros des bins Y de l'histos utilises pour le fit [j1,j2] 981 //| - L'erreur est celle associee a l'objet si existe, 1 sinon sauf si 982 //| E = erreur prise comme la racine de la valeur a fitter 983 //| Eaa.b = erreur prise aa.b*sqrt(val) 984 //| e = erreur prise egale a 1 pour toutes les valeurs 985 //| eaa.b = erreur prise egale a aa.b 986 //| xaa.b = demande de centrage: on fit x-aa.b au lieu de x) 987 //| x = demande de centrage: on fit x-xc au lieu de x 988 //| avec xc=abscisse du milieu de l'histogramme 989 //| Actif pour: exp+poly 1D, poly 1D 990 //| gauss+poly 1D (mais xc est le centre de la gaussienne) 991 //|2D yaa.b et y = idem "xaa.b et x" mais pour y 992 { 993 AnyDataObj* obj=mOmg->GetObj(nom); 994 if (obj == NULL) { 995 cout<<"Services2NObjMgr::Fit12D() Error , Pas d'objet de nom "<<nom<<endl; 996 return; 997 } 998 if (!mImgapp) return; 999 if(func.length()<=0) 1000 {cout<<"Services2NObjMgr::Fit12D() Donnez un nom de fonction a fitter."<<endl; 1001 return;} 1002 string ctyp = typeid(*obj).name(); 1003 1004 int ndim = 0, nbinx=0, nbiny=0, ndata = 0; 1005 Vector* v = NULL; Histo* h = NULL; 1006 Matrix* m = NULL; Histo2D* h2 = NULL; RzImage* im = NULL; 1007 GeneralFitData* g = NULL; 1008 1009 // 1D 1010 if (typeid(*obj) == typeid(Vector)) { 1011 ndim = 1; 1012 v = (Vector*) obj; nbinx = v->NElts(); nbiny = 1; 1013 } 1014 else if ( (typeid(*obj) == typeid(HProf)) || (typeid(*obj) == typeid(Histo)) ) { 1015 ndim = 1; 1016 h = (Histo*) obj; nbinx = h->NBins(); nbiny = 1; 1017 } 1018 else if (typeid(*obj) == typeid(Matrix)) { 1019 ndim = 2; 1020 m = (Matrix*) obj; nbinx = m->NCol(); nbiny = m->NRows(); 1021 } 1022 else if (typeid(*obj) == typeid(Histo2D)) { 1023 ndim = 2; 1024 h2 = (Histo2D*) obj; nbinx = h2->NBinX(); nbiny = h2->NBinY(); 1025 } 1026 else if (typeid(*obj) == typeid(GeneralFitData)) { 1027 g = (GeneralFitData*) obj; nbinx = g->NData(); nbiny = 1; 1028 if( g->NVar()==1) ndim = 1; 1029 else if(g->NVar()==2) ndim = 2; 1030 else { 1031 cout<<"GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: " 1032 <<((GeneralFitData*) obj)->NVar()<<endl; return; } 1033 } 1034 else if (dynamic_cast<RzImage*>(obj)) { 1035 ndim = 2; 1036 im = (RzImage*) obj; nbinx = im->XSize(); nbiny = im->YSize(); 1037 } 1038 else { 1039 cout<<"Services2NObjMgr::Fit12D() Error , Objet n'est pas un " 1040 <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl; 1041 return; 1042 } 1043 1044 ndata = nbinx*nbiny; 1045 if(ndata<=0) 1046 {cout<<"L'objet a "<<nbinx<<","<<nbiny<<" bins ("<<ndata<<")"<<endl; return;} 1047 1048 // Decodage des options et des parametres, mise en forme 1049 Vector Par(1); Vector Step(1); Vector Min(1); Vector Max(1); DFOPTIONS O; 1050 DecodeFitsOptions(par,step,min,max,opt,Par,Step,Min,Max,O); 1051 O.i1 = (O.i1<0||O.i1>=nbinx)? 0: O.i1; 1052 O.i2 = (O.i2<0||O.i2>=nbinx||O.i2<O.i1)? nbinx-1: O.i2; 1053 if(ndim>=2) { 1054 O.j1 = (O.j1<0||O.j1>=nbiny)? 0: O.j1; 1055 O.j2 = (O.j2<0||O.j2>=nbiny||O.j2<O.j1)? nbiny-1: O.j2; 1056 } else O.j2 = O.j1 = 0; 1057 if(O.polcx==2) { 1058 if(v||m) O.xc = (O.i2-O.i1+1)/2.; 1059 else if(h) O.xc = (h->XMin()+h->XMax())/2.; 1060 else if(h2) O.xc = (h2->XMin()+h2->XMax())/2.; 1061 else if(g) {double mini,maxi; g->GetMinMax(2,mini,maxi); O.xc=(mini+maxi)/2.;} 1062 else if(im) {O.xc = im->XOrg() * im->XPxSize()*(O.i2-O.i1+1)/2.;} 1063 } 1064 if(O.polcy==2 && ndim>=2) { 1065 if(m) O.yc = (O.j2-O.j1+1)/2.; 1066 if(h2) O.yc = (h2->YMin()+h2->YMax())/2.; 1067 if(g) {double mini,maxi; g->GetMinMax(12,mini,maxi); O.yc=(mini+maxi)/2.;} 1068 if(im) {O.yc = im->YOrg() * im->YPxSize()*(O.j2-O.j1+1)/2.;} 1069 } 1070 if(O.lp>0) 1071 cout<<"Fit["<<nbinx<<","<<nbiny<<"] ("<<ndata<<") dim="<<ndim<<":" 1072 <<" Int=["<<O.i1<<","<<O.i2<<"],["<<O.j1<<","<<O.j2<<"]"<<endl 1073 <<" Cent="<<O.polcx<<","<<O.polcy<<","<<O.xc<<"+x"<<","<<O.yc<<"+y" 1074 <<" TypE="<<O.err_e<<","<<O.err_E 1075 <<" StpX2="<<O.stc2 1076 <<" lp,lpg="<<O.lp<<","<<O.lpg<<endl; 1077 1078 /////////////////////////////////// 1079 // Remplissage de GeneralFitData // 1080 /////////////////////////////////// 1081 GeneralFitData mydata(ndim,ndata,0); 1082 {for(int i=O.i1;i<=O.i2;i++) for(int j=O.j1;j<=O.j2;j++) { 1083 double x,y,f,e; 1084 1085 if(v) 1086 {x= (double) i; f=(*v)(i); e=1.;} 1087 else if(h) 1088 {x=h->BinCenter(i); f=(*h)(i); e=(h->HasErrors())?h->Error(i):1.;} 1089 else if(m) 1090 {x=(double) i; y=(double) j; f=(*m)(j,i); e=1.;} 1091 else if(h2) 1092 {float xf,yf; h2->BinCenter(i,j,xf,yf); x=(double)xf; y=(double)yf; 1093 f=(*h2)(i,j); e=(h2->HasErrors())?h2->Error(i,j):1.;} 1094 else if(im) 1095 {x=im->XOrg()+(i+0.5)*im->XPxSize(); y=im->YOrg()+(j+0.5)*im->YPxSize(); 1096 f=im->DValue(i,j); e=1.;} 1097 else if(g&&ndim==1) {x= g->X(i); f=g->Val(i); e=g->EVal(i);} 1098 else if(g&&ndim==2) {x= g->X(i); y= g->Y(i); f=g->Val(i); e=g->EVal(i);} 1099 else x=y=f=e=0.; 1100 1101 // Gestion des erreurs a utiliser 1102 if(O.err_e>0.) e=O.err_e; 1103 else if(O.err_E>0.) {e=(y<-1.||y>1.)?O.err_E*sqrt(fabs(y)):O.err_E;} 1104 1105 // Remplissage de generalfit 1106 if(func[0]=='p') {x -= O.xc; if(ndim>=2) y -= O.yc;} 1107 if(ndim==1) mydata.AddData1(x,f,e); 1108 else if(ndim==2) mydata.AddData2(x,y,f,e); 1109 }} 1110 if(mydata.NData()<=0) 1111 {cout<<"Pas de donnees dans GeneralFitData: "<<mydata.NData()<<endl; 1112 return;} 1113 if(O.lpg>1) { 1114 mydata.PrintStatus(); 1115 mydata.PrintData(0); 1116 mydata.PrintData(mydata.NData()-1); 1117 } 1118 1119 //////////////////////////////////////////// 1120 // Identification de la fonction a fitter // 1121 //////////////////////////////////////////// 1122 GeneralFunction* myfunc = NULL; 1123 if(func[0]=='p' && ndim==1) { 1124 // Fit de polynome sans passer par les GeneralFit 1125 int degre = 0; 1126 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1127 cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl; 1128 Poly p1(0); 1129 double c2rl = mydata.PolFit(0,p1,degre); 1130 cout<<"C2r_lineaire = "<<c2rl<<endl; 1131 if(O.lp>0) cout<<p1<<endl; 1132 return; 1133 1134 } else if(func[0]=='P' && ndim==1) { 1135 // Fit de polynome 1136 int degre = 0; 1137 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1138 cout<<"Fit polynome 1D de degre "<<degre<<endl; 1139 Polyn1D* myf = new Polyn1D(degre,O.xc); 1140 myfunc = myf; 1141 1142 } else if(func[0]=='e' && ndim==1) { 1143 // Fit d'exponentielle 1144 int degre =-1; 1145 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1146 cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl; 1147 Exp1DPol* myf; 1148 if(degre>=0) myf = new Exp1DPol((unsigned int)degre,O.xc); 1149 else myf = new Exp1DPol(O.xc); 1150 myfunc = myf; 1151 1152 } else if(func[0]=='g' && ndim==1) { 1153 // Fit de gaussienne en hauteur 1154 int degre =-1; 1155 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1156 cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl; 1157 Gauss1DPol* myf; 1158 if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((O.polcx)?true:false)); 1159 else { bool bfg = (O.polcx)?true:false; myf = new Gauss1DPol(bfg); } 1160 myfunc = myf; 1161 1162 } else if(func[0]=='G' && ndim==1) { 1163 // Fit de gaussienne en volume 1164 int degre =-1; 1165 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1166 cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl; 1167 GaussN1DPol* myf; 1168 if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((O.polcx)?true:false)); 1169 else { bool bfg = (O.polcx)?true:false; myf = new GaussN1DPol(bfg); } 1170 myfunc = myf; 1171 1172 } else if(func[0]=='p' && ndim==2) { 1173 // Fit de polynome 2D sans passer par les GeneralFit 1174 int degre = 0; 1175 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1176 cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl; 1177 Poly2 p2(0); 1178 double c2rl = mydata.PolFit(0,1,p2,degre); 1179 cout<<"C2r_lineaire = "<<c2rl<<endl; 1180 if(O.lp>0) cout<<p2<<endl; 1181 return; 1182 1183 } else if(func[0]=='P' && ndim==2) { 1184 // Fit de polynome 2D 1185 int degre = 0; 1186 if(func.length()>1) sscanf(func.c_str()+1,"%d",°re); 1187 cout<<"Fit polynome 2D de degre "<<degre<<endl; 1188 Polyn2D* myf = new Polyn2D(degre,O.xc,O.yc); 1189 myfunc = myf; 1190 1191 } else if(func[0]=='G' && ndim==2) { 1192 // Fit de gaussienne+fond en volume 1193 int integ = 0; 1194 if(func.length()>1) if(func[1]=='i') integ=1; 1195 cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl; 1196 if(integ) {GauRhInt2D* myf = new GauRhInt2D; myfunc = myf;} 1197 else {GauRho2D* myf = new GauRho2D; myfunc = myf;} 1198 1199 } else if(func[0]=='d' && ndim==2) { 1200 // Fit de DL gaussienne+fond en volume 1201 int integ = 0; 1202 if(func.length()>1) if(func[1]=='i') integ=1; 1203 cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl; 1204 if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; myfunc = myf;} 1205 else {GdlRho2D* myf = new GdlRho2D; myfunc = myf;} 1206 1207 } else if(func[0]=='D' && ndim==2) { 1208 // Fit de DL gaussienne+fond avec coeff variable p6 en volume 1209 int integ = 0; 1210 if(func.length()>1) if(func[1]=='i') integ=1; 1211 cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl; 1212 if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; myfunc = myf;} 1213 else {Gdl1Rho2D* myf = new Gdl1Rho2D; myfunc = myf;} 1214 1215 } else if(func[0]=='M' && ndim==2) { 1216 // Fit de Moffat+fond (volume) 1217 int integ = 0; 1218 if(func.length()>1) if(func[1]=='i') integ=1; 1219 cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl; 1220 if(integ) {MofRhInt2D* myf = new MofRhInt2D; myfunc = myf;} 1221 else {MofRho2D* myf = new MofRho2D; myfunc = myf;} 1222 1223 } else { 1224 cout<<"Fonction "<<func<<" inconnue pour la dim "<<ndim<<endl; 1225 return; 1226 } 1227 1228 ///////////////////////// 1229 // Fit avec generalfit // 1230 ///////////////////////// 1231 if(myfunc->NPar()>Par.NElts()) 1232 {cout<<"Trop de parametres: "<<myfunc->NPar()<<">"<<Par.NElts()<<endl; 1233 if(myfunc) delete myfunc; return;} 1234 GeneralFit myfit(myfunc); 1235 myfit.SetDebug(O.lpg); 1236 myfit.SetData(&mydata); 1237 myfit.SetStopChi2(O.stc2); 1238 {for(int i=0;i<myfunc->NPar();i++) { 1239 char str[10]; 1240 sprintf(str,"P%d",i); 1241 myfit.SetParam(i,str,Par(i),Step(i),Min(i),Max(i)); 1242 }} 1243 if(O.lp>1) myfit.PrintFit(); 1244 double c2r = (double) myfit.Fit(); 1245 if(O.lp>0) myfit.PrintFit(); 1246 if(c2r>0.) { 1247 c2r = myfit.GetChi2Red(); 1248 cout<<"C2r_Reduit = "<<c2r<<" nstep="<<myfit.GetNStep()<<endl; 1249 Vector ParFit(myfunc->NPar()); 1250 for(int i=0;i<myfunc->NPar();i++) ParFit(i)=myfit.GetParm(i); 1251 } else { 1252 cout<<"echec Fit, rc = "<<c2r<<" nstep="<<myfit.GetNStep()<<endl; 1253 } 1254 1255 // Mise a disposition des resultats 1256 if(c2r>=0. && myfunc && (O.okres>0||O.okfun>0)) { 1257 string nomres = nom + "res"; 1258 string nomfun = nom + "fun"; 1259 if(v) { 1260 if(O.okres) {Vector* ob = v->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1261 if(O.okfun) {Vector* ob = v->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1262 } else if(h) { 1263 if(O.okres) {Histo* ob = h->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1264 if(O.okfun) {Histo* ob = h->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1265 } else if(m) { 1266 if(O.okres) {Matrix* ob = m->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1267 if(O.okfun) {Matrix* ob = m->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1268 } else if(h2) { 1269 if(O.okres) {Histo2D* ob = h2->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1270 if(O.okfun) {Histo2D* ob = h2->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1271 } else if(im) { 1272 if(O.okres) {RzImage* ob = im->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1273 if(O.okfun) {RzImage* ob = im->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1274 } else if(g) { 1275 if(O.okres) {GeneralFitData* ob = g->FitResidus(myfit); if(ob) mOmg->AddObj(ob,nomres);} 1276 if(O.okfun) {GeneralFitData* ob = g->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);} 1277 } 1278 } 1279 1280 // Nettoyage 1281 if(myfunc) delete myfunc; 1282 return; 1283 } 1284 1285 1286 /* --Methode-- */ 1287 void Services2NObjMgr::ComputeExpressions(NObjMgrAdapter* obja, string& expx, 1288 string& expy, string& expz, string& expt, string& expcut, 1289 NTuple* nt, Histo* h1, Histo2D* h2, HProf* hp) 75 1290 { 76 1291 if (obja == NULL) return; 77 1292 NTupleInterface* objnt = obja->GetNTupleInterface(); 78 1293 if (objnt == NULL) return; 79 string vardec = objnt->VarList_C("_zz6 1qq_");80 81 PlotExprFunc f = LinkExprFunc(vardec, expx, expy, expz, exp wt, expcut);1294 string vardec = objnt->VarList_C("_zz6qi_"); 1295 1296 PlotExprFunc f = LinkExprFunc(vardec, expx, expy, expz, expt, expcut); 82 1297 if (!f) { 83 cerr << " NamedObjMgr_ComputeExpressions() Error Creation PlotExprFunc " << endl;1298 cerr << "Services2NObjMgr::::ComputeExpressions() Error Creation PlotExprFunc " << endl; 84 1299 return; 85 1300 } … … 88 1303 float fxnt[10]; 89 1304 90 int i, j,k;1305 int i,k; 91 1306 for(i=0; i<10; i++) xnt[i] = 0.; 92 1307 … … 111 1326 cerr << endl; 112 1327 string es = PeidaExc(merr); 113 cerr << " NamedObjMgr_ComputeExpressions() Exception :" << merr << es;1328 cerr << "Services2NObjMgr::ComputeExpressions() Exception :" << merr << es; 114 1329 } ENDTRY; 115 1330 … … 122 1337 123 1338 /* --Methode-- */ 124 PlotExprFunc Services2NObjMgr::LinkExprFunc(string& vardec, string& expx, string& expy, string& expz,125 string& exp wt, string& cut)1339 PlotExprFunc Services2NObjMgr::LinkExprFunc(string& vardec, string& expx, string& expy, 1340 string& expz, string& expt, string& cut) 126 1341 { 127 1342 FILE *fip; … … 136 1351 if ((fip = fopen(fname.c_str(), "w")) == NULL) { 137 1352 string sn = fname; 138 cout << " NamedObjMgr/LinkExprFunc_Erreur: Pb. Ouverture " << sn << endl;1353 cout << "Services2NObjMgr/LinkExprFunc_Erreur: Pb. Ouverture " << sn << endl; 139 1354 return(NULL); 140 1355 } 141 1356 142 1357 // constitution du fichier a compiler 1358 fputs("#include <stdlib.h> \n", fip); 143 1359 fputs("#include <math.h> \n", fip); 144 fputs("int expf_pia_dl_func(double* _zz6 1qq_, double* _rx_61qq_, double* _ry_61qq_, double* _rz_61qq_, double* _wt_61qq_) \n{\n", fip);1360 fputs("int expf_pia_dl_func(double* _zz6qi_, double* _rx_6q_, double* _ry_6q_, double* _rz_6q_, double* _rt_6q_) \n{\n", fip); 145 1361 fprintf(fip,"%s \n", vardec.c_str()); 146 fprintf(fip, "if (!(%s)) { *_rx_6 1qq_ = *_ry_61qq_ = *_rz_61qq_ = *_wt_61qq_ = 0.; return(0); } \n", cut.c_str());147 fprintf(fip, "*_rx_6 1qq_ = %s ; \n", expx.c_str());148 fprintf(fip, "*_ry_6 1qq_ = %s ; \n", expy.c_str());149 fprintf(fip, "*_rz_6 1qq_ = %s ; \n", expz.c_str());150 fprintf(fip, "*_ wt_61qq_ = %s ; \n", expwt.c_str());1362 fprintf(fip, "if (!(%s)) { *_rx_6q_ = *_ry_6q_ = *_rz_6q_ = *_rt_6q_ = 0.; return(0); } \n", cut.c_str()); 1363 fprintf(fip, "*_rx_6q_ = %s ; \n", expx.c_str()); 1364 fprintf(fip, "*_ry_6q_ = %s ; \n", expy.c_str()); 1365 fprintf(fip, "*_rz_6q_ = %s ; \n", expz.c_str()); 1366 fprintf(fip, "*_rt_6q_ = %s ; \n", expt.c_str()); 151 1367 fputs("return(1); \n} \n", fip); 152 1368 fclose(fip); … … 163 1379 dynlink = PDynLinkMgr::BuildFromCFile(fname); 164 1380 if (dynlink == NULL) { 165 cerr << " NamedObjMgr/LinkFunctionFromFile_Erreur: Erreur creation/Ouverture SO " << endl;1381 cerr << "Services2NObjMgr/LinkFunctionFromFile_Erreur: Erreur creation/Ouverture SO " << endl; 166 1382 return(NULL); 167 1383 } … … 170 1386 if (retfunc == NULL) { 171 1387 string sn = funcname; 172 cerr << " NamedObjMgr/LinkExprFunc_Erreur: Erreur linking " << sn << endl;1388 cerr << "Services2NObjMgr/LinkExprFunc_Erreur: Erreur linking " << sn << endl; 173 1389 CloseDLL(); 174 1390 return(NULL); … … 183 1399 } 184 1400 185 /* --Methode-- */186 void Services2NObjMgr::PlotFunc(string const & expfunc, float xmin, float xmax, int np, string dopt)187 {188 FILE *fip;189 string fname = TmpDir + "func1_pia_dl.c";190 string cmd;191 int rc;192 193 if (!mImgapp) return;194 195 cmd = "rm -f " + fname;196 rc = system(cmd.c_str());197 // printf("PlotFunc_Do> %s (Rc=%d)\n", cmd.c_str(), rc);198 199 if ((fip = fopen(fname.c_str(), "w")) == NULL) {200 string sn = fname;201 cout << "NamedObjMgr/PlotFunc_Erreur: Pb. Ouverture " << sn << endl;202 return;203 }204 205 // constitution du fichier a compiler206 fputs("#include <math.h> \n", fip);207 fputs("double func1_pia_dl_func(double x) \n{\n", fip);208 fprintf(fip,"return(%s); \n}\n", expfunc.c_str());209 fclose(fip);210 211 string func = "func1_pia_dl_func";212 DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func);213 if (!f) return;214 PlotFunc(f, xmin, xmax, np, dopt);215 CloseDLL();216 return;217 }218 219 /* --Methode-- */220 void Services2NObjMgr::PlotFunc2D(string const & expfunc, float xmin, float xmax, float ymin, float ymax,221 int npx, int npy, string dopt)222 {223 FILE *fip;224 string fname = TmpDir + "func2_pia_dl.c";225 string cmd;226 int rc;227 228 if (!mImgapp) return;229 230 cmd = "rm " + fname;231 rc = system(cmd.c_str());232 // printf("PlotFunc2D_Do> %s (Rc=%d)\n", cmd.c_str(), rc);233 234 if ((fip = fopen(fname.c_str(), "w")) == NULL) {235 string sn = fname;236 cout << "NamedObjMgr/PlotFunc2D_Erreur: Pb. Ouverture " << sn << endl;237 return;238 }239 240 // constitution du fichier a compiler241 fputs("#include <math.h> \n", fip);242 fputs("double func2_pia_dl_func(double x, double y) \n{\n", fip);243 fprintf(fip,"return(%s); \n}\n", expfunc.c_str());244 fclose(fip);245 246 string func = "func2_pia_dl_func";247 DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func);248 if (!f) return;249 PlotFunc2D(f, xmin, xmax, ymin, ymax, npx, npy, dopt);250 CloseDLL();251 return;252 }253 254 /* --Methode-- */255 void Services2NObjMgr::PlotFuncFrCFile(string const & fname, string const & func, float xmin, float xmax,256 int np, string dopt)257 {258 DlFunctionOfX f = (DlFunctionOfX) LinkFunctionFromFile(fname, func);259 if (!f) return;260 PlotFunc(f, xmin, xmax, np, dopt);261 CloseDLL();262 return;263 }264 265 /* --Methode-- */266 void Services2NObjMgr::PlotFunc2DFrCFile(string const & fname, string const & func, float xmin, float xmax,267 float ymin, float ymax, int npx, int npy, string dopt)268 {269 DlFunctionOfXY f = (DlFunctionOfXY) LinkFunctionFromFile(fname, func);270 if (!f) return;271 PlotFunc2D(f, xmin, xmax, ymin, ymax, npx, npy, dopt);272 CloseDLL();273 return;274 }275 276 /* --Methode-- */277 void Services2NObjMgr::PlotFunc(DlFunctionOfX f, float xmin, float xmax, int np, string dopt)278 {279 if (!mImgapp) return;280 281 int k;282 if (np < 1) np = 1;283 if (xmax <= xmin) xmax = xmin+1.;284 Vector* vpy = new Vector(np);285 286 double xx;287 double dx = (xmax-xmin)/np;288 TRY {289 for(k=0; k<np; k++) { xx = xmin+dx*k; (*vpy)(k) = f(xx); }290 } CATCH(merr) {291 fflush(stdout);292 cout << endl;293 cerr << endl;294 string es = PeidaExc(merr);295 cerr << "Services2NObjMgr::PlotFunc() Exception :" << merr << es;296 delete vpy;297 vpy = NULL;298 } ENDTRY;299 300 301 if (vpy) {302 string nom = "Func";303 P1DArrayAdapter* vya = new POVectorAdapter(vpy, true);304 vya->DefineXCoordinate(xmin, (xmax-xmin)/np);305 PIYfXDrawer* dr = new PIYfXDrawer(vya, NULL, true) ;306 bool fgsr = true;307 dopt = "thinline," + dopt;308 int opt = DecodeDispOption(dopt, fgsr);309 int wrsid = mImgapp->DispScDrawer(dr, nom, opt);310 if (fgsr) mImgapp->RestoreGraphicAtt();311 }312 313 return;314 }315 316 /* --Methode-- */317 void Services2NObjMgr::PlotFunc2D(DlFunctionOfXY f, float xmin, float xmax, float ymin, float ymax,318 int npx, int npy, string dopt)319 {320 if (!mImgapp) return;321 322 if (npx < 3) npx = 3;323 if (npy < 3) npy = 3;324 if (npx > 250) npx = 250;325 if (npy > 250) npy = 250;326 if (xmax <= xmin) xmax = xmin+1.;327 if (ymax <= ymin) ymax = ymin+1.;328 329 Matrix* mtx = new Matrix(npy, npx);330 331 int i,j;332 double xx, yy;333 double dx = (xmax-xmin)/npx;334 double dy = (ymax-ymin)/npy;335 // printf(" -- DBG -- %d %d , %g %g , %g %g \n", npx, npy, xmin, xmax, ymin, ymax);336 TRY {337 for(j=0; j<npy; j++) {338 yy = ymin+dy*j;339 for(i=0; i<npx; i++) {340 xx = xmin+dx*i;341 (*mtx)(j, i) = f(xx, yy);342 }343 }344 } CATCH(merr) {345 fflush(stdout);346 cout << endl;347 cerr << endl;348 string es = PeidaExc(merr);349 cerr << "Services2NObjMgr::PlotFunc2D() Exception :" << merr << es;350 delete mtx; mtx = NULL;351 } ENDTRY;352 353 if (mtx) {354 string nom = "Func2";355 int wrsid = 0;356 bool fgsr = true;357 int opt = DecodeDispOption(dopt, fgsr);358 P2DArrayAdapter* arr = new POMatrixAdapter(mtx, true);359 arr->DefineXYCoordinates(xmin, ymin, dx, dy);360 PISurfaceDrawer* sdr = new PISurfaceDrawer(arr, true, true, true);361 wrsid = mImgapp->Disp3DDrawer(sdr, nom, opt);362 if (fgsr) mImgapp->RestoreGraphicAtt();363 }364 365 return;366 }367 1401 368 1402 /* --Methode-- */ … … 663 1697 664 1698 } 1699 1700 1701 // SANS_EVOLPLANCK Attention ! 1702 #include "pclassids.h" 1703 1704 /* --Methode-- */ 1705 char* Services2NObjMgr::PClassIdToClassName(int cid) 1706 { 1707 switch (cid) { 1708 case ClassId_Poly1 : 1709 return("Poly1"); 1710 case ClassId_Poly2 : 1711 return("Poly2"); 1712 case ClassId_Matrix : 1713 return("Matrix"); 1714 case ClassId_Vector : 1715 return("Vector"); 1716 1717 case ClassId_DVList : 1718 return("DVList"); 1719 1720 case ClassId_Histo1D : 1721 return("Histo1D"); 1722 case ClassId_Histo2D : 1723 return("Histo2D"); 1724 case ClassId_HProf : 1725 return("HProf"); 1726 case ClassId_NTuple : 1727 return("NTuple"); 1728 case ClassId_GeneralFitData : 1729 return("GeneralFitData"); 1730 1731 case ClassId_Image : 1732 return("RzImage"); 1733 case ClassId_Image + kuint_1 : 1734 return("ImageU1"); 1735 case ClassId_Image + kint_1 : 1736 return("ImageI1"); 1737 case ClassId_Image + kuint_2 : 1738 return("ImageU2"); 1739 case ClassId_Image + kint_2 : 1740 return("ImageI2"); 1741 case ClassId_Image + kuint_4 : 1742 return("ImageU4"); 1743 case ClassId_Image + kint_4 : 1744 return("ImageI4"); 1745 case ClassId_Image + kr_4 : 1746 return("ImageR4"); 1747 case ClassId_Image + kr_8 : 1748 return("ImageR8"); 1749 1750 case ClassId_ZFidu : 1751 return("ZFidu"); 1752 1753 case ClassId_StarList : 1754 return("StarList"); 1755 case ClassId_Transfo : 1756 return("Transfo"); 1757 case ClassId_PSF : 1758 return("PSF"); 1759 1760 1761 // - Ajout objet PPF 1762 default: 1763 return("AnyDataObj"); 1764 } 1765 } 1766
Note:
See TracChangeset
for help on using the changeset viewer.