Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 1963) +++ trunk/MagicSoft/Mars/Changelog (revision 1964) @@ -1,3 +1,13 @@ -*-*- END OF LINE -*-*- + + 2003/04/19: Abelardo Moralejo + + * mmontecarlo/MMcEnergyEst.[h,cc] + - Added. Contains routine for optimization of parameters of + energy estimator. + + * mmontecarlo/MMcEnergyMigration.[h,cc] + -Added. Task to fill the energy migration matrix histograms + contained in class MHMcEnergyMigration. 2003/04/17: Wolfgang Wittek Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 1964) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 1964) @@ -0,0 +1,326 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 1/2002 +! Author(s): Wolfgang Wittek 1/2002 +! Author(s): Abelardo Moralejo 4/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMcEnergyEst // +// // +// Class for otimizing the parameters of the energy estimator // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMcEnergyEst.h" + +#include +#include +#include + +#include "MParList.h" +#include "MTaskList.h" +#include "MGeomCamCT1.h" +#include "MFEventSelector.h" +#include "MReadTree.h" +#include "MFCT1SelFinal.h" +#include "MHMatrix.h" +#include "MEnergyEstParam.h" +#include "MMatrixLoop.h" +#include "MChisqEval.h" +#include "MEvtLoop.h" +#include "MDataElement.h" +#include "MDataMember.h" +#include "MLog.h" +#include "MLogManip.h" + +extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); +TMinuit *minuit; + +ClassImp(MMcEnergyEst); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMcEnergyEst::MMcEnergyEst(const char *name, const char *title) +{ + fName = name ? name : "MMcEnergyEst"; + fTitle = title ? title : "Optimizer of the energy estimator"; + + SetHillasName("MHillas"); + SetHillasSrcName("MHillasSrc"); + SetHadronnessName("MHadronness"); +} + +//------------------------------------------------------------------------ +// +// Optimization (via migrad minimization) of parameters of energy estimation. +// +void MMcEnergyEst::FindParams() +{ + MParList parlist; + + MFEventSelector selector; + selector.SetNumSelectEvts(fNevents); + cout << "Read events from file '" << fInFile << "'" << endl; + + MReadTree read("Events"); + read.AddFile(fInFile); + read.DisableAutoScheme(); + read.SetSelector(&selector); + + MFCT1SelFinal filterhadrons; + filterhadrons.SetCuts(fMaxHadronness, fMaxAlpha, fMaxDist); + filterhadrons.SetInverted(); + + cout << "Define columns of matrix" << endl; + MHMatrix matrix; + Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); + Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact"); + + if (colenergy < 0 || colimpact < 0) + { + cout << "colenergy, colimpact = " << colenergy << ", " + << colimpact << endl; + return; + } + + MEnergyEstParam eest(fHillasName); + eest.Add(fHillasSrcName); + eest.InitMapping(&matrix); + + cout << "--------------------------------------" << endl; + cout << "Fill events into the matrix" << endl; + if ( !matrix.Fill(&parlist, &read, &filterhadrons) ) + return; + cout << "Matrix was filled with " << matrix.GetNumRows() + << " events" << endl; + + //----------------------------------------------------------------------- + // + // Setup the eventloop which will be executed in the fcn of MINUIT + // + cout << "--------------------------------------" << endl; + cout << "Setup event loop to be used in 'fcn'" << endl; + + MTaskList tasklist; + + MMatrixLoop loop(&matrix); + + MChisqEval eval; + eval.SetY1(new MDataElement(&matrix, colenergy)); + eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); + eval.SetOwner(); + + // + // Entries in MParList + + parlist.AddToList(&tasklist); + + // + // Entries in MTaskList + + tasklist.AddToList(&loop); + tasklist.AddToList(&eest); + tasklist.AddToList(&eval); + + + cout << "Event loop was setup" << endl; + MEvtLoop evtloop; + evtloop.SetParList(&parlist); + + // + //---------- Start of minimization part -------------------- + // + // Do the minimization with MINUIT + // + // 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; + return; + } + + // + // Set initial values of the parameters (close enough to the final ones, taken + // from previous runs of the procedure). Parameter fA[4] is not used in the default + // energy estimation model (from D. Kranich). + // + fA.Set(5); + fB.Set(7); + + fA[0] = 21006.2; + fA[1] = -43.2648; + fA[2] = -690.403; + fA[3] = -0.0428544; + fA[4] = 1.; + fB[0] = -3391.05; + fB[1] = 136.58; + fB[2] = 0.253807; + fB[3] = 254.363; + fB[4] = 61.0386; + fB[5] = -0.0190984; + fB[6] = -421695; + + // + // Set starting values and step sizes for parameters + // + for (Int_t i=0; iDefineParameter(i, name, vinit, step, limlo, limup); + if (!rc) + continue; + + cout << "Error in defining parameter #" << i << endl; + return; + } + + for (Int_t i=0; iDefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup); + if (!rc) + continue; + + cout << "Error in defining parameter #" << i+fA.GetSize() << endl; + return; + } + + TStopwatch clock; + clock.Start(); + + // Now ready for minimization step: + + gLog.SetNullOutput(kTRUE); + Bool_t rc = minuit->Migrad(); + gLog.SetNullOutput(kFALSE); + + if (rc) + { + cout << "Migrad failed." << endl; + return; + } + + cout << endl; + clock.Stop(); + clock.Print(); + cout << endl; + + cout << "Resulting Chisq: " << minuit->fAmin << endl; + + // + // Update values of fA, fB: + // + for (Int_t i = 0; i < fA.GetSize(); i++) + { + Double_t x1, x2; + minuit->GetParameter(i,x1,x2); + fA[i] = x1; + } + for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++) + { + Double_t x1, x2; + minuit->GetParameter(i,x1,x2); + fB[i-fA.GetSize()] = x1; + } + + // eest.Print(); + eest.StopMapping(); + cout << "Minimization finished" << endl; + +} + +//------------------------------------------------------------------------ +// +// Print current values of parameters +// +void MMcEnergyEst::Print(Option_t *o) +{ + for (Int_t i = 0; i < fA.GetSize(); i++) + cout << "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; + + /* + // Print results + Double_t amin, edm, errdef; + Int_t nvpar, nparx, icstat; + minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat); + minuit->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 1964) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h (revision 1964) @@ -0,0 +1,56 @@ +#ifndef MARS_MMcEnergyEst +#define MARS_MMcEnergyEst + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MMcEnergyEst : public MParContainer +{ +private: + + Char_t *fInFile, *fOutFile; + Char_t *fHillasName, *fHillasSrcName, *fHadronnessName; + Int_t fNevents; + Float_t fMaxHadronness, fMaxAlpha, fMaxDist; + + TArrayD fA; + TArrayD fB; + +public: + 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;} + + 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;} + + Double_t GetCoeff(Int_t i) { return i +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMcEnergyMigration // +// // +// This task fills the histograms contained in the // +// class MHMcEnergyMigration // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MMcEnergyMigration.h" +#include "MHMcEnergyMigration.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +ClassImp(MMcEnergyMigration); + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MMcEnergyMigration::MMcEnergyMigration(const char *name, const char *title) +{ + fName = name ? name : "MMcEnergyMigration"; + fTitle = title ? title : "Task to fill histograms in MHMcEnergyMigration"; +} + +// -------------------------------------------------------------------------- +// +// Destructor +// +MMcEnergyMigration::~MMcEnergyMigration() +{ + if (fHistMigration) + delete fHistMigration; +} + +// -------------------------------------------------------------------------- +// +// Preprocess: look for the MHMcEnergyMigration object in the parameter list. +// +Bool_t MMcEnergyMigration::PreProcess(MParList *pList) +{ + fHistMigration = (MHMcEnergyMigration*)pList->FindObject("MHMcEnergyMigration"); + if (!fHistMigration) + { + *fLog << dbginf << "MHMcEnergyMigration not found... aborting." << endl; + return kFALSE; + } + + fHistMigration->SetupFill(pList); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms in the MHMcEnergyMigration object. +// +Bool_t MMcEnergyMigration::Process() +{ + fHistMigration->Fill(0); // Argument of Fill is dummy. + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Postprocessing does nothing at the moment. +// +Bool_t MMcEnergyMigration::PostProcess() +{ + return kTRUE; +} + + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyMigration.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyMigration.h (revision 1964) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyMigration.h (revision 1964) @@ -0,0 +1,26 @@ +#ifndef MARS_MMcEnergyMigration +#define MARS_MMcEnergyMigration + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MHMcEnergyMigration; + +class MMcEnergyMigration : public MTask +{ +private: + MHMcEnergyMigration *fHistMigration; + +public: + MMcEnergyMigration(const char *name=NULL, const char *title=NULL); + ~MMcEnergyMigration(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MMcEnergyMigration, 0) // Task to fill the Migration matrix histogram +}; + +#endif