Интеграция гистограмм CERN ROOT в Qt GUI Geant4

Собственно, возможность открывать гистограммы из GUI Geant4 я прикрутил уже давно. Реализацию кооперации ROOT и Geant4 можно найти в исходном коде ChargeExchangeMC, который поставляется вместе с Geant4 как один из advanced examples. Эта реализация основана на использовании бэкенда Qt CERN ROOT, поэтому она доступна только при использовании плагина визуализации Qt Geant4.

Базовая реализация интеграции гистограмм в ChargeExchangeMC (входящая в состав последней версии Geant4) была основана исключительно на вводе пользователем команд в командной строке Geant4. Но недавно я подумал, а почему бы не создать специальное меню, в котором будут перечислены все имеющиеся гистограммы в естественном древовидном порядке. При клике на имени гистограммы в меню она будет отображаться в канвасе: это намного проще и естественнее, чем набирать вручную команду листинга гистограмм, затем выбирать из полученного списка нужную гистограмму и вводить новую команду для ее отображения.

На следующей картинке изображена сессия Qt Geant4 c одной из гистограмм, изображающей пучок, прошедший сквозь входной мониторный счетчик. Гистограмма была открыта из меню Histo, находящемся в верхней левой части окна справа от меню Viewer. Все картинки ниже кликабельны.

А это меню Histo в развернутом виде:

Одна из поддиректорий развернута, внутри нее находятся десять гистограмм. Для удобства восприятия перед именами двух- и трехмерных гистограмм стоят соответствующие значки 2: и 3:.

Самой главной технической сложностью здесь было построение подменю. Интерфейсы GUI Geant4 их просто не поддерживают. Один из вариантов реализации поддержки подменю - субклассинг G4UIQt - мне не понравился, так как из пушки по воробьям. В итоге я решил добавить новую функцию AddSubmenu() прямо в класс CexmcHistoManager. Вот реализация функции BuildMenuTree(), которая рекурсивно строит меню гистограмм, и функции AddSubmenu():

void  CexmcHistoManager::BuildMenuTree( G4UIQt *  session,
                                        const G4String &  menu, TList *  ls )
{
    TIter      objs( ls );
    TObject *  obj( NULL );

    while ( ( obj = ( TObject * )objs() ) )
    {
        G4String  name( obj->GetName() );
        G4String  title( obj->GetTitle() );

        if ( obj->IsFolder() )
        {
            AddSubmenu( session, menu, name, title );
            BuildMenuTree( session, name, ( ( TDirectory * )obj )->GetList() );
        }
        else
        {
            G4String  options( name );

            do
            {
                if ( obj->InheritsFrom( TH3::Class() ) &&
                     ! drawOptions3D.empty() )
                {
                    title = G4String( "3: " ) + title;
                    options += G4String( " " ) + drawOptions3D;
                    break;
                }
                if ( obj->InheritsFrom( TH2::Class() ) &&
                     ! drawOptions2D.empty() )
                {
                    title = G4String( "2: " ) + title;
                    options += G4String( " " ) + drawOptions2D;
                    break;
                }
                if ( obj->InheritsFrom( TH1::Class() ) &&
                     ! drawOptions1D.empty() )
                {
                    options += G4String( " " ) + drawOptions1D;
                    break;
                }
            } while ( false );

            G4String  cmd( CexmcMessenger::histoDirName + "draw " + options );
            session->AddButton( menu, title.c_str(), cmd );
        }
    }
}

void  CexmcHistoManager::AddSubmenu( G4UIQt *  session,
                                     const G4String &  parent,
                                     const G4String &  name,
                                     const G4String &  label )
{
  QMenu *  menu( new QMenu( label.c_str() ) );
  QMenu *  parentMenu( ( QMenu * )session->GetInteractor( parent ) );

  parentMenu->addMenu( menu );
  session->AddInteractor( name, ( G4Interactor )menu );
}

Предполагается, что BuildMenuTree() будет вызвана в момент инициализации менеджера гистограмм примерно таким образом:

    G4UIQt *  qtSession( dynamic_cast< G4UIQt * >( session ) );

    if ( qtSession )
    {
        qtSession->AddMenu( histoMenuHandle, histoMenuLabel );
        BuildMenuTree( qtSession, histoMenuHandle, gDirectory->GetList() );
    }

Здесь session - указатель на объект типа G4UIsession, histoMenuHandle и histoMenuLabel - строки, первая задается произвольно - это дескриптор нового главного меню в строке меню (menu bar), вторая - имя этого меню, в нашем случае она равна Histo.

Этого кода еще нет в составе ChargeExchangeMC, но я думаю добавить его к следующему релизу Geant4.

CERN ROOT и надписи на русском языке

Как это ни странно, но проблема поддержки кириллицы в CERN ROOT действительно существует. Я не буду утверждать, что такая проблема присутствует в версии ROOT для Windows, но в Linux она актуальна. Однако, если учесть как эта поддержка реализована в нативной графической подсистеме ROOT, а также то, насколько глючно реализован Qt бекэнд, который поддерживает более широкий набор шрифтов, чем нативный вариант, такая проблема должна существовать и в Windows.

В качестве примера возьмем простой скрипт, с помощью которого можно создать следующее изображение:

Будем считать, что нашей целью является создание изображения в формате eps или pdf с целью дальнейшей вставки его в статью для публикации. Сам скрипт выглядит следующим образом:

void  createHisto( void )
{
    TCanvas *  c1 = new TCanvas();
    TH1 *      histo = new TH1F( "h1", "DCS vs. Energy", 10, 400., 600. );
    Double_t   content[] = { 0.0, 0.115, 0.145, 0.158, 0.170, 0.182, 0.180,
                             0.172, 0.154, 0.140, 0.132, 0.0 };
    Double_t   errors[] = { 0.0, 0.0015, 0.0012, 0.0008, 0.0016, 0.0023, 0.0018,
                            0.0013, 0.0014, 0.0006, 0.0002, 0.0 };
    histo->SetContent( content );
    histo->SetError( errors );
    histo->GetXaxis()->SetTitle( "E_{cm}, MeV" );
    histo->GetYaxis()->SetTitle( "d#sigma/d#Omega, mb/sr" );
    histo->SetStats( kFALSE );          // do not draw stats box
    histo->Draw();
    c1->SaveAs( "c1.eps" );
}

void histoTest( void )
{
    gROOT->SetStyle( "Plain" );
    gStyle->SetTitleBorderSize( 0 );    // do not draw frames around title
    gStyle->SetTitleOffset( 1.2, "xy" );// do not let axes' titles overlap ticks
    createHisto();
    if ( gROOT->IsBatch() )
        exit( 0 );
}

Обратите внимание, что в подписях к осям гистограммы используется вариант Latex для ROOT, стиль gStyle установлен в Plain, а рамка вокруг титула удалена. Отсутствие рамки вокруг титула упрощает решение задачи в первом и третьем подходах (см. ниже).

Для сохранения картинки наверху в файле c1.eps нужно создать скрипт histoTest.C с указанным содержимым и запустить root в пакетном режиме:

root -b histoTest.C

Пусть нашей задачей будет замена титула изображения на ДСР в зависимости от энергии и меток mb/sr, MeV и cm на мб/ср, МэВ и цм соответственно.

Наивная подстановка русских надписей в скрипт histoTest.C приведет к выводу в изображение c1.eps "кракозябр". Что же не так с ROOT? Нативный графический бэкенд ROOT не поддерживает кириллицу. Точка. Для решения проблемы ниже представлены три подхода, которые приводятся в порядке улучшения качества решения и полноты, хотя ни один из них не является идеальным. Идеи второго и третьего подходов взяты из обсуждений здесь и здесь.

1. Наиболее простое решение - сохранить исходное изображение в векторном формате svg и отредактировать его в каком-нибудь векторном редакторе, например в inkscape. Минусы здесь очевидны. Редактировать придется вручную. Кроме изменения текста придется изменять обрамляющие его элементы, например компоненты рамки вокруг титула (которую мы предусмотрительно убрали в исходном скрипте).
2. Использование Qt бэкенда (Qt layer) . Это был бы идеальный вариант, если бы не невообразимая глючность данного бэкенда. Из основных глюков стоит упомянуть: полупрозрачный или полностью прозрачный фон и шрифт абсолютно всех меню, невозможность использовать диалог сохранения файлов (!), невозможность сохранить валидное изображение в форматах eps и pdf (!). В данном подходе пришлось сохранять изображение как Save -> c1.png. Для использования Qt бэкенда нужно записать в файл ~/.rootrc следующие строки:
   

   Gui.Backend:                $(ROOTGUI)
   Gui.Factory:                $(ROOTGUI)
   Gui.DefaultFont:            -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-5
   Gui.MenuFont:               -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-5
   Gui.MenuHiFont:             -*-helvetica-bold-r-*-*-12-*-*-*-*-*-iso8859-5
   Gui.DocFixedFont:           -*-courier-medium-r-*-*-12-*-*-*-*-*-iso8859-5
   Gui.DocPropFont:            -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-5
   Gui.IconFont:               -*-helvetica-medium-r-*-*-10-*-*-*-*-*-iso8859-5
   Gui.StatusFont:             -*-helvetica-medium-r-*-*-10-*-*-*-*-*-iso8859-5
   

   Использование переменной среды ROOTGUI позволяет сохранить нативный бэкенд по умолчанию (однако придется добавить строку
   

   export ROOTGUI=native
   

   в ~/.bash_profile). Для того, чтобы запустить root с Qt бэкендом нужно ввести
   

   ROOTGUI=qt root histoTest_qt.C
   

   Скрипт histoTest_qt.C выглядит так:
   

   void  createHisto( void )
   {
       //TCanvas *  c1 = new TCanvas();
       TH1 *      histo = new TH1F( "h1", "´ÁÀ Ò ×ÐÒØáØÜÞáâØ Þâ íÝÕàÓØØ",
                                    10, 400., 600. );
       Double_t   content[] = { 0.0, 0.115, 0.145, 0.158, 0.170, 0.182, 0.180,
                                0.172, 0.154, 0.140, 0.132, 0.0 };
       Double_t   errors[] = { 0.0, 0.0015, 0.0012, 0.0008, 0.0016, 0.0023, 0.0018,
                               0.0013, 0.0014, 0.0006, 0.0002, 0.0 };
       histo->SetContent( content );
       histo->SetError( errors );
       histo->GetXaxis()->SetTitle( "E_{æÜ}, ¼í²" );
       histo->GetYaxis()->SetTitle( "d#sigma/d#Omega, ÜÑ/áà" );
       histo->SetStats( kFALSE );          // do not draw stats box
       histo->Draw();
       //c1->SaveAs( "c1_qt.eps" );
   }
   
   void histoTest_qt( void )
   {
       gROOT->SetStyle( "Plain" );
       gStyle->SetTitleBorderSize( 0 );    // do not draw frames around title
       gStyle->SetTitleOffset( 1.1, "xy" );// do not let axes' titles overlap ticks
       createHisto();
       if ( gROOT->IsBatch() )
           exit( 0 );
   }
   

   Заметили "кракозябры" внутри скрипта? Это еще одно неудобство. Скрипт выведен в общепринятой на сегодня кодировке UTF-8, а был записан в кодировке ISO-8859-5. В vim сохранение файла в данной кодировке обеспечивается командой
   

   :set fenc=iso-8859-5
   

   Однако при повторной загрузке этого файла в vim появятся всё те же "кракозябры" и новые русские записи не удастся добавить без полного переписывания русских символов заново. Возможно существуют редакторы, в которых проблема перевода в разные кодировки решается проще.
   
   Еще один интересный момент. Создание канваса c1 в данном варианте закомментировано. Если оставить его незакомментированным, то диапазон оси ординат изменится, и на месте гистограммы окажется чистое белое поле (!). Это еще раз подтверждает тот факт, что с текущей реализацией Qt бэкенда нужно быть очень осторожным.
   
   Изображение, полученное с помощью Qt бэкенда, выглядит так:
   
   
   
   
   Заметьте также, что Qt бэкенд применил преобразования локали: вместо точек в вещественных числах на координатных осях стоят запятые. 
3. Третий подход основан на Latex и двух латеховских модулях - standalone и psfrag. Оба модуля можно найти на ctan.org. В пакет texlive для Fedora 14 они не входят, однако существует отдельный репозиторий texlive, в котором эти модули присутствуют. Информация об этом репозитории здесь. Преимущество данного подхода в том, что он полностью автоматизирован. Скрипт histoTest_latex.C выглядит так:
   

   void  createHisto( void )
   {
       TCanvas *  c1 = new TCanvas();
       TH1 *      histo = new TH1F( "h1", "CT", 10, 400., 600. );
       Double_t   content[] = { 0.0, 0.115, 0.145, 0.158, 0.170, 0.182, 0.180,
                                0.172, 0.154, 0.140, 0.132, 0.0 };
       Double_t   errors[] = { 0.0, 0.0015, 0.0012, 0.0008, 0.0016, 0.0023, 0.0018,
                               0.0013, 0.0014, 0.0006, 0.0002, 0.0 };
       histo->SetContent( content );
       histo->SetError( errors );
       histo->GetXaxis()->SetTitle( "CX" );
       histo->GetYaxis()->SetTitle( "CY" );
       histo->SetStats( kFALSE );          // do not draw stats box
       histo->Draw();
       c1->SaveAs( "c1_latex.eps" );
   }
   
   void histoTest_latex( void )
   {
       gROOT->SetStyle( "Plain" );
       gStyle->SetTitleBorderSize( 0 );    // do not draw frames around title
       gStyle->SetTitleOffset( 1.2, "xy" );// do not let axes' titles overlap ticks
       createHisto();
       if ( gROOT->IsBatch() )
           exit( 0 );
   }
   

   Главным отличием от исходного скрипта является то, что оригинальные титул и подписи к осям заменены метками CT, CX и CY, которые будут обрабатываться в теховском скрипте с помощью psfrag. Идея psfrag весьма проста - заменить метки в исходном ps или eps файле указанными значениями. Отмечу сразу, что метки для подстановки следует выбирать короткие, иначе psfrag может сгенерировать ломаный файл dvi. Скрипт c1_latex.tex выглядит так:
   

   \documentclass{standalone}
   
   \usepackage[T2A]{fontenc}
   \usepackage[utf8x]{inputenc}
   
   \usepackage[dvips]{graphicx,color}
   \usepackage{psfrag}
   
   \begin{document}
   \psfrag{CT}[lt][lt][1.8]{\bf ДСР в зависимости от энергии}
   \psfrag{CX}[rt][rt][1.4]{\bf $\mathbf{E_{\mbox{цм}}}$, МэВ}
   \psfrag{CY}[rt][rt][1.4]{\bf $\mathbf{d\sigma/d\Omega}$, мб/ср}
   \includegraphics{c1_latex.eps}
   \end{document}
   

   О параметрах psfrag можно прочитать в документации к данному пакету; кратко, первый параметр в фигурных скобках является заменяемой меткой, последний параметр в фигурных скобках - подстановкой, первые два параметра в квадратных скобках после метки - позиционными якорями, третий параметр в квадратных скобках - коэффициентом масштабирования шрифта. Пакет standalone позволяет сохранить исходные размеры изображения.
   
   Русификация надписей может быть произведена с помощью отдельного шелл-скрипта mkrus:
   

   #!/bin/sh
   # source root macro
   root -b histoTest_latex.C
   # run psfrag
   latex c1_latex
   # make pdf
   dvips c1_latex
   ps2pdf c1_latex.ps
   # delete temporary files
   rm c1_latex.eps c1_latex.ps c1_latex.dvi c1_latex.aux c1_latex.log
   # make eps
   pdftops -eps c1_latex.pdf
   

   Этот скрипт создает два изображения: c1_latex.eps и c1_latex.pdf, которые выглядят так:
   
   
   
   
   В данном изображении используется шрифт Serif - просто в моем дистрибутиве texlive не нашлось русского Sans Serif шрифта c кодировкой T2A (которая используется в c1_latex.tex). Если необходим шрифт Sans Serif, нужно установить его из внешних источников.
   

Вышел Geant4 9.4 и ChargeExchangeMC (aka Cexmc)

17 декабря вышла новая версия библиотеки для моделирования прохождения элементарных частиц в веществе Geant4 9.4. Geant4 широко применяется при моделировании физических экспериментов, а также в медицине и космофизике.
Для меня данный релиз замечателен тем, что в него в качестве одного из advanced examples вошла программа ChargeExchangeMC, автором которой я являюсь. Программа применялась для моделирования реального эксперимента в Петербургском Институте Ядерной Физики. Первая версия, написанная на фортране, была в дальнейшем переписана с нуля на C++ с использованием Geant4.
Основные особенности программы:

1. Это достаточно большой проект, он включает 71 заголовочный файл и 63 файла исходников, общий размер программного кода превышает 20 тысяч строк
2. Геометрия установки полностью описана в формате GDML. Это позволит гибко изменять описание установки без изменения исходного кода
3. Для упрощения анализа данных используются гистограммы на основе библиотеки CERN ROOT
4. Реализована модель реконструкции данных с гибкой настройкой различных параметров
5. В программе реализован гибкий алгоритм на основе шаблонов C++ для возможности легкого изменения изучаемых взаимодействий
6. Персистентное хранение данных реализовано с помощью библиотеки boost::serialization
7. Пользовательский фильтр данных позволяет изменять полученные данные с помощью специальных скриптов. Грамматика скриптов описана с помощью библиотеки boost::spirit
8. Run manager имеет уникальный режим проигрывания событий (replay mode), который позволяет пересчитывать уже имеющиеся данные с параметрами, отличными от исходных

Документацию по установке и использованию программы можно найти здесь.

Как прикрутить вызов фортрановской функции из интерпретатора CINT

CERN ROOT - программа для анализа больших объемов данных, ориентированная прежде всего на использование в научной среде. Программа обладает широкими возможностями по визуализации данных и имеет встроенный интерпретатор языка C++ CINT. Как известно, в физике до сих пор имеет широкое применение язык программирования FORTRAN, на котором написано значительное количество интересных библиотек. Однажды мой коллега-физик попросил меня помочь разобраться в этом вопросе - ему нужно было вызывать функции фортрановской библиотеки MINUIT (разработанной тоже в CERN) непосредственно из интерпретатора CINT. Благодаря поискам решения этой задачи родился этот пример.
Итак, прежде всего создадим функцию на фортране, которая будет вызываться из CINT. Пусть это будет функция, складывающая два числа, также создадим в ней COMMON блок для того, чтобы иметь возможность поэкспериментировать с изменением глобальных переменных:

       INTEGER FUNCTION SUM( A, B )
       DOUBLE PRECISION GL1
       DOUBLE PRECISION GL2
       INTEGER GL3
       COMMON /GL/ GL1( 1, 2 ), GL2, GL3

       INTEGER A, B
       SUM = A + B
       GL1( 1, 1 ) = 0.1
       GL1( 1, 2 ) = 0.2
       GL2 = 10.89
       GL3 = 8
       RETURN
       END

Разместим этот код в файле Sum.f. Теперь нам нужно создать код на C++, который, во-первых, будет вызывать функцию SUM(), а во-вторых, сможет быть вызван из интерпретатора CINT. Пусть вызов SUM() будет организован в виде функтора (т.е. класса с определенным оператором ()) - это упростит семантику вызова из CINT. Разместим наш код в файлах Sum.h и Sum.c. Это содержимое Sum.h:

#include <iostream>
#include "TObject.h"


extern "C" { extern int  sum_( int*, int* ); }

extern struct  fort_struct { double  a[ 2 ][ 1 ]; double  b; int  c; }  gl_;


class  Sum : public  TObject
{
    public:
        int  operator()( int  a, int  b )
        {
            return sum_( &a, &b );
        }

        void  accessGl( void ) const
        {
            printGl();
            gl_.a[ 0 ][ 0 ] = 34.9;
            gl_.a[ 1 ][ 0 ] = 4.1;
            gl_.b = 0.88;
            gl_.c = 70;
            printGl();
        }

    private:
        void  printGl( void ) const
        {
            std::cout << "gl_ = " << gl_.a[ 0 ][ 0 ] << ", " <<
                         gl_.a[ 1 ][ 0 ] << ", " << gl_.b << ", " << gl_.c <<
                         std::endl;
        }

        ClassDef( Sum, 1 )
};

Это содержимое Sum.c:

#include "Sum.h"

ClassImp( Sum )

Сразу отметим наличие определений ClassDef() и ClassImp() - это то, что позволит сделать из нашего Sum объект, с которым интерпретатору CINT будет комфортно работать.
Теперь замечания по использованию фортрановского кода в коде C++:

1. Фортрановские функции должны связываться с C++ кодом как extern "C", при этом название функции записывается в нижнем регистре (несмотря на то, что в исходном коде оно может быть в верхнем регистре), а в конце добавляется символ _ - именно так компилятор фортрана записывает названия функций в объектном коде. В нашем случае функция SUM() трансформировалась в sum_().
2. Аргументы передаются в фортрановскую функцию через указатели, поэтому функция двух целочисленных аргументов SUM() превратилась в sum_( int*, int* ).
3. Глобальные переменные, объявленные в блоке COMMON, также записываются в нижнем регистре с добавлением символа _. Также обратите пристальное внимание на то, что столбцы и строки многомерных массивов в C++ и фортране меняются местами, так, переменная GL1( 1, 2 ) стала double a[ 2 ][ 1 ] в коде C++.

Исходный код готов, теперь нам нужно создать из него подключаемую библиотеку, ниже приводим содержимое Makefile:

all:
 rootcint -f dict_sum.cxx -c Sum.h
 g++ -fPIC -I$(ROOTSYS)/include/root -c dict_sum.cxx Sum.c
 g77 -fPIC -c Sum.f -o Sum_f.o
 g++ -shared -o Sum.so dict_sum.o Sum.o Sum_f.o

clean:
 rm *.o *.so dict_sum.cxx dict_sum.h

Обратите внимание на использование rootcint. Эта утилита создает вспомогательный код, необходимый для использования Sum в CINT, название файла dict_sum.cxx выбрано произвольно. Фортрановский код компилируется с помощью компилятора g77, но ничто не мешает заменить его на gfortran, если это необходимо. После запуска make создается объект Sum.so, который мы будем загружать в CINT, а затем использовать определенные в нем объекты. Пример сессии ROOT:

root [0] .L Sum.so
root [1] Sum s
root [2] int res = s(2,4)
root [3] .p res
(int)6
root [4] s.accessGl()
gl_ = 0.1, 0.2, 10.89, 8
gl_ = 34.9, 4.1, 0.88, 70
root [5]