00001
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015
00016 #include "FunctionController.h"
00017
00018 #include "Gammaq.h"
00019 #include "DisplayController.h"
00020
00021 #include "datareps/CompositeFunctionRep.h"
00022 #include "datareps/FunctionParameter.h"
00023 #include "datareps/FunctionRep1.h"
00024 #include "datareps/FunctionRep2.h"
00025
00026 #include "datasrcs/DataSourceController.h"
00027 #include "datasrcs/NTuple.h"
00028 #include "datasrcs/TupleCut.h"
00029
00030 #include "functions/FunctionBase.h"
00031 #include "functions/FunctionFactory.h"
00032
00033 #include "minimizers/Fitter.h"
00034 #include "minimizers/FitterFactory.h"
00035 #include "minimizers/NumLinAlg.h"
00036
00037 #include "pattern/string_convert.h"
00038 #include "plotters/PlotterBase.h"
00039 #include "plotters/CompositePlotter.h"
00040 #include "projectors/BinningProjector.h"
00041 #include "projectors/NTupleProjector.h"
00042
00043 #include <algorithm>
00044 #include <functional>
00045 #include <iterator>
00046 #include <stdexcept>
00047
00048 #include <cmath>
00049 #include <cassert>
00050
00051 #ifdef ITERATOR_MEMBER_DEFECT
00052 using namespace std;
00053 #else
00054 using std::distance;
00055 using std::find;
00056 using std::find_if;
00057 using std::list;
00058 using std::mem_fun;
00059 using std::string;
00060 using std::vector;
00061 using std::sqrt;
00062 using std::cos;
00063 using std::sin;
00064 using std::atan;
00065 using std::min;
00066 using std::max;
00067 using std::abs;
00068 #endif
00069
00070 using namespace hippodraw;
00071
00072 using namespace Numeric;
00073
00074 FunctionController * FunctionController::s_instance = 0;
00075
00076 FunctionController::FunctionController ( )
00077 : m_plotter ( 0 ),
00078 m_x ( 0 ),
00079 m_y ( 0 ),
00080 m_confid_count ( 0 )
00081 {
00082 m_deltaXSq.resize( 6 );
00083
00084 m_deltaXSq[ 0 ] = 15.1;
00085 m_deltaXSq[ 1 ] = 18.4;
00086 m_deltaXSq[ 2 ] = 21.1;
00087 m_deltaXSq[ 3 ] = 23.5;
00088 m_deltaXSq[ 4 ] = 25.7;
00089 m_deltaXSq[ 5 ] = 27.8;
00090 }
00091
00092 FunctionController::~FunctionController ( )
00093 {
00094 }
00095
00096 FunctionController * FunctionController::instance ( )
00097 {
00098 if ( s_instance == 0 ) {
00099 s_instance = new FunctionController ( );
00100 }
00101 return s_instance;
00102 }
00103
00104 const vector < string > &
00105 FunctionController::
00106 getFitterNames () const
00107 {
00108 FitterFactory * factory = FitterFactory::instance ();
00109
00110 return factory -> names ();
00111 }
00112
00113 const std::string &
00114 FunctionController::
00115 getDefaultFitter () const
00116 {
00117 FitterFactory * factory = FitterFactory::instance ();
00118
00119 return factory -> getDefault ();
00120 }
00121
00122 int FunctionController::
00123 getUniqueNonFunctionIndex ( const PlotterBase * plotter ) const
00124 {
00125 int index = -1;
00126 int count = 0;
00127
00128 int number = plotter->getNumDataReps ();
00129 for ( int i = 0; i < number; i++ ) {
00130 DataRep * r = plotter->getDataRep ( i );
00131 FunctionRep * fr = dynamic_cast < FunctionRep * > ( r );
00132 if ( fr != 0 ) continue;
00133 index = i;
00134 count++;
00135 }
00136 if ( count != 1 ) index = -1;
00137
00138 return index;
00139 }
00140
00141 void
00142 FunctionController::
00143 findFunctions ( const PlotterBase * plotter ) const
00144 {
00145 m_func_reps.clear ();
00146 int number = plotter->getNumDataReps ();
00147
00148 for ( int i = 0; i < number; i++ ) {
00149 DataRep * rep = plotter->getDataRep ( i );
00150 FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00151 if ( frep == 0 ) continue;
00152 m_func_reps.push_back ( frep );
00153 }
00154 }
00155
00156 FunctionRep *
00157 FunctionController::
00158 createFunctionRep ( const std::string & name,
00159 DataRep * rep )
00160 {
00161 FunctionFactory * factory = FunctionFactory::instance ();
00162 FunctionBase * function = factory -> create ( name );
00163 if ( function == 0 ) {
00164 string what ( "FunctionController: Unable to create function of type `" );
00165 what += name;
00166 what += "'";
00167 throw std::runtime_error ( what );
00168 }
00169
00170 return createFunctionRep ( function, rep );
00171 }
00172
00173 FunctionRep *
00174 FunctionController::
00175 createFunctionRep ( FunctionBase * function,
00176 DataRep * rep )
00177 {
00178 FunctionRep * frep = 0;
00179
00180 if ( function -> isComposite () ) {
00181 frep = new CompositeFunctionRep ( function, rep );
00182 }
00183 else {
00184 unsigned int dims = function -> dimensions ();
00185 if ( dims == 2 ) {
00186 frep = new FunctionRep2 ( function, rep );
00187 }
00188 else {
00189 frep = new FunctionRep1 ( function, rep );
00190 }
00191 }
00192
00193 if ( rep != 0 ) {
00194 FitterFactory * factory = FitterFactory::instance ();
00195 const string & name = factory -> getDefault ();
00196 bool ok = setFitter ( frep, name );
00197 if ( ok == false ) {
00198 delete frep;
00199 frep = 0;
00200 }
00201 }
00202
00203 return frep;
00204 }
00205
00209 FunctionBase *
00210 FunctionController::
00211 addFunction ( PlotterBase * plotter, const std::string & name )
00212 {
00213 DataRep * rep = plotter -> getDataRep ( 0 );
00214 FunctionRep * frep = getFunctionRep ( plotter );
00215 addFunction ( plotter, name, frep, rep );
00216
00217 frep = getFunctionRep ( plotter );
00218
00219 return frep -> getFunction ();
00220 }
00221
00225 FunctionRep *
00226 FunctionController::
00227 addFunction ( PlotterBase * plotter,
00228 const std::string & name,
00229 FunctionRep * frep,
00230 DataRep * rep )
00231 {
00232 FunctionRep * func_rep = createFunctionRep ( name, rep );
00233 func_rep -> initializeWith ( rep );
00234
00235 return addFunctionRep ( plotter, rep, frep, func_rep );
00236 }
00237
00238 void
00239 FunctionController::
00240 addFunction ( PlotterBase * plotter,
00241 FunctionRep * func_rep )
00242 {
00243 int index = plotter->activePlotIndex ();
00244
00245 if ( index < 0 ) {
00246 index = getUniqueNonFunctionIndex ( plotter );
00247 }
00248 if ( index >= 0 ) {
00249 plotter->setActivePlot ( index, false );
00250 DataRep * drep = plotter->getDataRep ( index );
00251 func_rep -> initializeWith ( drep );
00252 fillFunctionReps ( m_func_reps, plotter, drep );
00253 FunctionRep * frep = 0;
00254 if ( m_func_reps.empty () == false ) {
00255 frep = m_func_reps.front ();
00256 }
00257 addFunctionRep ( plotter, drep, frep, func_rep );
00258 }
00259 }
00260
00261
00262 void
00263 FunctionController::
00264 addDataRep ( PlotterBase * plotter, DataRep * rep )
00265 {
00266 FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00267 if ( frep == 0 ) {
00268 DisplayController * controller = DisplayController::instance ();
00269 controller -> addDataRep ( plotter, rep );
00270 }
00271 else {
00272 if ( frep -> isComposite () == false ) {
00273 DataRep * target = frep -> getTarget ();
00274 addFunctionRep ( plotter, target, 0, frep );
00275 }
00276 }
00277 }
00278
00279 void
00280 FunctionController::
00281 setTupleCut ( FunctionRep * rep )
00282 {
00283 rep -> setTupleCut ();
00284 rep -> setCutRange ( true );
00285 }
00286
00287 void
00288 FunctionController::
00289 setTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
00290 {
00291 data_rep -> addCut ();
00292 FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
00293 if ( func_rep != 0 ) {
00294 setTupleCut ( func_rep );
00295 }
00296 }
00297
00298 void
00299 FunctionController::
00300 removeTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
00301 {
00302 FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
00303 if ( func_rep != 0 ) {
00304 func_rep -> removeCut ();
00305 }
00306 }
00307
00308 FunctionRep *
00309 FunctionController::
00310 addFunctionRep ( PlotterBase * plotter,
00311 DataRep * drep,
00312 FunctionRep * frep,
00313 FunctionRep * func_rep )
00314 {
00315 FunctionRep * return_rep = frep;
00316
00317 plotter->addDataRep ( func_rep );
00318
00319 fillFunctionReps ( m_func_reps, plotter, drep );
00320
00321 FitterFactory * factory = FitterFactory::instance ();
00322 const string & name = factory -> getDefault ();
00323 setFitter ( func_rep, name );
00324
00325 if ( frep != 0 ) {
00326 if ( frep -> isComposite () ) {
00327 frep -> addToComposite ( func_rep );
00328 }
00329 else {
00330 FunctionRep * composite = createFunctionRep ( "Linear Sum", drep );
00331 composite -> initializeWith ( drep );
00332 plotter -> addDataRep ( composite );
00333 setFitter ( composite, name );
00334 setTupleCut ( func_rep );
00335 composite -> addToComposite ( frep );
00336 composite -> addToComposite ( func_rep );
00337 return_rep = composite;
00338 }
00339 }
00340 func_rep->notifyObservers ();
00341 setTupleCut ( func_rep );
00342 if ( return_rep == 0 ) return_rep = func_rep;
00343
00344 return return_rep;
00345 }
00346
00347 void
00348 FunctionController::
00349 removeFunction ( PlotterBase * plotter,
00350 FunctionRep * frep )
00351 {
00352 if ( frep -> isComposite () ) {
00353 CompositeFunctionRep * composite
00354 = dynamic_cast < CompositeFunctionRep * > ( frep );
00355 const vector < FunctionRep * > & freps = composite -> getFunctionReps ();
00356 unsigned int size = freps.size();
00357
00358
00359 for ( int i = size -1; i >= 0; i-- ) {
00360 FunctionRep * rep = freps[i];
00361 plotter -> removeDataRep ( rep );
00362 composite -> removeFromComposite ( rep );
00363 delete rep;
00364 }
00365 plotter -> removeDataRep ( frep );
00366 delete frep;
00367 }
00368 else {
00369 fillFunctionReps ( m_func_reps, plotter, 0 );
00370 vector < FunctionRep * > :: iterator it
00371 = find ( m_func_reps.begin(), m_func_reps.end (),
00372 frep );
00373 assert ( it != m_func_reps.end () );
00374 m_func_reps.erase ( it );
00375
00376 plotter -> removeDataRep ( frep );
00377 for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
00378 FunctionRep * rep = m_func_reps[i];
00379 if ( rep -> isComposite () ) {
00380 rep -> removeFromComposite ( frep );
00381 }
00382 }
00383
00384
00385 vector < FunctionRep * >::iterator iter = m_func_reps.begin ();
00386
00387
00388
00389 while ( iter != m_func_reps.end() ){
00390 FunctionRep * rep = *iter;
00391 if ( rep -> isComposite () ) {
00392 CompositeFunctionRep * composite
00393 = dynamic_cast < CompositeFunctionRep * > ( rep );
00394 if ( composite -> count() < 2 ) {
00395 plotter->removeDataRep ( composite );
00396 const CompositeFunctionRep::FunctionRepList_t & reps
00397 = composite -> getFunctionReps ();
00398 FunctionRep * crep = reps.front ();
00399 crep -> setInComposite ( false );
00400 delete composite;
00401 iter = m_func_reps.erase ( iter );
00402 }
00403 else {
00404 iter ++;
00405 }
00406 }
00407 else {
00408 iter ++;
00409 }
00410 }
00411 delete frep;
00412 }
00413 }
00414
00415 bool
00416 FunctionController::
00417 hasFunction ( const PlotterBase * plotter, const DataRep * rep )
00418 {
00419 bool yes = false;
00420 assert ( plotter != 0 );
00421
00422 fillFunctionReps ( m_func_reps, plotter, rep );
00423 yes = m_func_reps.empty () == false;
00424
00425 return yes;
00426 }
00427
00428 const vector< string > &
00429 FunctionController::
00430 functionNames ( PlotterBase * plotter, DataRep * rep )
00431 {
00432 fillFunctionReps ( m_func_reps, plotter, rep );
00433
00434 m_f_names.clear ();
00435 vector< FunctionRep * >:: const_iterator it = m_func_reps.begin ();
00436
00437 for ( ; it != m_func_reps.end (); ++it ) {
00438 FunctionRep * fp = *it;
00439 if ( rep != 0 && fp->getTarget() != rep ) continue;
00440 FunctionBase * function = fp->getFunction ();
00441 const string & name = function->name ();
00442 m_f_names.push_back ( name );
00443 }
00444
00445 return m_f_names;
00446 }
00447
00448 void
00449 FunctionController::
00450 fillFunctionReps ( std::vector < FunctionRep * > & freps,
00451 const PlotterBase * plotter,
00452 const DataRep * drep ) const
00453 {
00454 freps.clear ();
00455 int number = plotter -> getNumDataReps ();
00456 for ( int i = 0; i < number; i++ ) {
00457 DataRep * rep = plotter -> getDataRep ( i );
00458 FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00459 if ( frep != 0 ) {
00460 if ( drep != 0 ) {
00461 if ( frep -> getTarget () == drep ) {
00462 freps.push_back ( frep );
00463 }
00464 }
00465 else {
00466 freps.push_back ( frep );
00467 }
00468 }
00469 }
00470 }
00471
00472 void
00473 FunctionController::
00474 fillTopLevelFunctionReps ( std::vector < FunctionRep * > & freps,
00475 const PlotterBase * plotter,
00476 const DataRep * drep ) const
00477 {
00478 freps.clear ();
00479 fillFunctionReps ( m_func_reps, plotter, drep );
00480 unsigned int size = m_func_reps.size ();
00481 for ( unsigned int i = 0; i < size; i++ ) {
00482 FunctionRep * rep = m_func_reps[i];
00483 if ( rep -> isInComposite () == false ) {
00484 freps.push_back ( rep );
00485 }
00486 }
00487 }
00488
00489 void
00490 FunctionController::
00491 setErrorsFromComposite ( const PlotterBase * plotter,
00492 const FunctionRep * composite )
00493 {
00494 const vector < double > & errors = composite -> principleErrors();
00495 vector < double >::const_iterator begin = errors.begin();
00496
00497 DataRep * target = composite -> getTarget ();
00498 fillFunctionReps ( m_func_reps, plotter, target );
00499
00500 vector < FunctionRep * >::iterator first = m_func_reps.begin();
00501
00502 while ( first != m_func_reps.end() ) {
00503 FunctionRep * rep = *first++;
00504 if ( rep -> isComposite() ) continue;
00505 const vector < double > & t = rep -> principleErrors ();
00506 unsigned int number = t.size();
00507 rep->setPrincipleErrors ( begin,
00508 begin + number );
00509 begin += number;
00510 }
00511
00512 }
00513
00514 bool
00515 FunctionController::
00516 fitFunction ( PlotterBase * plotter, unsigned int )
00517 {
00518 FunctionRep * rep = getFunctionRep ( plotter );
00519
00520 return fitFunction ( plotter, rep );
00521 }
00522
00523 FunctionRep *
00524 FunctionController::
00525 getComposite ( const PlotterBase * plotter, FunctionRep * rep )
00526 {
00527 const DataRep * target = rep -> getTarget ();
00528 fillFunctionReps ( m_func_reps, plotter, target );
00529 for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
00530 FunctionRep * fr = m_func_reps[i];
00531 if ( fr -> isComposite () ) {
00532 CompositeFunctionRep * composite
00533 = dynamic_cast < CompositeFunctionRep * > ( fr );
00534 bool yes = composite -> isMember ( rep );
00535 if ( yes ) {
00536 rep = composite;
00537 break;
00538 }
00539 }
00540 }
00541
00542 return rep;
00543 }
00544
00545 bool
00546 FunctionController::
00547 fitFunction ( PlotterBase * plotter, FunctionRep * rep )
00548 {
00549 rep = getComposite ( plotter, rep );
00550 bool ok = rep -> fitFunction ();
00551 if ( rep -> isComposite () ) {
00552 setErrorsFromComposite ( plotter, rep );
00553 }
00554
00555 return ok;
00556 }
00557
00558 bool
00559 FunctionController::
00560 tryFitFunction ( PlotterBase * plotter, FunctionRep * func_rep )
00561 {
00562 saveParameters ( plotter );
00563 double chi_sq = func_rep -> objectiveValue ();
00564 bool ok = fitFunction ( plotter, func_rep );
00565 if ( ok ) {
00566 double new_chi_sq = func_rep -> objectiveValue ();
00567 if ( new_chi_sq > chi_sq ) {
00568 restoreParameters ( plotter );
00569 }
00570 } else {
00571 restoreParameters ( plotter );
00572 }
00573
00574 return ok;
00575 }
00576
00577 void
00578 FunctionController::
00579 setFitRange ( PlotterBase * plotter, const Range & range )
00580 {
00581 if ( hasFunction ( plotter, 0 ) ) {
00582 FunctionRep * frep = getFunctionRep ( plotter );
00583 frep -> setCutRange ( range );
00584 plotter -> update ();
00585 }
00586 }
00587
00588 void
00589 FunctionController::
00590 setFitRange ( PlotterBase * plotter, double low, double high )
00591 {
00592 const Range range ( low, high );
00593 setFitRange ( plotter, range );
00594 }
00595
00596 void FunctionController::saveParameters ( PlotterBase * plotter )
00597 {
00598 findFunctions ( plotter );
00599 vector< FunctionRep * >::iterator it = m_func_reps.begin ();
00600
00601 for ( ; it != m_func_reps.end (); ++it ) {
00602 FunctionRep * frep = *it;
00603 FunctionBase * function = frep->getFunction ();
00604 if ( function->isComposite () ) continue;
00605 frep->saveParameters ();
00606 }
00607 }
00608
00609 void FunctionController::restoreParameters ( PlotterBase * plotter )
00610 {
00611 findFunctions ( plotter );
00612 vector< FunctionRep * >::iterator it = m_func_reps.begin ();
00613
00614 for ( ; it != m_func_reps.end (); ++it ) {
00615 FunctionRep * frep = *it;
00616 FunctionBase * function = frep->getFunction ();
00617 if ( function->isComposite () ) continue;
00618 frep->restoreParameters ();
00619 }
00620
00621 }
00622
00623 FunctionRep *
00624 FunctionController::
00625 getFunctionRep ( const PlotterBase * plotter ) const
00626 {
00627 FunctionRep * frep = 0;
00628 DataRep * rep = 0;
00629 int index = plotter->activePlotIndex ();
00630
00631 if ( index >= 0 ) {
00632 rep = plotter -> getDataRep ( index );
00633 }
00634 else {
00635 rep = plotter -> getTarget ();
00636 }
00637
00638 FunctionRep * test = dynamic_cast < FunctionRep * > ( rep );
00639 if ( test != 0 ) {
00640 frep = test;
00641 }
00642 else {
00643 frep = getFunctionRep ( plotter, rep );
00644 }
00645
00646 return frep;
00647 }
00648
00652 FunctionRep *
00653 FunctionController::
00654 getFunctionRep ( const PlotterBase * plotter, const DataRep * datarep ) const
00655 {
00656 FunctionRep * frep = 0;
00657
00658 fillFunctionReps ( m_func_reps, plotter, datarep );
00659
00660 for ( unsigned int i = 0; i < m_func_reps.size(); i++ ) {
00661 FunctionRep * rep = m_func_reps[i];
00662 const DataRep * drep = rep -> getTarget ();
00663 if ( drep != datarep ) continue;
00664 frep = rep;
00665 if ( frep -> isComposite () ) break;
00666 }
00667
00668 return frep;
00669 }
00670
00671 ViewBase * FunctionController::
00672 createFuncView ( const ViewFactory * factory,
00673 PlotterBase * plotter,
00674 const std::string & name )
00675 {
00676 FunctionRep * frep = getFunctionRep ( plotter );
00677 assert ( frep != 0 );
00678
00679 DisplayController * controller = DisplayController::instance ();
00680 string nullstring ("");
00681
00682 return controller->createTextView ( factory, frep, name, nullstring );
00683 }
00684
00685 bool
00686 FunctionController::
00687 setFitter ( const PlotterBase * plotter, const std::string & name )
00688 {
00689 FunctionRep * frep = getFunctionRep ( plotter );
00690 bool ok = false;
00691 if ( frep != 0 ) {
00692 ok = setFitter ( frep, name );
00693 }
00694
00695 return ok;
00696 }
00697
00698 const string &
00699 FunctionController::
00700 getFitterName ( const PlotterBase * plotter )
00701 {
00702 FunctionRep * frep = getFunctionRep ( plotter );
00703 Fitter * fitter = frep -> getFitter ();
00704
00705 return fitter -> name ();
00706 }
00707
00708 Fitter *
00709 FunctionController::
00710 getFitter ( const PlotterBase * plotter )
00711 {
00712 FunctionRep * frep = getFunctionRep ( plotter );
00713 Fitter * fitter = frep -> getFitter ();
00714
00715 return fitter;
00716 }
00717
00718
00719 bool
00720 FunctionController::
00721 setFitter ( FunctionRep * frep, const std::string & name )
00722 {
00723 FitterFactory * factory = FitterFactory::instance ();
00724 Fitter * fitter = factory -> create ( name );
00725
00726 return frep -> setFitter ( fitter );
00727 }
00728
00729 void
00730 FunctionController::
00731 setDefaultFitter ( const std::string & name )
00732 {
00733 FitterFactory * factory = FitterFactory::instance ();
00734 factory -> setDefault ( name );
00735 }
00736
00737 bool
00738 FunctionController::
00739 changeFitter ( const PlotterBase * plotter,
00740 const DataRep * drep,
00741 const std::string & name )
00742 {
00743 FitterFactory * factory = FitterFactory::instance ();
00744 Fitter * proto = factory -> prototype ( name );
00745
00746 FunctionRep * frep = getFunctionRep ( plotter, drep );
00747 const FunctionBase * func = frep -> getFunction ();
00748 bool yes = proto -> isCompatible ( func );
00749 if ( yes ) {
00750 const Fitter * old_fitter = frep -> getFitter ();
00751 Fitter * new_fitter = factory -> create ( name );
00752 new_fitter -> copyFrom ( old_fitter );
00753 frep -> setFitter ( new_fitter );
00754 }
00755 return yes;
00756 }
00757
00758 double
00759 FunctionController::
00760 getObjectiveValue ( const PlotterBase * plotter, const DataRep * datarep )
00761 {
00762 FunctionRep * rep = getFunctionRep ( plotter, datarep );
00763
00764 return rep -> objectiveValue ();
00765 }
00766
00767 const vector < vector < double > > &
00768 FunctionController::
00769 getCovarianceMatrix ( const PlotterBase * plotter )
00770 {
00771 FunctionRep * rep = getFunctionRep ( plotter );
00772
00773 return rep -> covarianceMatrix ();
00774 }
00775
00780 double
00781 FunctionController::
00782 getChiSquared ( const PlotterBase * plotter )
00783 {
00784 const DataRep * rep = plotter -> getTarget ();
00785
00786 return getObjectiveValue ( plotter, rep );
00787 }
00788
00789 int
00790 FunctionController::
00791 getDegreesOfFreedom ( const PlotterBase * plotter )
00792 {
00793 FunctionRep * rep = getFunctionRep ( plotter );
00794
00795 return rep -> degreesOfFreedom ();
00796 }
00797
00798 NTuple *
00799 FunctionController::
00800 createNTuple ( const PlotterBase * plotter, const FunctionRep * rep )
00801 {
00802 NTuple * ntuple = 0;
00803 if ( rep == 0 ) {
00804 int old_index = plotter -> activePlotIndex ();
00805 int index = getUniqueNonFunctionIndex ( plotter );
00806 if ( index < 0 ) return 0;
00807
00808 int saved_index = plotter -> activePlotIndex ();
00809 PlotterBase * pb = const_cast < PlotterBase * > ( plotter );
00810 pb -> setActivePlot ( index, false );
00811 ntuple = plotter -> createNTuple ();
00812 pb -> setActivePlot ( saved_index, false );
00813 rep = getFunctionRep ( plotter );
00814 pb -> setActivePlot ( old_index, true );
00815 }
00816 else {
00817 const DataRep * target = rep -> getTarget ();
00818 ntuple = target -> createNTuple ();
00819 }
00820 if ( rep != 0 ) {
00821 const FunctionBase * function = rep -> getFunction ();
00822 vector < double > & x = ntuple -> getColumn ( 0 );
00823 vector < double > & y = ntuple -> getColumn ( 1 );
00824 unsigned int size = ntuple -> rows ();
00825 vector < double > values ( size );
00826 vector < double > residuals ( size );
00827
00828 for ( unsigned int i = 0; i < x.size(); i++ ) {
00829 values [i] = function -> operator () ( x[i] );
00830 residuals [i] = y[i] - values [i];
00831 }
00832
00833 ntuple -> addColumn ( function->name (), values );
00834 ntuple -> addColumn ( "Residuals", residuals );
00835 }
00836
00837 return ntuple;
00838 }
00839
00840 PlotterBase *
00841 FunctionController::
00842 createResidualsDisplay ( PlotterBase * plotter, const FunctionRep * frep )
00843 {
00844 NTuple * ntuple = createNTuple ( plotter, frep );
00845 ntuple -> setTitle ( plotter -> getTitle () );
00846 const string old_label ( "Residuals" );
00847 DataSourceController::instance () -> registerNTuple ( ntuple );
00848 vector < string > bindings ( 4 );
00849 bindings[0] = ntuple -> getLabelAt ( 0 );
00850 bindings[1] = old_label;
00851 bindings[2] = "nil";
00852 bindings[3] = ntuple -> getLabelAt ( 3 );
00853 bool axis_scaled = plotter -> isAxisScaled ( Axes::Y );
00854 if ( axis_scaled ) {
00855 int i = ntuple -> indexOf ( old_label );
00856 const string new_label ("Residuals (scaled)");
00857 ntuple -> setLabelAt ( new_label , i );
00858 bindings[1] = new_label;
00859 }
00860 DisplayController * controller = DisplayController::instance ();
00861
00862
00863 PlotterBase * new_plotter = controller -> createDisplay ( "XY Plot",
00864 *ntuple,
00865 bindings );
00866 if ( axis_scaled ) {
00867 double factor = plotter -> getScaleFactor ( Axes::Y );
00868 new_plotter -> setScaleFactor ( Axes::Y, factor );
00869 }
00870 controller->setLog( new_plotter, "x",
00871 controller -> getLog( plotter, "x" ) );
00872
00873 return new_plotter;
00874 }
00875
00876 PlotterBase *
00877 FunctionController::
00878 createNewEllipsoidDisplay ( PlotterBase * masterPlot, FunctionRep * rep )
00879 {
00880 int n = 100;
00881 double xmin, xmax, ymin, ymax;
00882
00883
00884 NTuple* ntuple = new( NTuple );
00885 ellipsoidNTuple ( masterPlot, rep, ntuple, n, xmin, xmax, ymin, ymax );
00886 ntuple -> setTitle ( masterPlot -> getTitle () );
00887
00888 vector < string > bindings ( 1 );
00889 bindings[0] = ntuple -> getLabelAt ( 0 );
00890
00891 DisplayController * dcontroller = DisplayController::instance ();
00892
00893
00894 PlotterBase * ellipsePlot = dcontroller ->
00895 createDisplay ( "Image", *ntuple, bindings );
00896
00897
00898
00899
00900 ellipsePlot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
00901 ellipsePlot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
00902
00903 ellipsePlot -> setOffsets ( xmin, ymin );
00904
00905 ellipsePlot -> setRange ( string( "X" ), xmin, xmax );
00906 ellipsePlot -> setRange ( string( "Y" ), ymin, ymax );
00907
00908
00909
00910
00911 int index = dcontroller -> activeDataRepIndex( masterPlot );
00912 ellipsePlot -> setParentDataRepIndex( index );
00913 ellipsePlot -> setParentPlotter( masterPlot );
00914
00915 return ellipsePlot;
00916 }
00917
00918 void
00919 FunctionController::
00920 ellipsoidNTuple ( PlotterBase* masterPlot, FunctionRep * frep,
00921 NTuple* nt, int n,
00922 double & xmin, double & xmax,
00923 double & ymin, double & ymax )
00924 {
00925 string text ( "Confidence ellipse " );
00926 text += String::convert ( ++m_confid_count );
00927 nt -> setName ( text );
00928
00929 DataSourceController * controller = DataSourceController::instance ();
00930 controller -> registerNTuple ( nt );
00931
00932 assert( frep );
00933
00934
00935
00936
00937
00938 vector< vector< double > > Sigma( 2 );
00939 Sigma[0].resize( 2, 0.0 );
00940 Sigma[1].resize( 2, 0.0 );
00941
00942 const vector < vector < double > > & covariance
00943 = frep -> covarianceMatrix ();
00944
00945 Sigma[0][0] = covariance [m_x][m_x];
00946 Sigma[0][1] = covariance [m_x][m_y];
00947 Sigma[1][0] = covariance [m_y][m_x];
00948 Sigma[1][1] = covariance [m_y][m_y];
00949
00950
00951 vector< vector< double > > SigmaInv;
00952 invertMatrix( Sigma, SigmaInv );
00953
00954
00955
00956
00957 vector< double > xbar( 2 );
00958
00959 const Fitter * fitter = getFitter ( masterPlot );
00960 vector < double > free_parms;
00961 fitter -> fillFreeParameters ( free_parms );
00962 xbar[ 0 ] = free_parms [ m_x ];
00963 xbar[ 1 ] = free_parms [ m_y ];
00964
00965
00966
00967
00968 unsigned int mu = free_parms.size();
00969 NTuple * boundingTuple = ellipse( xbar, SigmaInv, m_deltaXSq[ mu - 1 ] );
00970
00971 xmin = boundingTuple -> minElement ( 0 );
00972 xmax = boundingTuple -> maxElement ( 0 );
00973 ymin = boundingTuple -> minElement ( 1 );
00974 ymax = boundingTuple -> maxElement ( 1 );
00975
00976 delete boundingTuple;
00977
00978
00979
00980
00981 vector< double > p( n * n ), a( 2 );
00982
00983 double dx = ( xmax - xmin ) / ( n - 1.0 );
00984 double dy = ( ymax - ymin ) / ( n - 1.0 );
00985 double delta = 0.0;
00986
00987 for( int i = 0; i < n; i++ )
00988 for( int j = 0; j < n; j++ )
00989 {
00990 a[ 0 ] = xmin + i * dx - xbar[ 0 ] ;
00991 a[ 1 ] = ymin + j * dy - xbar[ 1 ];
00992 delta = quadraticProduct( SigmaInv, a );
00993 p[ n * i + j ] = 1 - gammq( mu / 2.0, delta / 2.0 );
00994 }
00995
00996 nt -> clear();
00997 nt -> addColumn( "Percent", 100 * p );
00998
00999 return;
01000 }
01001
01002 PlotterBase *
01003 FunctionController::
01004 refreshEllipsoidDisplay ( PlotterBase * plot, FunctionRep * frep )
01005 {
01006 int n = 100;
01007 double xmin, xmax, ymin, ymax;
01008
01009 PlotterBase* masterPlot = plot->getParentPlotter();
01010
01011
01012 NTuple* ntuple = new( NTuple );
01013 ellipsoidNTuple ( masterPlot, frep, ntuple, n, xmin, xmax, ymin, ymax );
01014 ntuple -> setTitle ( masterPlot -> getTitle () );
01015
01016
01017
01018
01019 DataRep * drep = plot -> selectedDataRep();
01020 NTupleProjector * ntProjector =
01021 dynamic_cast < NTupleProjector * > ( drep -> getProjector() );
01022 ntProjector -> setNTuple ( ntuple );
01023
01024 NTuple * nt = const_cast < NTuple * > ( ntuple );
01025 nt -> addObserver ( ntProjector );
01026
01027
01028
01029
01030 plot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
01031 plot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
01032
01033 plot -> setOffsets ( xmin, ymin );
01034
01035 plot -> setRange ( string( "X" ), xmin, xmax );
01036 plot -> setRange ( string( "Y" ), ymin, ymax );
01037
01038 return plot;
01039 }
01040
01041 NTuple *
01042 FunctionController::
01043 ellipse ( const std::vector< double > & xbar,
01044 std::vector< std::vector < double > > & SigmaInv,
01045 double Csquare )
01046 {
01047
01048 assert( xbar.size() == 2 );
01049
01050
01051 assert( SigmaInv.size() == 2 );
01052 assert( SigmaInv[0].size() == 2 || SigmaInv[1].size() == 2 );
01053
01054
01055
01056
01057 double temp = ( SigmaInv[0][1] + SigmaInv[1][0] ) / 2;
01058 SigmaInv[0][1] = temp;
01059 SigmaInv[1][0] = temp;
01060
01061
01062 assert( Csquare > DBL_EPSILON );
01063
01064
01065 double b = -( SigmaInv[0][0] + SigmaInv[1][1] );
01066 double c = SigmaInv[0][0] * SigmaInv[1][1] - SigmaInv[0][1] * SigmaInv[1][0];
01067
01068 double lambda1 = ( -b + sqrt( b * b - 4 * c ) ) / 2;
01069 double lambda2 = ( -b - sqrt( b * b - 4 * c ) ) / 2;
01070
01071
01072 double alpha1 = atan( (SigmaInv[0][0] - lambda1) / SigmaInv[0][1] );
01073 double a = sqrt( Csquare / lambda1 );
01074 b = sqrt( Csquare / lambda2 );
01075
01076
01077 double Rot00 = cos( alpha1 );
01078 double Rot11 = cos( alpha1 );
01079 double Rot01 = -sin( alpha1 );
01080 double Rot10 = sin( alpha1 );
01081
01082
01083 int N = 100;
01084 vector< double > x, y;
01085
01086 x.resize( N, 0.0 );
01087 y.resize( N, 0.0 );
01088
01089 for( int i = 0; i < N; i++)
01090 {
01091
01092 double theta = ( 2 * M_PI * i ) / ( (double) (N-1) );
01093 double x0 = a * cos( theta );
01094 double x1 = b * sin( theta );
01095
01096
01097 double xrot0 = Rot00 * x0 + Rot01 * x1;
01098 double xrot1 = Rot10 * x0 + Rot11 * x1;
01099
01100
01101 x[i] = xrot0 + xbar[0];
01102 y[i] = xrot1 + xbar[1];
01103 }
01104
01105
01106 NTuple* ntuple = new( NTuple );
01107 ntuple -> addColumn( "X", x );
01108 ntuple -> addColumn( "Y", y );
01109
01110 return ntuple;
01111 }
01112
01113 int FunctionController::setEllpsoidParamIndex( hippodraw::Axes::Type axes,
01114 int index )
01115 {
01116 assert( axes == Axes::X || axes == Axes::Y );
01117
01118 if( axes == Axes::X )
01119 m_x = index;
01120
01121 if( axes == Axes::Y )
01122 m_y = index;
01123
01124 return EXIT_SUCCESS;
01125 }
01126
01127 bool
01128 FunctionController::
01129 functionExists ( const std::string & name )
01130 {
01131 FunctionFactory * factory = FunctionFactory::instance ();
01132
01133 return factory -> exists ( name );
01134 }
01135
01136 bool
01137 FunctionController::
01138 isCompatible ( const std::string & function,
01139 const std::string & fitter )
01140 {
01141 bool yes = true;
01142
01143 FunctionFactory * fun_fac = FunctionFactory::instance ();
01144 FunctionBase * fun_proto = fun_fac -> prototype ( function );
01145
01146 FitterFactory * fit_fac = FitterFactory::instance ();
01147 Fitter * fit_proto = fit_fac -> prototype ( fitter );
01148
01149 yes = fit_proto -> isCompatible ( fun_proto );
01150
01151 return yes;
01152 }
01153
01154 const vector < string > &
01155 FunctionController::
01156 getFunctionNames () const
01157 {
01158 return FunctionFactory::instance () -> names ();
01159 }