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