Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2030) +++ trunk/MagicSoft/Mars/Changelog (revision 2031) @@ -1,3 +1,12 @@ -*-*- END OF LINE -*-*- + 2003/04/28: Abelardo Moralejo + + * mmontecarlo/MMcEnergyEst.[h,cc] + - Lots of fixes after Thomas suggestions. Now cuts are not part + of the class, but introduced via a new MFilter* member. Changed + all Char_t* for TString. Changed own TMiniut pointer by gMinuit. + Removed couts and used fLog instead. Function fcn is no longer + declared external. + 2003/04/27: Thomas Bretz Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 2030) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 2031) @@ -30,4 +30,9 @@ // // // Class for otimizing the parameters of the energy estimator // +// // +// FIXME: the class must be made more flexible, allowing for different // +// parametrizations to be used. Also a new class is needed, a container // +// for the parameters of the energy estimator. // +// FIXME: the starting values of the parameters are now fixed. // // // // // @@ -55,6 +60,24 @@ #include "MLogManip.h" -extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); -TMinuit *minuit; + +//------------------------------------------------------------------------ +// +// fcn calculates the function to be minimized (using TMinuit::Migrad) +// +void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit(); + + MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); + + MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); + MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); + + eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); + + evtloop->Eventloop(); + + f = eval->GetChisq(); +} ClassImp(MMcEnergyEst); @@ -71,5 +94,4 @@ SetHillasName("MHillas"); SetHillasSrcName("MHillasSrc"); - SetHadronnessName("MHadronness"); } @@ -84,5 +106,5 @@ MFEventSelector selector; selector.SetNumSelectEvts(fNevents); - cout << "Read events from file '" << fInFile << "'" << endl; + *fLog << inf << "Read events from file '" << fInFile << "'" << endl; MReadTree read("Events"); @@ -91,10 +113,5 @@ read.SetSelector(&selector); - MFCT1SelFinal filterhadrons; - filterhadrons.SetHadronnessName(fHadronnessName); - filterhadrons.SetCuts(fMaxHadronness, fMaxAlpha, fMaxDist); - filterhadrons.SetInverted(); - - cout << "Define columns of matrix" << endl; + *fLog << inf << "Define columns of matrix" << endl; MHMatrix matrix; Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); @@ -103,6 +120,6 @@ if (colenergy < 0 || colimpact < 0) { - cout << "colenergy, colimpact = " << colenergy << ", " - << colimpact << endl; + *fLog << err << dbginf << "colenergy, colimpact = " << colenergy << ", " + << colimpact << endl; return; } @@ -112,10 +129,10 @@ eest.InitMapping(&matrix); - cout << "--------------------------------------" << endl; - cout << "Fill events into the matrix" << endl; - if ( !matrix.Fill(&parlist, &read, &filterhadrons) ) + *fLog << inf << "--------------------------------------" << endl; + *fLog << inf << "Fill events into the matrix" << endl; + if ( !matrix.Fill(&parlist, &read, fEventFilter) ) return; - cout << "Matrix was filled with " << matrix.GetNumRows() - << " events" << endl; + *fLog << inf << "Matrix was filled with " << matrix.GetNumRows() + << inf << " events" << endl; //----------------------------------------------------------------------- @@ -123,6 +140,6 @@ // Setup the eventloop which will be executed in the fcn of MINUIT // - cout << "--------------------------------------" << endl; - cout << "Setup event loop to be used in 'fcn'" << endl; + *fLog << inf << "--------------------------------------" << endl; + *fLog << inf << "Setup event loop to be used in 'fcn'" << endl; MTaskList tasklist; @@ -148,5 +165,5 @@ - cout << "Event loop was setup" << endl; + *fLog << inf << "Event loop was setup" << endl; MEvtLoop evtloop; evtloop.SetParList(&parlist); @@ -159,18 +176,18 @@ // Be careful: This is not thread safe // - cout << "--------------------------------------" << endl; - cout << "Start minimization" << endl; - - minuit = new TMinuit(12); - minuit->SetPrintLevel(-1); - - minuit->SetFCN(fcn); - minuit->SetObjectFit(&evtloop); - - // Ready for: minuit->mnexcm("SET ERR", arglist, 1, ierflg) - - if (minuit->SetErrorDef(1)) - { - cout << "SetErrorDef failed." << endl; + *fLog << inf << "--------------------------------------" << endl; + *fLog << inf << "Start minimization" << endl; + + gMinuit = new TMinuit(12); + gMinuit->SetPrintLevel(-1); + + gMinuit->SetFCN(fcn); + gMinuit->SetObjectFit(&evtloop); + + // Ready for: gMinuit->mnexcm("SET ERR", arglist, 1, ierflg) + + if (gMinuit->SetErrorDef(1)) + { + *fLog << err << dbginf << "SetErrorDef failed." << endl; return; } @@ -211,9 +228,9 @@ Double_t limup = 0; - Bool_t rc = minuit->DefineParameter(i, name, vinit, step, limlo, limup); + Bool_t rc = gMinuit->DefineParameter(i, name, vinit, step, limlo, limup); if (!rc) continue; - cout << "Error in defining parameter #" << i << endl; + *fLog << err << dbginf << "Error in defining parameter #" << i << endl; return; } @@ -230,9 +247,9 @@ Double_t limup = 0; - Bool_t rc = minuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup); + Bool_t rc = gMinuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup); if (!rc) continue; - cout << "Error in defining parameter #" << i+fA.GetSize() << endl; + *fLog << err << dbginf << "Error in defining parameter #" << i+fA.GetSize() << endl; return; } @@ -244,19 +261,19 @@ gLog.SetNullOutput(kTRUE); - Bool_t rc = minuit->Migrad(); + Bool_t rc = gMinuit->Migrad(); gLog.SetNullOutput(kFALSE); if (rc) { - cout << "Migrad failed." << endl; + *fLog << err << dbginf << "Migrad failed." << endl; return; } - cout << endl; + *fLog << inf << endl; clock.Stop(); clock.Print(); - cout << endl; - - cout << "Resulting Chisq: " << minuit->fAmin << endl; + *fLog << inf << endl; + + *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; // @@ -266,5 +283,5 @@ { Double_t x1, x2; - minuit->GetParameter(i,x1,x2); + gMinuit->GetParameter(i,x1,x2); fA[i] = x1; } @@ -272,5 +289,5 @@ { Double_t x1, x2; - minuit->GetParameter(i,x1,x2); + gMinuit->GetParameter(i,x1,x2); fB[i-fA.GetSize()] = x1; } @@ -278,5 +295,5 @@ // eest.Print(); eest.StopMapping(); - cout << "Minimization finished" << endl; + *fLog << inf << "Minimization finished" << endl; } @@ -289,8 +306,8 @@ { for (Int_t i = 0; i < fA.GetSize(); i++) - cout << "Parameter " << i << ": " << fA[i] << endl; + *fLog << inf << "Parameter " << i << ": " << fA[i] << endl; for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++) - cout << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl; + *fLog << inf << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl; /* @@ -298,30 +315,9 @@ Double_t amin, edm, errdef; Int_t nvpar, nparx, icstat; - minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat); - minuit->mnprin(3, amin); + gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat); + gMinuit->mnprin(3, amin); */ } -//------------------------------------------------------------------------ -// -// fcn calculates the function to be minimized (using TMinuit::Migrad) -// -void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - MEvtLoop *evtloop = (MEvtLoop*)minuit->GetObjectFit(); - - MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); - - MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); - MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); - - eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); - - evtloop->Eventloop(); - - f = eval->GetChisq(); -} - - - + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h (revision 2030) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h (revision 2031) @@ -10,12 +10,16 @@ #endif +#include "MFilter.h" + class MMcEnergyEst : public MParContainer { private: - Char_t *fInFile, *fOutFile; - Char_t *fHillasName, *fHillasSrcName, *fHadronnessName; + TString fInFile, fOutFile; + TString fHillasName; + TString fHillasSrcName; Int_t fNevents; - Float_t fMaxHadronness, fMaxAlpha, fMaxDist; + + MFilter *fEventFilter; //! TArrayD fA; @@ -25,23 +29,16 @@ MMcEnergyEst(const char *name=NULL, const char *title=NULL); - void SetInFile(Char_t *name) {fInFile = name;} - void SetOutFile(Char_t *name) {fOutFile = name;} - void SetHillasName(Char_t *name) {fHillasName = name;} - void SetHillasSrcName(Char_t *name) {fHillasSrcName = name;} - void SetHadronnessName(Char_t *name) {fHadronnessName = name;} - void SetNevents(Int_t n) {fNevents = n;} - void SetMaxHadronness(Float_t x) {fMaxHadronness = x;} - void SetMaxAlpha(Float_t x) {fMaxAlpha = x;} - void SetMaxDist(Float_t x) {fMaxDist = x;} + void SetInFile(const TString &name) {fInFile = name;} + void SetOutFile(const TString &name) {fOutFile = name;} + void SetHillasName(const TString &name) {fHillasName = name;} + void SetHillasSrcName(const TString &name) {fHillasSrcName = name;} + void SetEventFilter(MFilter *filter) {fEventFilter = filter;} + void SetNevents(Int_t n) {fNevents = n;} - Char_t* GetInFile() {return fInFile;} - Char_t* GetOutFile() {return fOutFile;} - Char_t* GetHillasName() {return fHillasName;} - Char_t* GetHillasSrcName() {return fHillasSrcName;} - Char_t* GetHadronnessName() {return fHadronnessName;} - Int_t GetNevents() {return fNevents;} - Float_t GetMaxHadronness() {return fMaxHadronness;} - Float_t GetMaxAlpha() {return fMaxAlpha;} - Float_t GetMaxDist() {return fMaxDist;} + TString GetInFile() const {return fInFile;} + TString GetOutFile() const {return fOutFile;} + TString GetHillasName() const {return fHillasName;} + TString GetHillasSrcName() const {return fHillasSrcName;} + Int_t GetNevents() const {return fNevents;} Double_t GetCoeff(Int_t i) { return i