TestRootTree example

App/App.pro

@@ -31,7 +31,9 @@ SOURCES += \
     src/TestIsGISAXS9.cpp \
     src/TestIsGISAXS10.cpp \
     src/TestMesoCrystal.cpp \
+    src/TestRootTree.cpp \
     src/TestRoughness.cpp \
+    src/TestEventStructure.cpp
 
 HEADERS += \
     inc/App.h \
@@ -53,7 +55,9 @@ HEADERS += \
     inc/TestIsGISAXS9.h \
     inc/TestIsGISAXS10.h \
     inc/TestMesoCrystal.h \
+    inc/TestRootTree.h \
     inc/TestRoughness.h \
+    inc/TestEventStructure.h
 
 INCLUDEPATH += ./inc ../Core/Algorithms/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc
 DEPENDPATH  += ./inc ../Core/Algorithms/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc
@@ -103,7 +107,7 @@ CONFIG(JCNS) {
 exists($$(ROOTSYS)/bin/root-config){
   INCLUDEPATH += $$system($ROOTSYS/bin/root-config --incdir)
   #LIBS += $$system($ROOTSYS/bin/root-config --glibs)
-  LIBS += -L$$system($ROOTSYS/bin/root-config --libdir ) -lGui -lCore -lCint -lRIO -lHist -lGraf -lGraf3d -lGpad  -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -lm -ldl
+  LIBS += -L$$system($ROOTSYS/bin/root-config --libdir ) -lGui -lCore -lCint -lRIO -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -lm -ldl
 
   MYROOTCINT = ${ROOTSYS}/bin/rootcint
 }

App/inc/App.h

@@ -16,5 +16,7 @@
 
 #include "ISingleton.h"
 #include "DrawHelper.h"
+//#include "OutputData.h"
+#include "TestEventStructure.h"
 
 #endif // APP_H

App/inc/AppLinkDef.h

@@ -21,4 +21,8 @@
 #pragma link C++ class ISingleton<DrawHelper>+;
 #pragma link C++ class DrawHelper+;
 
+//#pragma link C++ class MultiIndex+;
+//#pragma link C++ class OutputData<double>+;
+#pragma link C++ class TestEventStructure+;
+
 #endif

App/inc/TestEventStructure.h

+#ifndef TESTEVENTSTRUCTURE_H
+#define TESTEVENTSTRUCTURE_H
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   TestEventStructure.h
+//! @brief  Definition of TestEventStructure class for writing in root files
+//! @author Scientific Computing Group at FRM II
+//! @date   19.07.2012
+
+
+#include "TObject.h"
+#include <vector>
+
+//- -------------------------------------------------------------------
+//! @class TestEventStructure
+//! @brief test structure to save in the root file
+//- -------------------------------------------------------------------
+class TestEventStructure
+{
+public:
+    TestEventStructure();
+    ~TestEventStructure(){}
+
+    double alpha_i;
+    double alpha_f;
+    std::vector<double > vec;
+    //OutputData<double> m_data;
+
+ClassDef(TestEventStructure,1)
+};
+
+
+
+#endif // TESTEVENTSTRUCTURE_H

App/inc/TestRootTree.h

+#ifndef TESTROOTTREE_H
+#define TESTROOTTREE_H
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   TestRootTree.h
+//! @brief  Definition of TestRootTree class for testing ROOT trees
+//! @author Scientific Computing Group at FRM II
+//! @date   18.07.2012
+
+#include "IFunctionalTest.h"
+#include "ISample.h"
+#include "OutputData.h"
+#include "GISASExperiment.h"
+
+
+//- -------------------------------------------------------------------
+//! @class TestRootTree.h
+//! @brief using ROOT trees to read/write data from/to disk
+//- -------------------------------------------------------------------
+class TestRootTree : public IFunctionalTest
+{
+public:
+    TestRootTree();
+    virtual ~TestRootTree();
+
+    virtual void execute();
+
+private:
+    //! example showing writing in the tree simple data structures
+    void simple_write();
+
+    //! example showing reading from the tree simple data structures
+    void simple_read();
+
+    //! example showing writing in the tree complex data structures
+    void complex_write();
+
+    //! example showing reading from the tree complex data structures
+    void complex_read();
+
+    //! prepare for calculations
+    void prepare_experiment();
+
+    ISample *m_sample;
+    GISASExperiment *m_experiment;
+    OutputData<double> *m_data;
+};
+
+
+
+#endif // TESTROOTTREE_H

App/src/CommandLine.cpp

@@ -48,6 +48,7 @@ CommandLine::CommandLine(int argc, char **argv) : m_argc(argc), m_argv(argv)
     m_description["diffuse"]      = "functional test: diffuse scattering from multi layer with roughness";
     m_description["formfactor"]   = "functional test: some formfactor";
     m_description["roughness"]    = "functional test: roughness parameters";
+    m_description["roottree"]     = "functional test: using root trees to read/write data from/to disk";
 
 }
 

App/src/IsGISAXSTools.cpp

@@ -140,11 +140,11 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
         double x = left_clone->next();
         if(x!=0) {
             x = log10(fabs(x));
-            h1_spect.Fill( x );
         } else {
             // lets put the cases then the difference is exactly 0 to underflow bin
             x = -21.;
         }
+        h1_spect.Fill( x );
     }
 
     gPad->SetLogy();

App/src/TestEventStructure.cpp

+#include "TestEventStructure.h"
+
+TestEventStructure::TestEventStructure()
+{
+}

App/src/TestFactory.cpp

@@ -10,6 +10,7 @@
 #include "TestConvolution.h"
 #include "TestDetectorResolution.h"
 #include "TestMesoCrystal.h"
+#include "TestRootTree.h"
 
 #include "TBenchmark.h"
 
@@ -30,6 +31,7 @@ TestFactory::TestFactory() : m_benchmark(0)
     registerItem("convolution", IFactoryCreateFunction<TestConvolution, IFunctionalTest> );
     registerItem("detectorresolution", IFactoryCreateFunction<TestDetectorResolution, IFunctionalTest> );
     registerItem("mesocrystal", IFactoryCreateFunction<TestMesoCrystal, IFunctionalTest> );
+    registerItem("roottree",    IFactoryCreateFunction<TestRootTree, IFunctionalTest> );
 
     m_benchmark = new TBenchmark();
 }

App/src/TestIsGISAXS9.cpp

@@ -18,12 +18,6 @@
 #include "TVector3.h"
 #include "TRandom.h"
 
-// сравнение 1d и 2d of two output data classes
-// случай с qx,qy=0
-// 45 degree
-// output data: const index
-// output data: export to *double, export to numpy array
-// trree example simple, complex
 
 
 TestIsGISAXS9::TestIsGISAXS9()
@@ -93,7 +87,7 @@ void TestIsGISAXS9::finalise()
         c1->cd(2); gPad->SetLogz();
         IsGISAXSTools::drawOutputDataInPad(*isgi_data, "CONT4 Z", "IsGisaxs pyramid FF");
 
-        IsGISAXSTools::resetMinimum(); IsGISAXSTools::resetMaximum();
+        //IsGISAXSTools::resetMinimum(); IsGISAXSTools::resetMaximum();
         IsGISAXSTools::setMinimum(-0.0001);
         IsGISAXSTools::setMaximum(0.0001);
 

App/src/TestRootTree.cpp

+#include "TestRootTree.h"
+#include "MultiLayer.h"
+#include "MaterialManager.h"
+#include "NanoParticleDecoration.h"
+#include "NanoParticle.h"
+#include "LayerDecorator.h"
+#include "GISASExperiment.h"
+#include "FormFactors.h"
+#include "Units.h"
+#include "InterferenceFunctionNone.h"
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "IsGISAXSTools.h"
+#include "TestEventStructure.h"
+
+
+TestRootTree::TestRootTree() : m_sample(0), m_experiment(0), m_data(0)
+{
+
+}
+
+
+TestRootTree::~TestRootTree()
+{
+    delete m_sample;
+    delete m_experiment;
+    delete m_data;
+}
+
+
+void TestRootTree::execute()
+{
+
+    // preparing some real data for tree playing
+    prepare_experiment();
+
+    // example showing handling of tree with simple data structures
+    simple_write();
+    simple_read();
+
+    // example showing handling of tree with complex data structures
+//    complex_write();
+//    complex_read();
+}
+
+
+
+/* ************************************************************************* */
+// example showing writing in the tree simple data structures
+/* ************************************************************************* */
+void TestRootTree::simple_write()
+{
+    // variables below will be written in the tree
+    double intens1(0), intens2(0), alpha_i(0), phi_i(0), alpha_f(0), phi_f(0);
+    int nev;
+    //std::vector<double > *somevector;
+
+    std::string root_file_name = "mydata.root";
+    std::string tree_name = "mytree";
+
+    // preparing root file for writing
+    TFile *top = new TFile(root_file_name.c_str(),"RECREATE");
+    if( !top->IsOpen() ) {
+        throw RuntimeErrorException("TestRootTree::simple_write() -> Can't open file "+root_file_name+" for writing");
+    }
+
+    // creating new tree
+    TTree *tree = new TTree(tree_name.c_str(),"Oh, my data");
+
+    // pointing tree to local variables
+    tree->Branch("intens1",&intens1, "intens1/D");
+    tree->Branch("intens2",&intens2, "intens2/D");
+    tree->Branch("alpha_i",&alpha_i, "alpha_i/D");
+    tree->Branch("alpha_f",&alpha_f, "alpha_f/D");
+    tree->Branch("phi_f",&phi_f, "phi_f/D");
+    tree->Branch("nev",&nev, "nev/I");
+    //somevector = new std::vector<double>;
+    //tree->Branch("somevector","vector<double >", &somevector);
+
+    TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
+
+    // our experiment begins here
+    const int nTotalEvents = 10;
+    TRandom mr(0);
+    for(int i_ev=0; i_ev<nTotalEvents; i_ev++) {
+        if(i_ev%10 ==0 ) std::cout << "nevent:" << i_ev << std::endl;
+
+        alpha_i = -0.3 + 0.1*mr.Rndm(); // generating random alpha_i in the interval
+        //phi_i = M_PI*2.*mr.Rndm(); // generating random phi_i in the interval
+        phi_i = 0;
+        nev = i_ev;
+
+        m_experiment->setBeamParameters(1.0*Units::angstrom, alpha_i*Units::degree, phi_i);
+        m_experiment->runSimulation();
+
+        m_data = m_experiment->getOutputData();
+        // accessing to scattering data
+        NamedVector<double> *axis0 = dynamic_cast<NamedVector<double>*>(m_data->getAxes()[0]);
+        NamedVector<double> *axis1 = dynamic_cast<NamedVector<double>*>(m_data->getAxes()[1]);
+        std::string axis0_name = axis0->getName();
+        std::string axis1_name = axis1->getName();
+
+        c1->cd(); gPad->SetLogz();
+        c1->Clear();
+        IsGISAXSTools::setMinimum(1.);
+        IsGISAXSTools::drawOutputDataInPad(*m_data, "CONT4 Z", "IsGisaxs pyramid FF");
+        c1->Modified();
+        c1->Update();
+        c1->Write();
+
+        m_data->resetIndex();
+        while (m_data->hasNext())
+        {
+            size_t index_phi_f =  m_data->getCurrentIndexOfAxis(axis0_name.c_str());
+            size_t index_alpha_f = m_data->getCurrentIndexOfAxis(axis1_name.c_str());
+            phi_f = Units::rad2deg( (*axis0)[index_phi_f]);
+            alpha_f = Units::rad2deg( (*axis1)[index_alpha_f] );
+            //std::cout << phi_f << " " << alpha_f << std::endl;
+            intens1 =m_data->next();
+            tree->Fill();
+
+        }
+
+        delete m_data;
+    }
+
+    tree->Write();
+    top->Close();
+    delete top;
+
+    //delete somevector;
+
+    std::string info;
+    info += "Root file '" + root_file_name + "' has been successfully created. \n";
+    info += "Hits for root session: \n";
+    info += "root -l " + root_file_name + "\n";
+    info += "TBrowser b; \n";
+    info += "TTree *tree = (TTree *)_file0->Get(\"" + tree_name + "\"); \n";
+    info += "gPad->SetLogz(); \n";
+    info += "tree->Draw(\"alpha_i\"); \n";
+    info += "tree->Draw(\"alpha_f:phi_f\"); \n";
+    info += "tree->Draw(\"intens1:alpha_f:phi_f\",\"alpha_i>-0.2&&alpha_i<-0.1 \",\"prof CONT4 Z\"); \n";
+
+    std::cout << info << std::endl;
+
+}
+
+
+/* ************************************************************************* */
+// example showing reading from the tree simple data structures
+/* ************************************************************************* */
+void TestRootTree::simple_read()
+{
+    std::cout << "TestRootTree::simple_read() -> going to red tree back from file" << std::endl;
+
+    std::string root_file_name = "mydata.root";
+    std::string tree_name = "mytree";
+
+    // preparing root file for reading
+    TFile *top = new TFile(root_file_name.c_str(),"READ");
+    if( !top->IsOpen() ) {
+        throw RuntimeErrorException("TestRootTree::simple_read() -> Can't open file "+root_file_name+" for reading");
+    }
+
+    // getting existing tree
+    TTree *tree = (TTree *)top->Get(tree_name.c_str());
+    if( !tree ) {
+        throw RuntimeErrorException("TestRootTree::simple_read() -> Can't get tree with name '" + tree_name + "' from root file");
+    }
+
+    TCanvas *c2 = new TCanvas("c2","reading tree back",1024, 768);
+    c2->cd(); gPad->SetLogz();
+    tree->Print();
+    tree->Draw("intens1:alpha_f:phi_f","intens1>0","prof CONT4 Z");
+
+    // not reading the tree in another way
+    // variables below will be red from the tree
+    double intens1(0), intens2(0), alpha_i(0), phi_i(0), alpha_f(0), phi_f(0);
+    int nev;
+
+    // pointing tree to local variables
+    tree->SetBranchAddress("intens1",&intens1);
+    tree->SetBranchAddress("intens2",&intens2);
+    tree->SetBranchAddress("alpha_i",&alpha_i);
+    tree->SetBranchAddress("alpha_f",&alpha_f);
+    tree->SetBranchAddress("phi_f",&phi_f);
+    tree->SetBranchAddress("nev",&nev);
+
+    // loop over all records stored in the tree
+    for(int i=0; i<tree->GetEntries(); i++) {
+        tree->GetEntry(i);
+        // at this point all local variables are filled with values from the tree
+        std::cout << "alpha_i:" << alpha_i << " phi_i:" << phi_i  << " alpha_f:" << alpha_f << " phi_f:" << phi_f << " intens:" << intens1 << std::endl;
+        if(i>10) break;
+    }
+
+    top->Close();
+
+}
+
+
+
+/* ************************************************************************* */
+// example showing writing in the tree simple data structures
+/* ************************************************************************* */
+void TestRootTree::complex_write()
+{
+
+    std::string root_file_name = "mydata2.root";
+    std::string tree_name = "mytree2";
+
+    // preparing root file for writing
+    TFile *top = new TFile(root_file_name.c_str(),"RECREATE");
+    if( !top->IsOpen() ) {
+        throw RuntimeErrorException("TestRootTree::complex_write() -> Can't open file "+root_file_name+" for writing");
+    }
+
+    // creating new tree
+    TTree *tree = new TTree(tree_name.c_str(),"Oh, my data");
+
+    TestEventStructure *event = new TestEventStructure();
+    tree->Branch("Event",&event,16000,2);
+
+    // writing 10 events
+    for(int i=0; i<10; i++) {
+        event->alpha_i = i+1;
+        event->alpha_f = i+10;
+        event->vec.resize(10, 99.9);
+        tree->Fill();
+    }
+
+    tree->Write();
+    top->Close();
+}
+
+
+/* ************************************************************************* */
+// example showing reading from the tree simple data structures
+/* ************************************************************************* */
+void TestRootTree::complex_read()
+{
+    std::cout << "TestRootTree::complex_read() -> going to red tree back from file" << std::endl;
+
+    std::string root_file_name = "mydata2.root";
+    std::string tree_name = "mytree2";
+
+    // preparing root file for reading
+    TFile *top = new TFile(root_file_name.c_str(),"READ");
+    if( !top->IsOpen() ) {
+        throw RuntimeErrorException("TestRootTree::complex_read() -> Can't open file "+root_file_name+" for reading");
+    }
+
+    // getting existing tree
+    TTree *tree = (TTree *)top->Get(tree_name.c_str());
+    if( !tree ) {
+        throw RuntimeErrorException("TestRootTree::complex_read() -> Can't get tree with name '" + tree_name + "' from root file");
+    }
+
+    TestEventStructure *event = 0;
+    tree->SetBranchAddress("Event", &event);
+
+    // reading data from the tree
+    for(int i=0; i<tree->GetEntries(); i++) {
+        tree->GetEntry(i);
+        std::cout << event->alpha_i << " " << event->alpha_f << " " << std::endl;
+    }
+
+    top->Close();
+
+}
+
+
+
+/* ************************************************************************* */
+// prepare for calculations
+/* ************************************************************************* */
+void TestRootTree::prepare_experiment()
+{
+    if(m_sample) delete m_sample;
+    m_sample=0;
+
+    if(m_experiment) delete m_experiment;
+    m_experiment = 0;
+
+    if(m_data) delete m_data;
+    m_data = 0;
+
+    // making sample
+    MultiLayer *multi_layer = new MultiLayer();
+    complex_t n_air(1.0, 0.0);
+    complex_t n_substrate(1.0-6e-6, 2e-8);
+    complex_t n_particle(1.0-6e-4, 2e-8);
+    const IMaterial *p_air_material = MaterialManager::instance().addHomogeneousMaterial("Air", n_air);
+    const IMaterial *p_substrate_material = MaterialManager::instance().addHomogeneousMaterial("Substrate", n_substrate);
+    Layer air_layer;
+    air_layer.setMaterial(p_air_material);
+    Layer substrate_layer;
+    substrate_layer.setMaterial(p_substrate_material);
+    NanoParticleDecoration particle_decoration(
+                new NanoParticle(n_particle, new FormFactorPyramid(5*Units::nanometer, 5*Units::nanometer, Units::deg2rad(54.73 ) ) ) );
+    particle_decoration.addInterferenceFunction(new InterferenceFunctionNone());
+    LayerDecorator air_layer_decorator(air_layer, particle_decoration);
+
+    multi_layer->addLayer(air_layer_decorator);
+    multi_layer->addLayer(substrate_layer);
+
+    m_sample = multi_layer;
+
+    // setting experiment
+    m_experiment = new GISASExperiment();
+    m_experiment->setDetectorParameters(0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, 100, true);
+    m_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+
+    m_experiment->setSample(m_sample);
+
+}
+
+