Commit 5e56589a authored by andreasgad's avatar andreasgad
Browse files

Initial

parents
cmake_minimum_required(VERSION 2.8)
project(I257)
# Tell cmake where to look for cmake modules.
set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${CMAKE_SOURCE_DIR}/cmake")
#Lets write modern C++
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp -std=c++11")
#----------------------------------------------------------------------------
# Dependencies
#
# Variable for holding dependant libraries
set(LIBRARIES "")
# Include AUSAlib. Look in cmake/FindAUSALIB.cmake for how the specify location
find_package(AUSALIB REQUIRED)
include_directories(${AUSALIB_INCLUDES})
list(APPEND LIBRARIES ${AUSALIB_LIBRARIES})
# Include ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
link_directories(${ROOT_LIBRARY_DIR})
list(APPEND LIBRARIES ${ROOT_LIBRARIES} Spectrum MathMore)
find_package(CURL)
if (CURL_FOUND)
include_directories(${CURL_INCLUDE_DIRS})
list(APPEND LIBRARIES ${CURL_LIBRARIES})
add_definitions("-DAUSA_HAS_LIBCURL")
endif(CURL_FOUND)
link_libraries(${LIBRARIES})
include_directories(external/jsoncons/include)
add_executable(Singles SinglesAnalysis.cpp)
#include <TVector3.h>
#include <TLorentzVector.h>
#include "ParticleType.h"
struct Hit {
double deposited;
double paddeposited;
double E, FE, BE, Edssd, dE, EBeta, fbdiff;
double TF, TB, TPad, T;
double tarTrav, detTrav;
double Ectarget;
TVector3 direction, position, origin;
double theta;
double angle;
double phi;
int numInDet, detectorMul;
short fseg, bseg;
std::vector<ParticleType> type;
bool canBeAlpha;
bool canBeBeta;
size_t index;
TLorentzVector lVector, lVectorBeta, lVectorAlpha;
};
#ifndef PARTICLETYPE_H
#define PARTICLETYPE_H
#include <ausa/eloss/Ion.h>
class ParticleType {
public:
explicit ParticleType(const std::string &name);
static ParticleType constructUnknown() {
auto ret = ParticleType{};
ret.unknown = true;
return ret;
}
double mass{};
std::string name{};
int A{}, Z{};
bool unknown{};
ParticleType() = default;
};
ParticleType::ParticleType(const std::string &name) {
AUSA::EnergyLoss::Ion ion{name};
mass = ion.getMass();
this->name = ion.getName();
Z = ion.getZ();
A = ion.getA();
unknown = false;
}
#endif //PARTICLETYPE_H
#include <iostream>
#include <string>
#include <ausa/json/IO.h>
#include <ausa/geometry/Box.h>
#include <ausa/sort/analyzer/AbstractSortedAnalyzer.h>
#include <ausa/sort/SortedReader.h>
#include <ausa/util/StringUtil.h>
#include <ausa/util/FileUtil.h>
#include <ausa/setup/DoubleSidedSiliconDetector.h>
#include <ausa/setup/SingleSidedSiliconDetector.h>
#include <ausa/util/DynamicBranchVector.h>
#include <ausa/eloss/Ion.h>
#include <ausa/eloss/Default.h>
#include <ausa/util/memory>
#include <ausa/constants/Mass.h>
#include <ausa/output/OutputConvenience.h>
#include <Math/Vector3D.h>
#include "Hit.h"
#include <TROOT.h>
#include <ctime>
using namespace std;
using namespace AUSA;
using namespace ROOT::Math;
using namespace AUSA::Sort;
using namespace AUSA::EnergyLoss;
using namespace AUSA::Constants;
class SinglesAnalysis : public AbstractSortedAnalyzer {
public:
SinglesAnalysis(Target &target, TFile *output) : target(target) {
NUM = 0;
ALPHA = new ParticleType{"He4"};
t = new TTree("a", "a");
t->Branch("mul", &mul);
t->Branch("clock", &CLOCK);
t->Branch("num", &NUM);
v_dir = make_unique<DynamicBranchVector<TVector3>>(*t, "dir");
v_pos = make_unique<DynamicBranchVector<TVector3>>(*t, "pos");
v_theta = make_unique<DynamicBranchVector<double>>(*t, "theta", "mul");
v_ang = make_unique<DynamicBranchVector<double>>(*t, "angle", "mul");
v_E = make_unique<DynamicBranchVector<double>>(*t, "E", "mul");
v_BE = make_unique<DynamicBranchVector<double>>(*t, "BE", "mul");
v_FE = make_unique<DynamicBranchVector<double>>(*t, "FE", "mul");
v_FT = make_unique<DynamicBranchVector<double>>(*t, "FT", "mul");
v_BT = make_unique<DynamicBranchVector<double>>(*t, "BT", "mul");
v_dE = make_unique<DynamicBranchVector<double>>(*t, "dE", "mul");
v_i = make_unique<DynamicBranchVector<short>>(*t, "id", "mul");
v_F = make_unique<DynamicBranchVector<short>>(*t, "FI", "mul");
v_B = make_unique<DynamicBranchVector<short>>(*t, "BI", "mul");
SiCalc = defaultRangeInverter("He4", "Silicon");
for (auto &layer : target.getLayers()) {
targetCalcs.push_back(defaultRangeInverter(Ion::predefined("He4"), layer.getMaterial()));
}
}
void setup(const SortedSetupOutput &output) override {
AbstractSortedAnalyzer::setup(output);
for (size_t i = 0; i < output.dssdCount(); ++i) {
auto dl = getFrontDeadLayer(output.getDssdOutput(i).detector());
auto dlB = getBackDeadLayer(output.getDssdOutput(i).detector());
auto dlP = getFrontDeadLayer(output.getSingleOutput(i).detector());
deadlayerF.push_back(dl);
deadlayerB.push_back(dlB);
deadlayerP.push_back(dlP);
}
}
void analyze() override {
clear();
CLOCK = output.getScalerOutput("CLOCK").getValue();
findHits();
doAnalysis();
if (mul > 0) { t->Fill(); }
hits.clear();
NUM++;
}
void findHits() {
for (size_t i = 0; i < output.dssdCount(); i++) {
auto &o = output.getDssdOutput(i);
auto &p = output.getSingleOutput(i);
auto &d = o.detector();
auto m = AUSA::mul(o);
for (UInt_t j = 0; j < m; j++) {
Hit hit;
auto dE = fEnergy(o, j) - bEnergy(o, j);
hit.dE = dE;
auto eDssd = energy(o, j);
auto eFDssd = fEnergy(o, j);
auto eBDssd = bEnergy(o, j);
auto ePad = p.energy(0);
auto BI = bSeg(o, j);
auto FI = fSeg(o, j);
hit.fseg = short(FI);
hit.bseg = short(BI);
auto position = o.detector().getUniformPixelPosition(FI, BI);
auto origin = target.getCenter();
hit.position = position;
auto direction = (position - origin).Unit();
hit.direction = direction;
hit.theta = hit.direction.Theta();
if (!simulation) {
hit.TF = fTime(o, j);
hit.TB = bTime(o, j);
hit.TPad = p.time(0);
} else {
hit.TF = 42;
hit.TB = 42;
hit.TPad = 42;
}
auto angle = hit.direction.Angle(-d.getNormal());
hit.angle = angle;
auto tF = deadlayerF[i] / abs(cos(angle));
auto tB = deadlayerB[i] / abs(cos(angle));
auto tP = deadlayerP[i] / abs(cos(angle));
double E = 0.0;
double FE = 0.0;
double BE = 0.0;
E += eDssd;
E += SiCalc->getTotalEnergyCorrection(E, tF);
FE += eFDssd;
FE += SiCalc->getTotalEnergyCorrection(FE, tF);
BE += eBDssd;
BE += SiCalc->getTotalEnergyCorrection(BE, tF);
hit.E = E;
hit.BE = BE;
hit.FE = FE;
auto &from = position;
for (auto &intersection : target.getIntersections(from, target.getCenter())) {
auto &calc = targetCalcs[intersection.index];
hit.E += calc->getTotalEnergyCorrection(hit.E, intersection.transversed);
hit.FE += calc->getTotalEnergyCorrection(hit.FE, intersection.transversed);
hit.BE += calc->getTotalEnergyCorrection(hit.BE, intersection.transversed);
}
hit.index = i;
hit.lVector = {sqrt(2 * hit.E * ALPHA_MASS) * hit.direction, hit.E + ALPHA_MASS};
hits.emplace_back(move(hit));
}
}
}
void doAnalysis() {
if (hits.empty()) return;
mul = hits.size();
for (auto &hit : hits) {
v_pos->add(hit.position);
v_dir->add(hit.direction);
v_theta->add(hit.theta * TMath::RadToDeg());
v_E->add(hit.E);
v_BE->add(hit.BE);
v_FE->add(hit.FE);
v_dE->add(hit.dE);
v_ang->add(hit.angle * TMath::RadToDeg());
v_i->add(static_cast<short>(hit.index));
v_F->add(hit.fseg);
v_B->add(hit.bseg);
v_FT->add(hit.TF);
v_BT->add(hit.TB);
}
}
void terminate() override {
AbstractSortedAnalyzer::terminate();
gDirectory->WriteTObject(t);
}
void clear() {
mul = 0;
AUSA::clear(
*v_E, *v_theta,
*v_i, *v_FE, *v_BE,
*v_F, *v_B,
*v_ang, *v_pos, *v_dir,
*v_dE, *v_FT, *v_BT
);
}
int NUM;
TTree *t;
unique_ptr<DynamicBranchVector<TVector3>> v_dir, v_pos;
unique_ptr<DynamicBranchVector<double>> v_E, v_BE, v_FE, v_theta, v_dE;
unique_ptr<DynamicBranchVector<short>> v_i;
unique_ptr<DynamicBranchVector<short>> v_F, v_B;
unique_ptr<DynamicBranchVector<double>> v_ang;
unique_ptr<DynamicBranchVector<double>> v_FT, v_BT;
UInt_t mul{};
vector<Hit> hits;
unique_ptr<EnergyLossRangeInverter> SiCalc;
vector<unique_ptr<EnergyLossRangeInverter>> targetCalcs;
vector<double> deadlayerF, deadlayerB, deadlayerP;
Target &target;
bool simulation = false;
ParticleType *ALPHA;
UInt_t CLOCK{};
double IonRange, targetThickness;
};
void print(vector<string> const &input) {
for (auto &in : input) {
cout << in << ' ';
}
}
int main(int argc, char *argv[]) {
vector<string> input;
bool sim = false;
int run = 0;
auto setup = JSON::readSetupFromJSON("../setup/setup.json");
auto target = JSON::readTargetFromJSON("../setup/targets/target.json");
for (int i = 1; i < argc; i++) {
try {
run = stoi(argv[i]);
findFilesMatchingWildcard(Form("../data/sorted/%d*", run), input);
} catch (invalid_argument &ia) {
cerr << "SIM!" << endl;
findFilesMatchingWildcard(argv[i], input);
sim = true;
}
}
for (auto &in : input) {
clock_t start = clock();
SortedReader reader{*setup};
reader.add(in);
reader.setVerbose(true);
auto base = stripFileExtension(extractFileName(in));
TString outfile = ("../data/singles/" + base + "lio.root").c_str();
cout << endl << endl << "Reading from: " << in << ' ' << endl;
cout << "Printing to: " << outfile << endl;
TFile output(outfile, "RECREATE");
auto analysis = make_shared<SinglesAnalysis>(target, &output);
analysis->simulation = sim;
reader.attach(analysis);
reader.run();
clock_t stop = clock();
double elapsed = (double) (stop - start) / CLOCKS_PER_SEC;
printf("\nTime elapsed: %.5f\n", elapsed);
}
}
# Tries to determine the location of AUSAlib.
# It will try to look in the systems include/lib paths
#
# You have two options if you want to have a local build of AUSAlib.
#
# 1) Set these two environment variables
# AUSALIB_INC_DIR Include directory
# AUSALIB_LIB_DIR Directory containing the build library
#
# 2) Set this single environment variable
# AUSALIBPATH Path that should contain libAUSA.a and ausa/AUSA.h directly or inside include
if (AUSALIB_INCLUDES)
# Already in cache, be silent
set (AUSALIB_FIND_QUIETLY TRUE)
endif (AUSALIB_INCLUDES)
find_path (AUSALIB_INCLUDES ausa/AUSA.h
HINTS $ENV{AUSALIBPATH}/include ENV AUSALIB_INC_DIR ENV AUSALIBPATH)
find_library(AUSALIB_LIBRARIES AUSA
HINTS ENV AUSALIB_LIB_DIR ENV AUSALIBPATH)
# Handle the QUIETLY and REQUIRED arguments and set AUSALIB_FOUND to TRUE if
# all listed variables are TRUE
include (FindPackageHandleStandardArgs)
find_package_handle_standard_args (AUSALIB DEFAULT_MSG AUSALIB_LIBRARIES AUSALIB_INCLUDES)
\ No newline at end of file
# - Find ROOT instalation
# This module tries to find the ROOT installation on your system.
# It tries to find the root-config script which gives you all the needed information.
# If the system variable ROOTSYS is set this is straight forward.
# If not the module uses the pathes given in ROOT_CONFIG_SEARCHPATH.
# If you need an other path you should add this path to this varaible.
# The root-config script is then used to detect basically everything else.
# This module defines a number of key variables and macros.
# F.Uhlig@gsi.de (fairroot.gsi.de)
MESSAGE(STATUS "Looking for Root...")
SET(ROOT_CONFIG_SEARCHPATH
${SIMPATH}/tools/rootcd /bin
$ENV{ROOTSYS}/bin
)
SET(ROOT_DEFINITIONS "")
SET(ROOT_INSTALLED_VERSION_TOO_OLD FALSE)
SET(ROOT_CONFIG_EXECUTABLE ROOT_CONFIG_EXECUTABLE-NOTFOUND)
FIND_PROGRAM(ROOT_CONFIG_EXECUTABLE NAMES root-config PATHS
${ROOT_CONFIG_SEARCHPATH})
IF (${ROOT_CONFIG_EXECUTABLE} MATCHES "ROOT_CONFIG_EXECUTABLE-NOTFOUND")
MESSAGE( FATAL_ERROR "ROOT not installed in the searchpath and ROOTSYS is not set. Please
set ROOTSYS or add the path to your ROOT installation in the Macro FindROOT.cmake in the
subdirectory cmake/modules.")
ELSE (${ROOT_CONFIG_EXECUTABLE} MATCHES "ROOT_CONFIG_EXECUTABLE-NOTFOUND")
STRING(REGEX REPLACE "(^.*)/bin/root-config" "\\1" ../test ${ROOT_CONFIG_EXECUTABLE})
SET( ENV{ROOTSYS} ${test})
set( ROOTSYS ${test})
ENDIF (${ROOT_CONFIG_EXECUTABLE} MATCHES "ROOT_CONFIG_EXECUTABLE-NOTFOUND")
IF (ROOT_CONFIG_EXECUTABLE)
SET(ROOT_FOUND FALSE)
EXEC_PROGRAM(${ROOT_CONFIG_EXECUTABLE} ARGS "--version" OUTPUT_VARIABLE ROOTVERSION)
MESSAGE(STATUS "Looking for Root... - found $ENV{ROOTSYS}/bin/root")
MESSAGE(STATUS "Looking for Root... - version ${ROOTVERSION} ")
# we need at least version 5.00/00
IF (NOT ROOT_MIN_VERSION)
SET(ROOT_MIN_VERSION "5.00/00")
ENDIF (NOT ROOT_MIN_VERSION)
# now parse the parts of the user given version string into variables
STRING(REGEX REPLACE "^([0-9]+)\\.[0-9][0-9]+\\/[0-9][0-9]+" "\\1" req_root_major_vers "${ROOT_MIN_VERSION}")
STRING(REGEX REPLACE "^[0-9]+\\.([0-9][0-9])+\\/[0-9][0-9]+.*" "\\1" req_root_minor_vers "${ROOT_MIN_VERSION}")
STRING(REGEX REPLACE "^[0-9]+\\.[0-9][0-9]+\\/([0-9][0-9]+)" "\\1" req_root_patch_vers "${ROOT_MIN_VERSION}")
# and now the version string given by qmake
STRING(REGEX REPLACE "^([0-9]+)\\.[0-9][0-9]+\\/[0-9][0-9]+.*" "\\1" found_root_major_vers "${ROOTVERSION}")
STRING(REGEX REPLACE "^[0-9]+\\.([0-9][0-9])+\\/[0-9][0-9]+.*" "\\1" found_root_minor_vers "${ROOTVERSION}")
STRING(REGEX REPLACE "^[0-9]+\\.[0-9][0-9]+\\/([0-9][0-9]+).*" "\\1" found_root_patch_vers "${ROOTVERSION}")
IF (found_root_major_vers LESS 5)
MESSAGE( FATAL_ERROR "Invalid ROOT version \"${ROOTERSION}\", at least major version 4 is required, e.g. \"5.00/00\"")
ENDIF (found_root_major_vers LESS 5)
# compute an overall version number which can be compared at once
MATH(EXPR req_vers "${req_root_major_vers}*10000 + ${req_root_minor_vers}*100 + ${req_root_patch_vers}")
MATH(EXPR found_vers "${found_root_major_vers}*10000 + ${found_root_minor_vers}*100 + ${found_root_patch_vers}")
IF (found_vers LESS req_vers)
SET(ROOT_FOUND FALSE)
SET(ROOT_INSTALLED_VERSION_TOO_OLD TRUE)
ELSE (found_vers LESS req_vers)
SET(ROOT_FOUND TRUE)
ENDIF (found_vers LESS req_vers)
ENDIF (ROOT_CONFIG_EXECUTABLE)
IF (ROOT_FOUND)
# ask root-config for the library dir
# Set ROOT_LIBRARY_DIR
EXEC_PROGRAM( ${ROOT_CONFIG_EXECUTABLE}
ARGS "--libdir"
OUTPUT_VARIABLE ROOT_LIBRARY_DIR_TMP )
IF(EXISTS "${ROOT_LIBRARY_DIR_TMP}")
SET(ROOT_LIBRARY_DIR ${ROOT_LIBRARY_DIR_TMP} )
ELSE(EXISTS "${ROOT_LIBRARY_DIR_TMP}")
MESSAGE("Warning: ROOT_CONFIG_EXECUTABLE reported ${ROOT_LIBRARY_DIR_TMP} as library path,")
MESSAGE("Warning: but ${ROOT_LIBRARY_DIR_TMP} does NOT exist, ROOT must NOT be installed correctly.")
ENDIF(EXISTS "${ROOT_LIBRARY_DIR_TMP}")
# ask root-config for the binary dir
EXEC_PROGRAM(${ROOT_CONFIG_EXECUTABLE}
ARGS "--bindir"
OUTPUT_VARIABLE root_bins )
SET(ROOT_BINARY_DIR ${root_bins})
# ask root-config for the include dir
EXEC_PROGRAM( ${ROOT_CONFIG_EXECUTABLE}
ARGS "--incdir"
OUTPUT_VARIABLE root_headers )
SET(ROOT_INCLUDE_DIR ${root_headers})
# CACHE INTERNAL "")
# ask root-config for the library varaibles
EXEC_PROGRAM( ${ROOT_CONFIG_EXECUTABLE}
# ARGS "--noldflags --noauxlibs --libs"
ARGS "--glibs"
OUTPUT_VARIABLE root_flags )
# STRING(REGEX MATCHALL "([^ ])+" root_libs_all ${root_flags})
# STRING(REGEX MATCHALL "-L([^ ])+" root_library ${root_flags})
# REMOVE_FROM_LIST(root_flags "${root_libs_all}" "${root_library}")
SET(ROOT_LIBRARIES ${root_flags})
# Make variables changeble to the advanced user
MARK_AS_ADVANCED( ROOT_LIBRARY_DIR ROOT_INCLUDE_DIR ROOT_DEFINITIONS)
# Set ROOT_INCLUDES
SET( ROOT_INCLUDES ${ROOT_INCLUDE_DIR})
SET(LD_LIBRARY_PATH ${LD_LIBRARY_PATH} ${ROOT_LIBRARY_DIR})
#######################################
#
# Check the executables of ROOT
# ( rootcint )
#
#######################################
FIND_PROGRAM(ROOT_CINT_EXECUTABLE
NAMES rootcint
PATHS ${ROOT_BINARY_DIR}
NO_DEFAULT_PATH
)
ENDIF (ROOT_FOUND)
###########################################
#
# Macros for building ROOT dictionary
#
###########################################
MACRO (ROOT_GENERATE_DICTIONARY_OLD )
set(INFILES "")
foreach (_current_FILE ${ARGN})
IF (${_current_FILE} MATCHES "^.*\\.h$")
IF (${_current_FILE} MATCHES "^.*Link.*$")
set(LINKDEF_FILE ${_current_FILE})
ELSE (${_current_FILE} MATCHES "^.*Link.*$")
set(INFILES ${INFILES} ${_current_FILE})
ENDIF (${_current_FILE} MATCHES "^.*Link.*$")
ELSE (${_current_FILE} MATCHES "^.*\\.h$")
IF (${_current_FILE} MATCHES "^.*\\.cxx$")
set(OUTFILE ${_current_FILE})
ELSE (${_current_FILE} MATCHES "^.*\\.cxx$")
set(INCLUDE_DIRS ${INCLUDE_DIRS} -I${_current_FILE})
ENDIF (${_current_FILE} MATCHES "^.*\\.cxx$")
ENDIF (${_current_FILE} MATCHES "^.*\\.h$")
endforeach (_current_FILE ${ARGN})
# MESSAGE("INFILES: ${INFILES}")
# MESSAGE("OutFILE: ${OUTFILE}")
# MESSAGE("LINKDEF_FILE: ${LINKDEF_FILE}")
# MESSAGE("INCLUDE_DIRS: ${INCLUDE_DIRS}")
STRING(REGEX REPLACE "(^.*).cxx" "\\1.h" bla "${OUTFILE}")
# MESSAGE("BLA: ${bla}")
SET (OUTFILES ${OUTFILE} ${bla})
ADD_CUSTOM_COMMAND(OUTPUT ${OUTFILES}
COMMAND ${ROOT_CINT_EXECUTABLE}
ARGS -f ${OUTFILE} -c -DHAVE_CONFIG_H ${INCLUDE_DIRS} ${INFILES} ${LINKDEF_FILE} DEPENDS ${INFILES})
# MESSAGE("ROOT_CINT_EXECUTABLE has created the dictionary ${OUTFILE}")
ENDMACRO (ROOT_GENERATE_DICTIONARY_OLD)
###########################################
#
# Macros for building ROOT dictionary