00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "TGo4FitMinuit.h"
00017
00018 #include "Riostream.h"
00019
00020 #include "TMinuit.h"
00021 #include "TArrayD.h"
00022 #include "TMatrixD.h"
00023 #include "TObjString.h"
00024
00025 #include "TGo4FitterAbstract.h"
00026 #include "TGo4FitMinuitResult.h"
00027
00028 class TMinuitEx : public TMinuit {
00029 public:
00030
00031 TMinuitEx(Int_t NumPars, TGo4FitterAbstract* fitter);
00032 virtual ~TMinuitEx();
00033
00034 virtual Int_t Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *pars, Int_t iflag);
00035
00036 protected:
00037
00038 TGo4FitterAbstract* fxFitter;
00039 };
00040
00041 TMinuitEx::TMinuitEx(Int_t NumPars, TGo4FitterAbstract* fitter) :
00042 TMinuit(NumPars), fxFitter(fitter) {
00043 }
00044
00045 TMinuitEx::~TMinuitEx() {
00046 }
00047
00048 Int_t TMinuitEx::Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *pars, Int_t iflag) {
00049 if (fxFitter) fval = fxFitter->CalculateFitFunction(pars);
00050 else fval = 0.;
00051 return 0;
00052 }
00053
00054
00055
00056 TGo4FitMinuit::TGo4FitMinuit() : TGo4FitterAction(), fxCommands(), fxResults() {
00057 }
00058
00059 TGo4FitMinuit::TGo4FitMinuit(const char* Name) :
00060 TGo4FitterAction(Name,"Fitter minimization using TMinuit object"), fxCommands(), fxResults() {
00061 fxCommands.SetOwner(kTRUE);
00062 fxResults.SetOwner(kTRUE);
00063 }
00064
00065 TGo4FitMinuit::~TGo4FitMinuit() {
00066 }
00067
00068 void TGo4FitMinuit::AddCommand(const char* iCommand)
00069 {
00070 fxCommands.Add( new TObjString(iCommand) );
00071 }
00072
00073 const char* TGo4FitMinuit::GetCommand(Int_t n)
00074 {
00075 return ((TObjString*) fxCommands[n])->GetString().Data();
00076 }
00077
00078 void TGo4FitMinuit::DoAction(TGo4FitterAbstract* Fitter) {
00079
00080 TMinuitEx fMinuit(Fitter->NumPars(), Fitter);
00081 fMinuit.SetPrintLevel(-1);
00082
00083 for(Int_t n=0;n<Fitter->NumPars();n++) {
00084 const char* FullName = Fitter->GetParFullName(n);
00085 Int_t ierflg = 0;
00086 Double_t epsilon = 0, RangeMin = 0, RangeMax = 0;
00087 if (!Fitter->GetParEpsilon(FullName,epsilon)) epsilon = 1.;
00088
00089 if (!Fitter->GetParRange(FullName,RangeMin,RangeMax)) { RangeMin = 0; RangeMax = 0; }
00090
00091 fMinuit.mnparm(n, FullName, Fitter->GetParValue(FullName),
00092 epsilon, RangeMin, RangeMax, ierflg);
00093
00094 if (Fitter->GetParFixed(FullName)) fMinuit.FixParameter(n);
00095 }
00096
00097 if (fxCommands.GetLast()<0) fMinuit.Command("MIGRAD 500 1");
00098 else
00099 for(Int_t n=0;n<=fxCommands.GetLast();n++) {
00100 TString cmd ( ((TObjString*) fxCommands[n])->GetString() );
00101 if (cmd[0] == 'r') {
00102 if (cmd.Index("result",6,0,TString::kIgnoreCase) == 0) {
00103 cmd.Remove(0,6);
00104 while ((cmd.Length()>0) && (cmd[0]==' ')) cmd.Remove(0,1);
00105 if (cmd.Length()==0) cmd = "1000";
00106 if (cmd.Length()<4) { cout << "invalid result command syntax" << endl; break; }
00107 Bool_t getpar = (cmd[0]=='1');
00108 Bool_t geterr = (cmd[1]=='1');
00109 Bool_t getmatr = (cmd[2]=='1');
00110 Bool_t getcontr = (cmd[3]=='1');
00111 cmd.Remove(0,4);
00112 while ((cmd.Length()>0) && (cmd[0]==' ')) cmd.Remove(0,1);
00113 if (cmd.Length()==0) cmd="Result";
00114
00115 TGo4FitMinuitResult *res = new TGo4FitMinuitResult(cmd,"TMinuit result object");
00116 fxResults.Add(res);
00117 res->CallMNSTAT(&fMinuit);
00118 if (getpar) res->CallMNPOUT(&fMinuit,Fitter->NumPars());
00119 if (geterr) res->CallMNERRS(&fMinuit,Fitter->NumPars());
00120 if (getmatr) res->CallMNEMAT(&fMinuit,Fitter->NumPars(), kTRUE);
00121 if (getcontr) res->GetContourPlot(&fMinuit);
00122 }
00123
00124 } else fMinuit.Command(cmd);
00125 }
00126
00127 for(Int_t n=0;n<Fitter->NumPars();n++) {
00128 Double_t value,error;
00129 fMinuit.GetParameter(n, value, error);
00130 Fitter->SetParValue(Fitter->GetParFullName(n),value);
00131 Fitter->SetParError(Fitter->GetParFullName(n),error);
00132 }
00133 }
00134
00135 TGo4FitMinuitResult* TGo4FitMinuit::GetResult(Int_t indx)
00136 {
00137 return (indx>=0) && (indx<=fxResults.GetLast()) ? (TGo4FitMinuitResult*) fxResults.At(indx) : 0;
00138 }
00139
00140 TGo4FitMinuitResult* TGo4FitMinuit::FindResult(const char* ResName)
00141 {
00142 return (TGo4FitMinuitResult*) fxResults.FindObject(ResName);
00143 }
00144
00145 void TGo4FitMinuit::AddResult(TGo4FitMinuitResult* res)
00146 {
00147 fxResults.Add(res);
00148 }
00149
00150 void TGo4FitMinuit::RemoveResult(TGo4FitMinuitResult* res)
00151 {
00152 fxResults.Remove(res);
00153 fxResults.Compress();
00154 }
00155
00156 void TGo4FitMinuit::Print(Option_t* option) const {
00157 TGo4FitterAction::Print(option);
00158 if (fxCommands.GetLast()>=0)
00159 cout << "List of commands:" << endl;
00160 for(Int_t n=0;n<=fxCommands.GetLast();n++)
00161 cout << " " << ((TObjString*) fxCommands[n])->String().Data() << endl;
00162 if (fxResults.GetLast()>=0) {
00163 cout << "List of stored results:" << endl;
00164 fxResults.Print(option);
00165 }
00166 }
00167
00168