diff --git a/App/App.pro b/App/App.pro index e402ee3bfcb89c8aba98a376bf6d09d54a988019..7c46ed6a8a4685d46ce0f477b3e81fee6f3904cb 100644 --- a/App/App.pro +++ b/App/App.pro @@ -16,8 +16,7 @@ SOURCES += \ src/AppOptionsDescription.cpp \ src/CommandLine.cpp \ src/DrawHelper.cpp \ - src/EventFrame.cpp \ - src/FittingHelper.cpp \ + src/FitSuiteHelper.cpp \ src/FunctionalTestFactory.cpp \ src/IFunctionalTest.cpp \ src/IsGISAXSTools.cpp \ @@ -43,7 +42,8 @@ SOURCES += \ src/TestMultiLayerRoughness.cpp \ src/TestPerformance.cpp \ src/TestRootTree.cpp \ - src/TestRoughness.cpp + src/TestRoughness.cpp \ + src/TreeEventStructure.cpp HEADERS += \ inc/App.h \ @@ -51,9 +51,8 @@ HEADERS += \ inc/AppOptionsDescription.h \ inc/CommandLine.h \ inc/DrawHelper.h \ - inc/EventFrame.h \ + inc/FitSuiteHelper.h \ inc/FunctionalTestFactory.h \ - inc/FittingHelper.h \ inc/IFunctionalTest.h \ inc/IsGISAXSTools.h \ inc/ROOTMinimizer.h \ @@ -77,7 +76,8 @@ HEADERS += \ inc/TestMultiLayerRoughness.h \ inc/TestPerformance.h \ inc/TestRootTree.h \ - inc/TestRoughness.h + inc/TestRoughness.h \ + inc/TreeEventStructure.h INCLUDEPATH += ./inc ../Core/Algorithms/inc ../Core/FormFactors/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc DEPENDPATH += ./inc ../Core/Algorithms/inc ../Core/FormFactors/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc diff --git a/App/inc/App.h b/App/inc/App.h index e6244f69a1a236013d57642763b7336aa78dabb6..dee1d9fdea64c5d952db712bca0237ad9570da79 100644 --- a/App/inc/App.h +++ b/App/inc/App.h @@ -16,6 +16,6 @@ #include "ISingleton.h" #include "DrawHelper.h" -#include "EventFrame.h" +#include "TreeEventStructure.h" #endif // APP_H diff --git a/App/inc/AppLinkDef.h b/App/inc/AppLinkDef.h index b297e139e7776c636a28331a0f57e29fb4cf504a..e93d4c2ad0f48604f0cb03fa7cac90b3ba5e5ca4 100644 --- a/App/inc/AppLinkDef.h +++ b/App/inc/AppLinkDef.h @@ -20,7 +20,7 @@ #pragma link C++ class ISingleton<DrawHelper>+; #pragma link C++ class DrawHelper+; - -#pragma link C++ class EventFrame+; +#pragma link C++ class TreeEventOutputData+; +#pragma link C++ class TreeEventFitData+; #endif diff --git a/App/inc/FittingHelper.h b/App/inc/FitSuiteHelper.h similarity index 55% rename from App/inc/FittingHelper.h rename to App/inc/FitSuiteHelper.h index 81a3b2eb809a2074d60c35fcb4ac005af9f3597a..7000380a897ee34dea4a40310a259fd0ef758e19 100644 --- a/App/inc/FittingHelper.h +++ b/App/inc/FitSuiteHelper.h @@ -1,5 +1,5 @@ -#ifndef FITTINGHELPER_H -#define FITTINGHELPER_H +#ifndef FITSUITEHELPER_H +#define FITSUITEHELPER_H // ******************************************************************** // * The BornAgain project * // * Simulation of neutron and x-ray scattering at grazing incidence * @@ -9,7 +9,7 @@ // * eget quam orci. Quisque porta varius dui, quis posuere nibh * // * mollis quis. Mauris commodo rhoncus porttitor. * // ******************************************************************** -//! @file FittingHelper.h +//! @file FitSuiteHelper.h //! @brief Collection of helper classes for fitting from App //! @author Scientific Computing Group at FRM II //! @date 08.10.2012 @@ -18,17 +18,15 @@ #include "IObserver.h" #include <string> -namespace FittingHelper { - //- ------------------------------------------------------------------- -//! @class ObserveAndDraw -//! @brief Draw fit progress at the end of each FitSuite's minimization function call +//! @class FitSuiteObserverDraw +//! @brief Draw fit progress at the end of each FitSuite's iteration //- ------------------------------------------------------------------- -class ObserveAndDraw : public IObserver +class FitSuiteObserverDraw : public IObserver { public: - ObserveAndDraw(std::string canvas_name) : m_ncall(0), m_canvas_name(canvas_name) {} + FitSuiteObserverDraw(std::string canvas_name) : m_ncall(0), m_canvas_name(canvas_name) {} void update(IObservable *subject); private: int m_ncall; //! number of call of given observer @@ -36,7 +34,22 @@ private: }; -} +//- ------------------------------------------------------------------- +//! @class FitSuiteObserverWrite +//! @brief Save results of each fit iteration to the disk in the form of ROOT tree +//! If tree exist, data will be appended to it +//- ------------------------------------------------------------------- +class FitSuiteObserverWriteTree : public IObserver +{ +public: + FitSuiteObserverWriteTree(std::string file_name) : m_ncall(0), m_file_name(file_name) {} + void update(IObservable *subject); +private: + int m_ncall; + std::string m_file_name; //! canvas name were to draw +}; + + -#endif // FITTINGHELPER_H +#endif // FITSUITEHELPER_H diff --git a/App/inc/ROOTMinimizer.h b/App/inc/ROOTMinimizer.h index 94ea46df617b34b8947d33e8bbbd7099408898e6..3d7d86717eda6940771cda97cc8c4f92616fad9b 100644 --- a/App/inc/ROOTMinimizer.h +++ b/App/inc/ROOTMinimizer.h @@ -44,7 +44,7 @@ public: //! return minimum function value virtual double getMinValue() { return m_root_minimizer->MinValue(); } - //! return pointer to the parameters values at the minimum + //! return value of variable corresponding the minimum of the function virtual double getValueOfVariableAtMinimum(size_t i) { if(i >= m_root_minimizer->NDim() ) throw OutOfBoundsException("ROOTMinimizer::getVariableAtMinimum() -> Wrong number of the variable"); return m_root_minimizer->X()[i]; diff --git a/App/inc/TestMiscellaneous.h b/App/inc/TestMiscellaneous.h index 6cc8beba26a11c9df48451fd6f1e81490c6b9617..fc1f0a769f72e7dd7ecd7664147507862dbf2abc 100644 --- a/App/inc/TestMiscellaneous.h +++ b/App/inc/TestMiscellaneous.h @@ -44,6 +44,9 @@ public: //! test of reading of OutputData from ASCII file void test_OutputDataIOFactory(); + //! test kvector container + void test_KVectorContainer(); + }; #endif // TESTMISCELLANEOUS_H diff --git a/App/inc/EventFrame.h b/App/inc/TreeEventStructure.h similarity index 51% rename from App/inc/EventFrame.h rename to App/inc/TreeEventStructure.h index 25f582eedfd7befbcc6cedd07b0211f9086b6980..8782679774af94bc3bb00f4f29cf5d33ad1ddb65 100644 --- a/App/inc/EventFrame.h +++ b/App/inc/TreeEventStructure.h @@ -1,5 +1,5 @@ -#ifndef EVENTFRAME_H -#define EVENTFRAME_H +#ifndef TREEEVENTSTRUCTURE_H +#define TREEEVENTSTRUCTURE_H // ******************************************************************** // * The BornAgain project * // * Simulation of neutron and x-ray scattering at grazing incidence * @@ -9,24 +9,27 @@ // * eget quam orci. Quisque porta varius dui, quis posuere nibh * // * mollis quis. Mauris commodo rhoncus porttitor. * // ******************************************************************** -//! @file EventFrame.h -//! @brief Definition of EventFrame class for writing in root files +//! @file TreeEventStructure.h +//! @brief Collection of classes to write output data in ROOT trees //! @author Scientific Computing Group at FRM II -//! @date 19.07.2012 - +//! @date 09.10.2012 #include "TObject.h" #include <vector> +#include <string> + + //- ------------------------------------------------------------------- -//! @class EventFrame -//! @brief event structure to save in the root file +//! @class TreeEventOutputData +//! @brief Data structure respresenting OutputData and mesocrystal settings +//! for writing/reading in/from the ROOT tree. //- ------------------------------------------------------------------- -class EventFrame +class TreeEventOutputData { public: - EventFrame(); - virtual ~EventFrame(){} + TreeEventOutputData(); + virtual ~TreeEventOutputData() { } void clear(); @@ -49,9 +52,35 @@ public: std::vector<std::vector<double > > valpha_f; // values of alpha_f for the frame std::vector<std::vector<double > > vphi_f; // values of phi_f for the frames - ClassDef(EventFrame,1) + ClassDef(TreeEventOutputData,1) +}; + + +//- ------------------------------------------------------------------- +//! @class TreeEventFitData +//! @brief Represent fit results after each iteration for writing/reading +//! in/from the ROOT tree +//- ------------------------------------------------------------------- +class TreeEventFitData +{ +public: + TreeEventFitData(); + virtual ~TreeEventFitData() { } + void clear(); + + int niter; // number of iteration + std::vector<std::vector<double > > real_data; // real data + std::vector<std::vector<double > > fit_data; // current fit iteration data + std::vector<std::vector<double > > diff; // chi2 difference between real and fit data + std::vector<std::vector<double > > axis0; // values of phi_f (made in 2D for convenient drawing) + std::vector<std::vector<double > > axis1; // values of alpha_f (made in 2D for convenient drawing) + double chi2; // current value of the function to minimize + std::vector<double> parvalues; // vector of minimization parameters + std::vector<std::string> parnames; // names of parameters + + ClassDef(TreeEventFitData,1) }; -#endif // EventFrame_H +#endif // TREEEVENTSTRUCTURE_H diff --git a/App/src/DrawHelper.cpp b/App/src/DrawHelper.cpp index 7b7018b69e6a989af1bb09118ef245aea2779761..eb55ce722e3d43c9836c493a814833b9ba201503 100644 --- a/App/src/DrawHelper.cpp +++ b/App/src/DrawHelper.cpp @@ -319,7 +319,7 @@ void DrawHelper::DrawMultilayer(const MultiLayer *sample) TCanvas *DrawHelper::createAndRegisterCanvas(std::string name, std::string title, int xsize, int ysize) { if(xsize==0) xsize = m_default_canvas_xsize; - if(ysize==0) xsize = m_default_canvas_ysize; + if(ysize==0) ysize = m_default_canvas_ysize; TCanvas *c1 = new TCanvas(name.c_str(), title.c_str(), xsize, ysize); SetMagnifier(c1); diff --git a/App/src/EventFrame.cpp b/App/src/EventFrame.cpp deleted file mode 100644 index e291a06a9b0b09d63ba0cfd5d17f217a2097f48d..0000000000000000000000000000000000000000 --- a/App/src/EventFrame.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include "EventFrame.h" - -EventFrame::EventFrame() -{ - clear(); -} - -void EventFrame::clear() -{ - nframe = 0; - alpha_i = 0; - phi_i = 0; - - nphi_f = 0; - phi_f_min = 0; - phi_f_max = 0; - - nalpha_f = 0; - alpha_f_min = 0; - alpha_f_max = 0; - - malpha = 0; - mphi = 0; - npR = 0; - - valpha_f.clear(); - vphi_f.clear(); - vi.clear(); -} diff --git a/App/src/FitSuiteHelper.cpp b/App/src/FitSuiteHelper.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ae42c0d5233f26a9595eccf0d71022fe79b727ab --- /dev/null +++ b/App/src/FitSuiteHelper.cpp @@ -0,0 +1,114 @@ +#include "FitSuite.h" +#include "FitSuiteHelper.h" +#include "TreeEventStructure.h" +#include "IsGISAXSTools.h" + +#include "TCanvas.h" +#include "TPaveText.h" +#include "ChiSquaredModule.h" +#include "TROOT.h" +#include "TFile.h" +#include "TTree.h" + + +/* ************************************************************************* */ +// draw results of fit iteration in ROOT's canvas +/* ************************************************************************* */ +void FitSuiteObserverDraw::update(IObservable *subject) +{ + FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject); + if( !fitSuite ) throw NullPointerException("FitSuiteObserverDraw::update() -> Error! Can't cast FitSuite"); + + std::cout << "FitObserver: " << " ncall" << m_ncall << " chi2:" << fitSuite->getChiSquaredModule()->getValue(); + // printing parameter values + for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { + std::cout << " " << (*it)->getName() << " " << (*it)->getValue(); + // std::cout << *(*it); + } + std::cout << std::endl; + + TCanvas *c1 = dynamic_cast<TCanvas *>( gROOT->FindObject(m_canvas_name.c_str()) ); + if(!c1) throw NullPointerException("FitSuiteObserverDraw::update() -> No access to canvas"); + c1->cd(3); + gPad->SetLogz(); + IsGISAXSTools::drawOutputDataInPad(*fitSuite->getChiSquaredModule()->getSimulationData(), "CONT4 Z", "current simulated data"); + c1->cd(4); + IsGISAXSTools::drawOutputDataDifference2D(*fitSuite->getChiSquaredModule()->getSimulationData(), *fitSuite->getChiSquaredModule()->getRealData(), "CONT4 Z", "difference"); + c1->cd(5); + + TPaveText *pt = new TPaveText(.05,.1,.95,.8); + char str[256]; + sprintf(str,"Iteration %d",m_ncall); + pt->AddText(str); + sprintf(str,"chi2 %e",fitSuite->getChiSquaredModule()->getValue()); + pt->AddText(str); + //pt->AddLine(.0,.5,1.,.5); + for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { + sprintf(str,"%s %f", (*it)->getName().c_str(), (*it)->getValue()); + pt->AddText(str); + } + pt->Draw(); + + m_ncall++; +} + + +/* ************************************************************************* */ +// Save results of each fit iteration to the disk in the form of ROOT tree +/* ************************************************************************* */ +void FitSuiteObserverWriteTree::update(IObservable *subject) +{ + std::string tree_name("FitSuiteTree"); + + FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject); + if( !fitSuite ) throw NullPointerException("FitSuiteObserverWriteTree::update() -> Error! Can't cast FitSuite"); + + // preparing root file for writing + // if it is first call the file will be opened in 'recreate' mode, otherwise in 'update' mode + TFile *top(0); + if(m_ncall == 0) { + top = new TFile(m_file_name.c_str(),"RECREATE"); + } else { + top = new TFile(m_file_name.c_str(),"UPDATE"); + } + if( !top->IsOpen() ) { + throw RuntimeErrorException("FitSuiteObserverWriteTree::update() -> Can't open file "+ m_file_name + " for writing"); + } + + + // data object to write in the tree + TreeEventFitData *event = new TreeEventFitData(); + + // creating new tree + TTree *tree = dynamic_cast<TTree *>(top->Get(tree_name.c_str())); + if( tree == 0) { + // tree doesn't exist due to new file, making new tree + tree = new TTree(tree_name.c_str(),"Oh, my data"); + tree->Branch("Event",&event,16000,2); + } else { + // tree exists, pointing it to the new data + tree->SetBranchAddress("Event", &event); + } + + // filling data object with data from FitSuite + const OutputData<double > *real_data = fitSuite->getChiSquaredModule()->getRealData(); + const OutputData<double > *simu_data = fitSuite->getChiSquaredModule()->getSimulationData(); + IsGISAXSTools::exportOutputDataInVectors2D(*real_data, event->real_data, event->axis0, event->axis1); + IsGISAXSTools::exportOutputDataInVectors2D(*simu_data, event->fit_data, event->axis0, event->axis1); + event->chi2 = fitSuite->getChiSquaredModule()->getValue(); + for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { + event->parvalues.push_back( (*it)->getValue() ); + event->parnames.push_back( (*it)->getName().c_str() ); + } + event->niter = m_ncall; + + // appending data to the tree + tree->Fill(); + tree->Write(tree_name.c_str(), TObject::kOverwrite); + top->Close(); + delete top; // there is no need to delete tree since ROOT file takes care about all objects opened afterwards + delete event; + + m_ncall++; +} + diff --git a/App/src/FittingHelper.cpp b/App/src/FittingHelper.cpp deleted file mode 100644 index 320daa8cadfbec306d314d94f2caf94a941dc780..0000000000000000000000000000000000000000 --- a/App/src/FittingHelper.cpp +++ /dev/null @@ -1,48 +0,0 @@ -#include "FittingHelper.h" -#include "FitSuite.h" -#include "IsGISAXSTools.h" - -#include "TCanvas.h" -#include "TPaveText.h" -#include "ChiSquaredModule.h" -#include "TROOT.h" - -void FittingHelper::ObserveAndDraw::update(IObservable *subject) -{ - FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject); - if( !fitSuite ) throw NullPointerException("LocalFitObserver::IObserver() -> Error! Can't cast FitSuite"); - - std::cout << "FitObserver: " << " ncall" << m_ncall << " chi2:" << fitSuite->getChiSquaredModule()->getValue(); - // printing parameter values - for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { - std::cout << " " << (*it)->getName() << " " << (*it)->getValue(); - // std::cout << *(*it); - } - std::cout << std::endl; - - TCanvas *c1 = dynamic_cast<TCanvas *>( gROOT->FindObject(m_canvas_name.c_str()) ); - if(!c1) throw NullPointerException("LocalFitIbserver::update() -> No access to canvas"); - c1->cd(3); - gPad->SetLogz(); - IsGISAXSTools::drawOutputDataInPad(*fitSuite->getChiSquaredModule()->getSimulationData(), "CONT4 Z", "current simulated data"); - c1->cd(4); - IsGISAXSTools::drawOutputDataDifference2D(*fitSuite->getChiSquaredModule()->getSimulationData(), *fitSuite->getChiSquaredModule()->getRealData(), "CONT4 Z", "difference"); - c1->cd(5); - - TPaveText *pt = new TPaveText(.05,.1,.95,.8); - char str[256]; - sprintf(str,"Iteration %d",m_ncall); - pt->AddText(str); - sprintf(str,"chi2 %e",fitSuite->getChiSquaredModule()->getValue()); - pt->AddText(str); - //pt->AddLine(.0,.5,1.,.5); - for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { - sprintf(str,"%s %f", (*it)->getName().c_str(), (*it)->getValue()); - pt->AddText(str); - } - pt->Draw(); - - m_ncall++; -} - - diff --git a/App/src/TestFittingModule.cpp b/App/src/TestFittingModule.cpp index 51ffb8e4e14f213949b19bea0508a3b06171542c..60f0b8eeadb804176e4b6cc8e2c116a98b9b248b 100644 --- a/App/src/TestFittingModule.cpp +++ b/App/src/TestFittingModule.cpp @@ -12,7 +12,7 @@ #include "FormFactors.h" #include "Exceptions.h" #include "DrawHelper.h" -#include "FittingHelper.h" +#include "FitSuiteHelper.h" #include "IObserver.h" #include "FitSuite.h" @@ -80,10 +80,14 @@ void TestFittingModule::execute() fitSuite->addFitParameter("*/MultiLayer/Layer0/thickness", 12*Units::nanometer, 1*Units::nanometer, TRange<double>(1.0, 20.0) ); fitSuite->addFitParameter("*/FormFactorCylinder/radius", 2*Units::nanometer, 1*Units::nanometer, TRange<double>(1.0, 20.0) ); - FittingHelper::ObserveAndDraw *observer = new FittingHelper::ObserveAndDraw(canvas_name); - fitSuite->attachObserver(observer); + FitSuiteObserverDraw *drawObserver = new FitSuiteObserverDraw(canvas_name); + fitSuite->attachObserver(drawObserver); + FitSuiteObserverWriteTree *writeObserver = new FitSuiteObserverWriteTree("fitsuite.root"); + fitSuite->attachObserver(writeObserver); fitSuite->runFit(); + delete drawObserver; + delete writeObserver; std::cout << "------ RESULTS ---------" << std::endl; std::cout << "FitSuite > MinValue:" << fitSuite->getMinimizer()->getMinValue() << " " << fitSuite->getMinimizer()->getValueOfVariableAtMinimum(0) << std::endl; diff --git a/App/src/TestMesoCrystal2.cpp b/App/src/TestMesoCrystal2.cpp index b37c8b8653e4abe372106936aa6a1566479c4a6b..bd3e78e4002eabcf9ba00e8baaa055fc7b964b0a 100644 --- a/App/src/TestMesoCrystal2.cpp +++ b/App/src/TestMesoCrystal2.cpp @@ -24,7 +24,7 @@ #include "SampleFactory.h" #include "FitMultiParameter.h" #include "TRange.h" -#include "FittingHelper.h" +#include "FitSuiteHelper.h" #include "TCanvas.h" #include "TH2D.h" @@ -83,9 +83,9 @@ void TestMesoCrystal2::execute() fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") ); fitSuite->addFitParameter("/MultiLayer/LayerInterface1/roughness/sigma", 0.1, 0.05, TRange<double>(0.0, 10.0) ); - FittingHelper::ObserveAndDraw *observer = new FittingHelper::ObserveAndDraw(canvas_name); - fitSuite->attachObserver(observer); - + FitSuiteObserverDraw *drawObserver = new FitSuiteObserverDraw(canvas_name); + fitSuite->attachObserver(drawObserver); +\ fitSuite->runFit(); for(FitSuite::fitparameters_t::iterator it = fitSuite->fitparams_begin(); it!=fitSuite->fitparams_end(); ++it) { diff --git a/App/src/TestMiscellaneous.cpp b/App/src/TestMiscellaneous.cpp index 63214676990e5bc05dfa6a9113bf8999de34564a..79adc2c0d17a639a1fa26cb31021c4431c7b6ece 100644 --- a/App/src/TestMiscellaneous.cpp +++ b/App/src/TestMiscellaneous.cpp @@ -16,6 +16,7 @@ #include "MathFunctions.h" #include "OutputDataIOFactory.h" #include "Utils.h" +#include "Types.h" #include "TGraph.h" #include "TH2D.h" @@ -34,7 +35,8 @@ TestMiscellaneous::TestMiscellaneous() void TestMiscellaneous::execute() { - test_OutputDataIOFactory(); + test_KVectorContainer(); + //test_OutputDataIOFactory(); //test_FastSin(); //test_DoubleToComplexInterpolatingFunction(); //test_FormFactor(); @@ -42,6 +44,20 @@ void TestMiscellaneous::execute() } +/* ************************************************************************* */ +// test of reading of OutputData from ASCII file +/* ************************************************************************* */ +void TestMiscellaneous::test_KVectorContainer() +{ + KVectorContainer cc; + + for(size_t i=0; i<100; ++i) { + cc.push_back(kvector_t(i,0,0)); + cc.print(); + } +} + + /* ************************************************************************* */ // test of reading of OutputData from ASCII file diff --git a/App/src/TestRootTree.cpp b/App/src/TestRootTree.cpp index 496025111cb30c512907bee26db4acc0f42420f9..732d8b178231acd8a042b76edbaa320d1df0639e 100644 --- a/App/src/TestRootTree.cpp +++ b/App/src/TestRootTree.cpp @@ -16,7 +16,7 @@ #include "TRandom.h" #include "TCanvas.h" #include "IsGISAXSTools.h" -#include "EventFrame.h" +#include "TreeEventStructure.h" #include "TestMesoCrystal1.h" #include <vector> @@ -71,7 +71,7 @@ void TestRootTree::complex_write() // creating new tree TTree *tree = new TTree(tree_name.c_str(),"Oh, my data"); - EventFrame *event = new EventFrame(); + TreeEventOutputData *event = new TreeEventOutputData(); tree->Branch("Event",&event,16000,2); // preparing set of meso parameters to simulate @@ -181,7 +181,8 @@ void TestRootTree::complex_read() throw RuntimeErrorException("TestRootTree::complex_read() -> Can't get tree with name '" + tree_name + "' from root file"); } - EventFrame *event = 0; + TreeEventOutputData *event = 0; + tree->SetBranchAddress("Event", &event); // reading data from the tree diff --git a/App/src/TreeEventStructure.cpp b/App/src/TreeEventStructure.cpp new file mode 100644 index 0000000000000000000000000000000000000000..15fcebc5e99c62ca3a90f310c9d7b7bdc259afe6 --- /dev/null +++ b/App/src/TreeEventStructure.cpp @@ -0,0 +1,57 @@ +#include "TreeEventStructure.h" + + +/* ************************************************************************* */ +// Output data +/* ************************************************************************* */ +TreeEventOutputData::TreeEventOutputData() +{ + clear(); +} + + +void TreeEventOutputData::clear() +{ + nframe = 0; + alpha_i = 0; + phi_i = 0; + + nphi_f = 0; + phi_f_min = 0; + phi_f_max = 0; + + nalpha_f = 0; + alpha_f_min = 0; + alpha_f_max = 0; + + malpha = 0; + mphi = 0; + npR = 0; + + valpha_f.clear(); + vphi_f.clear(); + vi.clear(); +} + + +/* ************************************************************************* */ +// Fit data +/* ************************************************************************* */ +TreeEventFitData::TreeEventFitData() +{ + clear(); +} + +void TreeEventFitData::clear() +{ + niter = 0; + real_data.clear(); + fit_data.clear(); + diff.clear(); + axis0.clear(); + axis1.clear(); + chi2=0; + parvalues.clear(); + parnames.clear(); +} + diff --git a/App/src/main.cpp b/App/src/main.cpp index f7edbe472b3dbe49a33fde3e7a534b095744fa84..3f253b35fdf25e9a58e6318c3a52e3d86bd6cc84 100644 --- a/App/src/main.cpp +++ b/App/src/main.cpp @@ -1,16 +1,18 @@ #include "FunctionalTestFactory.h" #include "DrawHelper.h" -#include "CommandLine.h" -#include "Utils.h" +//#include "CommandLine.h" +//#include "Utils.h" #include <iostream> #include <string> #include "TROOT.h" #include "TApplication.h" #include "ProgramOptions.h" -#include "AppOptionsDescription.h" -#include "Macros.h" -#include "GISASExperiment.h" +//#include "AppOptionsDescription.h" +//#include "Macros.h" +//#include "GISASExperiment.h" + + int main(int argc, char **argv) { diff --git a/Core/Core.pro b/Core/Core.pro index aa45ce48d430bed30c2a0034341d209e866219a7..a7b49e3b8663d025c93e35da29e1cc51d033bb7b 100644 --- a/Core/Core.pro +++ b/Core/Core.pro @@ -4,9 +4,8 @@ TARGET = ScattCore TEMPLATE = lib CONFIG += plugin # to remove versions from file name -#CONFIG += debug QT -= core gui -CONFIG += BUILD_PYTHON_BOOST_MODULE # to generate python interface +#CONFIG += BUILD_PYTHON_BOOST_MODULE # to generate python interface # including common project properties include($$PWD/../shared.pri) @@ -111,7 +110,8 @@ SOURCES += \ PythonAPI/src/PythonListConverter.cpp \ PythonAPI/src/PythonModule.cpp \ PythonAPI/src/PythonPlusplusHelper.cpp \ - PythonAPI/src/PythonOutputData.cpp + PythonAPI/src/PythonOutputData.cpp \ + Tools/src/IObserver.cpp HEADERS += \ Algorithms/inc/Beam.h \ diff --git a/Core/FormFactors/src/FormFactorCrystal.cpp b/Core/FormFactors/src/FormFactorCrystal.cpp index 0d633911a5e56ec7ce524fd8ec8ac7cbe0d33178..5d0f3b2da6b699232723aa6a9dccf061cd2c4f68 100644 --- a/Core/FormFactors/src/FormFactorCrystal.cpp +++ b/Core/FormFactors/src/FormFactorCrystal.cpp @@ -46,12 +46,16 @@ complex_t FormFactorCrystal::evaluate_for_q(const cvector_t &q) const // calculate the used radius in function of the reciprocal lattice scale double radius = 1.1*m_max_rec_length; // retrieve nearest reciprocal lattice vectors - std::vector<kvector_t> rec_vectors = - m_lattice.getReciprocalLatticeVectorsWithinRadius(q_real, radius); +// std::vector<kvector_t> rec_vectors = +// m_lattice.getReciprocalLatticeVectorsWithinRadius(q_real, radius); + + m_lattice.computeReciprocalLatticeVectorsWithinRadius(q_real, radius); + const KVectorContainer &rec_vectors = m_lattice.getKVectorContainer(); // perform convolution on these lattice vectors complex_t result(0.0, 0.0); - for (std::vector<kvector_t>::const_iterator it = rec_vectors.begin(); it != rec_vectors.end(); ++it) { + //for (std::vector<kvector_t>::const_iterator it = rec_vectors.begin(); it != rec_vectors.end(); ++it) { + for (KVectorContainer::const_iterator it = rec_vectors.begin(); it != rec_vectors.end(); ++it) { cvector_t q_i((*it).x(), (*it).y(), (*it).z()); cvector_t q_min_q_i = q - q_i; complex_t basis_factor = mp_basis_form_factor->evaluate(q_i, k_zero, 0.0, 0.0); diff --git a/Core/Samples/inc/Lattice.h b/Core/Samples/inc/Lattice.h index 05f415bf284316384409e426af72dff6a51b1220..64a9bd14473a071feb45db88ac987cfe133fb02c 100644 --- a/Core/Samples/inc/Lattice.h +++ b/Core/Samples/inc/Lattice.h @@ -58,14 +58,18 @@ public: Coordinate3D<int> getNearestReciprocalLatticeVectorCoordinates(const kvector_t &vector_in) const; //! get a list of lattice vectors within a specified distance of a given vector - std::vector<kvector_t> getLatticeVectorsWithinRadius( - const kvector_t &input_vector, double radius) const; +// std::vector<kvector_t> getLatticeVectorsWithinRadius( +// const kvector_t &input_vector, double radius) const; //! get a list of reciprocal lattice vectors within a specified distance of a given vector - std::vector<kvector_t> getReciprocalLatticeVectorsWithinRadius( +// std::vector<kvector_t> getReciprocalLatticeVectorsWithinRadius( +// const kvector_t &input_vector, double radius) const; + void computeReciprocalLatticeVectorsWithinRadius( const kvector_t &input_vector, double radius) const; + //! get a list of rotation angles within a specified range that hit a maximal set of small Bragg peaks + //std::vector<double> collectBraggAngles(size_t size, double max_radius, const TRange<double> &phi_range, const TRange<double> &z_range) const; std::vector<double> collectBraggAngles(size_t size, double max_radius, const TRange<double> &phi_range, const TRange<double> &z_range) const; //! set a selection rule for the reciprocal vectors @@ -80,12 +84,20 @@ public: static Lattice createTrigonalLattice(double a, double c); + const KVectorContainer &getKVectorContainer() const { return m_kvector_container; } private: Lattice &operator=(const Lattice &lattice); - std::vector<kvector_t> getVectorsWithinRadius(const kvector_t &input_vector, +// std::vector<kvector_t> getVectorsWithinRadius(const kvector_t &input_vector, +// const Coordinate3D<int> &nearest_coords, double radius, +// const kvector_t &v1, const kvector_t &v2, const kvector_t &v3, +// const kvector_t &rec1, const kvector_t &rec2, const kvector_t &rec3) const; + + void computeVectorsWithinRadius(const kvector_t &input_vector, const Coordinate3D<int> &nearest_coords, double radius, const kvector_t &v1, const kvector_t &v2, const kvector_t &v3, const kvector_t &rec1, const kvector_t &rec2, const kvector_t &rec3) const; + + void computeReciprocalVectors() const; void computeInverseLatticeVectors() const; void computeInverseReciprocalLatticeVectors() const; @@ -95,6 +107,8 @@ private: mutable kvector_t m_b1, m_b2, m_b3; //!< Cache of basis vectors in reciprocal space mutable kvector_t m_amin1, m_amin2, m_amin3, m_bmin1, m_bmin2, m_bmin3; //!< Cache of inverse vectors for calculation of coordinates mutable bool m_cache_ok, m_is_zero; //!< Boolean indicating if the reciprocal vectors are already initialized in the cache + + mutable KVectorContainer m_kvector_container; }; diff --git a/Core/Samples/src/Lattice.cpp b/Core/Samples/src/Lattice.cpp index 1183acd68d4e5e3f664c971cf1f8d84d17c570d9..281da36ea0f00f7ea35eec0f5795b37b31bb6030 100644 --- a/Core/Samples/src/Lattice.cpp +++ b/Core/Samples/src/Lattice.cpp @@ -85,26 +85,61 @@ Coordinate3D<int> Lattice::getNearestReciprocalLatticeVectorCoordinates(const kv return Coordinate3D<int>(c1, c2, c3); } -std::vector<kvector_t> Lattice::getLatticeVectorsWithinRadius( +//std::vector<kvector_t> Lattice::getLatticeVectorsWithinRadius( +// const kvector_t& input_vector, double radius) const +//{ +// if (!m_cache_ok) { +// initialize(); +// } +// Coordinate3D<int> nearest_coords = getNearestLatticeVectorCoordinates(input_vector); +// return getVectorsWithinRadius(input_vector, nearest_coords, radius, m_a1, m_a2, m_a3, m_b1, m_b2, m_b3); +//} + +//std::vector<kvector_t> Lattice::getReciprocalLatticeVectorsWithinRadius( +// const kvector_t& input_vector, double radius) const +//{ +// if (!m_cache_ok) { +// initialize(); +// } +// Coordinate3D<int> nearest_coords = getNearestReciprocalLatticeVectorCoordinates(input_vector); +// return getVectorsWithinRadius(input_vector, nearest_coords, radius, m_b1, m_b2, m_b3, m_a1, m_a2, m_a3); + +//} + +void Lattice::computeReciprocalLatticeVectorsWithinRadius( const kvector_t& input_vector, double radius) const { if (!m_cache_ok) { initialize(); } - Coordinate3D<int> nearest_coords = getNearestLatticeVectorCoordinates(input_vector); - return getVectorsWithinRadius(input_vector, nearest_coords, radius, m_a1, m_a2, m_a3, m_b1, m_b2, m_b3); + Coordinate3D<int> nearest_coords = getNearestReciprocalLatticeVectorCoordinates(input_vector); + computeVectorsWithinRadius(input_vector, nearest_coords, radius, m_b1, m_b2, m_b3, m_a1, m_a2, m_a3); } -std::vector<kvector_t> Lattice::getReciprocalLatticeVectorsWithinRadius( - const kvector_t& input_vector, double radius) const -{ - if (!m_cache_ok) { - initialize(); - } - Coordinate3D<int> nearest_coords = getNearestReciprocalLatticeVectorCoordinates(input_vector); - return getVectorsWithinRadius(input_vector, nearest_coords, radius, m_b1, m_b2, m_b3, m_a1, m_a2, m_a3); -} +//std::vector<double> Lattice::collectBraggAngles(size_t size, double max_radius, +// const TRange<double>& phi_range, const TRange<double>& z_range) const +//{ +// std::vector<double> result; +//// int granularity = std::max(1000, (int)size); // +// double brillouin_volume = 8*M_PI*M_PI*M_PI/getVolume(); +// double max_volume = max_radius*max_radius*phi_range.getDifference()*z_range.getDifference()/2.0; +// int max_nbr_angles = (int)(max_volume/brillouin_volume); +// if (size < (size_t)max_nbr_angles) { +// max_radius *= (double)size/max_nbr_angles; +// } +// double radius = std::max(max_radius, z_range.getMax()); +// std::vector<kvector_t> rec_vectors = getReciprocalLatticeVectorsWithinRadius(kvector_t(0.0, 0.0, 0.0), radius); +// for (size_t i=0; i<rec_vectors.size(); ++i) { +// kvector_t rvec = rec_vectors[i]; +// double phi = rvec.phi(); +// if (rvec.rho()<max_radius && phi_range.inRange(phi) && z_range.inRange(rvec.z())) { +// result.push_back(phi); +// } +// } +// std::cout << "Returning " << result.size() << " angles" << std::endl; +// return result; +//} std::vector<double> Lattice::collectBraggAngles(size_t size, double max_radius, @@ -119,9 +154,11 @@ std::vector<double> Lattice::collectBraggAngles(size_t size, double max_radius, max_radius *= (double)size/max_nbr_angles; } double radius = std::max(max_radius, z_range.getMax()); - std::vector<kvector_t> rec_vectors = getReciprocalLatticeVectorsWithinRadius(kvector_t(0.0, 0.0, 0.0), radius); - for (size_t i=0; i<rec_vectors.size(); ++i) { - kvector_t rvec = rec_vectors[i]; + + computeReciprocalLatticeVectorsWithinRadius(kvector_t(0.0, 0.0, 0.0), radius); + const KVectorContainer &rec_vectors = getKVectorContainer(); + for(KVectorContainer::const_iterator it = rec_vectors.begin(); it!= rec_vectors.end(); ++it) { + const kvector_t &rvec = (*it); double phi = rvec.phi(); if (rvec.rho()<max_radius && phi_range.inRange(phi) && z_range.inRange(rvec.z())) { result.push_back(phi); @@ -131,6 +168,8 @@ std::vector<double> Lattice::collectBraggAngles(size_t size, double max_radius, return result; } + + Lattice Lattice::createFCCLattice(double a) { double b = a/2.0; @@ -158,7 +197,38 @@ void Lattice::computeReciprocalVectors() const m_b3 = 2*M_PI/DotProduct(m_a3, a12)*a12; } -std::vector<kvector_t> Lattice::getVectorsWithinRadius(const kvector_t &input_vector, +//std::vector<kvector_t> Lattice::getVectorsWithinRadius(const kvector_t &input_vector, +// const Coordinate3D<int> &nearest_coords, double radius, const kvector_t& v1, +// const kvector_t& v2, const kvector_t& v3, const kvector_t& rec1, +// const kvector_t& rec2, const kvector_t& rec3) const +//{ +// int max_X = (int)std::floor( rec1.mag()*radius/(2*M_PI) ); +// int max_Y = (int)std::floor( rec2.mag()*radius/(2*M_PI) ); +// int max_Z = (int)std::floor( rec3.mag()*radius/(2*M_PI) ); + +// std::vector<kvector_t> result; +// for (int index_X = -max_X; index_X <= max_X; ++index_X) +// { +// for (int index_Y = -max_Y; index_Y <= max_Y; ++index_Y) +// { +// for (int index_Z = -max_Z; index_Z <= max_Z; ++index_Z) +// { +// Coordinate3D<int> coords(index_X + nearest_coords[0], +// index_Y + nearest_coords[1], index_Z + nearest_coords[2]); +// if (mp_selection_rule && !mp_selection_rule->coordinateSelected(coords)) continue; +// kvector_t latticePoint = coords[0]*v1 + coords[1]*v2 + coords[2]*v3; +// if ((latticePoint - input_vector).mag() <= radius) +// { +// result.push_back(latticePoint); +// } +// } +// } +// } +// return result; +//} + + +void Lattice::computeVectorsWithinRadius(const kvector_t &input_vector, const Coordinate3D<int> &nearest_coords, double radius, const kvector_t& v1, const kvector_t& v2, const kvector_t& v3, const kvector_t& rec1, const kvector_t& rec2, const kvector_t& rec3) const @@ -167,7 +237,7 @@ std::vector<kvector_t> Lattice::getVectorsWithinRadius(const kvector_t &input_ve int max_Y = (int)std::floor( rec2.mag()*radius/(2*M_PI) ); int max_Z = (int)std::floor( rec3.mag()*radius/(2*M_PI) ); - std::vector<kvector_t> result; + m_kvector_container.clear(); for (int index_X = -max_X; index_X <= max_X; ++index_X) { for (int index_Y = -max_Y; index_Y <= max_Y; ++index_Y) @@ -180,12 +250,11 @@ std::vector<kvector_t> Lattice::getVectorsWithinRadius(const kvector_t &input_ve kvector_t latticePoint = coords[0]*v1 + coords[1]*v2 + coords[2]*v3; if ((latticePoint - input_vector).mag() <= radius) { - result.push_back(latticePoint); + m_kvector_container.push_back(latticePoint); } } } } - return result; } diff --git a/Core/Tools/inc/IObserver.h b/Core/Tools/inc/IObserver.h index 509ee6849ee97c5e03100f51e7367dc50a712895..5d4abff69a503a7015fb3bd8c58ae263bfa90ebf 100644 --- a/Core/Tools/inc/IObserver.h +++ b/Core/Tools/inc/IObserver.h @@ -25,11 +25,16 @@ class IObservable; //- ------------------------------------------------------------------- class IObserver { public: - virtual ~IObserver(){} - //! method which is used by observable object to notify change in status + //! destructor detach observer from observed subject + virtual ~IObserver(); + //! method which is used by observable subject to notify change in status virtual void update (IObservable *subject) = 0; + //! set pointer to observed subject + virtual void setObservedSubject(IObservable *subject); protected: - IObserver(){} + IObserver() : m_observed_subject(0) {} +private: + IObservable *m_observed_subject; }; @@ -44,18 +49,13 @@ public: virtual ~IObservable(){} //! attach observer to the list of observers - virtual void attachObserver(IObserver *obj) { m_observers.push_back(obj); } + virtual void attachObserver(IObserver *obj); //! detach observer from observers list - virtual void detachObserver(IObserver *obj) { m_observers.remove(obj); } + virtual void detachObserver(IObserver *obj); //! notify observers about change in status - virtual void notifyObservers() - { - for(observers_t::iterator it = m_observers.begin(); it!=m_observers.end(); ++it) { - (*it)->update(this); - } - } + virtual void notifyObservers(); protected: IObservable(){} diff --git a/Core/Tools/inc/Types.h b/Core/Tools/inc/Types.h index a1d842b54b4d2dc43d6233576feaddcf33a0014a..2a70e81250ec178a35873662638d0cee31b41dba 100644 --- a/Core/Tools/inc/Types.h +++ b/Core/Tools/inc/Types.h @@ -26,6 +26,49 @@ typedef std::vector<double > vdouble1d_t; typedef std::vector<vdouble1d_t > vdouble2d_t; +// container for holding kvectors +class KVectorContainer { +public: + typedef std::vector<kvector_t > container_t; + typedef container_t::const_iterator const_iterator; + KVectorContainer(int buff_size = 3) : m_current_position(0), m_max_buff_size(buff_size) + { + m_buffer.reserve(m_max_buff_size); + for(size_t i=0; i<m_max_buff_size; ++i) m_buffer.push_back(kvector_t(0,0,0)); + } + + inline void push_back(const kvector_t &k) { + if(m_current_position == m_max_buff_size) { + std::cout << "KVectorContainer::push_nack() -> Info. Increasing size of the buffer from " << m_max_buff_size; + m_max_buff_size *=2; + std::cout << " to " << m_max_buff_size << std::endl; + m_buffer.resize(m_max_buff_size); + } + m_buffer[m_current_position][0] = k[0]; + m_buffer[m_current_position][1] = k[1]; + m_buffer[m_current_position][2] = k[2]; + m_current_position++; + } + + inline void clear() { m_current_position = 0; } + + inline size_t size() { return m_current_position; } + + void print() { + for(size_t i=0; i<m_max_buff_size; ++i) std::cout << i << " " << m_buffer[i] << std::endl; + std::cout << "size:" << size() << std::endl; + } + + const_iterator begin() const { return m_buffer.begin(); } + const_iterator end() const { return m_buffer.begin()+m_current_position; } + + +private: + size_t m_current_position; + size_t m_max_buff_size; + container_t m_buffer; +}; + //// we need forward declaration here to be able to redefine ostream as friend //template<typename T> class KVector; //template<typename T> std::ostream& operator<< (std::ostream& o, const KVector<T>& ); diff --git a/Core/Tools/src/IObserver.cpp b/Core/Tools/src/IObserver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..42c874f96cda01f81728499c9c775d96b48cc544 --- /dev/null +++ b/Core/Tools/src/IObserver.cpp @@ -0,0 +1,43 @@ +#include "IObserver.h" + + + +/* ************************************************************************* */ +// Observer +/* ************************************************************************* */ +IObserver::~IObserver() +{ + if(m_observed_subject) m_observed_subject->detachObserver(this); +} + + +void IObserver::setObservedSubject(IObservable *subject) +{ + m_observed_subject = subject; +} + + + +/* ************************************************************************* */ +// Observable +/* ************************************************************************* */ +void IObservable::attachObserver(IObserver *obj) +{ + obj->setObservedSubject(this); + m_observers.push_back(obj); +} + + +void IObservable::detachObserver(IObserver *obj) +{ + m_observers.remove(obj); +} + + +void IObservable::notifyObservers() +{ + for(observers_t::iterator it = m_observers.begin(); it!=m_observers.end(); ++it) { + (*it)->update(this); + } +} + diff --git a/Examples/Analysis/RootMacros/README b/Examples/Analysis/RootMacros/README new file mode 100644 index 0000000000000000000000000000000000000000..d5c6dc996c4333d18c1dd23b8eac05327f8869ed --- /dev/null +++ b/Examples/Analysis/RootMacros/README @@ -0,0 +1,9 @@ +Directory contains exmpales with ROOT analysis macros + +ex01_MesoGif - analysis of the data produced by TestRootTree +ex02_ReadFitData - analysis of the data produced by TestFittingModule + + + +To clean remainings from ROOT session +rm -f */*.so */*.d */*~ diff --git a/Examples/Analysis/MesoData/MesoData.C b/Examples/Analysis/RootMacros/ex01_MesoGif/MesoData.C similarity index 100% rename from Examples/Analysis/MesoData/MesoData.C rename to Examples/Analysis/RootMacros/ex01_MesoGif/MesoData.C diff --git a/Examples/Analysis/MesoData/MesoData.h b/Examples/Analysis/RootMacros/ex01_MesoGif/MesoData.h similarity index 100% rename from Examples/Analysis/MesoData/MesoData.h rename to Examples/Analysis/RootMacros/ex01_MesoGif/MesoData.h diff --git a/Examples/Analysis/MesoData/analysis1.C b/Examples/Analysis/RootMacros/ex01_MesoGif/analysis1.C similarity index 100% rename from Examples/Analysis/MesoData/analysis1.C rename to Examples/Analysis/RootMacros/ex01_MesoGif/analysis1.C diff --git a/Examples/Analysis/MesoData/analysis2.C b/Examples/Analysis/RootMacros/ex01_MesoGif/analysis2.C similarity index 100% rename from Examples/Analysis/MesoData/analysis2.C rename to Examples/Analysis/RootMacros/ex01_MesoGif/analysis2.C diff --git a/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.C b/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.C new file mode 100644 index 0000000000000000000000000000000000000000..84e6e5e8fbb413f2cdac35bc12bb151affcd0e0f --- /dev/null +++ b/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.C @@ -0,0 +1,43 @@ +#define FitData_cxx +#include "FitData.h" +#include <TH2.h> +#include <TStyle.h> +#include <TCanvas.h> + +void FitData::Loop() +{ +// In a ROOT session, you can do: +// Root > .L FitData.C +// Root > FitData t +// Root > t.GetEntry(12); // Fill t data members with entry number 12 +// Root > t.Show(); // Show values of entry 12 +// Root > t.Show(16); // Read and show values of entry 16 +// Root > t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentry<nentries;jentry++) { + Long64_t ientry = LoadTree(jentry); + if (ientry < 0) break; + nb = fChain->GetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + } +} diff --git a/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.h b/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.h new file mode 100644 index 0000000000000000000000000000000000000000..f8875b1b44c67934d1a65d7e00d4f9a0f221b8b8 --- /dev/null +++ b/Examples/Analysis/RootMacros/ex02_ReadFitData/FitData.h @@ -0,0 +1,154 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Wed Oct 10 15:14:35 2012 by ROOT version 5.34/01 +// from TTree FitSuiteTree/Oh, my data +// found on file: ../../../../fitsuite.root +////////////////////////////////////////////////////////// + +#ifndef FitData_h +#define FitData_h + +#include <TROOT.h> +#include <TChain.h> +#include <TFile.h> + +// Header file for the classes stored in the TTree if any. + +// Fixed size dimensions of array or collections stored in the TTree if any. + +class FitData { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + + // Declaration of leaf types + //TreeEventFitData *Event; + Int_t niter; + vector<vector<double> > real_data; + vector<vector<double> > fit_data; + vector<vector<double> > diff; + vector<vector<double> > axis0; + vector<vector<double> > axis1; + Double_t chi2; + vector<double> params; + vector<string> parnames; + + // List of branches + TBranch *b_Event_niter; //! + TBranch *b_Event_real_data; //! + TBranch *b_Event_fit_data; //! + TBranch *b_Event_diff; //! + TBranch *b_Event_axis0; //! + TBranch *b_Event_axis1; //! + TBranch *b_Event_chi2; //! + TBranch *b_Event_params; //! + TBranch *b_Event_parnames; //! + + FitData(TTree *tree=0); + virtual ~FitData(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); +}; + +#endif + +#ifdef FitData_cxx +FitData::FitData(TTree *tree) : fChain(0) +{ +// if parameter tree is not specified (or zero), connect the file +// used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("../../../../fitsuite.root"); + if (!f || !f->IsOpen()) { + f = new TFile("../../../../fitsuite.root"); + } + f->GetObject("FitSuiteTree",tree); + + } + Init(tree); +} + +FitData::~FitData() +{ + if (!fChain) return; + delete fChain->GetCurrentFile(); +} + +Int_t FitData::GetEntry(Long64_t entry) +{ +// Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} +Long64_t FitData::LoadTree(Long64_t entry) +{ +// Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (fChain->GetTreeNumber() != fCurrent) { + fCurrent = fChain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void FitData::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("niter", &niter, &b_Event_niter); + fChain->SetBranchAddress("real_data", &real_data, &b_Event_real_data); + fChain->SetBranchAddress("fit_data", &fit_data, &b_Event_fit_data); + fChain->SetBranchAddress("diff", &diff, &b_Event_diff); + fChain->SetBranchAddress("axis0", &axis0, &b_Event_axis0); + fChain->SetBranchAddress("axis1", &axis1, &b_Event_axis1); + fChain->SetBranchAddress("chi2", &chi2, &b_Event_chi2); + fChain->SetBranchAddress("params", ¶ms, &b_Event_params); + fChain->SetBranchAddress("parnames", &parnames, &b_Event_parnames); + Notify(); +} + +Bool_t FitData::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void FitData::Show(Long64_t entry) +{ +// Print contents of entry. +// If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} +Int_t FitData::Cut(Long64_t entry) +{ +// This function may be called from Loop. +// returns 1 if entry is accepted. +// returns -1 otherwise. + return 1; +} +#endif // #ifdef FitData_cxx diff --git a/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata1.C b/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata1.C new file mode 100644 index 0000000000000000000000000000000000000000..ac1c306b6dcb08e578738149e7657308be198b94 --- /dev/null +++ b/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata1.C @@ -0,0 +1,61 @@ +// Simple drawing of fit data stored in the FitSuiteTree using TTree::Draw() method +// (data produced by FitSuiteObserverWriteTree un TestFittingModule) +// Usage: +// +// Way #1 (script mode) +// > root -l analyse_fitdata1.C +// +// Way #2 (precompiled mode) +// > root -l +// > .L analyse_fitdata1.C++ +// analyse_fitdata1() + +#include <iostream> +#include <string> +#include <vector> +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" + +#include "FitData.h" + + +void analyse_fitdata1() +{ + const char *file_name ="../../../../fitsuite.root"; + const char *tree_name = "FitSuiteTree"; + + // opening root file + TFile *top = new TFile(file_name,"READ"); + if( !top->IsOpen() ) { + std::cout << "Can't open file '" << file_name << "'" << std::endl; + return; + } + + // access tree stored in the file + TTree *tree = (TTree *) top->Get(tree_name); + if( !tree ) { + std::cout << "Can't get tree '" << tree_name << "'" << std::endl; + return; + } + + TCanvas *c1 = new TCanvas("c1","c1",1024, 768); + c1->Divide(2,2); + + // loop over tree entries and drawing + for(int i=0; i<tree->GetEntries(); i++){ + // defining criteria for Draw method -> it will draw only data corresponding to the event with given number + char criteria[128]; + sprintf(criteria,"Entry$==%d",i); + c1->cd(1); gPad->SetLogz(); + tree->Draw("real_data:axis1:axis0",criteria,"prof COLZ"); + c1->cd(2); gPad->SetLogz(); + tree->Draw("fit_data:axis1:axis0",criteria,"prof COLZ"); + c1->Update(); + } + + + +} + + diff --git a/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata2.C b/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata2.C new file mode 100644 index 0000000000000000000000000000000000000000..1456163035643bb9af41f3614545bb8270cabc14 --- /dev/null +++ b/Examples/Analysis/RootMacros/ex02_ReadFitData/analyse_fitdata2.C @@ -0,0 +1,103 @@ +// Simple drawing of fit data stored in the FitSuiteTree (using FitData class) +// (data produced by FitSuiteObserverWriteTree un TestFittingModule) +// Usage: +// +// Way #1 (script mode) +// > root -l analyse_fitdata1.C +// +// Way #2 (precompiled mode) +// > root -l +// > .L analyse_fitdata1.C++ +// analyse_fitdata1() + +#include <iostream> +#include <string> +#include <vector> +#include "TFile.h" +#include "TChain.h" +#include "TCanvas.h" +#include "TH2D.h" +#include "TSystem.h" +#include "TPaveText.h" + +using namespace std; +#include "FitData.C" + + +void analyse_fitdata2() +{ + const char *file_name ="../../../../fitsuite.root"; + const char *tree_name = "FitSuiteTree"; + + TChain *chain = new TChain(tree_name); + chain->Add(file_name); + + FitData *event = new FitData(chain); + + // reading first event + chain->GetEntry(0); + + // finding histogram limits + std::vector<double> axis0; + for(size_t i=0; i<event->axis0.size(); ++i) axis0.push_back(event->axis0[i][0]); + std::vector<double> axis1 = event->axis1[0]; + double dphi = (axis0.back()-axis0.front())/double(axis0.size()-1); + double phi1 = axis0.front() - dphi/2.; + double phi2 = axis0.back() + dphi/2.; + double dalpha = (axis1.back()-axis1.front())/double(axis1.size()-1); + double alpha1 = axis1.front() - dalpha/2.; + double alpha2 = axis1.back() + dalpha/2.; + std::cout << "dphi: " << dphi << " phi1:" << phi1 << " phi2:" << phi2 << std::endl; + std::cout << "dalpha: " << dalpha << " alpha1:" << alpha1 << " alpha2:" << alpha2 << std::endl; + + TH2D *myHisto = new TH2D("myHisto","myHisto", axis0.size(), phi1, phi2, axis0.size(), alpha1, alpha2); + myHisto->SetContour(50); + myHisto->SetTitle(""); + myHisto->SetMinimum(1); + myHisto->SetStats(0); + + TCanvas *c1 = new TCanvas("c1","c1",1024, 768); + c1->Divide(2,2); + // loop over entries + int i_entry(0); + while(chain->GetEntry(i_entry++)){ + for(size_t i=0; i<axis0.size(); i++){ + for(size_t j=0; j<axis1.size(); j++) { + myHisto->Fill(event->axis0[i][j], event->axis1[i][j], event->real_data[i][j]); + } + } + c1->cd(1); gPad->SetLogz(); + myHisto->DrawCopy("colz"); + myHisto->Reset(); + for(size_t i=0; i<axis0.size(); i++){ + for(size_t j=0; j<axis1.size(); j++) { + myHisto->Fill(event->axis0[i][j], event->axis1[i][j], event->fit_data[i][j]); + } + } + c1->cd(2); gPad->SetLogz(); + myHisto->DrawCopy("colz"); + myHisto->Reset(); + + // printing parameters + c1->cd(3); + TPaveText *pt = new TPaveText(.05,.1,.95,.8); + char str[256]; + sprintf(str,"Iteration %d",event->niter); + pt->AddText(str); + sprintf(str,"chi2 %e",event->chi2); + pt->AddText(str); + for(size_t i =0; i<event->params.size(); ++i) { + sprintf(str,"%s %f", event->parnames[i].c_str(), event->params[i]); + pt->AddText(str); + } + pt->Draw(); + + c1->Update(); + gSystem->Sleep(5); // wait a bit for smooth animation + } // loop over recorder events + + + +} + + diff --git a/Examples/Analysis/RootMacros/ex02_ReadFitData/make_class.C b/Examples/Analysis/RootMacros/ex02_ReadFitData/make_class.C new file mode 100644 index 0000000000000000000000000000000000000000..1334911927542d664a73fa2456eb255e33eb3ad8 --- /dev/null +++ b/Examples/Analysis/RootMacros/ex02_ReadFitData/make_class.C @@ -0,0 +1,37 @@ +// Generates classes to read ROOT tree from existing ROOT file +// +// Usage: +// Way 1 (script mode): +// > root -l -q .x make_class.C +// Way 2 (precompiled mode): +// > root -l +// > .L make_class.C++ +// > make_class() + + +#include "TFile.h" +#include "TTree.h" +#include <iostream> + +void make_class(const char *file_name ="../../../../fitsuite.root" + ,const char *tree_name = "FitSuiteTree" + ,const char *class_name = "FitData") +{ + + // open root file + TFile *top = new TFile(file_name,"READ"); + + // access to the tree stored in the file + if( !top->IsOpen() ) { + std::cout << "File is absent " << file_name << std::endl; + return; + } + TTree *tree = (TTree *)top->Get(tree_name); + if( !tree ) { + std::cout << "Tree is absent " << tree_name << std::endl; + return; + } + + // generating class files (*.h, *.C) + tree->MakeClass(class_name); +} diff --git a/Examples/Performance/perf_history.txt b/Examples/Performance/perf_history.txt index 78c703f8950592cf82990598551194ecfdb7f175..e4dd0dab31e2c010aaf2a0ee7c8078b83a4840e8 100644 --- a/Examples/Performance/perf_history.txt +++ b/Examples/Performance/perf_history.txt @@ -120,3 +120,8 @@ # current status of the art (not possible to trace what was changed) 2012-10-08 17:36:40 | jcnsopc73 | macosx64, 2800 MHz | 270270 | 13.986 | 13.8889 | 4.7619 | 2012-10-08 17:37:28 | jcnsopc73 | macosx64, 2800 MHz | 266667 | 13.986 | 13.986 | 4.7619 | +2012-10-09 14:00:42 | jcnsopc73 | macosx64, 2800 MHz | 266667 | 14.1844 | 13.986 | 4.7619 | + +# start rearanging ISample (Iparametrized) +2012-10-10 15:38:43 | jcnsopc73 | macosx64, 2800 MHz | 266667 | 14.1844 | 13.986 | 4.65116 | +2012-10-10 15:39:04 | jcnsopc73 | macosx64, 2800 MHz | 259740 | 14.0845 | 13.6986 | 4.65116 | diff --git a/Macros/GitUtils/gisasfw_loc.png b/Macros/GitUtils/gisasfw_loc.png index 5cccb7fa5df8baea7f69757132fd19cfd486030a..4965f53a19e2660f1b665fb0010496306c7d3056 100644 Binary files a/Macros/GitUtils/gisasfw_loc.png and b/Macros/GitUtils/gisasfw_loc.png differ