/* ======================================================================== *\ ! ! * ! * 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): Eva Domingo, 12/2004 ! Wolfgang Wittek, 12/2004 ! ! Copyright: MAGIC Software Development, 2000-2005 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MFindDisp // // // // Class for otimizing the estimation of Disp // // // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MFindDisp.h" #include // fabs #include #include #include #include #include #include #include #include #include "MBinning.h" #include "MContinue.h" #include "MDispParameters.h" #include "MDispCalc.h" #include "MDataElement.h" #include "MDataMember.h" #include "MEvtLoop.h" #include "MFSelFinal.h" #include "MF.h" #include "MFDisp.h" #include "MFEventSelector.h" #include "MFEventSelector2.h" #include "MFillH.h" #include "MFEventSelector.h" #include "MGeomCamMagic.h" #include "MH3.h" #include "MHDisp.h" #include "MHFindSignificance.h" #include "MHMatrix.h" #include "MHOnSubtraction.h" #include "MLog.h" #include "MLogManip.h" #include "MMatrixLoop.h" #include "MMinuitInterface.h" #include "MParameters.h" #include "MParList.h" #include "MProgressBar.h" #include "MReadMarsFile.h" #include "MReadTree.h" #include "MStatusDisplay.h" #include "MTaskList.h" ClassImp(MFindDisp); using namespace std; //------------------------------------------------------------------------ // fcnDisp //------------------------------------------------------------------------ // // - calculates the quantity to be minimized (by TMinuit) // // - the quantity to be minimized is defined in // MHDisp::Fill() and MHDisp::Finalize() // // - the parameters to be varied in the minimization are the parameters // appearing in the parametrization of Disp (MDispCalc::Calc()) // //------------------------------------------------------------------------ static void fcnDisp(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { cout << "entry fcnDisp" << endl; // save pointer to the MINUIT object (for optimizing Disp) for the case // that it is overwritten by the pointer to another Minuit object // TMinuit *savePointer = gMinuit; //------------------------------------------------------------- MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); // get needed Disp containers from the existing parList MDispParameters *dispparams = (MDispParameters*)plistfcn->FindObject("MDispParameters"); if (!dispparams) { gLog << "fcnDisp : MDispParameters object '" << "MDispParameters" << "' not found... aborting" << endl; return; } MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp"); if (!hdisp) { gLog << "fcnDisp : MHDisp object '" << "MHDisp" << "' not found... aborting" << endl; return; } MParameterD *minpar = (MParameterD*)plistfcn->FindObject("MinimizationParameter"); // // transfer current parameter values to MDispParameters // // Attention : npar is the number of variable parameters // not the total number of parameters // Double_t fMin, fEdm, fErrdef; Int_t fNpari, fNparx, fIstat; gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); dispparams->SetParameters(TArrayD(fNparx, par)); // checking that parameters have been properly set TArrayD checkparameters = dispparams->GetParameters(); gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ", " << fNparx << endl; gLog << "fcnDisp : i, par, checkparameters =" << endl; for (Int_t i=0; iEventloop(); tasklistfcn->PrintStatistics(0, kTRUE); //------------------------------------------- // get the Minimization parameter value for the current values // of the Disp parameters if (minpar) f = minpar->GetVal(); else { gLog << "fcnDisp : MParameterD object '" << "MinimizationParameter" << "' not found... aborting" << endl; return; } } // -------------------------------------------------------------------------- // // Default constructor. // MFindDisp::MFindDisp(MFDisp *fdisp, const char *name, const char *title) : fHowManyTrain(10000), fHowManyTest(10000) { fName = name ? name : "MFindDisp"; fTitle = title ? title : "Optimizer of Disp"; //--------------------------- // camera geometry is needed for conversion mm ==> degree fCam = new MGeomCamMagic; // initialize MFDisp filter cuts to default (no cutting) values // (done at MFDisp constructor) fDispFilter = fdisp; // matrices to contain the training/test samples fMatrixTrainCalc = new MHMatrix("MatrixTrainCalc"); fMatrixTrainHists = new MHMatrix("MatrixTrainHists"); fMatrixTestCalc = new MHMatrix("MatrixTestCalc"); fMatrixTestHists = new MHMatrix("MatrixTestHists"); // objects of MDispCalc where the first part of the matrices mapping is defined fDispCalcTrain = new MDispCalc("DispTrain","MDispParameters"); fDispCalcTest = new MDispCalc("DispTest","MDispParametersTest"); // objects of MHDisp where the second part of the matrices mapping is defined fHDispTrain = new MHDisp("DispTrain"); fHDispTest = new MHDisp("DispTest"); // define columns of matrices // Train matrix fDispCalcTrain->InitMapping(fMatrixTrainCalc); fHDispTrain->InitMapping(fMatrixTrainHists); // Test matrix fDispCalcTest->InitMapping(fMatrixTestCalc); fHDispTest->InitMapping(fMatrixTestHists); } // -------------------------------------------------------------------------- // // Default destructor. // MFindDisp::~MFindDisp() { delete fCam; delete fMatrixTrainCalc; delete fMatrixTrainHists; delete fMatrixTestCalc; delete fMatrixTestHists; delete fDispCalcTrain; delete fDispCalcTest; delete fHDispTrain; delete fHDispTest; } // -------------------------------------------------------------------------- // // Define the matrices 'fMatrixTrainCalc' and 'fMatrixTrainHists' // for the TRAINING sample // // alltogether 'howmanytrain' events are read from file 'nametrain'; // the events are selected according to a target distribution 'hreftrain' // // Bool_t MFindDisp::DefineTrainMatrix( const TString &nametrain, MH3 &hreftrain, const Int_t howmanytrain, const TString &filetrain, Int_t iflag) { if (nametrain.IsNull() || howmanytrain <= 0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "Fill TRAINING Matrices from file '" << nametrain << endl; *fLog << " select " << howmanytrain << " events " << endl; if (!hreftrain.GetHist().GetEntries()==0) { *fLog << " according to a distribution given by the MH3 object '" << hreftrain.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } *fLog << "=============================================" << endl; MParList plist; MTaskList tlist; // initialize display to check the reference distribution // for the event selection (just if iflag==1) MStatusDisplay *display = NULL; if (iflag) display = new MStatusDisplay; MReadMarsFile read("Events", nametrain); read.DisableAutoScheme(); // apply cuts for event selection MContinue contdisp(fDispFilter); contdisp.SetName("ContFilterSelector2"); // choose a reference distribution MFEventSelector2 seltrain(hreftrain); seltrain.SetNumMax(howmanytrain); seltrain.SetName("selectTrain"); MFillH filltraincalc(fMatrixTrainCalc); filltraincalc.SetFilter(&seltrain); filltraincalc.SetName("fillMatrixTrainCalc"); MFillH filltrainhists(fMatrixTrainHists); filltrainhists.SetFilter(&seltrain); filltrainhists.SetName("fillMatrixTrainHists"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrainCalc); plist.AddToList(fMatrixTrainHists); //****************************** // entries in MTaskList tlist.AddToList(&read); if (fDispFilter != NULL) tlist.AddToList(&contdisp); tlist.AddToList(&seltrain); tlist.AddToList(&filltraincalc); tlist.AddToList(&filltrainhists); //****************************** MProgressBar bar; MEvtLoop evtloop; if (display != NULL) evtloop.SetDisplay(display); evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrain"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) return kFALSE; tlist.PrintStatistics(0, kTRUE); // print the filled Training Matrices fMatrixTrainCalc->Print("SizeCols"); fMatrixTrainHists->Print("SizeCols"); // check that number of generated events is compatible with the resquested Int_t howmanygeneratedcalc = fMatrixTrainCalc->GetM().GetNrows(); if (TMath::Abs(howmanygeneratedcalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" << howmanygeneratedcalc << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } Int_t howmanygeneratedhists = fMatrixTrainHists->GetM().GetNrows(); if (TMath::Abs(howmanygeneratedhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" << howmanygeneratedhists << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } *fLog << "TRAINING Matrices were filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out training matrices if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrainCalc->Write(); fMatrixTrainHists->Write(); filetr.Close(); *fLog << "MFindDisp::DefineTrainMatrix; Training matrices were written onto file '" << filetrain << "'" << endl; } if (display != NULL) delete display; return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices 'fMatrixTestCalc' and 'fMatrixTestHists' // for the TEST sample // // alltogether 'howmanytest' events are read from file 'nametest' // the events are selected according to a target distribution 'hreftest' // // Bool_t MFindDisp::DefineTestMatrix( const TString &nametest, MH3 &hreftest, const Int_t howmanytest, const TString &filetest, Int_t iflag) { if (nametest.IsNull() || howmanytest<=0) return kFALSE; *fLog << "=============================================" << endl; *fLog << "Fill TEST Matrices from file '" << nametest << endl; *fLog << " select " << howmanytest << " events " << endl; if (!hreftest.GetHist().GetEntries()==0) { *fLog << " according to a distribution given by the MH3 object '" << hreftest.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } MParList plist; MTaskList tlist; // initialize display to check the reference distribution // for the event selection (just if iflag==1) MStatusDisplay *display = NULL; if (iflag) display = new MStatusDisplay; MReadMarsFile read("Events", nametest); read.DisableAutoScheme(); // apply cuts for event selection MContinue contdisp(fDispFilter); contdisp.SetName("ContFilterSelector2"); // choose a reference distribution MFEventSelector2 seltest(hreftest); seltest.SetNumMax(howmanytest); seltest.SetName("selectTest"); MFillH filltestcalc(fMatrixTestCalc); filltestcalc.SetFilter(&seltest); filltestcalc.SetName("fillMatrixTestCalc"); MFillH filltesthists(fMatrixTestHists); filltesthists.SetFilter(&seltest); filltesthists.SetName("fillMatrixTestHists"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTestCalc); plist.AddToList(fMatrixTestHists); //****************************** // entries in MTaskList tlist.AddToList(&read); if (fDispFilter != NULL) tlist.AddToList(&contdisp); tlist.AddToList(&seltest); tlist.AddToList(&filltestcalc); tlist.AddToList(&filltesthists); //****************************** MProgressBar bar; MEvtLoop evtloop; if (display != NULL) evtloop.SetDisplay(display); evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTest"); evtloop.SetProgressBar(&bar); if (!evtloop.Eventloop()) return kFALSE; tlist.PrintStatistics(0, kTRUE); // print the filled Test Matrices fMatrixTestCalc->Print("SizeCols"); fMatrixTestHists->Print("SizeCols"); // check that number of generated events is compatible with the resquested const Int_t howmanygeneratedcalc = fMatrixTestCalc->GetM().GetNrows(); if (TMath::Abs(howmanygeneratedcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" << howmanygeneratedcalc << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } const Int_t howmanygeneratedhists = fMatrixTestHists->GetM().GetNrows(); if (TMath::Abs(howmanygeneratedhists-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" << howmanygeneratedhists << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } *fLog << "TEST Matrices were filled" << endl; *fLog << "=============================================" << endl; //-------------------------- // write out test matrices if (filetest != "") { TFile filete(filetest, "RECREATE"); fMatrixTestCalc->Write(); fMatrixTestHists->Write(); filete.Close(); *fLog << "MFindDisp::DefineTestMatrix; Test matrices were written onto file '" << filetest << "'" << endl; } if (display != NULL) delete display; return kTRUE; } // -------------------------------------------------------------------------- // // Define the matrices for the TRAINING and TEST sample respectively // // - this is for the case that training and test samples are generated // from the same input file (which has the name 'name') // Bool_t MFindDisp::DefineTrainTestMatrix( const TString &name, MH3 &href, const Int_t howmanytrain, const Int_t howmanytest, const TString &filetrain, const TString &filetest, Int_t iflag) { *fLog << "=============================================" << endl; *fLog << "Fill TRAINING and TEST Matrices from file '" << name << endl; *fLog << " select " << howmanytrain << " training and " << howmanytest << " test events " << endl; if (!href.GetHist().GetEntries()==0) { *fLog << " according to a distribution given by the MH3 object '" << href.GetName() << "'" << endl; } else { *fLog << " randomly" << endl; } MParList plist; MTaskList tlist; // initialize display to check the reference distribution // for the event selection (just if iflag==1) MStatusDisplay *display = NULL; if (iflag) display = new MStatusDisplay; MReadMarsFile read("Events", name); read.DisableAutoScheme(); // apply cuts for event selection MContinue contdisp(fDispFilter); contdisp.SetName("ContFilterSelector2"); // choose a reference distribution MFEventSelector2 selector(href); selector.SetNumMax(howmanytrain+howmanytest); selector.SetName("selectTrainTest"); selector.SetInverted(); MContinue cont(&selector); cont.SetName("ContTrainTest"); // choose randomly the events to fill the Training Matrix Double_t prob = ( (Double_t) howmanytrain ) / ( (Double_t)(howmanytrain+howmanytest) ); MFEventSelector split; split.SetSelectionRatio(prob); MFillH filltraincalc(fMatrixTrainCalc); filltraincalc.SetFilter(&split); filltraincalc.SetName("fillMatrixTrainCalc"); MFillH filltrainhists(fMatrixTrainHists); filltrainhists.SetFilter(&split); filltrainhists.SetName("fillMatrixTrainHists"); // consider this event as candidate for a test event // only if event was not accepted as a training event MContinue conttrain(&split); conttrain.SetName("ContTrain"); MFillH filltestcalc(fMatrixTestCalc); filltestcalc.SetName("fillMatrixTestCalc"); MFillH filltesthists(fMatrixTestHists); filltesthists.SetName("fillMatrixTestHists"); //****************************** // entries in MParList plist.AddToList(&tlist); plist.AddToList(fCam); plist.AddToList(fMatrixTrainCalc); plist.AddToList(fMatrixTrainHists); plist.AddToList(fMatrixTestCalc); plist.AddToList(fMatrixTestHists); //****************************** // entries in MTaskList tlist.AddToList(&read); if (fDispFilter != NULL) tlist.AddToList(&contdisp); tlist.AddToList(&cont); tlist.AddToList(&filltraincalc); tlist.AddToList(&filltrainhists); tlist.AddToList(&conttrain); tlist.AddToList(&filltestcalc); tlist.AddToList(&filltesthists); //****************************** MProgressBar bar; MEvtLoop evtloop; if (display != NULL) evtloop.SetDisplay(display); evtloop.SetParList(&plist); evtloop.SetName("EvtLoopMatrixTrainTest"); evtloop.SetProgressBar(&bar); Int_t maxev = -1; if (!evtloop.Eventloop(maxev)) return kFALSE; tlist.PrintStatistics(0, kTRUE); // print the filled Train Matrices fMatrixTrainCalc->Print("SizeCols"); fMatrixTrainHists->Print("SizeCols"); // check that number of generated events is compatible with the resquested const Int_t generatedtraincalc = fMatrixTrainCalc->GetM().GetNrows(); if (TMath::Abs(generatedtraincalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" << generatedtraincalc << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } const Int_t generatedtrainhists = fMatrixTrainHists->GetM().GetNrows(); if (TMath::Abs(generatedtrainhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" << generatedtrainhists << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } // print the filled Test Matrices fMatrixTestCalc->Print("SizeCols"); fMatrixTestHists->Print("SizeCols"); // check that number of generated events is compatible with the resquested const Int_t generatedtestcalc = fMatrixTestCalc->GetM().GetNrows(); if (TMath::Abs(generatedtestcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" << generatedtestcalc << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } const Int_t generatedtesthists = fMatrixTestHists->GetM().GetNrows(); if (TMath::Abs(generatedtesthists-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" << generatedtesthists << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } *fLog << "TRAINING and TEST Matrices were filled" << endl; *fLog << "=============================================" << endl; //---------------------------- // write out training matrices if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); fMatrixTrainCalc->Write(); fMatrixTrainHists->Write(); filetr.Close(); *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrices were written onto file '" << filetrain << "'" << endl; } //-------------------------- // write out test matrices if (filetest != "") { TFile filete(filetest, "RECREATE"); fMatrixTestCalc->Write(); fMatrixTestHists->Write(); filete.Close(); *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrices were written onto file '" << filetest << "'" << endl; } if (display != NULL) delete display; return kTRUE; } // -------------------------------------------------------------------------- // // Read training and test matrices from files // // Bool_t MFindDisp::ReadMatrix(const TString &filetrain, const TString &filetest) { //-------------------------- // read in training matrices TFile filetr(filetrain); fMatrixTrainCalc->Read("MatrixTrainCalc"); fMatrixTrainHists->Read("MatrixTrainHists"); fMatrixTrainCalc->Print("SizeCols"); fMatrixTrainHists->Print("SizeCols"); *fLog << "MFindDisp::ReadMatrix; Training matrices were read in from file '" << filetrain << "'" << endl; filetr.Close(); //-------------------------- // read in test matrices TFile filete(filetest); fMatrixTestCalc->Read("MatrixTestCalc"); fMatrixTestHists->Read("MatrixTestHists"); fMatrixTestCalc->Print("SizeCols"); fMatrixTestHists->Print("SizeCols"); *fLog << "MFindDisp::ReadMatrix; Test matrices were read in from file '" << filetest << "'" << endl; filete.Close(); return kTRUE; } //------------------------------------------------------------------------ // // Steering program for optimizing Disp // ------------------------------------ // // the criterion for the 'optimum' is defined in // MHDisp::Fill() and MHDisp::Finalize(); // for example : smallest sum (over all events) of d^2, where d is the // distance between the estimated source position // (calculated using the current value of Disp) and // the true source position // // The various steps are : // // - setup the event loop to be executed for each call to fcnDisp // // - call TMinuit to do the minimization : // the fcnDisp function calculates the parameter to minimize // for the current Disp parameter values; // for this - Disp is calculated in the event loop by calling // MDispCalc::Process() ==> MDispCalc::Calc() // - the Minimization parameter contributions are summed up // in the event loop by calling MHDisp::Fill() // - after the event loop the final value of the Minimization // parameter is calculated by calling MHDisp::Finalize() // // Needed as input : (to be set by the Set functions) // // - fFilenameParam name of file to which optimum values of the // parameters are written // // - for the minimization, the starting values of the parameters are taken // - from the file parDispInit (if it is != "") // - or from the arrays params and/or steps // - or from the DispParameters constructor // //---------------------------------------------------------------------- Bool_t MFindDisp::FindParams(TString parDispInit, TArrayD ¶ms, TArrayD &steps) { // Setup the event loop which will be executed in the // fcnDisp function of MINUIT // // parDispInit is the name of the file containing the initial values // of the parameters; // if parDispInit = "" 'params' and 'steps' are taken as initial values // *fLog << "=============================================" << endl; *fLog << "Setup event loop for fcnDisp" << endl; if (fMatrixTrainCalc == NULL || fMatrixTrainHists == NULL) { *fLog << "MFindDisp::FindParams; training matrices are not defined... aborting" << endl; return kFALSE; } if (fMatrixTrainCalc->GetM().GetNrows() <= 0 || fMatrixTrainHists->GetM().GetNrows() <= 0) { *fLog << "MFindDisp::FindParams; training matrices have no entries" << endl; return kFALSE; } //--------------------------------------------------------- MParList parlistfcn; MTaskList tasklistfcn; // loop over rows of matrix MMatrixLoop loopcalc(fMatrixTrainCalc); MMatrixLoop loophists(fMatrixTrainHists); //-------------------------------- // create container for the Disp parameters // and set them to their initial values MDispParameters *dispparams = fDispCalcTrain->GetDispParameters(); // take initial values from file parDispInit if (parDispInit != "") { TFile inparam(parDispInit); dispparams->Read("MDispParameters"); inparam.Close(); *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file " << parDispInit << endl; } // take initial values from 'params' and/or 'steps' else if (params.GetSize() != 0 || steps.GetSize() != 0 ) { if (params.GetSize() != 0) { *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'" << endl; dispparams->SetParameters(params); } if (steps.GetSize() != 0) { *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'" << endl; dispparams->SetStepsizes(steps); } } else { *fLog << "MFindDisp::FindParams; initial values and step sizes are taken " << "from the MDispParameters constructor" << endl; } // fill the plots for Disp and sum up the Minimization parameter contributions MFillH filldispplots("MHDisp", ""); //****************************** // entries in MParList parlistfcn.AddToList(&tasklistfcn); parlistfcn.AddToList(dispparams); parlistfcn.AddToList(fHDispTrain); parlistfcn.AddToList(fCam); parlistfcn.AddToList(fMatrixTrainCalc); parlistfcn.AddToList(fMatrixTrainHists); //****************************** // entries in MTaskList tasklistfcn.AddToList(&loopcalc); tasklistfcn.AddToList(&loophists); tasklistfcn.AddToList(fDispCalcTrain); tasklistfcn.AddToList(&filldispplots); //****************************** MEvtLoop evtloopfcn("EvtLoopFCN"); evtloopfcn.SetParList(&parlistfcn); *fLog << "Event loop for fcnDisp has been setup" << endl; // address of evtloopfcn is used in CallMinuit() //----------------------------------------------------------------------- // //---------- Start of minimization part -------------------- // // Do the minimization with MINUIT // // Be careful: This is not thread safe // *fLog << "========================================================" << endl; *fLog << "Start minimization for Disp" << endl; // ------------------------------------------- // prepare call to MINUIT // // get initial values of parameters fVinit = dispparams->GetParameters(); fStep = dispparams->GetStepsizes(); TString name[fVinit.GetSize()]; fStep.Set(fVinit.GetSize()); fLimlo.Set(fVinit.GetSize()); fLimup.Set(fVinit.GetSize()); fFix.Set(fVinit.GetSize()); fNpar = fVinit.GetSize(); // define names, step sizes, lower and upper limits of the parameters // fFix[] = 0; means the parameter is not fixed for (UInt_t i=0; iWrite(); outparam.Close(); *fLog << "Optimum parameter values for Disp were written onto file '" << fFilenameParam << "' :" << endl; const TArrayD &check = dispparams->GetParameters(); for (Int_t i=0; iDraw(); *fLog << "End of Disp optimization part" << endl; *fLog << "======================================================" << endl; return kTRUE; } // ----------------------------------------------------------------------- // // Test the result of the Disp optimization on the test sample // Bool_t MFindDisp::TestParams() { if (fMatrixTestCalc == NULL || fMatrixTestHists == NULL) { *fLog << "MFindDisp::TestParams; test matrices are not defined... aborting" << endl; return kFALSE; } if (fMatrixTestCalc->GetM().GetNrows() <= 0 || fMatrixTestHists->GetM().GetNrows() <= 0) { *fLog << "MFindDisp::TestParams; test matrices have no entries" << endl; return kFALSE; } // ------------- TEST the Disp optimization ------------------- // *fLog << "Test the Disp optimization on the test sample" << endl; // ----------------------------------------------------------------- // read optimum parameter values from file fFilenameParam // into array 'dispPar' TFile inparam(fFilenameParam); MDispParameters dispin; dispin.Read("MDispParameters"); inparam.Close(); *fLog << "Optimum parameter values for Disp were read from file '"; *fLog << fFilenameParam << "' :" << endl; const TArrayD &dispPar = dispin.GetParameters(); for (Int_t i=0; iGetDispParameters(); dispparams->SetParameters(dispPar); MMatrixLoop loopcalc(fMatrixTestCalc); MMatrixLoop loophists(fMatrixTestHists); MHMatrix *imatrix = NULL; TArrayI imap; fHDispTest->GetMatrixMap(imatrix,imap); // attention: argument of MHDisp is the name of MImageParDisp container, that should // be the same than the name given to it when creating MDispCalc object at the MFindDisp // constructor: fDispCalcTrain = new MDispCalc("DispTest"); // fill the plots for Disp and sum up the Minimization parameter contributions MHDisp hdisp1("DispTest"); hdisp1.SetName("MHDispCorr"); hdisp1.SetSelectedPos(1); hdisp1.SetMatrixMap(imatrix,imap); MFillH filldispplots1("MHDispCorr[MHDisp]", ""); MHDisp hdisp2("DispTest"); hdisp2.SetName("MHDispWrong"); hdisp2.SetSelectedPos(2); hdisp2.SetMatrixMap(imatrix,imap); MFillH filldispplots2("MHDispWrong[MHDisp]", ""); MHDisp hdisp3("DispTest"); hdisp3.SetName("MHDispM3Long"); hdisp3.SetSelectedPos(3); hdisp3.SetMatrixMap(imatrix,imap); MFillH filldispplots3("MHDispM3Long[MHDisp]", ""); MHDisp hdisp4("DispTest"); hdisp4.SetName("MHDispAsym"); hdisp4.SetSelectedPos(4); hdisp4.SetMatrixMap(imatrix,imap); MFillH filldispplots4("MHDispAsym[MHDisp]", ""); //****************************** // entries in MParList parlist2.AddToList(&tasklist2); parlist2.AddToList(dispparams); parlist2.AddToList(&hdisp1); parlist2.AddToList(&hdisp2); parlist2.AddToList(&hdisp3); parlist2.AddToList(&hdisp4); parlist2.AddToList(fCam); parlist2.AddToList(fMatrixTestCalc); parlist2.AddToList(fMatrixTestHists); //****************************** // entries in MTaskList tasklist2.AddToList(&loopcalc); tasklist2.AddToList(&loophists); tasklist2.AddToList(fDispCalcTest); tasklist2.AddToList(&filldispplots1); tasklist2.AddToList(&filldispplots2); tasklist2.AddToList(&filldispplots3); tasklist2.AddToList(&filldispplots4); //****************************** MProgressBar bar2; MEvtLoop evtloop2; evtloop2.SetParList(&parlist2); evtloop2.SetName("EvtLoopTestParams"); evtloop2.SetProgressBar(&bar2); Int_t maxevents2 = -1; if (!evtloop2.Eventloop(maxevents2)) return kFALSE; tasklist2.PrintStatistics(0, kTRUE); //------------------------------------------- // draw the plots hdisp1.Draw(); hdisp2.Draw(); hdisp3.Draw(); hdisp4.Draw(); //------------------------------------------- *fLog << "Test of Disp otimization finished" << endl; *fLog << "======================================================" << endl; return kTRUE; }