fox/BluetAnalyzer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Init

Browse Source
This commit is contained in:
liuyihui 2024-11-23 16:56:50 +08:00
commit 58ee1a2a2a
9 changed files with 3003 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.clang-format Normal file
View File

@ -0,0 +1,4 @@
---
Language: Cpp
BasedOnStyle: llvm
...

9
.gitignore vendored Normal file
View File

@ -0,0 +1,9 @@
# Config
.vscode
cfg
# Data
data
# build
build/

78
CMakeLists.txt Normal file
View File

@ -0,0 +1,78 @@
cmake_minimum_required(VERSION 3.8...3.18)
project(BluetAnalyzer)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_FLAGS_DEBUG "$ENV{CXXFLAGS} -O0 -Wall -g -ggdb")
set(CMAKE_CXX_FLAGS_RELEASE "$ENV{CXXFLAGS} -O2 -Wall")
set(CMAKE_BUILD_TYPE "Debug")
# Find fmt
find_package(fmt)
# Find Eigen
find_package(Eigen3)
include_directories(/usr/include/eigen3)
# Find Geant4
find_package(Geant4 REQUIRED ui_all vis_all)
# Find ROOT
find_package(ROOT REQUIRED COMPONENTS RIO Net Spectrum Gui)
# Find Garfield
find_package(Garfield REQUIRED)
include_directories($ENV{GARFIELD_INSTALL}/include/Garfield)
# Setup include directory for this project
include(${Geant4_USE_FILE})
include_directories(${PROJECT_SOURCE_DIR}/include)
# Include bluet file
set(bluet_base "/home/liuyihui/bluet/")
# include(${bluet_base}/sources/config/include/BluetConfig.hh)
# include(${bluet_base}/sources/modules/include/BluetDataModel.hh)
include_directories(${bluet_base}/sources/config/include)
include_directories(${bluet_base}/sources/modules/include)
include_directories(${bluet_base}/sources/runner/include)
include_directories(${bluet_base}/utils/include)
include_directories($ENV{GARFIELD_INSTALL}/include/Garfield)
# Locate sources and headers for this project
file(GLOB analyzer ${PROJECT_SOURCE_DIR}/include/* ${PROJECT_SOURCE_DIR}/src/*)
# file(GLOB config ${bluet_base}/sources/config/include/* ${bluet_base}/sources/config/src/*)
# file(GLOB modules ${bluet_base}/sources/modules/include/* ${bluet_base}/sources/modules/src/*)
# file(GLOB runner ${bluet_base}/sources/runner/include/* ${bluet_base}/sources/runner/src/*)
# file(GLOB utils ${bluet_base}/utils/include/* ${bluet_base}/utils/src/*)
list(FILTER analyzer EXCLUDE REGEX "BluetDrawLinkDef")
# ROOT_GENERATE_DICTIONARY(D_${PROJECT_NAME} /home/fox/bluet/sources/modules/include/DrawEvent.hh LINKDEF /home/fox/bluet/sources/modules/include/BluetDrawLinkDef.hh)
add_library(SO_${PROJECT_NAME} SHARED
${analyzer}
# Bluet
${bluet_base}/sources/config/include/BluetDataModel.hh
${bluet_base}/utils/include/pugixml.hpp
${bluet_base}/utils/src/pugixml.cpp
${bluet_base}/utils/include/xmlparse.hh
${bluet_base}/utils/src/xmlparse.cc
${bluet_base}/utils/include/stringhandle.hh
${bluet_base}/utils/src/stringhandle.cc
# D_${PROJECT_NAME}.cxx
)
target_link_libraries(SO_${PROJECT_NAME} PUBLIC
${Geant4_LIBRARIES}
${ROOT_LIBRARIES} ROOT::RIO ROOT::Net ROOT::Spectrum ROOT::Gui
Garfield::Garfield
fmt::fmt
gsl gslcblas
)
# Add the executable, and link it to the libraries
add_executable(${PROJECT_NAME} main.cc)
target_link_libraries(${PROJECT_NAME} SO_${PROJECT_NAME})
#--- Set default installation prefix ------------------------------------------
set(CMAKE_INSTALL_PREFIX
${PROJECT_SOURCE_DIR}/install
CACHE PATH "Install path prefix, prepended onto install directories." FORCE)
INSTALL(TARGETS ${PROJECT_NAME} SO_${PROJECT_NAME})

151
include/BluetAnalyzer.hh Normal file
View File

@ -0,0 +1,151 @@
//---------------------------------------------------------------------------
// version: v1.0
// code log:
// 1) read and study on the BluetFactory output data in BluetDataModel format
// 2) get the phi30 beam profile from 6Li(n,t) in the low energy region
// 3) get drift velocity
// author: yih@ihep.ac.cn
// date: 2023-04-04
//--------------------------------------------------------------------------
// #include "../../sources/config/include/BluetDataModel.hh"
#include "ROOT/TDataFrame.hxx"
#include "TApplication.h"
#include "TChain.h"
#include "TF1.h"
#include "TF2.h"
#include "TGTab.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TInterpreter.h"
#include "fmt/core.h"
#include "string"
#include "BluetDataModel.hh"
const int sizeX = 900;
const int sizeY = 900;
const double rBeam = 35. / 2;
const double rTarget = 73. / 2;
const double rSubstrate = 89. / 2;
const double rPCB = rSubstrate + 10.;
class BluetAnalyzer {
public:
BluetAnalyzer() {
// Declare main frame
int argc = 0;
char **argv = nullptr;
app = new TApplication("app", &argc, argv);
mFrame = new TGMainFrame(gClient->GetRoot(), sizeX, sizeY);
mTab = new TGTab(mFrame, sizeX, sizeY);
mFrame->AddFrame(mTab, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
// Declare Data Models
const string bluetBase = getenv("bluet_base");
const string datamodel = fmt::format(
"#include \"{}/sources/config/include/BluetDataModel.hh\"", bluetBase);
gInterpreter->Declare(datamodel.c_str());
};
~BluetAnalyzer();
private:
TApplication *app;
TGMainFrame *mFrame;
TGTab *mTab;
double sidelength;
int Nlayer;
double driftL;
double sampleT;
double EcutMin;
double EcutMax;
double tofCutMin;
double tofCutMax;
double stampCutMin;
double stampCutMax;
double cathAmpCutMin;
double cathAmpCutMax;
double cathTimeCutMax;
double cathTimeCutMin;
double meshAmpCutMin;
double meshAmpCutMax;
double meshTimeCutMax;
double meshTimeCutMin;
int NtrkCutMin;
int NtrkCutMax;
double lengthCutMin;
double lengthCutMax;
bool rCut;
double rCutMin;
double rCutMax;
double cosCutMin;
double CenterPos[2];
double sliceW;
TF1 *fbeam1;
TF2 *fbeam2;
TF1 *fbeamr;
TF1 *fdgaus;
TChain *eventTree = new TChain("eventTree");
BluetEvent *fEvent = nullptr;
int Nevent;
int Nbin;
double tmin;
double tmax;
int Nrbin;
double rmax;
double *Sbin;
TH2F *hprofile;
TH1F *hslicex;
TH1F *hslicey;
TH1F *hbeamr;
TH1F *hrinte;
TH2F *htofE;
TH1F *htof;
TH1F *htofr0;
TH1F *htofr1;
TH1F *hpull;
TH1F *hratio;
TH2F *hLengthE;
TH1F *htcath;
TH1F *htpad;
TH1F *hv0;
TH1F *hv1;
TH1F *hcos;
TH1F *hphi;
TH2F *hdeltacos1;
TH2F *hdeltacos2;
TH2F *hts;
TH2F *htn;
TH2F *hamp;
TH2F *han;
TH2F *hta;
TH2F *htt;
public:
void readCutParameters(string cfgfile);
void defineFitFunctions();
void defineHistograms();
bool isGoodEvent(BluetEvent *event);
bool isGoodTrack(BluetTrack *track);
void readTreeData(vector<string> files);
void FillHistograms();
void postOperations();
void drawHistograms();
private:
void drawVertex();
void drawTOF();
void drawDrift();
void drawAngle();
void draw1DFeatures();
void draw2DFeatures();
void drawPadFeatures();
};

1414
include/pgbar/pgbar.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

346
include/pgbar/range.hpp Normal file
View File

@ -0,0 +1,346 @@
// This code is licensed under the MIT License.
// Please see the LICENSE file in the root of the repository for the full license text.
// Copyright (c) 2023-2024 Konvt
#pragma once
#ifndef __PGBAR_RANGE_HPP__
# define __PGBAR_RANGE_HPP__
#include "pgbar.hpp" // Other required header files have been included in pgbar.hpp.
#include <algorithm> // std::distance
#if defined(__GNUC__) || defined(__clang__)
# define __PGBAR_NODISCARD__ __attribute__((warn_unused_result))
#elif defined(_MSC_VER)
# define __PGBAR_NODISCARD__ _Check_return_
#else
# define __PGBAR_NODISCARD__
#endif
#if defined(_MSVC_VER) && defined(_MSVC_LANG) // for msvc
# define __PGBAR_CMP_V__ _MSVC_LANG
#else
# define __PGBAR_CMP_V__ __cplusplus
#endif
#if __PGBAR_CMP_V__ >= 202002L
# include <ranges> // std::ranges::begin、std::ranges::end
# define __PGBAR_CXX20__ 1
#else
# define __PGBAR_CXX20__ 0
#endif // __cplusplus >= 202002L
namespace pgbar {
namespace __detail {
template<typename EleT, typename BarT>
class range_iterator final { // for number
static_assert(
is_pgbar<BarT>::value,
"pgbar::__detail::range_iterator: Only available for `pgbar::pgbar`"
);
numeric_iterator<EleT> iter;
BarT* bar;
public:
using iterator_category = std::output_iterator_tag;
using value_type = EleT;
using difference_type = void;
using pointer = void;
using reference = value_type;
explicit range_iterator( EleT _startpoint, EleT _endpoint, EleT _step, BarT& _bar )
: iter { _startpoint, _endpoint, _step } {
if ( _endpoint < _startpoint && _step > 0 )
throw bad_pgbar { "pgbar::__detail::range_iterator: invalid iteration range" };
bar = &_bar;
bar->reset().set_task( iter.extent() ).set_step( 1 );
}
__PGBAR_NODISCARD__ range_iterator begin() const noexcept {
return *this;
} // invokes copy constructor
__PGBAR_NODISCARD__ range_iterator end() const noexcept {
range_iterator ed_pnt = *this;
ed_pnt.iter = ed_pnt.iter.end();
return ed_pnt;
}
__PGBAR_NODISCARD__ reference operator*() const noexcept {
return *iter;
}
__PGBAR_NODISCARD__ bool operator==( const range_iterator& lhs ) const noexcept {
return iter == lhs.iter;
}
__PGBAR_NODISCARD__ bool operator!=( const range_iterator& lhs ) const noexcept {
return !(operator==( lhs ));
}
range_iterator& operator++() {
++iter; bar->update(); return *this;
}
range_iterator operator++( int ) {
range_iterator before = *this; operator++(); return before;
}
};
template<typename IterT, typename BarT>
class container_iterator final { // for container
static_assert( // `IterT` means iterator type
!std::is_arithmetic<IterT>::value &&
is_pgbar<BarT>::value,
"pgbar::__detail::container_iterator: Only available for container types"
);
using SizeT = size_t;
using EleT = typename std::iterator_traits<IterT>::value_type;
BarT* bar;
IterT start, terminus, current;
SizeT extent;
bool is_reversed;
public:
using iterator_category = std::output_iterator_tag;
using value_type = EleT;
using difference_type = void;
using pointer = EleT*;
using reference = EleT&;
explicit container_iterator( IterT _begin, IterT _endpoint, BarT& _bar ) {
static_assert(
!std::is_same<typename std::iterator_traits<IterT>::difference_type, void>::value,
"pgbar::__detail::container_iterator: the difference_type of the iterator shouldn't be 'void'"
);
auto dist = std::distance( _begin, _endpoint );
if ( dist > 0 ) {
is_reversed = false;
extent = static_cast<SizeT>(dist);
} else {
is_reversed = true;
extent = static_cast<SizeT>(-dist);
}
start = _begin; terminus = _endpoint;
current = start; bar = &_bar;
bar->set_task( extent ).set_step( 1 );
}
container_iterator( const container_iterator& _other ) {
bar = _other.bar;
start = _other.start; terminus = _other.terminus;
current = start; extent = _other.extent;
is_reversed = _other.is_reversed;
}
container_iterator( container_iterator&& _rhs ) {
bar = _rhs.bar; // same as copy constructor
start = std::move( _rhs.start ); terminus = std::move( _rhs.terminus );
current = std::move( _rhs.current ); extent = _rhs.extent;
is_reversed = _rhs.is_reversed;
_rhs.bar = nullptr; // clear up
_rhs.start = {}; _rhs.terminus = {};
_rhs.terminus = {}; _rhs.extent = 0;
_rhs.is_reversed = false;
}
~container_iterator() {
bar = nullptr;
}
__PGBAR_NODISCARD__ container_iterator begin() const {
return container_iterator( *this );
} // invokes copy constructor
__PGBAR_NODISCARD__ container_iterator end() const {
auto ed_pnt = container_iterator( *this );
ed_pnt.start = terminus;
ed_pnt.current = terminus;
return ed_pnt;
}
__PGBAR_NODISCARD__ reference operator*() noexcept {
return *current;
}
__PGBAR_NODISCARD__ pointer operator->() noexcept {
return std::addressof( current );
}
__PGBAR_NODISCARD__ bool operator==( const container_iterator& _other ) const noexcept {
return current == _other.current;
}
__PGBAR_NODISCARD__ bool operator!=( const container_iterator& _other ) const noexcept {
return !(operator==( _other ));
}
container_iterator& operator++() {
if ( is_reversed ) --current; else ++current; bar->update(); return *this;
}
container_iterator operator++( int ) {
container_iterator before = *this; operator++(); return before;
}
};
#if __PGBAR_CXX20__
template<typename T>
concept ArithmeticType =
std::is_arithmetic_v<std::decay_t<T>>;
#endif // __PGBAR_CXX20__
} // namespace __detail
/// @brief Update the progress bar based on the range specified by the parameters.
/// @tparam EleT The type of generated elements.
/// @tparam BarT The type of the progress bar.
/// @param _bar The progress bar that will be updated.
/// @return Return an iterator that moves unidirectionally within the range `[_startpoint, _endpoint-1]`.
template<
#if __PGBAR_CXX20__
__detail::ArithmeticType _EleT, __detail::PgbarType BarT
, typename EleT = std::decay_t<_EleT>
> __PGBAR_NODISCARD__ __detail::range_iterator<EleT, BarT> // return type
#else
typename _EleT, typename BarT
, typename EleT = typename std::decay<_EleT>::type
> __PGBAR_NODISCARD__ typename std::enable_if<
std::is_arithmetic<EleT>::value &&
is_pgbar<BarT>::value
, __detail::range_iterator<EleT, BarT>>::type
#endif // __PGBAR_CXX20__
inline range( _EleT&& _startpoint, _EleT&& _endpoint, _EleT&& _step, BarT& _bar ) {
return __detail::range_iterator<EleT, BarT>( _startpoint, _endpoint, _step, _bar );
}
#if __PGBAR_CXX20__
template<typename _EleT, __detail::PgbarType BarT
, typename EleT = typename std::decay_t<_EleT>
> requires std::is_arithmetic_v<EleT>
__PGBAR_NODISCARD__ __detail::range_iterator<EleT, BarT> // return type
#else
template<
typename _EleT, typename BarT
, typename EleT = typename std::decay<_EleT>::type
> __PGBAR_NODISCARD__ typename std::enable_if<
std::is_arithmetic<EleT>::value &&
is_pgbar<BarT>::value
, __detail::range_iterator<EleT, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( _EleT&& _endpoint, _EleT&& _step, BarT& _bar ) {
return __detail::range_iterator<EleT, BarT>( {}, _endpoint, _step, _bar );
}
#if __PGBAR_CXX20__
template<typename _EleT, __detail::PgbarType BarT
, typename EleT = typename std::decay_t<_EleT>
> requires std::integral<EleT>
__PGBAR_NODISCARD__ __detail::range_iterator<EleT, BarT> // return type
#else
template<
typename _EleT, typename BarT
, typename EleT = typename std::decay<_EleT>::type
> __PGBAR_NODISCARD__ typename std::enable_if<
std::is_integral<EleT>::value &&
is_pgbar<BarT>::value
, __detail::range_iterator<EleT, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( _EleT&& _startpoint, _EleT&& _endpoint, BarT& _bar ) {
return __detail::range_iterator<EleT, BarT>( _startpoint, _endpoint, 1, _bar );
}
#if __PGBAR_CXX20__
template<typename _EleT, __detail::PgbarType BarT
, typename EleT = typename std::decay_t<_EleT>
> requires std::integral<EleT>
__PGBAR_NODISCARD__ __detail::range_iterator<EleT, BarT> // return type
#else
template<
typename _EleT, typename BarT
, typename EleT = typename std::decay<_EleT>::type
> __PGBAR_NODISCARD__ typename std::enable_if<
std::is_integral<EleT>::value &&
is_pgbar<BarT>::value
, __detail::range_iterator<EleT, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( _EleT&& _endpoint, BarT& _bar ) {
return __detail::range_iterator<EleT, BarT>( {}, _endpoint, 1, _bar );
}
/// @brief Accepts the beginning and end iterators,
/// @brief and updates the passed `_bar` based on the range defined by these two iterators.
/// @tparam IterT The type of the iterators.
template<
typename _BeginT, typename _endpointT, typename BarT
, typename IterT = typename std::decay<_BeginT>::type
#if __PGBAR_CXX20__
> requires std::conjunction_v<
std::negation<std::is_arithmetic<IterT>>,
std::is_same<IterT, std::decay_t<_endpointT>>,
is_pgbar<BarT>
> __PGBAR_NODISCARD__ __detail::container_iterator<IterT, BarT>
#else
> __PGBAR_NODISCARD__ typename std::enable_if<
!std::is_arithmetic<IterT>::value &&
std::is_same<IterT, typename std::decay<_endpointT>::type>::value &&
is_pgbar<BarT>::value
, __detail::container_iterator<IterT, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( _BeginT&& _startpoint, _endpointT&& _endpoint, BarT& _bar ) {
return __detail::container_iterator<IterT, BarT>( _startpoint, _endpoint, _bar );
}
/// @brief Accepts a iterable container,
/// @brief and updates the passed `_bar` based on the elements in the container.
/// @param container A lvalue container whose iterator type name is iterator.
/// @param _bar A pgbar object.
/// @tparam ConT The type of the container.
template<typename ConT, typename BarT>
#if __PGBAR_CXX20__
requires std::conjunction_v<
std::is_lvalue_reference<ConT>,
std::negation<std::is_array<std::remove_reference_t<ConT>>>,
is_pgbar<BarT>
> __PGBAR_NODISCARD__ __detail::container_iterator<typename std::decay_t<ConT>::iterator, BarT>
#else
__PGBAR_NODISCARD__ typename std::enable_if<
std::is_lvalue_reference<ConT>::value &&
!std::is_array<typename std::remove_reference<ConT>>::value &&
is_pgbar<BarT>::value
, __detail::container_iterator<typename std::decay<ConT>::type::iterator, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( ConT&& container, BarT& _bar ) {
#if __PGBAR_CXX20__
return range( std::ranges::begin( std::forward<ConT>( container ) ), std::ranges::end( std::forward<ConT>( container ) ), _bar );
#else
using std::begin; using std::end; // ADL
return range( begin( std::forward<ConT>( container ) ), end( std::forward<ConT>( container ) ), _bar );
#endif // __PGBAR_CXX20__
}
/// @brief Accepts a original array,
/// @brief and updates the passed `_bar` based on the elements in the container.
/// @param container An original array.
/// @param _bar A pgbar object.
/// @tparam ArrT The type of the array.
template<typename ArrT, typename BarT>
#if __PGBAR_CXX20__
requires std::conjunction_v<
std::is_lvalue_reference<ArrT>,
std::is_array<std::remove_reference_t<ArrT>>,
is_pgbar<BarT>
> __PGBAR_NODISCARD__ __detail::container_iterator<std::decay_t<ArrT>, BarT> // here's the difference between array type and container type
#else
__PGBAR_NODISCARD__ typename std::enable_if<
std::is_lvalue_reference<ArrT>::value &&
std::is_array<typename std::remove_reference<ArrT>::type>::value &&
is_pgbar<BarT>::value
, __detail::container_iterator<typename std::decay<ArrT>::type, BarT>
>::type
#endif // __PGBAR_CXX20__
inline range( ArrT&& container, BarT& _bar ) {// for original arrays
#if __PGBAR_CXX20__
return range( std::ranges::begin( std::forward<ArrT>( container ) ), std::ranges::end( std::forward<ArrT>( container ) ), _bar );
#else
using std::begin; using std::end; // ADL
return range( begin( std::forward<ArrT>( container ) ), end( std::forward<ArrT>( container ) ), _bar );
#endif // __PGBAR_CXX20__
}
} // namespace pgbar
#undef __PGBAR_CMP_V__
#undef __PGBAR_CXX20__
#endif

28
main.cc Normal file
View File

@ -0,0 +1,28 @@
#include "BluetAnalyzer.hh"
#include "iostream"
using namespace std;
int main(int argc, char **argv) {
if (argc < 3) {
cout << "Usage: " << argv[0] << " <config file> <data file>" << endl;
return 1;
}
string cfgfile = argv[1];
string datafile = argv[2];
vector<string> runfiles;
runfiles.push_back(datafile);
BluetAnalyzer *bluet = new BluetAnalyzer();
bluet->readCutParameters(cfgfile);
bluet->defineFitFunctions();
bluet->defineHistograms();
bluet->readTreeData(runfiles);
bluet->FillHistograms();
bluet->postOperations();
bluet->drawHistograms();
return 0;
}

965
src/BluetAnalyzer.cc Normal file
View File

@ -0,0 +1,965 @@
#include "BluetAnalyzer.hh"
#include "stringhandle.hh"
#include "xmlparse.hh"
#include "TCanvas.h"
#include "TEllipse.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMath.h"
#include "TRootEmbeddedCanvas.h"
#include "TString.h"
#include "TStyle.h"
#include "pgbar/pgbar.hpp"
#include "stdio.h"
using namespace TMath;
/**
* @brief Computes the sum of two Gaussian functions.
*
* This function calculates the sum of two Gaussian functions with the given
* parameters.
*
* @param par Array of parameters for the Gaussian functions:
* - par[0]: Mean of both two Gaussian.
* - par[1]: Amplitude of the first Gaussian.
* - par[2]: Standard deviation of the first Gaussian.
* - par[3]: Amplitude of the second Gaussian.
* - par[4]: Standard deviation of the second Gaussian.
* @return The sum of the two Gaussian functions evaluated at the input value.
*/
double dgaus(double *x, double *par) {
double ret = par[1] * TMath::Exp(-(x[0] - par[0]) * (x[0] - par[0]) / 2. /
par[2] / par[2]) +
par[3] * TMath::Exp(-(x[0] - par[0]) * (x[0] - par[0]) / 2. /
par[4] / par[4]);
return ret;
}
static double dv0;
static int nInner = 0;
static int nRing = 0;
void fitdgaus(TH1F *hist, TF1 *func, double *&par, double m0) {
double mean = m0 >= 0 ? m0 : hist->GetMean();
double sigma = hist->GetRMS();
double amp = hist->GetMaximum();
double xl = mean - 3 * sigma;
double xr = mean + 3 * sigma;
func->SetParameters(mean, amp / 2, sigma, amp / 2, sigma);
hist->Fit("fdgaus", "LQ", "", xl, xr);
auto pars = func->GetParameters();
mean = m0 >= 0 ? m0 : pars[0];
double amp1 = pars[1];
double sigma1 = pars[2];
double amp2 = pars[3];
double sigma2 = pars[4];
sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
xl = mean - 3. * sigma;
xr = mean + 3. * sigma;
par = new double[7]{mean, amp1, sigma1, amp2, sigma2, xl, xr};
}
void adjustAxisRange(TH2F *hist) {
int nBinsX = hist->GetNbinsX();
int nBinsY = hist->GetNbinsY();
double minX = hist->GetXaxis()->GetXmax();
double maxX = hist->GetXaxis()->GetXmin();
double minY = hist->GetYaxis()->GetXmax();
double maxY = hist->GetYaxis()->GetXmin();
for (int i = 1; i <= nBinsX; ++i) {
for (int j = 1; j <= nBinsY; ++j) {
if (hist->GetBinContent(i, j) != 0) {
double x = hist->GetXaxis()->GetBinCenter(i);
double y = hist->GetYaxis()->GetBinCenter(j);
if (x < minX)
minX = x;
if (x > maxX)
maxX = x;
if (y < minY)
minY = y;
if (y > maxY)
maxY = y;
}
}
}
hist->GetXaxis()->SetRangeUser(minX * 0.9, maxX * 1.1);
hist->GetYaxis()->SetRangeUser(minY * 0.9, maxY * 1.1);
}
void adjustAxisRange(TH1F *hist) {
int nBins = hist->GetNbinsX();
double minX = hist->GetXaxis()->GetXmax();
double maxX = hist->GetXaxis()->GetXmin();
double minY = hist->GetYaxis()->GetXmax();
double maxY = hist->GetYaxis()->GetXmin();
for (int i = 1; i <= nBins; ++i) {
double y = hist->GetBinContent(i);
if (y != 0) {
double x = hist->GetXaxis()->GetBinCenter(i);
if (x < minX)
minX = x;
if (x > maxX)
maxX = x;
if (y < minY)
minY = y;
if (y > maxY)
maxY = y;
}
}
hist->GetXaxis()->SetRangeUser(minX * 0.9, maxX * 1.1);
hist->GetYaxis()->SetRangeUser(minY * 0.9, maxY * 1.1);
}
void BluetAnalyzer::readCutParameters(string cfgfile) {
pugi::xml_document document;
document.load_file(cfgfile.c_str());
pugi::xml_node root = document.child("Bluet");
pugi::xml_node detector = root.child("DetectorParameters");
sidelength = bluet::XMLToQuantity(detector.child("PadSideLength")); // mm
Nlayer = bluet::XMLToQuantity(detector.child("PadArrayNlayer"));
driftL = bluet::XMLToQuantity(detector.child("DriftLength"));
pugi::xml_node electronics = root.child("ElectronicsParameters");
sampleT = bluet::XMLToQuantity(electronics.child("WaveformSamplingPeriod"));
pugi::xml_node gas = root.child("GasParameters");
dv0 = bluet::XMLToQuantity(gas.child("DriftVelocity")) * 1e4;
pugi::xml_node cut = root.child("CutParameters");
EcutMin = bluet::XMLToQuantity(cut.child("ECutMin"));
EcutMax = bluet::XMLToQuantity(cut.child("ECutMax"));
tofCutMin = bluet::XMLToQuantity(cut.child("TOFCutMin"));
tofCutMax = bluet::XMLToQuantity(cut.child("TOFCutMax"));
stampCutMin = bluet::XMLToQuantity(cut.child("TimeStampCutMin"));
stampCutMax = bluet::XMLToQuantity(cut.child("TimeStampCutMax"));
cathAmpCutMin = bluet::XMLToQuantity(cut.child("CathodeAmpCutMin"));
cathAmpCutMax = bluet::XMLToQuantity(cut.child("CathodeAmpCutMax"));
cathTimeCutMin = bluet::XMLToQuantity(cut.child("CathodeTimeCutMin"));
cathTimeCutMax = bluet::XMLToQuantity(cut.child("CathodeTimeCutMax"));
meshAmpCutMin = bluet::XMLToQuantity(cut.child("MeshAmpCutMin"));
meshAmpCutMax = bluet::XMLToQuantity(cut.child("MeshAmpCutMax"));
meshTimeCutMin = bluet::XMLToQuantity(cut.child("MeshTimeCutMin"));
meshTimeCutMax = bluet::XMLToQuantity(cut.child("MeshTimeCutMax"));
NtrkCutMin = bluet::XMLToQuantity(cut.child("NtrackCutMin"));
NtrkCutMax = bluet::XMLToQuantity(cut.child("NtrackCutMax"));
lengthCutMin = bluet::XMLToQuantity(cut.child("LengthCutMin"));
lengthCutMax = bluet::XMLToQuantity(cut.child("LengthCutMax"));
cosCutMin = bluet::XMLToQuantity(cut.child("cosCutMin"));
CenterPos[0] = -3.03631;
CenterPos[1] = 0.634384;
sliceW = 2.5; // [mm]
rCut = bluet::XMLToQuantity(cut.child("radiusCutEnable"));
rCutMin = bluet::XMLToQuantity(cut.child("radiusCutMin"));
rCutMax = bluet::XMLToQuantity(cut.child("radiusCutMax"));
Nbin = 100;
tmin = 2.e2; // [ns] gamma
tmax = 5.e6; // [ns] 1 eV
Nrbin = 300;
rmax = 80.;
Sbin = new double[Nrbin];
double rstep = rmax / Nrbin;
for (int ibin = 0; ibin < Nrbin; ++ibin) {
Sbin[ibin] = ((ibin + 1) * rstep) * ((ibin + 1) * rstep) -
(ibin * rstep) * (ibin * rstep);
}
}
void BluetAnalyzer::defineFitFunctions() {
fbeam1 = new TF1("fbeam1",
"[0]+0.5*[1]*(TMath::Erf((x+[2])/(TMath::Sqrt(2.)*[3])) - "
"TMath::Erf((x-[2])/(TMath::Sqrt(2.)*[3])) )",
-80., 80.);
fbeam1->SetParameters(1., 250., 15., 5.);
fbeam1->SetLineWidth(2.);
fbeam2 = new TF2(
"fbeam2",
"[0]+0.5*[1]*("
"TMath::Erf((TMath::Sqrt(TMath::Power((x-[4]),2)+TMath::Power((y-[5])"
",2))+[2])/(TMath::Sqrt(2.)*[3])) - "
"TMath::Erf((TMath::Sqrt(TMath::Power((x-[4]),2)+TMath::Power((y-[5])"
",2))-[2])/(TMath::Sqrt(2.)*[3])) )",
-80., 80., -80., 80.);
fbeam2->SetParameters(1., 60., 15., 2., -4., 1.);
fbeam2->SetParLimits(0, 0., 5.);
fbeam2->SetParLimits(1, 10., 500.);
fbeam2->SetParLimits(2, 10., 20.);
fbeam2->SetParLimits(3, 1., 10.);
fbeam2->SetParLimits(4, -5., 5.);
fbeam2->SetParLimits(5, -5., 5.);
fbeam2->SetLineWidth(short(1.5));
fbeam2->SetLineColor(kGreen);
fbeamr = new TF1("fbeamr", "[0]+[1]/(1.+TMath::Exp((x-[2])/[3] ) )", 0., 80.);
fbeamr->SetParameters(1., 50., 15., 1.);
fbeamr->SetParLimits(3, 1., 5.);
fbeamr->SetLineWidth(2.);
fbeamr->SetLineColor(kYellow);
fdgaus = new TF1("fdgaus", dgaus, 0., 50., 5);
fdgaus->SetParNames("#mu", "amp_1", "#sigma_1", "amp_2", "#sigma_2");
fdgaus->SetParLimits(0, 0., 1.e6);
fdgaus->SetParLimits(1, 0., 1.e8);
fdgaus->SetParLimits(2, 0., 1.e6);
fdgaus->SetParLimits(3, 0., 1.e8);
fdgaus->SetParLimits(4, 0., 1.e6);
}
void BluetAnalyzer::defineHistograms() {
hprofile = new TH2F("hprofile", "hprofile;xpos [mm];ypos [mm];Count", 80,
-80., 80., 80, -80., 80.);
hprofile->GetXaxis()->CenterTitle();
hprofile->GetYaxis()->CenterTitle();
hprofile->GetZaxis()->CenterTitle("Count");
hprofile->GetZaxis()->SetTitleOffset(1.4);
hprofile->GetZaxis()->SetMaxDigits(1);
hslicex = new TH1F("hslicex", "hslicex;pos [mm];Count", 50, -80., 80.);
hslicex->GetXaxis()->CenterTitle();
hslicex->GetYaxis()->CenterTitle();
hslicex->SetMarkerColor(kRed);
hslicex->SetMarkerStyle(20);
hslicex->GetYaxis()->SetMaxDigits(2);
hslicey = new TH1F("hslicey", "hslicey;pos [mm];Count", 50, -80., 80.);
hslicey->GetXaxis()->CenterTitle();
hslicey->GetYaxis()->CenterTitle();
hslicey->SetMarkerColor(kBlue);
hslicey->SetMarkerStyle(20);
hslicey->GetYaxis()->SetMaxDigits(2);
hbeamr = new TH1F("hbeamr", "hbeamr;Radius [mm];Count", Nrbin, 0., rmax);
hbeamr->GetYaxis()->SetMaxDigits(2);
hbeamr->GetXaxis()->CenterTitle();
hbeamr->GetYaxis()->CenterTitle();
hbeamr->SetLineWidth(2);
hrinte = new TH1F("hrinte", "hrinte;Radius [mm];Accumulate", Nrbin, 0., rmax);
hrinte->GetYaxis()->SetMaxDigits(2);
hrinte->SetLineColor(kRed);
hrinte->SetLineWidth(2);
const double tstep = (Log10(tmax) - Log10(tmin)) / Nbin;
double tpos[Nbin + 1];
for (int it = 0; it < (Nbin + 1); ++it)
tpos[it] = tmin * TMath::Power(10., it * tstep);
htofE = new TH2F("htofE", "htofE;TOF [ns];Energy [ADC]", Nbin, tpos, 250, 0.,
5000.);
htofE->GetXaxis()->CenterTitle();
htofE->GetYaxis()->CenterTitle();
htofE->GetYaxis()->SetMaxDigits(2);
htofr0 = new TH1F("htofr0", "htofr0;TOF [ns];Count", Nbin, tpos);
htofr0->GetXaxis()->CenterTitle();
htofr0->GetYaxis()->CenterTitle();
htofr0->GetYaxis()->SetMaxDigits(2);
htofr0->SetLineColor(kRed);
htofr0->SetLineWidth(2.);
htofr1 = new TH1F("htofr1", "htofr1;TOF [ns];Count", Nbin, tpos);
htofr1->GetXaxis()->CenterTitle();
htofr1->GetYaxis()->CenterTitle();
htofr1->GetYaxis()->SetMaxDigits(2);
htofr1->SetLineColor(kBlue);
htofr1->SetLineWidth(2.);
hpull = new TH1F("hpull", "hpull;TOF [ns];Pull Value", Nbin, tpos);
hpull->GetXaxis()->CenterTitle();
hpull->GetYaxis()->CenterTitle();
hpull->GetYaxis()->SetMaxDigits(2);
hpull->SetMarkerStyle(21);
hpull->SetMarkerColor(9);
hpull->SetLineColor(9);
hratio = new TH1F("hratio", "hratio;TOF [ns];Ratio", Nbin, tpos);
hratio->GetXaxis()->CenterTitle();
hratio->GetYaxis()->CenterTitle();
hratio->SetMarkerStyle(21);
hratio->SetMarkerColor(9);
hratio->SetLineColor(9);
hLengthE = new TH2F("hLengthE", "hLengthE;Length [mm];Energy [ADC]", 200, 0.,
200., 250, 0., 5000.);
hLengthE->GetXaxis()->CenterTitle();
hLengthE->GetYaxis()->CenterTitle();
hLengthE->GetYaxis()->SetMaxDigits(2);
htcath = new TH1F("htcath", "htcath;T_{max}-T_{cathode} [ns];Count", 100, 0.,
10000.);
htcath->GetXaxis()->CenterTitle();
htcath->GetYaxis()->CenterTitle();
htcath->SetMarkerColor(kBlack);
htcath->SetMarkerStyle(20);
htpad =
new TH1F("htpad", "htpad;T_{max}-T_{min} [ns];Count", 100, 0., 10000.);
htpad->GetXaxis()->CenterTitle();
htpad->GetYaxis()->CenterTitle();
htpad->SetMarkerColor(kBlack);
htpad->SetMarkerStyle(20);
hv0 = new TH1F("hv0", "hv0;Drift Veloctiy [mm/#mus];Count", 300, 0., 300.);
hv0->GetXaxis()->CenterTitle();
hv0->GetYaxis()->CenterTitle();
hv0->SetMarkerColor(kBlack);
hv0->SetMarkerStyle(20);
hv1 = new TH1F("hv1", "hv1;Drift Veloctiy [mm/#mus];Count", 300, 0., 300.);
hv1->GetXaxis()->CenterTitle();
hv1->GetYaxis()->CenterTitle();
hv1->SetMarkerColor(kBlack);
hv1->SetMarkerStyle(20);
hcos = new TH1F("hcos", "hcos;cos#theta;Count", 100, 0., 1.);
hcos->GetXaxis()->CenterTitle();
hcos->GetYaxis()->CenterTitle();
hcos->SetMarkerColor(kBlack);
hcos->SetLineWidth(2.);
hcos->SetMarkerStyle(20);
hphi = new TH1F("hphi", "hphi;#phi [deg];Count", 72, -180., 180.);
hphi->GetXaxis()->CenterTitle();
hphi->GetYaxis()->CenterTitle();
hphi->SetMarkerColor(kBlack);
hphi->SetLineWidth(2.);
hphi->SetMarkerStyle(20);
hdeltacos1 =
new TH2F("hdeltacos1", "hdeltacos1;cos#theta;T_{max}-T_{cathode} [ns]",
100, 0., 1., 100, 0., 10000.);
hdeltacos1->GetXaxis()->CenterTitle();
hdeltacos1->GetYaxis()->CenterTitle();
hdeltacos1->GetYaxis()->SetMaxDigits(2);
hdeltacos1->GetYaxis()->SetTitleOffset(1.5);
hdeltacos2 =
new TH2F("hdeltacos2", "hdeltacos2;cos#theta;T_{max}-T_{min} [ns]", 100,
0., 1., 110, -1000., 10000.);
hdeltacos2->GetXaxis()->CenterTitle();
hdeltacos2->GetYaxis()->CenterTitle();
hdeltacos2->GetYaxis()->SetMaxDigits(2);
hdeltacos2->GetYaxis()->SetTitleOffset(1.5);
hts =
new TH2F("hts", "hts;TimeStamp;TOF0 [ns]", 1000, 0., 1.e6, 2000, 0, 2.e6);
hts->GetXaxis()->CenterTitle();
hts->GetYaxis()->CenterTitle();
hts->GetYaxis()->SetTitleOffset(1.5);
htn = new TH2F("htn", "htn;TimeStamp;Nhit", 200, 1.5e3, 5.e3, 300, 0, 300);
htn->GetXaxis()->CenterTitle();
htn->GetYaxis()->CenterTitle();
htn->GetYaxis()->SetTitleOffset(1.5);
hamp = new TH2F("hamp", "hamp;MeshAmp [ADC];CathodeAmp [ADC]", 200, 0, 4e3,
200, 0, 4e3);
hamp->GetXaxis()->CenterTitle();
hamp->GetYaxis()->CenterTitle();
hamp->GetYaxis()->SetTitleOffset(1.5);
han = new TH2F("han", "han;Nhit;CathodeAmp [ADC]", 300, 0, 300, 200, 0, 4e3);
han->GetXaxis()->CenterTitle();
han->GetYaxis()->CenterTitle();
han->GetYaxis()->SetTitleOffset(1.5);
hta =
new TH2F("hta", "hta;#tau [us];Amplitude [ADC]", 100, 0, 5, 300, 0,
3e3);
hta->GetXaxis()->CenterTitle();
hta->GetYaxis()->CenterTitle();
hta->GetYaxis()->SetTitleOffset(1.5);
htt = new TH2F("htt", "htt;T_0 [us];#tau [us]", 100, 0, 10, 100, 0, 5);
htt->GetXaxis()->CenterTitle();
htt->GetYaxis()->CenterTitle();
htt->GetYaxis()->SetTitleOffset(1.5);
}
bool BluetAnalyzer::isGoodEvent(BluetEvent *event) {
if (event->TOF0 < tofCutMin || event->TOF0 > tofCutMax)
return false;
if (event->TimeStamp < stampCutMin || event->TimeStamp > stampCutMax)
return false;
if (event->CathodeAmp < cathAmpCutMin || event->CathodeAmp > cathAmpCutMax)
return false;
if (event->CathodeTime < cathTimeCutMin ||
event->CathodeTime > cathTimeCutMax)
return false;
if (event->MeshAmp < meshAmpCutMin || event->MeshAmp > meshAmpCutMax)
return false;
if (event->MeshTime < meshTimeCutMin || event->MeshTime > meshTimeCutMax)
return false;
if (event->Tracks.size() < NtrkCutMin || event->Tracks.size() > NtrkCutMax)
return false;
return true;
}
bool BluetAnalyzer::isGoodTrack(BluetTrack *track) {
if (track->TrSumEdep < EcutMin || track->TrSumEdep > EcutMax)
return false;
if (track->TrRawLength < lengthCutMin || track->TrRawLength > lengthCutMax)
return false;
if (track->TrCosTheta < cosCutMin)
return false;
double *xypos = track->ZmaxHit_xy;
double radius =
TMath::Sqrt((xypos[0] - CenterPos[0]) * (xypos[0] - CenterPos[0]) +
(xypos[1] - CenterPos[1]) * (xypos[1] - CenterPos[1]));
if (rCut && (radius < rCutMin || radius > rCutMax))
return false;
return true;
}
void BluetAnalyzer::readTreeData(vector<string> files) {
for (int i = 0; i < files.size(); ++i) {
string infile = files[i];
eventTree->Add(infile.c_str());
}
eventTree->SetBranchAddress("BluetEvent", &fEvent);
Nevent = eventTree->GetEntries();
printf("===> Total event = %d\n", Nevent);
}
void BluetAnalyzer::FillHistograms() {
printf("===> Start reading data ...\n");
pgbar::pgbar<> bar{Nevent};
bar.set_step(500);
for (Long_t ievent = 0; ievent < Nevent; ievent++) {
if (ievent % 500 == 0)
bar.update();
eventTree->GetEntry(ievent);
if (!isGoodEvent(fEvent))
continue;
for (size_t ihit = 0; ihit < fEvent->Hits.size(); ++ihit) {
BluetHit *hit = &(fEvent->Hits[ihit]);
hta->Fill(hit->HitTau, hit->HitAmp);
htt->Fill(hit->HitT * sampleT / 1e3, hit->HitTau);
}
for (size_t itrk = 0; itrk < fEvent->Tracks.size(); ++itrk) {
if (!isGoodTrack(&(fEvent->Tracks[itrk])))
continue;
double ts = fEvent->TimeStamp;
double tof = fEvent->TOF0;
double sumE = fEvent->Tracks[itrk].TrSumEdep;
double length = fEvent->Tracks[itrk].TrRawLength;
double tmax = fEvent->Tracks[itrk].Tmax;
double tmin = fEvent->Tracks[itrk].Tmin;
double tcath = fEvent->CathodeTime;
double *xypos = fEvent->Tracks[itrk].ZmaxHit_xy;
double radius =
TMath::Sqrt((xypos[0] - CenterPos[0]) * (xypos[0] - CenterPos[0]) +
(xypos[1] - CenterPos[1]) * (xypos[1] - CenterPos[1]));
if (radius <= rBeam)
nInner++;
if (radius >= rTarget && rSubstrate <= rSubstrate + 5)
nRing++;
double v0 = driftL / (tmax - tcath) / sampleT;
double v1 = driftL / (tmax - tmin) / sampleT;
hv0->Fill(v0 * 1000.); // [mm/ns]->[mm/us]
hv1->Fill(v1 * 1000.); // [mm/ns]->[mm/us]
htpad->Fill((tmax - tmin) * sampleT);
htcath->Fill((tmax - tcath) * sampleT);
hcos->Fill(fabs(fEvent->Tracks[itrk].TrCosTheta));
hphi->Fill(fEvent->Tracks[itrk].TrPhi);
hdeltacos1->Fill(fabs(fEvent->Tracks[itrk].TrCosTheta),
(tmax - tcath) * sampleT);
hdeltacos2->Fill(fabs(fEvent->Tracks[itrk].TrCosTheta),
(tmax - tmin) * sampleT);
hLengthE->Fill(length, sumE);
hprofile->Fill(xypos[0], xypos[1]);
hbeamr->Fill(radius);
hts->Fill(ts, tof);
htn->Fill(ts, fEvent->Tracks[itrk].NTrHits);
hamp->Fill(fEvent->MeshAmp, fEvent->CathodeAmp);
han->Fill(fEvent->Tracks[itrk].NTrHits, fEvent->CathodeAmp);
if (fabs(xypos[1] - CenterPos[1]) < sliceW)
hslicex->Fill(xypos[0] - CenterPos[0]);
if (fabs(xypos[0] - CenterPos[0]) < sliceW)
hslicey->Fill(xypos[1] - CenterPos[1]);
if (radius < rCutMin)
htofr0->Fill(tof);
if ((radius > rCutMin) && (radius < rCutMax))
htofr1->Fill(tof);
if (radius < rCutMax) {
htofE->Fill(tof, sumE);
}
}
}
printf("\n===> Finish reading data ...\n");
}
void BluetAnalyzer::postOperations() {
printf("===> After Filtering = %f\n", hv0->Integral());
double *aCount = hbeamr->GetIntegral();
for (int ibin = 0; ibin < Nrbin; ++ibin)
hrinte->SetBinContent(ibin + 1, aCount[ibin]);
for (int ibin = 0; ibin < Nrbin; ++ibin) {
double count = hbeamr->GetBinContent(ibin + 1);
count = count / Sbin[ibin];
hbeamr->SetBinContent(ibin + 1, count);
}
hrinte->Scale(hbeamr->GetBinContent(hbeamr->GetMaximumBin()) /
hrinte->GetBinContent(hrinte->GetMaximumBin()));
// htof, htofr0, htofr1
htof = (TH1F *)(htofE->ProjectionX());
htof->SetLineWidth(2.);
htof->SetLineColor(kBlack);
htofr0->Scale(htof->Integral() / htofr0->Integral());
htofr1->Scale(htof->Integral() / htofr1->Integral());
// hpull, hratio
for (int ibin = 0; ibin < Nbin; ++ibin) {
double n0 = htofr0->GetBinContent(ibin + 1);
double n1 = htofr1->GetBinContent(ibin + 1);
double n = htof->GetBinContent(ibin + 1);
double e0 = htofr0->GetBinError(ibin + 1);
double e1 = htofr1->GetBinError(ibin + 1);
// double e = htof->GetBinError(ibin + 1);
if ((n0 > 0) && (n1 > 0) && (n > 0)) {
double pull = (n0 - n1) / TMath::Sqrt(e0 * e0 + e1 * e1);
double ratio = n0 / n1;
double rErr = TMath::Sqrt(e0 * e0 / (n1 * n1) +
e1 * e1 * n0 * n0 / (n1 * n1 * n1 * n1));
hpull->SetBinContent(ibin + 1, pull);
hratio->SetBinContent(ibin + 1, ratio);
hratio->SetBinError(ibin + 1, rErr);
}
}
}
void BluetAnalyzer::drawVertex() {
TGCompositeFrame *tab = mTab->AddTab("Vertex");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva Vertex", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("v2d", "Vertex 2D", 0.02, 0.02, 0.48, 0.48);
pad->SetTopMargin(0.15);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
gPad->SetGridx();
gPad->SetGridy();
hprofile->Draw("colz");
hprofile->Fit("fbeam2", "QN");
TEllipse *cPCB = new TEllipse(CenterPos[0], CenterPos[1], rPCB);
cPCB->SetLineColor(kMagenta);
cPCB->SetLineWidth(2);
cPCB->SetFillColorAlpha(kWhite, 0);
cPCB->Draw();
TEllipse *cSubstrate = new TEllipse(CenterPos[0], CenterPos[1], rSubstrate);
cSubstrate->SetLineColor(kBlue);
cSubstrate->SetLineWidth(2);
cSubstrate->SetFillColorAlpha(kWhite, 0);
cSubstrate->Draw();
TEllipse *cTarget = new TEllipse(CenterPos[0], CenterPos[1], rTarget);
cTarget->SetLineColor(kGreen);
cTarget->SetLineWidth(2);
cTarget->SetFillColorAlpha(kWhite, 0);
cTarget->Draw();
TEllipse *cBeam = new TEllipse(CenterPos[0], CenterPos[1], rBeam);
cBeam->SetLineColor(kRed);
cBeam->SetLineWidth(2);
cBeam->SetFillColorAlpha(kWhite, 0);
cBeam->Draw();
TLatex *ratio =
new TLatex(-75, 70, Form("Ratio: %.2f", nInner * 1.0 / nRing));
ratio->SetTextSize(0.05);
ratio->Draw();
c0->cd();
pad = new TPad("v3d", "Vertex 3D", 0.52, 0.02, 0.98, 0.48);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
hprofile->SetFillColor(kYellow);
hprofile->Draw("lego2z 0");
fbeam2->Draw("same surf");
c0->cd();
pad = new TPad("xyd", "X/Y Distribution", 0.02, 0.52, 0.48, 0.98);
pad->Draw();
pad->cd();
hslicex->Draw("E1");
hslicex->Fit("fbeam1", "QNL");
fbeam1->SetLineColor(kRed);
((TF1 *)fbeam1->Clone())->Draw("same");
double top = fbeam1->GetParameter(1);
double r10 = fbeam1->GetX(0.1 * top, 0., 80.);
printf("===> x-direction R10 = %f\n", r10);
hslicey->Draw("same E1");
hslicey->Fit("fbeam1", "QNL");
fbeam1->SetLineColor(kBlue);
((TF1 *)fbeam1->Clone())->Draw("same");
top = fbeam1->GetParameter(1);
r10 = fbeam1->GetX(0.1 * top, 0., 80.);
printf("===> y-direction R10 = %f\n", r10);
TLegend *leg = new TLegend(0.2, 0.8, 0.35, 0.9);
leg->AddEntry(hslicex, "x-slice");
leg->AddEntry(hslicey, "y-slice");
leg->SetTextSize(0.03);
leg->Draw();
c0->cd();
pad = new TPad("rd", "R Distribution", 0.52, 0.52, 0.98, 0.98);
pad->Draw();
pad->cd();
gPad->SetGridx();
gPad->SetGridy();
hbeamr->Draw();
hrinte->Draw("same HIST");
hbeamr->Fit("fbeamr", "Q");
TLegend *leg1 = new TLegend(0.68, 0.65, 0.93, 0.75);
leg1->AddEntry(hbeamr, "r-distribution");
leg1->AddEntry(hrinte, "accumulative");
leg1->SetTextSize(0.03);
leg1->Draw();
c0->cd();
c0->Update();
}
void BluetAnalyzer::drawTOF() {
TGCompositeFrame *tab = mTab->AddTab("TOF");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva TOF", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("tof", "", 0.02, 0.02, 0.48, 0.48);
pad->Draw();
pad->cd();
gPad->SetLogx();
gPad->SetLogy();
htofr0->Draw("HIST");
htofr1->Draw("same HIST");
htof->Draw("same HIST");
TLegend *leg = new TLegend(0.45, 0.75, 0.8, 0.9);
leg->AddEntry(htofr0, Form("r < %.1f mm", rCutMin));
leg->AddEntry(htofr1, Form("%.1f < r < %.1f mm", rCutMin, rCutMax));
leg->AddEntry(htof, Form("r < %.1f mm", rCutMax));
leg->SetTextSize(0.03);
leg->Draw();
c0->cd();
pad = new TPad("pull", "", 0.52, 0.02, 0.98, 0.48);
pad->Draw();
pad->cd();
hpull->GetYaxis()->SetRangeUser(-8., 8.);
hpull->Draw("P");
TLine *line = new TLine();
line->DrawLine(tmin, 0., tmax, 0.);
c0->cd();
pad = new TPad("ratio", "", 0.02, 0.52, 0.48, 0.98);
pad->Draw();
pad->cd();
hratio->Draw("E1");
line->DrawLine(tmin, 1., tmax, 1.);
c0->cd();
pad = new TPad("tof.ts", "", 0.52, 0.52, 0.98, 0.98);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(hts);
hts->Draw("colz");
c0->cd();
c0->Update();
}
void BluetAnalyzer::drawDrift() {
double *par = nullptr;
TGCompositeFrame *tab = mTab->AddTab("Drift Time & Velocity");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva dtv", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("tcath", "", 0.02, 0.02, 0.48, 0.48);
pad->Draw();
pad->cd();
adjustAxisRange(htcath);
htcath->Draw("E1");
fitdgaus(htcath, fdgaus, par, 70 / dv0 * 1000);
fdgaus->SetParameters(par[0], par[1], par[2], par[3], par[4]);
htcath->Fit("fdgaus", "LQ", "", par[5], par[6]);
htcath->GetFunction("fdgaus")->SetLineColor(kRed);
c0->cd();
pad = new TPad("dv0", "", 0.52, 0.02, 0.98, 0.48);
pad->Draw();
pad->cd();
adjustAxisRange(hv0);
hv0->Draw("E1");
fitdgaus(hv0, fdgaus, par, dv0);
fdgaus->SetParameters(par[0], par[1], par[2], par[3], par[4]);
hv0->Fit("fdgaus", "LQ", "", par[5], par[6]);
hv0->GetFunction("fdgaus")->SetLineColor(kRed);
c0->cd();
pad = new TPad("tpad", "", 0.02, 0.52, 0.48, 0.98);
pad->Draw();
pad->cd();
adjustAxisRange(htpad);
htpad->Draw("E1");
fitdgaus(htpad, fdgaus, par, 70 / dv0 * 1000);
fdgaus->SetParameters(par[0], par[1], par[2], par[3], par[4]);
htpad->Fit("fdgaus", "LQ", "", par[5], par[6]);
htpad->GetFunction("fdgaus")->SetLineColor(kRed);
c0->cd();
pad = new TPad("dv1", "", 0.52, 0.52, 0.98, 0.98);
pad->Draw();
pad->cd();
adjustAxisRange(hv1);
fitdgaus(hv1, fdgaus, par, dv0);
fdgaus->SetParameters(par[0], par[1], par[2], par[3], par[4]);
hv1->Draw("E1");
hv1->Fit("fdgaus", "LQ", "", par[5], par[6]);
hv1->GetFunction("fdgaus")->SetLineColor(kRed);
c0->cd();
c0->Update();
}
void BluetAnalyzer::drawAngle() {
TGCompositeFrame *tab = mTab->AddTab("Angle");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva angle", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("cos", "", 0.02, 0.02, 0.48, 0.48);
pad->Draw();
pad->cd();
hcos->GetYaxis()->SetRangeUser(0., hcos->GetMaximum() * 1.1);
hcos->Draw("HIST");
c0->cd();
pad = new TPad("phi", "", 0.52, 0.02, 0.98, 0.48);
pad->Draw();
pad->cd();
hphi->GetYaxis()->SetRangeUser(0., hphi->GetMaximum() * 1.1);
hphi->Draw("HIST");
c0->cd();
pad = new TPad("tcath.cos", "", 0.02, 0.52, 0.48, 0.98);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
gPad->SetGridy();
hdeltacos1->Draw("colz");
c0->cd();
pad = new TPad("tpad.cos", "", 0.52, 0.52, 0.98, 0.98);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
gPad->SetGridy();
hdeltacos2->Draw("colz");
c0->cd();
c0->Update();
}
void BluetAnalyzer::draw1DFeatures() {
TGCompositeFrame *tab = mTab->AddTab("1D Features");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva 1d", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
TH1F *htmp;
c0->cd();
pad = new TPad("meshamp", "", 0.02, 0.02, 0.48, 0.48);
pad->Draw();
pad->cd();
htmp = (TH1F *)(hamp->ProjectionX());
adjustAxisRange(htmp);
htmp->Draw("HIST");
c0->cd();
pad = new TPad("cathamp", "", 0.52, 0.02, 0.98, 0.48);
pad->Draw();
pad->cd();
htmp = (TH1F *)(hamp->ProjectionY());
adjustAxisRange(htmp);
htmp->Draw("HIST");
c0->cd();
pad = new TPad("length", "", 0.02, 0.52, 0.48, 0.98);
pad->Draw();
pad->cd();
htmp = (TH1F *)(hLengthE->ProjectionX());
adjustAxisRange(htmp);
htmp->Draw("HIST");
c0->cd();
pad = new TPad("energy", "", 0.52, 0.52, 0.98, 0.98);
pad->Draw();
pad->cd();
htmp = (TH1F *)(hLengthE->ProjectionY());
adjustAxisRange(htmp);
htmp->Draw("HIST");
c0->cd();
c0->Update();
}
void BluetAnalyzer::draw2DFeatures() {
TGCompositeFrame *tab = mTab->AddTab("2D Features");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva 2d", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("le", "", 0.02, 0.02, 0.48, 0.48);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(hLengthE);
hLengthE->Draw("colz");
c0->cd();
pad = new TPad("tn", "", 0.52, 0.02, 0.98, 0.48);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
gPad->SetLogx();
adjustAxisRange(htn);
htn->Draw("colz");
c0->cd();
pad = new TPad("amp", "", 0.02, 0.52, 0.48, 0.98);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(hamp);
hamp->Draw("colz");
c0->cd();
pad = new TPad("an", "", 0.52, 0.52, 0.98, 0.98);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(han);
han->Draw("colz");
c0->cd();
c0->Update();
}
void BluetAnalyzer::drawPadFeatures() {
TGCompositeFrame *tab = mTab->AddTab("Pad Features");
TRootEmbeddedCanvas *ec =
new TRootEmbeddedCanvas("Canva Pad", tab, sizeX, sizeY);
tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
TCanvas *c0 = ec->GetCanvas();
TPad *pad;
c0->cd();
pad = new TPad("le", "", 0.02, 0.02, 0.48, 0.48);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(hta);
hta->Draw("colz");
c0->cd();
pad = new TPad("tn", "", 0.52, 0.02, 0.98, 0.48);
pad->SetRightMargin(0.15);
pad->Draw();
pad->cd();
adjustAxisRange(htt);
htt->Draw("colz");
c0->cd();
pad = new TPad("amp", "", 0.02, 0.52, 0.48, 0.98);
pad->Draw();
pad->cd();
c0->cd();
pad = new TPad("an", "", 0.52, 0.52, 0.98, 0.98);
pad->Draw();
pad->cd();
c0->cd();
c0->Update();
}
void BluetAnalyzer::drawHistograms() {
gROOT->SetStyle("ATLAS");
gStyle->SetOptFit(1);
gStyle->SetStatFontSize(0.03);
gStyle->SetCanvasPreferGL();
drawVertex();
drawTOF();
drawDrift();
drawAngle();
draw1DFeatures();
draw2DFeatures();
drawPadFeatures();
mFrame->SetWindowName("Bluet Analyzer");
mFrame->MapSubwindows();
mFrame->Resize(mFrame->GetDefaultSize());
mFrame->MapWindow();
app->Run(kTRUE);
return;
}

8
temp_env.xml Normal file
View File

@ -0,0 +1,8 @@
<?xml version="1.0"?>
<root>
<bluet_config value="/home/fox/bluet/BluetConfig/BluetConfigFile_np.xml" />
<bluet_rawdata value="/home/fox/bluet/myBluetData/RawData/" />
<bluet_eventdata value="/home/fox/bluet/myBluetData/EventData/" />
<bluet_rootdata value="/home/fox/bluet/myBluetData/RootData/" />
<bluet_simdata value="/home/fox/bluet/myBluetData/SimData/" />
</root>