diff --git a/analysis/CMakeLists.txt b/analysis/CMakeLists.txt index 0dc5c11d8e2eb411d395ea55b33955a3ed38a08a..b854f7910bef10c173bab30d3a502f7f9b0cc614 100644 --- a/analysis/CMakeLists.txt +++ b/analysis/CMakeLists.txt @@ -12,6 +12,8 @@ add_subdirectory (PWGCHA/jpsi) add_subdirectory (PWGHAD/hadron) +add_subdirectory (common/analysis_tree_converter) + add_subdirectory (detectors) diff --git a/analysis/common/analysis_tree_converter/CMakeLists.txt b/analysis/common/analysis_tree_converter/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e6ffe26db114ebc2aa32a6fcbfaabe59fd2bbc07 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CMakeLists.txt @@ -0,0 +1,95 @@ +Set(LIBRARY_NAME CbmAnalysisTreeInterface) + +set(SRCS + CbmConverterManager.cxx + CbmConverterTask.cxx + CbmSimEventHeaderConverter.cxx + CbmRecEventHeaderConverter.cxx + CbmStsTracksConverter.cxx + CbmSimTracksConverter.cxx + CbmPsdModulesConverter.cxx + CbmTofHitsConverter.cxx + ) + +Set(INCLUDE_DIRECTORIES + ${CMAKE_CURRENT_SOURCE_DIR} + + ${CBMROOT_SOURCE_DIR}/core/base + ${CBMROOT_SOURCE_DIR}/core/data + ${CBMROOT_SOURCE_DIR}/core/data/base + ${CBMROOT_SOURCE_DIR}/core/data/global + ${CBMROOT_SOURCE_DIR}/core/data/mvd + ${CBMROOT_SOURCE_DIR}/core/data/sts + ${CBMROOT_SOURCE_DIR}/core/data/tof + ${CBMROOT_SOURCE_DIR}/core/data/psd + + ${CBMROOT_SOURCE_DIR}/sim/transport/generators + ${CBMROOT_SOURCE_DIR}/sim/transport/steer + ${CBMROOT_SOURCE_DIR}/sim/transport/geosetup + + ${CBMROOT_SOURCE_DIR}/reco/KF + ${CBMROOT_SOURCE_DIR}/reco/KF/Interface + ${CBMROOT_SOURCE_DIR}/reco/KF/KFQA + ${CBMROOT_SOURCE_DIR}/reco/L1 + ${CBMROOT_SOURCE_DIR}/reco/L1/L1Algo + ${CBMROOT_SOURCE_DIR}/reco/L1/ParticleFinder + ${KFParticle_INCLUDE_DIR} + ${AnalysisTree_INCLUDE_DIR} +) + +Include_Directories (${INCLUDE_DIRECTORIES}) + +Set(SYSTEM_INCLUDE_DIRECTORIES + ${VC_INCLUDE_DIRS} + ${BASE_INCLUDE_DIRECTORIES} + ${Boost_INCLUDE_DIR} +) + +Include_Directories (SYSTEM ${SYSTEM_INCLUDE_DIRECTORIES}) + +set (LINK_DIRECTORIES + ${ROOT_LIBRARY_DIR} + ${FAIRROOT_LIBRARY_DIR} + ${Boost_LIBRARY_DIRS} + ${Vc_LIB_DIR} + ${KFParticle_LIB_DIR} + ${AnalysisTree_LIBRARY_DIR} +) + +link_directories(${LINK_DIRECTORIES}) + +IF (SSE_FOUND) + Message(STATUS "${LIBRARY_NAME} will be compiled with SSE support") + ADD_DEFINITIONS(-DHAVE_SSE) + SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-msse -O3") +ELSE (SSE_FOUND) + MESSAGE(STATUS "${LIBRARY_NAME} will be compiled without SSE support") + SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-msse -O3") +ENDIF (SSE_FOUND) + +Set(LINKDEF CbmAnalysisTreeInterfaceLinkDef.h) + +If(UNIX AND NOT APPLE) + Set(_AnalysisTree_LIB AnalysisTreeBase.so AnalysisTreeInfra.so) + Else() + Set(_AnalysisTree_LIB AnalysisTreeBase.dylib AnalysisTreeInfra.dylib) +EndIf() + +Set(DEPENDENCIES + ${_AnalysisTree_LIB} + CbmData + CbmBase + KF + L1 + Vc.a) +Set(DEFINITIONS -DDO_TPCCATRACKER_EFF_PERFORMANCE -DNonhomogeneousField -DCBM -DUSE_TIMERS) + +ADD_DEFINITIONS(${DEFINITIONS}) + +GENERATE_LIBRARY() + +Install(FILES ../../../macro/analysis/common/analysis_tree_converter/run_treemaker.C + DESTINATION share/cbmroot/macro/analysis_tree + ) + +Add_Dependencies(CbmAnalysisTreeInterface ANALYSISTREE) \ No newline at end of file diff --git a/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h new file mode 100644 index 0000000000000000000000000000000000000000..0c8265b512cac8d7a2faf7f9f3d05fe5b6756c7f --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h @@ -0,0 +1,20 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class CbmConverterManager + ; +#pragma link C++ class CbmConverterTask + ; +#pragma link C++ class CbmSimEventHeaderConverter + ; +#pragma link C++ class CbmRecEventHeaderConverter + ; +#pragma link C++ class CbmSimTracksConverter + ; +#pragma link C++ class CbmStsTracksConverter + ; +#pragma link C++ class CbmTofHitsConverter + ; +#pragma link C++ class CbmPsdModulesConverter + ; + +//#pragma link C++ class CbmMuchTracksConverter+; +//#pragma link C++ class CbmTrdTracksConverter+; +//#pragma link C++ class CbmRichRingsConverter+; + +#endif diff --git a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b30fab670bc5878e996a02c6e1f61861e6173db1 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx @@ -0,0 +1,127 @@ +#include <iostream> + +#include "AnalysisTree/DataHeader.hpp" +#include "AnalysisTree/TaskManager.hpp" + +#include "TGeoBBox.h" +#include "TGeoManager.h" + +#include "CbmConverterManager.h" +#include "CbmConverterTask.h" + +ClassImp(CbmConverterManager) + + InitStatus CbmConverterManager::Init() { + assert(!out_file_name_.empty() && !out_tree_name_.empty()); + + out_file_ = TFile::Open(out_file_name_.c_str(), "RECREATE"); + out_file_->cd(); + out_tree_ = new TTree(out_tree_name_.c_str(), "Analysis Tree"); + + FillDataHeader(); + + std::map<std::string, void*> unused_map {}; + + out_config_ = new AnalysisTree::Configuration; + for (auto* task : tasks_) { + task->SetOutTree(out_tree_); + task->SetOutConfiguration(out_config_); + task->SetDataHeader(data_header_); + task->Init(unused_map); + } + + out_config_->Write("Configuration"); + + return kSUCCESS; +} + +void CbmConverterManager::Exec(Option_t* opt) { + index_map_.clear(); + for (auto* task : tasks_) { + task->SetIndexesMap(&index_map_); + task->Exec(); + // if(!task->GetOutIndexesMap().empty()) { + index_map_.insert( + std::make_pair(task->GetOutputBranchName(), task->GetOutIndexesMap())); + // } + } + out_tree_->Fill(); +} + +void CbmConverterManager::Finish() { + TDirectory* old_dir = gDirectory; + + out_file_->cd(); + for (auto* task : tasks_) { + task->Finish(); + } + out_tree_->Write(); + out_file_->Close(); + + old_dir->cd(); +} + +void CbmConverterManager::FillDataHeader() { + + assert( + !system_.empty() + && beam_mom_); // Force user to write data info //TODO is there a way to + // read it from a file automatically? + + data_header_ = new AnalysisTree::DataHeader; + + std::cout << "ReadDataHeader" << std::endl; + data_header_->SetSystem(system_); + data_header_->SetBeamMomentum(beam_mom_); + + auto& psd_mod_pos = data_header_->AddDetector(); + const int psd_node_id = 6; + const char* module_name_prefix = "module"; + + auto* geoMan = gGeoManager; + auto* caveNode = geoMan->GetTopNode(); + auto* psdNode = caveNode->GetDaughter(psd_node_id); + std::cout << "-I- " << psdNode->GetName() << std::endl; + + auto psdGeoMatrix = psdNode->GetMatrix(); + auto psdBox = (TGeoBBox*) psdNode->GetVolume()->GetShape(); + TVector3 frontFaceLocal(0, 0, -psdBox->GetDZ()); + + TVector3 frontFaceGlobal; + psdGeoMatrix->LocalToMaster(&frontFaceLocal[0], &frontFaceGlobal[0]); + + for (int i_d = 0; i_d < psdNode->GetNdaughters(); ++i_d) { + auto* daughter = psdNode->GetDaughter(i_d); + TString daughterName(daughter->GetName()); + if (daughterName.BeginsWith(module_name_prefix)) { + + auto geoMatrix = daughter->GetMatrix(); + TVector3 translation(geoMatrix->GetTranslation()); + + int modID = daughter->GetNumber(); + double x = translation.X(); + double y = translation.Y(); + + std::cout << "mod" << modID << " : " << Form("(%.3f, %3f)", x, y) + << std::endl; + + auto* module = psd_mod_pos.AddChannel(); + module->SetPosition(x, y, frontFaceGlobal[2]); + } + } + + TDirectory* curr = gDirectory; // TODO check why this is needed + TFile* currentFile = gFile; + + out_file_->cd(); + data_header_->Write("DataHeader"); + + gFile = currentFile; + gDirectory = curr; +} +CbmConverterManager::~CbmConverterManager() { + delete out_config_; + delete data_header_; + delete out_tree_; + delete out_file_; +} diff --git a/analysis/common/analysis_tree_converter/CbmConverterManager.h b/analysis/common/analysis_tree_converter/CbmConverterManager.h new file mode 100644 index 0000000000000000000000000000000000000000..0eb217ff98dffb9076a8a27a73f84e87057739d6 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmConverterManager.h @@ -0,0 +1,51 @@ +#ifndef ANALYSIS_TREE_CONVERTERMANAGER_H_ +#define ANALYSIS_TREE_CONVERTERMANAGER_H_ + +#include "FairTask.h" + +namespace AnalysisTree { + class Configuration; + class DataHeader; +} // namespace AnalysisTree + +class CbmConverterTask; + +class CbmConverterManager : public FairTask { + +public: + CbmConverterManager() = default; + ~CbmConverterManager() override; + + InitStatus Init() override; + void Exec(Option_t* opt) override; + void Finish() override; + + void AddTask(CbmConverterTask* task) { tasks_.emplace_back(task); } + void SetOutFileName(std::string name) { out_file_name_ = std::move(name); } + void SetOutTreeName(std::string name) { out_tree_name_ = std::move(name); } + + void SetSystem(const std::string& system) { system_ = system; } + void SetBeamMomentum(float beam_mom) { beam_mom_ = beam_mom; } + +private: + void FillDataHeader(); + + TFile* out_file_ {nullptr}; + TTree* out_tree_ {nullptr}; + std::string out_file_name_ {""}; + std::string out_tree_name_ {""}; + + std::string system_ {""}; + float beam_mom_ {0.}; + + AnalysisTree::Configuration* out_config_ {nullptr}; + AnalysisTree::DataHeader* data_header_ {nullptr}; + std::vector<CbmConverterTask*> tasks_ {}; + + std::map<std::string, std::map<int, int>> + index_map_ {}; ///< map CbmRoot to AT of indexes for a given branch + + ClassDefOverride(CbmConverterManager, 1) +}; + +#endif // ANALYSIS_TREE_CONVERTERMANAGER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmConverterTask.cxx b/analysis/common/analysis_tree_converter/CbmConverterTask.cxx new file mode 100644 index 0000000000000000000000000000000000000000..59dc5bb2fbd6daa0f6e272570e51db62b0b5a3ea --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmConverterTask.cxx @@ -0,0 +1 @@ +#include "CbmConverterTask.h" diff --git a/analysis/common/analysis_tree_converter/CbmConverterTask.h b/analysis/common/analysis_tree_converter/CbmConverterTask.h new file mode 100644 index 0000000000000000000000000000000000000000..e37d0cf74282056953452ddde288f56dbeb54b66 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmConverterTask.h @@ -0,0 +1,39 @@ +#ifndef ANALYSIS_TREE_CONVERTERTASK_H_ +#define ANALYSIS_TREE_CONVERTERTASK_H_ + +#include "AnalysisTree/FillTask.hpp" + +class FairRootManager; + +class CbmConverterTask : public AnalysisTree::FillTask { + +public: + CbmConverterTask() = delete; + explicit CbmConverterTask(std::string out_branch_name, + std::string match_to = "") { + out_branch_ = std::move(out_branch_name); + match_to_ = std::move(match_to); + }; + + ~CbmConverterTask() override = default; + + // void SetIoMan(FairRootManager* man) { ioman_ = man; } + const std::map<int, int>& GetOutIndexesMap() const { + return out_indexes_map_; + } + + void SetIndexesMap(std::map<std::string, std::map<int, int>>* indexes_map) { + indexes_map_ = indexes_map; + } + +protected: + // FairRootManager* ioman_{nullptr}; + std::map<int, int> out_indexes_map_ {}; ///< CbmRoot to AnalysisTree indexes + ///< map for output branch + std::map<std::string, std::map<int, int>>* + indexes_map_ {}; ///< CbmRoot to AnalysisTree indexes map for branches + ///< from other tasks + std::string match_to_ {}; ///< AT branch to match +}; + +#endif // ANALYSIS_TREE_CONVERTERTASK_H_ diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..02dbf1f8eeb67434f2766fc06ec762a8711c165f --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx @@ -0,0 +1,66 @@ +#include <cassert> +#include <vector> + +#include "TClonesArray.h" + +#include "FairRootManager.h" + +#include "AnalysisTree/DataHeader.hpp" +#include "AnalysisTree/Detector.hpp" + +#include "CbmPsdHit.h" +#include "CbmPsdModulesConverter.h" + +ClassImp(CbmPsdModulesConverter) + + void CbmPsdModulesConverter::Init(std::map<std::string, void*>&) { + assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ + && out_tree_); + auto* ioman = FairRootManager::Instance(); + assert(ioman != nullptr); + cbm_psd_hits_ = (TClonesArray*) ioman->GetObject("PsdHit"); + + AnalysisTree::BranchConfig psd_branch(out_branch_, + AnalysisTree::DetType::kModule); + out_config_->AddBranchConfig(std::move(psd_branch)); + psd_modules_ = new AnalysisTree::ModuleDetector(out_config_->GetLastId()); + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::ModuleDetector", &psd_modules_); +} + +void CbmPsdModulesConverter::Exec() { + assert(cbm_psd_hits_); + psd_modules_->ClearChannels(); + + CbmPsdHit* hit {nullptr}; + Float_t psd_energy {0.}; + + // //TODO fix this logic. At the moment I don't know how to access number of + // PSD modules from her hit = (CbmPsdHit*) cbm_psd_hits_->At(nPsdHits-1); // + // last hit const int last_module_id = hit->GetModuleID()-1; + const int n_psd_modules = + data_header_->GetModulePositions(0).GetNumberOfChannels(); + psd_modules_->Reserve(n_psd_modules); + for (int i = 0; i < n_psd_modules; ++i) { + auto* module = psd_modules_->AddChannel(); + module->SetSignal(0.f); + } + + const int nPsdHits = cbm_psd_hits_->GetEntriesFast(); + for (int i = 0; i < nPsdHits; ++i) { + hit = (CbmPsdHit*) cbm_psd_hits_->At(i); + if (hit == nullptr) continue; + auto& module = psd_modules_->GetChannel(hit->GetModuleID() - 1); + module.SetNumber(i + 1); + module.SetSignal(hit->GetEdep()); + psd_energy += hit->GetEdep(); + } + // rec_event_header_ -> SetField(psd_energy, + // config_.GetBranchConfig(rec_event_header_->GetId() ).GetFieldId("Epsd")); +} + +void CbmPsdModulesConverter::Finish() { + delete psd_modules_; + delete cbm_psd_hits_; +} +CbmPsdModulesConverter::~CbmPsdModulesConverter() { delete psd_modules_; }; diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..a96028c6756d0830e3e07737bf64f9fd557a0b05 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h @@ -0,0 +1,31 @@ +#ifndef ANALYSIS_TREE_PSDMODULESCONVERTER_H_ +#define ANALYSIS_TREE_PSDMODULESCONVERTER_H_ + +#include "AnalysisTree/Detector.hpp" +#include "CbmConverterTask.h" + +class TClonesArray; + +class CbmPsdModulesConverter final : public CbmConverterTask { +public: + explicit CbmPsdModulesConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) { + out_branch_ = "PsdModules"; + in_branches_.emplace_back("PsdHits"); + }; + + ~CbmPsdModulesConverter() final; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final; + +private: + AnalysisTree::ModuleDetector* psd_modules_ {nullptr}; + TClonesArray* cbm_psd_hits_ {nullptr}; + + ClassDef(CbmPsdModulesConverter, 1) +}; + +#endif // ANALYSIS_TREE_PSDMODULESCONVERTER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e5bdcc73a94347829be6252acc894602f7123db4 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx @@ -0,0 +1,53 @@ +#include "CbmVertex.h" +#include "FairMCEventHeader.h" +#include "FairRootManager.h" +#include "TClonesArray.h" +#include "cassert" + +#include "CbmRecEventHeaderConverter.h" + +ClassImp(CbmRecEventHeaderConverter) + + void CbmRecEventHeaderConverter::Init(std::map<std::string, void*>&) { + assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ + && out_tree_); + auto* ioman = FairRootManager::Instance(); + assert(ioman != nullptr); + cbm_header_ = (FairMCEventHeader*) ioman->GetObject(in_branches_[0].c_str()); + cbm_prim_vertex_ = (CbmVertex*) ioman->GetObject(in_branches_[1].c_str()); + + // ***** RecEventHeader ******* + AnalysisTree::BranchConfig RecEventHeaderBranch( + "RecEventHeader", AnalysisTree::DetType::kEventHeader); + RecEventHeaderBranch.AddField<float>("psi_EP"); + RecEventHeaderBranch.AddField<float>("vtx_chi2"); + RecEventHeaderBranch.AddField<float>("Epsd"); + RecEventHeaderBranch.AddField<int>("M"); + RecEventHeaderBranch.AddField<int>("evt_id"); + out_config_->AddBranchConfig(RecEventHeaderBranch); + rec_event_header_ = new AnalysisTree::EventHeader(out_config_->GetLastId()); + rec_event_header_->Init(RecEventHeaderBranch); + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::EventHeader", &rec_event_header_); +} + +void CbmRecEventHeaderConverter::Exec() { + const auto& conf = out_config_->GetBranchConfig(rec_event_header_->GetId()); + if (!cbm_prim_vertex_) { + std::cout << "WARNING! No fPrimVtx!" << std::endl; + + rec_event_header_->SetVertexPosition3({-999., -999., -999.}); + rec_event_header_->SetField(float(-999.), conf.GetFieldId("vtx_chi2")); + return; + } + rec_event_header_->SetVertexPosition3({cbm_prim_vertex_->GetX(), + cbm_prim_vertex_->GetY(), + cbm_prim_vertex_->GetZ()}); + // rec_event_header_ -> SetField(int(cbm_sts_tracks_->GetEntries()), + // config_.GetBranchConfig(rec_event_header_->GetId() ).GetFieldId("M")); + rec_event_header_->SetField(int(cbm_header_->GetEventID()), + conf.GetFieldId("evt_id")); + rec_event_header_->SetField( + float(cbm_prim_vertex_->GetChi2() / cbm_prim_vertex_->GetNDF()), + conf.GetFieldId("vtx_chi2")); +} diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..1208bf7361f0f2123f47fb00daf223967da0dd06 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h @@ -0,0 +1,33 @@ +#ifndef ANALYSIS_TREE_RECEVENTHEADERCONVERTER_H_ +#define ANALYSIS_TREE_RECEVENTHEADERCONVERTER_H_ + +#include "AnalysisTree/EventHeader.hpp" + +#include "CbmConverterTask.h" + +class FairMCEventHeader; +class CbmVertex; + +class CbmRecEventHeaderConverter final : public CbmConverterTask { +public: + explicit CbmRecEventHeaderConverter(std::string out_branch_name) + : CbmConverterTask(std::move(out_branch_name)) { + in_branches_.emplace_back("MCEventHeader."); + in_branches_.emplace_back("PrimaryVertex."); + }; + ~CbmRecEventHeaderConverter() final = default; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final { delete rec_event_header_; }; + +private: + AnalysisTree::EventHeader* rec_event_header_ {nullptr}; + + FairMCEventHeader* cbm_header_ {nullptr}; + CbmVertex* cbm_prim_vertex_ {nullptr}; + + ClassDef(CbmRecEventHeaderConverter, 1) +}; + +#endif // ANALYSIS_TREE_RECEVENTHEADERCONVERTER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..77f7bdd263550976143d8d17667f4702021b531b --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx @@ -0,0 +1,46 @@ +#include "FairMCEventHeader.h" +#include "FairRootManager.h" +#include "TClonesArray.h" +#include "cassert" + +#include "CbmSimEventHeaderConverter.h" + +ClassImp(CbmSimEventHeaderConverter) + + void CbmSimEventHeaderConverter::Init(std::map<std::string, void*>&) { + assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ + && out_tree_); + auto* ioman = FairRootManager::Instance(); + assert(ioman != nullptr); + cbm_header_ = (FairMCEventHeader*) ioman->GetObject(in_branches_[0].c_str()); + + // ***** SimEventHeader ******* + AnalysisTree::BranchConfig SimEventHeaderBranch( + "SimEventHeader", AnalysisTree::DetType::kEventHeader); + SimEventHeaderBranch.AddField<float>("psi_RP"); + SimEventHeaderBranch.AddField<float>("b"); + out_config_->AddBranchConfig(SimEventHeaderBranch); + sim_event_header_ = new AnalysisTree::EventHeader(out_config_->GetLastId()); + sim_event_header_->Init(SimEventHeaderBranch); + + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::EventHeader", &sim_event_header_); +} + +void CbmSimEventHeaderConverter::Exec() { + if (!cbm_header_) { + throw std::runtime_error( + "CbmSimEventHeaderConverter::Exec - ERROR! No fHeader!"); + } + + TVector3 pos = { + cbm_header_->GetX(), cbm_header_->GetY(), cbm_header_->GetZ()}; + sim_event_header_->SetVertexPosition3(pos); + sim_event_header_->SetField( + float(cbm_header_->GetRotZ()), + out_config_->GetBranchConfig(sim_event_header_->GetId()) + .GetFieldId("psi_RP")); + sim_event_header_->SetField( + float(cbm_header_->GetB()), + out_config_->GetBranchConfig(sim_event_header_->GetId()).GetFieldId("b")); +} diff --git a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..db81811f9f81bae9c2c6ec7f3219748eab1dcc3c --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h @@ -0,0 +1,30 @@ +#ifndef ANALYSIS_TREE_SIMEVENTHEADERCONVERTER_H_ +#define ANALYSIS_TREE_SIMEVENTHEADERCONVERTER_H_ + +#include "AnalysisTree/EventHeader.hpp" + +#include "CbmConverterTask.h" + +class FairMCEventHeader; +class CbmVertex; + +class CbmSimEventHeaderConverter final : public CbmConverterTask { +public: + explicit CbmSimEventHeaderConverter(std::string out_branch_name) + : CbmConverterTask(std::move(out_branch_name)) { + in_branches_.emplace_back("MCEventHeader."); + }; + ~CbmSimEventHeaderConverter() final = default; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final { delete sim_event_header_; }; + +private: + AnalysisTree::EventHeader* sim_event_header_ {nullptr}; + FairMCEventHeader* cbm_header_ {nullptr}; + + ClassDef(CbmSimEventHeaderConverter, 1) +}; + +#endif // ANALYSIS_TREE_SIMEVENTHEADERCONVERTER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..20122644f0e6bf25c1487140b3e8cd5490d3c239 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx @@ -0,0 +1,68 @@ +#include <cassert> +#include <vector> + +#include "TClonesArray.h" + +#include "FairRootManager.h" + +#include "CbmMCTrack.h" +#include "CbmSimTracksConverter.h" + +ClassImp(CbmSimTracksConverter) + + void CbmSimTracksConverter::Init(std::map<std::string, void*>&) { + + assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ + && out_tree_); + auto* ioman = FairRootManager::Instance(); + cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject(in_branches_[0].c_str()); + + AnalysisTree::BranchConfig SimTracksBranch(out_branch_, + AnalysisTree::DetType::kParticle); + SimTracksBranch.AddField<int>("mother_id"); + out_config_->AddBranchConfig(std::move(SimTracksBranch)); + sim_tracks_ = new AnalysisTree::Particles(out_config_->GetLastId()); + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::Particles", &sim_tracks_); +} + +void CbmSimTracksConverter::Exec() { + assert(cbm_mc_tracks_); + out_indexes_map_.clear(); + + std::cout << "ReadMcTracks" << std::endl; + sim_tracks_->ClearChannels(); + + const int nMcTracks = cbm_mc_tracks_->GetEntries(); + const int imother_id = + out_config_->GetBranchConfig(sim_tracks_->GetId()).GetFieldId("mother_id"); + sim_tracks_->Reserve(nMcTracks); + + for (int iMcTrack = 0; iMcTrack < nMcTracks; ++iMcTrack) { + const auto* mctrack = (CbmMCTrack*) cbm_mc_tracks_->At(iMcTrack); + if (mctrack->GetStartZ() > 200.) + continue; // to exclude shower in PSD etc //TODO check if better cut can + // be applied. + + auto* Track = sim_tracks_->AddChannel(); + Track->Init(out_config_->GetBranchConfig(sim_tracks_->GetId())); + + out_indexes_map_.insert(std::make_pair(iMcTrack, Track->GetId())); + + Track->SetMomentum(mctrack->GetPx(), mctrack->GetPy(), mctrack->GetPz()); + Track->SetMass(float(mctrack->GetMass())); + Track->SetPid(int(mctrack->GetPdgCode())); + + if (mctrack->GetMotherId() + == -1) { // mother id should < than track id, so we can do it here + Track->SetField(int(-1), imother_id); + } else { + auto p = out_indexes_map_.find(mctrack->GetMotherId()); + if (p == out_indexes_map_.end()) // match is not found + Track->SetField(int(-999), imother_id); + else + Track->SetField(int(p->second), imother_id); + } + } +} +CbmSimTracksConverter::~CbmSimTracksConverter() { delete sim_tracks_; }; diff --git a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..dd3c92a4bd601db7308d93e02b1b9b56f15b3a01 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h @@ -0,0 +1,31 @@ +#ifndef ANALYSIS_TREE_SIMTRACKSCONVERTER_H_ +#define ANALYSIS_TREE_SIMTRACKSCONVERTER_H_ + +#include "AnalysisTree/Detector.hpp" + +#include "CbmConverterTask.h" + +class TClonesArray; + +class CbmSimTracksConverter final : public CbmConverterTask { +public: + explicit CbmSimTracksConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) { + in_branches_.emplace_back("MCTrack"); + }; + + ~CbmSimTracksConverter() final; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final {}; + +private: + AnalysisTree::Particles* sim_tracks_ {nullptr}; + TClonesArray* cbm_mc_tracks_ {nullptr}; + + ClassDef(CbmSimTracksConverter, 1) +}; + +#endif diff --git a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..91f2dc4e1a4a9c134cf0104a153916f5745ffd00 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx @@ -0,0 +1,341 @@ +#include <CbmMCTrack.h> +#include <cassert> + +#include "TClonesArray.h" + +#include "FairRootManager.h" + +#include "AnalysisTree/Matching.hpp" + +#include "CbmStsTrack.h" +#include "CbmTrackMatchNew.h" +#include "CbmVertex.h" +#include "Interface/CbmKFVertex.h" +#include "L1Field.h" +#include "ParticleFinder/CbmL1PFFitter.h" + +#include "CbmStsTracksConverter.h" + +ClassImp(CbmStsTracksConverter) + + void CbmStsTracksConverter::Exec() { + assert(cbm_sts_tracks_ != nullptr && cbm_mc_tracks_ != nullptr); + + vtx_tracks_2_sim_->Clear(); + out_indexes_map_.clear(); + + ReadVertexTracks(); + MapTracks(); +} + +CbmStsTracksConverter::~CbmStsTracksConverter() { + delete vtx_tracks_; + delete vtx_tracks_2_sim_; +}; + +// TODO misleading name, move field filling somewhere else? +float CbmStsTracksConverter::ExtrapolateToVertex(CbmStsTrack* sts_track, + AnalysisTree::Track* track, + int pdg) { + + vector<CbmStsTrack> tracks = {*sts_track}; + CbmL1PFFitter fitter; + vector<float> chi2_to_vtx; + vector<L1FieldRegion> field; + CbmKFVertex kfVertex = CbmKFVertex(*cbm_prim_vertex_); + if (is_reproduce_cbmkfpf_) { + std::vector<int> pdgVector = {pdg}; + fitter.Fit(tracks, pdgVector); + } + fitter.GetChiToVertex(tracks, field, chi2_to_vtx, kfVertex, 3.); + *sts_track = tracks[0]; + + if (is_write_kfinfo_) { + track->SetField(float(field.at(0).cx0[0]), imf_); + track->SetField(float(field.at(0).cx1[0]), imf_ + 1); + track->SetField(float(field.at(0).cx2[0]), imf_ + 2); + track->SetField(float(field.at(0).cy0[0]), imf_ + 3); + track->SetField(float(field.at(0).cy1[0]), imf_ + 4); + track->SetField(float(field.at(0).cy2[0]), imf_ + 5); + track->SetField(float(field.at(0).cz0[0]), imf_ + 6); + track->SetField(float(field.at(0).cz1[0]), imf_ + 7); + track->SetField(float(field.at(0).cz2[0]), imf_ + 8); + track->SetField(float(field.at(0).z0[0]), imf_ + 9); + } + + return chi2_to_vtx[0]; +} + +void CbmStsTracksConverter::ReadVertexTracks() { + assert(cbm_prim_vertex_ && cbm_sts_tracks_); + + vtx_tracks_->ClearChannels(); + + const int iq = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("q"); + const int indf = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("ndf"); + const int ichi2 = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("chi2"); + const int inhits = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("nhits"); + const int inhits_mvd = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("nhits_mvd"); + const int idcax = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("dcax"); + const int ivtx_chi2 = + out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("vtx_chi2"); + + const int n_sts_tracks = cbm_sts_tracks_->GetEntries(); + vtx_tracks_->Reserve(n_sts_tracks); + + for (short i_track = 0; i_track < n_sts_tracks; ++i_track) { + auto* sts_track = (CbmStsTrack*) cbm_sts_tracks_->At(i_track); + if (!sts_track) { throw std::runtime_error("empty track!"); } + + auto* track = vtx_tracks_->AddChannel(); + track->Init(out_config_->GetBranchConfig(vtx_tracks_->GetId())); + + int pdg = GetMcPid( + sts_track, (CbmTrackMatchNew*) cbm_sts_match_->At(i_track), track); + + float chi2_vertex = ExtrapolateToVertex(sts_track, track, pdg); + const FairTrackParam* trackParamFirst = sts_track->GetParamFirst(); + TVector3 momRec; + trackParamFirst->Momentum(momRec); + const Int_t q = trackParamFirst->GetQp() > 0 ? 1 : -1; + + track->SetMomentum3(momRec); + track->SetField(int(q), iq); + track->SetField(int(sts_track->GetNDF()), indf); + track->SetField(float(sts_track->GetChiSq()), ichi2); + track->SetField(int(sts_track->GetNofHits()), inhits); + track->SetField(float(trackParamFirst->GetX() - cbm_prim_vertex_->GetX()), + idcax); + track->SetField(float(trackParamFirst->GetY() - cbm_prim_vertex_->GetY()), + idcax + 1); + track->SetField(float(trackParamFirst->GetZ() - cbm_prim_vertex_->GetZ()), + idcax + 2); + track->SetField(int(sts_track->GetNofMvdHits()), inhits_mvd); + track->SetField(float(chi2_vertex), ivtx_chi2); + + out_indexes_map_.insert(std::make_pair(i_track, track->GetId())); + + if (is_write_kfinfo_) { WriteKFInfo(track, sts_track); } + } +} + +void CbmStsTracksConverter::WriteKFInfo(AnalysisTree::Track* track, + const CbmStsTrack* sts_track) const { + assert(track && sts_track); + const FairTrackParam* trackParamFirst = sts_track->GetParamFirst(); + + track->SetField(float(trackParamFirst->GetX()), ipar_); + track->SetField(float(trackParamFirst->GetY()), ipar_ + 1); + track->SetField(float(trackParamFirst->GetZ()), ipar_ + 2); + track->SetField(float(trackParamFirst->GetTx()), ipar_ + 3); + track->SetField(float(trackParamFirst->GetTy()), ipar_ + 4); + track->SetField(float(trackParamFirst->GetQp()), ipar_ + 5); + + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { + track->SetField(float(trackParamFirst->GetCovariance(i, j)), + icov_ + iCov); + } + } + track->SetField(IsGoodCovMatrix(sts_track), ipasscuts_); +} + +bool CbmStsTracksConverter::IsGoodCovMatrix( + const CbmStsTrack* sts_track) const { + + if (!is_reproduce_cbmkfpf_) return true; + + assert(sts_track); + const FairTrackParam* trackParamFirst = sts_track->GetParamFirst(); + + Double_t cov_matrix[15] = {0.f}; + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { + cov_matrix[iCov] = trackParamFirst->GetCovariance(i, j); + } + } + // Cuts, coded in MZ's CbmKFParticleFinder.cxx + bool ok = true; + ok = ok && finite(sts_track->GetParamFirst()->GetX()); + ok = ok && finite(sts_track->GetParamFirst()->GetY()); + ok = ok && finite(sts_track->GetParamFirst()->GetZ()); + ok = ok && finite(sts_track->GetParamFirst()->GetTx()); + ok = ok && finite(sts_track->GetParamFirst()->GetTy()); + ok = ok && finite(sts_track->GetParamFirst()->GetQp()); + + for (auto element : cov_matrix) { + ok = ok && finite(element); + } + ok = ok && (cov_matrix[0] < 1. && cov_matrix[0] > 0.) + && (cov_matrix[2] < 1. && cov_matrix[2] > 0.) + && (cov_matrix[5] < 1. && cov_matrix[5] > 0.) + && (cov_matrix[9] < 1. && cov_matrix[9] > 0.) + && (cov_matrix[14] < 1. && cov_matrix[14] > 0.); + + ok = ok && sts_track->GetChiSq() < 10. * sts_track->GetNDF(); + + return ok; +} + +int CbmStsTracksConverter::GetMcPid(const CbmStsTrack* sts_track, + const CbmTrackMatchNew* match, + AnalysisTree::Track* track) const { + + if (!is_write_kfinfo_) { return -2; } + //-----------PID as in MZ's + //CbmKFPF---------------------------------------------------------- + Int_t nMCTracks = cbm_mc_tracks_->GetEntriesFast(); + + int pdgCode = -2; + if (match->GetNofLinks() > 0) { // there is at least one matched MC-track + Float_t bestWeight = 0.f; + Float_t totalWeight = 0.f; + Int_t mcTrackId = -1; + + for (int iLink = 0; iLink < match->GetNofLinks(); iLink++) { + totalWeight += match->GetLink(iLink).GetWeight(); + if (match->GetLink(iLink).GetWeight() > bestWeight) { + bestWeight = match->GetLink(iLink).GetWeight(); + mcTrackId = match->GetLink(iLink).GetIndex(); + } + } + + if (!((bestWeight / totalWeight < 0.7) + || (mcTrackId >= nMCTracks || mcTrackId < 0))) { + auto* mctrack = (CbmMCTrack*) cbm_mc_tracks_->At(mcTrackId); + + if (!(TMath::Abs(mctrack->GetPdgCode()) == 11 + || TMath::Abs(mctrack->GetPdgCode()) == 13 + || TMath::Abs(mctrack->GetPdgCode()) == 211 + || TMath::Abs(mctrack->GetPdgCode()) == 321 + || TMath::Abs(mctrack->GetPdgCode()) == 2212 + || TMath::Abs(mctrack->GetPdgCode()) == 1000010020 + || TMath::Abs(mctrack->GetPdgCode()) == 1000010030 + || TMath::Abs(mctrack->GetPdgCode()) == 1000020030 + || TMath::Abs(mctrack->GetPdgCode()) == 1000020040)) { + pdgCode = -1; + } else { + pdgCode = mctrack->GetPdgCode(); + } + if (mctrack->GetMotherId() > -1) { + track->SetField( + int(((CbmMCTrack*) cbm_mc_tracks_->At(mctrack->GetMotherId())) + ->GetPdgCode()), + imother_pdg_); + } + } + } + track->SetField(pdgCode, imc_pdg_); + + return pdgCode; +} + +void CbmStsTracksConverter::InitInput() { + auto* ioman = FairRootManager::Instance(); + + cbm_prim_vertex_ = + (CbmVertex*) ioman->GetObject(in_branches_.at(ePrimiryVertex).c_str()); + cbm_sts_tracks_ = + (TClonesArray*) ioman->GetObject(in_branches_.at(eStsTracks).c_str()); + cbm_mc_tracks_ = + (TClonesArray*) ioman->GetObject(in_branches_.at(eSimTracks).c_str()); + cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch"); +} + +void CbmStsTracksConverter::Init(map<std::string, void*>& branches) { + InitInput(); + + AnalysisTree::BranchConfig vtx_tracks_config(out_branch_, + AnalysisTree::DetType::kTrack); + vtx_tracks_config.AddField<float>("chi2"); + vtx_tracks_config.AddField<float>("vtx_chi2"); + vtx_tracks_config.AddFields<float>({"dcax", "dcay", "dcaz"}); + vtx_tracks_config.AddField<int>("ndf"); + vtx_tracks_config.AddField<int>("q"); + vtx_tracks_config.AddField<int>("nhits"); + vtx_tracks_config.AddField<int>("nhits_mvd"); + + if (is_write_kfinfo_) { + vtx_tracks_config.AddFields<float>({"x", "y", "z", "tx", "ty", "qp"}); + vtx_tracks_config.AddFields<float>( + {"cx0", "cx1", "cx2", "cy0", "cy1", "cy2", "cz0", "cz1", "cz2", "z0"}); + vtx_tracks_config.AddFields<float>({"cov1", + "cov2", + "cov3", + "cov4", + "cov5", + "cov6", + "cov7", + "cov8", + "cov9", + "cov10", + "cov11", + "cov12", + "cov13", + "cov14", + "cov15"}); + + vtx_tracks_config.AddField<int>("mother_pdg"); + vtx_tracks_config.AddField<int>("mc_pdg"); + vtx_tracks_config.AddField<bool>("pass_cuts"); + + ipar_ = vtx_tracks_config.GetFieldId("x"); + imf_ = vtx_tracks_config.GetFieldId("cx0"); + icov_ = vtx_tracks_config.GetFieldId("cov1"); + imc_pdg_ = vtx_tracks_config.GetFieldId("mc_pdg"); + imother_pdg_ = vtx_tracks_config.GetFieldId("mother_pdg"); + ipasscuts_ = vtx_tracks_config.GetFieldId("pass_cuts"); + } + + out_config_->AddBranchConfig(std::move(vtx_tracks_config)); + vtx_tracks_ = new AnalysisTree::TrackDetector(out_config_->GetLastId()); + vtx_tracks_2_sim_ = new AnalysisTree::Matching( + out_config_->GetLastId(), out_config_->GetBranchConfig(match_to_).GetId()); + + out_config_->AddMatch(vtx_tracks_2_sim_); + + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::TrackDetector", &vtx_tracks_); + out_tree_->Branch((out_branch_ + "2" + match_to_).c_str(), + "AnalysisTree::Matching", + &vtx_tracks_2_sim_); + + branches.emplace(out_branch_, vtx_tracks_); +} + +void CbmStsTracksConverter::MapTracks() { + const auto it = indexes_map_->find(match_to_); + if (it == indexes_map_->end()) { + throw std::runtime_error(match_to_ + + " is not found to match with vertex tracks"); + } + auto sim_tracks_map = it->second; + + assert(!out_indexes_map_.empty() || !sim_tracks_map.empty()); + CbmTrackMatchNew* match {nullptr}; + + for (const auto& track_id : out_indexes_map_) { + const int cbm_id = track_id.first; + const int out_id = track_id.second; + match = (CbmTrackMatchNew*) cbm_sts_match_->At(cbm_id); + if (match->GetNofLinks() > 0) { // there is at least one matched MC-track + const int mc_track_id = match->GetMatchedLink().GetIndex(); + + if (match->GetTrueOverAllHitsRatio() + < 0.5) // not enough hits to be a match + continue; + + auto p = sim_tracks_map.find(mc_track_id); + if (p == sim_tracks_map.end()) // match is not found + continue; + + vtx_tracks_2_sim_->AddMatch(out_id, p->second); + } + } +} diff --git a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..24dd66aa8f70962419cb5f43ea96fe4a6b4af6b2 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h @@ -0,0 +1,82 @@ +#ifndef ANALYSIS_TREE_STSTRACKSCONVERTER_H_ +#define ANALYSIS_TREE_STSTRACKSCONVERTER_H_ + +#include "CbmConverterTask.h" + +#include "AnalysisTree/Detector.hpp" + +class TClonesArray; +class CbmVertex; +class CbmStsTrack; +class CbmTrackMatchNew; + +namespace AnalysisTree { + class Matching; +} + +class CbmStsTracksConverter final : public CbmConverterTask { + +public: + enum kInBranches { + eStsTracks, + ePrimiryVertex, + eSimTracks, + eNumberOfInputBranches + }; + + explicit CbmStsTracksConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) { + in_branches_.resize(eNumberOfInputBranches); + in_branches_.at(eStsTracks) = "StsTrack"; + in_branches_.at(ePrimiryVertex) = "PrimaryVertex."; + in_branches_.at(eSimTracks) = "MCTrack"; + } + + ~CbmStsTracksConverter() final; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final {} + + void SetIsWriteKFInfo(bool is = true) { is_write_kfinfo_ = is; } + void SetIsReproduceCbmKFPF(bool is = true) { is_reproduce_cbmkfpf_ = is; } + +private: + void ReadVertexTracks(); + void MapTracks(); + void InitInput(); + float ExtrapolateToVertex(CbmStsTrack* sts_track, + AnalysisTree::Track* track, + int pdg); + void WriteKFInfo(AnalysisTree::Track* track, + const CbmStsTrack* sts_track) const; + bool IsGoodCovMatrix(const CbmStsTrack* sts_track) const; + int GetMcPid(const CbmStsTrack* sts_track, + const CbmTrackMatchNew* match, + AnalysisTree::Track* track) const; + + AnalysisTree::TrackDetector* vtx_tracks_ { + nullptr}; ///< raw pointers are needed for TTree::Branch + AnalysisTree::Matching* vtx_tracks_2_sim_ { + nullptr}; ///< raw pointers are needed for TTree::Branch + + CbmVertex* cbm_prim_vertex_ {nullptr}; ///< non-owning pointer + TClonesArray* cbm_mc_tracks_ {nullptr}; ///< non-owning pointer + TClonesArray* cbm_sts_tracks_ {nullptr}; ///< non-owning pointer + TClonesArray* cbm_sts_match_ {nullptr}; ///< non-owning pointer + + bool is_write_kfinfo_ {true}; + bool is_reproduce_cbmkfpf_ {false}; + + int ipar_ {-1}; + int imf_ {-1}; + int icov_ {-1}; + int imc_pdg_ {-1}; + int imother_pdg_ {-1}; + int ipasscuts_ {-1}; + + ClassDef(CbmStsTracksConverter, 1) +}; + +#endif // ANALYSIS_TREE_STSTRACKSCONVERTER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9819cea643e1292cbcb8dda5b0741d04de4291d5 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx @@ -0,0 +1,174 @@ +#include <cassert> +#include <vector> + +#include "TClonesArray.h" + +#include <FairMCPoint.h> +#include <FairRootManager.h> + +#include "AnalysisTree/Matching.hpp" + +#include <CbmGlobalTrack.h> +#include <CbmMCTrack.h> +#include <CbmTofHit.h> +#include <CbmTrackMatchNew.h> + +#include "CbmTofHitsConverter.h" + +ClassImp(CbmTofHitsConverter) + + void CbmTofHitsConverter::Init(std::map<std::string, void*>&) { + + assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ + && out_tree_); + auto* ioman = FairRootManager::Instance(); + + cbm_tof_hits_ = (TClonesArray*) ioman->GetObject("TofHit"); + cbm_global_tracks_ = (TClonesArray*) ioman->GetObject("GlobalTrack"); + // cbm_tof_match_ = (TClonesArray*) ioman->GetObject("TofHitMatch"); + // cbm_tof_points_ = (TClonesArray*) ioman->GetObject("TofPoint"); + // cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); + // cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch"); + + AnalysisTree::BranchConfig tof_branch("TofHits", AnalysisTree::DetType::kHit); + tof_branch.AddField<float>("mass2"); + tof_branch.AddField<float>("l"); + tof_branch.AddField<float>("t"); + tof_branch.AddField<float>("qp_tof"); + tof_branch.AddFields<float>({"dx", "dy", "dz"}); + tof_branch.AddField<int>("mc_pdg"); + tof_branch.AddField<bool>("is_correct_match"); + out_config_->AddBranchConfig(std::move(tof_branch)); + tof_hits_ = new AnalysisTree::HitDetector(out_config_->GetLastId()); + + vtx_tracks_2_tof_ = new AnalysisTree::Matching( + out_config_->GetBranchConfig(match_to_).GetId(), out_config_->GetLastId()); + out_config_->AddMatch(vtx_tracks_2_tof_); + + out_tree_->Branch( + out_branch_.c_str(), "AnalysisTree::HitDetector", &tof_hits_); + out_tree_->Branch((match_to_ + "2" + out_branch_).c_str(), + "AnalysisTree::Matching", + &vtx_tracks_2_tof_); +} + +void CbmTofHitsConverter::ExtrapolateStraightLine(FairTrackParam* params, + float z) { + + const Float_t Tx = params->GetTx(); + const Float_t Ty = params->GetTy(); + const Float_t old_z = params->GetZ(); + const Float_t dz = z - old_z; + + const Float_t x = params->GetX() + Tx * dz; + const Float_t y = params->GetY() + Ty * dz; + + params->SetPosition({x, y, z}); +} + +void CbmTofHitsConverter::FillTofHits() { + assert(cbm_tof_hits_); + tof_hits_->ClearChannels(); + vtx_tracks_2_tof_->Clear(); + + const int i_mass2 = + out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("mass2"); + const int i_qp = + out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("qp_tof"); + const int i_dx = + out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("dx"); + const int i_t = + out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("t"); + const int i_l = + out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("l"); + // const int i_is_correct = + // out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("is_correct_match"); + // const int i_pdg = + // out_config_->GetBranchConfig(tof_hits_->GetId()).GetFieldId("mc_pdg"); + + const auto it = indexes_map_->find(match_to_); + if (it == indexes_map_->end()) { + throw std::runtime_error(match_to_ + + " is not found to match with TOF hits"); + } + auto rec_tracks_map = it->second; + + tof_hits_->Reserve(cbm_global_tracks_->GetEntries()); + + for (Int_t igt = 0; igt < cbm_global_tracks_->GetEntries(); igt++) { + const auto* globalTrack = + static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(igt)); + const Int_t tofHitIndex = globalTrack->GetTofHitIndex(); + if (tofHitIndex < 0) continue; + + const auto* tofHit = + static_cast<const CbmTofHit*>(cbm_tof_hits_->At(tofHitIndex)); + + FairTrackParam param_last = *(globalTrack->GetParamLast()); + TVector3 p_tof; + param_last.Momentum(p_tof); + + const Float_t p = p_tof.Mag(); + const Int_t q = param_last.GetQp() > 0 ? 1 : -1; + const Float_t l = + globalTrack->GetLength(); // l is calculated by global tracking + const Float_t time = tofHit->GetTime(); + const Float_t beta = l / (time * 29.9792458); + const Float_t m2 = p * p * (1. / (beta * beta) - 1.); + + const Float_t hitX = tofHit->GetX(); + const Float_t hitY = tofHit->GetY(); + const Float_t hitZ = tofHit->GetZ(); + + ExtrapolateStraightLine(¶m_last, tofHit->GetZ()); + + auto* Hit = tof_hits_->AddChannel(); + Hit->Init(out_config_->GetBranchConfig(tof_hits_->GetId())); + Hit->SetPosition(hitX, hitY, hitZ); + Hit->SetSignal(time); + Hit->SetField(m2, i_mass2); + Hit->SetField(float(q * p_tof.Mag()), i_qp); + Hit->SetField(float(param_last.GetX() - hitX), i_dx); + Hit->SetField(float(param_last.GetY() - hitY), i_dx + 1); + Hit->SetField(float(param_last.GetZ() - hitZ), i_dx + 2); + Hit->SetField(l, i_l); + Hit->SetField(time, i_t); + + if (rec_tracks_map.empty()) { return; } + const Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); + if (rec_tracks_map.find(stsTrackIndex) != rec_tracks_map.end()) { + vtx_tracks_2_tof_->AddMatch(rec_tracks_map.find(stsTrackIndex)->second, + Hit->GetId()); + } + + // const auto* tofMatch = + // dynamic_cast<CbmMatch*>(cbm_tof_match_->At(tofHitIndex)); if(tofMatch + // != nullptr && tofMatch->GetNofLinks()>0) { + // + // const auto* tofPoint = dynamic_cast<FairMCPoint*>( + // cbm_tof_points_->At(tofMatch->GetMatchedLink().GetIndex()) ); + // + // Int_t itofMC = (tofPoint ? tofPoint->GetTrackID() : -1 ); + // if(itofMC >= 0){ + // const auto* mc_track = dynamic_cast<const + // CbmMCTrack*>(cbm_mc_tracks_->At(itofMC)); + // Hit->SetField(mc_track->GetPdgCode(), i_pdg); + // + // const Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); + // if(stsTrackIndex<0) return; + // + // auto* match = (CbmTrackMatchNew*) cbm_sts_match_->At(stsTrackIndex); + // if(match == nullptr || match->GetNofLinks() == 0) continue; + // const int mc_id_sts = match->GetMatchedLink().GetIndex(); + // Hit->SetField( bool(mc_id_sts == itofMC), i_is_correct); + // } + // } + } +} + +void CbmTofHitsConverter::Exec() { FillTofHits(); } + +CbmTofHitsConverter::~CbmTofHitsConverter() { + delete tof_hits_; + delete vtx_tracks_2_tof_; +}; diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..990102074b51d67f84ede97113b9ea58dbdb6f57 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h @@ -0,0 +1,47 @@ +#ifndef ANALYSIS_TREE_TOFHITSCONVERTER_H +#define ANALYSIS_TREE_TOFHITSCONVERTER_H + +#include "CbmConverterTask.h" + +#include "AnalysisTree/Detector.hpp" + +class TClonesArray; +class FairTrackParam; + +namespace AnalysisTree { + class Matching; +} + +class CbmTofHitsConverter final : public CbmConverterTask { +public: + explicit CbmTofHitsConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) { + in_branches_.emplace_back("TofHits"); + }; + + ~CbmTofHitsConverter() final; + ; + + void Init(std::map<std::string, void*>&) final; + void Exec() final; + void Finish() final {}; + +private: + void FillTofHits(); + static void ExtrapolateStraightLine(FairTrackParam* params, float z); + + TClonesArray* cbm_global_tracks_ {nullptr}; + TClonesArray* cbm_tof_hits_ {nullptr}; + // TClonesArray *cbm_tof_points_{nullptr}; + // TClonesArray *cbm_tof_match_{nullptr}; + // TClonesArray* cbm_mc_tracks_{nullptr}; + // TClonesArray* cbm_sts_match_{nullptr}; + + AnalysisTree::HitDetector* tof_hits_ {nullptr}; + AnalysisTree::Matching* vtx_tracks_2_tof_ {nullptr}; + + ClassDef(CbmTofHitsConverter, 1) +}; + +#endif // ANALYSIS_TREE_TOFHITSCONVERTER_H diff --git a/analysis/common/analysis_tree_converter/README.md b/analysis/common/analysis_tree_converter/README.md new file mode 100644 index 0000000000000000000000000000000000000000..907e0ec47f3b271a06537d7ea217c7c28c4357ff --- /dev/null +++ b/analysis/common/analysis_tree_converter/README.md @@ -0,0 +1,15 @@ +# Converter from CBMROOT to AnalysisTree format + +## AnalysisTree format + +https://github.com/HeavyIonAnalysis/AnalysisTree + +## Usage + +Macro location: + +macro/analysis/common/analysis_tree_converter/run_treemaker.C + +## Known problems + + * no TRD, RICh, MuCh information \ No newline at end of file diff --git a/external/.gitignore b/external/.gitignore index 2e92e649953c8e312036737f561d4b637f09f8b0..3fc963cf47ce80ef2bd9e803607b8d62060af100 100644 --- a/external/.gitignore +++ b/external/.gitignore @@ -1,4 +1,6 @@ DataTree +AnalysisTree +AnalysisTreeQA DataTreeQA KFParticle NicaFemto diff --git a/external/CMakeLists.txt b/external/CMakeLists.txt index d567fd5d097ec782334433022e10192cd37344f1..e92f1b7f6dfce6be597e6dff27e807f31b9eedc7 100644 --- a/external/CMakeLists.txt +++ b/external/CMakeLists.txt @@ -28,6 +28,8 @@ Include(InstallKFParticle.cmake) Include(InstallDataTree.cmake) Include(InstallDataTreeQA.cmake) Include(InstallNicaFemto.cmake) +Include(InstallAnalysisTree.cmake) +Include(InstallAnalysisTreeQA.cmake) Include(InstallParameter.cmake) Include(InstallInput.cmake) diff --git a/external/InstallAnalysisTree.cmake b/external/InstallAnalysisTree.cmake new file mode 100644 index 0000000000000000000000000000000000000000..3179a253b610136b6a25144cce33b1cbdd370ea1 --- /dev/null +++ b/external/InstallAnalysisTree.cmake @@ -0,0 +1,65 @@ +set(ANALYSISTREE_VERSION v1.0) + +set(ANALYSISTREE_SRC_URL "https://github.com/HeavyIonAnalysis/AnalysisTree.git") +set(ANALYSISTREE_DESTDIR "${CMAKE_BINARY_DIR}/external/ANALYSISTREE-prefix") + +list(APPEND ANALYSISTREE_LIBNAME "AnalysisTreeBase" "AnalysisTreeInfra") + +#set(ANALYSISTREE_LIBNAME "${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTreeBase${CMAKE_SHARED_LIBRARY_SUFFIX}" "${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTreeInfra${CMAKE_SHARED_LIBRARY_SUFFIX}") + +download_project_if_needed(PROJECT AnalysisTree_source + GIT_REPOSITORY ${ANALYSISTREE_SRC_URL} + GIT_TAG ${ANALYSISTREE_VERSION} + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/AnalysisTree + TEST_FILE CMakeLists.txt + ) + +If(ProjectUpdated) + File(REMOVE_RECURSE ${ANALYSISTREE_DESTDIR}) + Message("AnalysisTree source directory was changed so build directory was deleted") +EndIf() + +ExternalProject_Add(ANALYSISTREE + BUILD_IN_SOURCE 0 + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/AnalysisTree + BUILD_BYPRODUCTS ${ANALYSISTREE_LIBRARY} + LOG_DOWNLOAD 1 LOG_CONFIGURE 1 LOG_BUILD 1 LOG_INSTALL 1 + CMAKE_ARGS -G ${CMAKE_GENERATOR} + -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + -DCMAKE_C_FLAGS=${CMAKE_C_FLAGS} + -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} + -DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS} + -DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR} + -DCMAKE_CXX_STANDARD=11 + -DROOTSYS=${SIMPATH} + -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON + INSTALL_COMMAND ${CMAKE_COMMAND} --build . --target install + ) + +add_library(AnalysisTreeBase SHARED IMPORTED) +set_target_properties(AnalysisTreeBase PROPERTIES IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/lib) +#add_dependencies(AnalysisTreeBase ANALYSISTREE) + +add_library(AnalysisTreeInfra SHARED IMPORTED) +set_target_properties(AnalysisTreeInfra PROPERTIES IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/lib) +#add_dependencies(AnalysisTreeInfra ANALYSISTREE) + +set(AnalysisTree_LIB_DIR ${CMAKE_BINARY_DIR}/lib) +set(AnalysisTree_LIBRARIES AnalysisTreeBase AnalysisTreeInfra) +set(AnalysisTree_INCLUDE_DIR "${CMAKE_BINARY_DIR}/include") +set(AnalysisTree_FOUND TRUE) + +foreach(LIB_NAME ${ANALYSISTREE_LIBNAME}) + Install( FILES + "${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${LIB_NAME}${CMAKE_SHARED_LIBRARY_SUFFIX}" + "${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${LIB_NAME}_rdict.pcm" + "${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${LIB_NAME}.rootmap" + DESTINATION lib) +endforeach(LIB_NAME) + + +#Install(FILES ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTree.rootmap +# ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTree_rdict.pcm +# DESTINATION lib +# ) diff --git a/external/InstallAnalysisTreeQA.cmake b/external/InstallAnalysisTreeQA.cmake new file mode 100644 index 0000000000000000000000000000000000000000..17a44edeb2a1a03416a91caf35c636e60002859e --- /dev/null +++ b/external/InstallAnalysisTreeQA.cmake @@ -0,0 +1,51 @@ +set(ANALYSISTREEQA_VERSION v1.0) + +set(ANALYSISTREEQA_SRC_URL "https://github.com/HeavyIonAnalysis/AnalysisTreeQA.git") +set(ANALYSISTREEQA_DESTDIR "${CMAKE_BINARY_DIR}/external/ANALYSISTREEQA-prefix") +set(ANALYSISTREEQA_LIBNAME "${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTreeQA${CMAKE_STATIC_LIBRARY_SUFFIX}") + +download_project_if_needed(PROJECT AnalysisTreeQA_source + GIT_REPOSITORY ${ANALYSISTREEQA_SRC_URL} + GIT_TAG ${ANALYSISTREEQA_VERSION} + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/AnalysisTreeQA + TEST_FILE CMakeLists.txt + ) + +If(ProjectUpdated) + File(REMOVE_RECURSE ${ANALYSISTREEQA_DESTDIR}) + Message("AnalysisTreeQA source directory was changed so build directory was deleted") +EndIf() + +ExternalProject_Add(ANALYSISTREEQA + BUILD_IN_SOURCE 0 + SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/AnalysisTreeQA + BUILD_BYPRODUCTS ${ANALYSISTREEQA_LIBRARY} + LOG_DOWNLOAD 1 LOG_CONFIGURE 1 LOG_BUILD 1 LOG_INSTALL 1 + CMAKE_ARGS -G ${CMAKE_GENERATOR} + -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + -DCMAKE_C_FLAGS=${CMAKE_C_FLAGS} + -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} + -DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS} + -DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR} + -DCMAKE_CXX_STANDARD=11 + -DROOTSYS=${SIMPATH} + -DUSEBOOST=TRUE + -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON + INSTALL_COMMAND ${CMAKE_COMMAND} --build . --target install + ) + +add_library(AnalysisTreeQA SHARED IMPORTED) +set_target_properties(AnalysisTreeQA PROPERTIES IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/lib) +add_dependencies(AnalysisTreeQA ANALYSISTREE) + +set(AnalysisTreeQA_LIB_DIR ${CMAKE_BINARY_DIR}/lib) +set(AnalysisTreeQA_LIBRARIES AnalysisTreeQA) +set(AnalysisTreeQA_INCLUDE_DIR "${CMAKE_BINARY_DIR}/include") +set(AnalysisTreeQA_FOUND TRUE) + +Install(FILES ${CMAKE_BINARY_DIR}/lib/${ANALYSISTREEQA_LIBNAME} + ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTreeQA.rootmap + ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}AnalysisTreeQA_rdict.pcm + DESTINATION lib + ) diff --git a/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh b/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh new file mode 100755 index 0000000000000000000000000000000000000000..733836b97d2b7fbb541a6cf4ad926053862548e9 --- /dev/null +++ b/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh @@ -0,0 +1,39 @@ +#!/bin/bash +#SBATCH -J CbmAnalysisTree +#SBATCH -o out/%j.out.log +#SBATCH -e error/%j.err.log +#SBATCH --time=2:00:00 +#SBATCH --array=1-50 + +N_EVENTS=100 +SETUP="sis100_electron" + +CBMROOT_DIR=/lustre/cbm/users/klochkov/soft/cbmroot/c2f_fork/ + +source $CBMROOT_DIR/build/config.sh + +INDEX=$SLURM_ARRAY_TASK_ID + +# INPUT_DIR=/lustre/cbm/users/ogolosov/mc/cbmsim/test_prod_AT/dcmqgsm_smm_pluto/auau/12agev/mbias/psd44_hole20_pipe0/TGeant3/ +OUTPUT_DIR=/lustre/nyx/cbm/users/klochkov/cbm/test_prod_AT/dcmqgsm_smm_pluto/auau/12agev/mbias/psd44_hole20_pipe0/TGeant3/$INDEX +INPUT_FILE=/lustre/cbm/users/ogolosov/mc/generators/urqmd/v3.4/auau/pbeam12agev_eos0/mbias/root/urqmd_$INDEX.root + +DATA_SET=$INDEX +MACRO_DIR=$CBMROOT_DIR/macro/analysis/common/analysis_tree_converter/ + +mkdir -p $OUTPUT_DIR +cd $OUTPUT_DIR + +INPUT_DIR=$INPUT_DIR/$DATA_SET/ + +cp $CBMROOT_DIR/macro/run/run_transport.C ./ +cp $CBMROOT_DIR/macro/run/run_digi.C ./ +cp $CBMROOT_DIR/macro/run/run_reco_event.C ./ +cp $CBMROOT_DIR/macro/include/.rootrc ./ + +cp $MACRO_DIR/run_treemaker.C ./ + +root -l -q -b "run_transport.C($N_EVENTS, \"$SETUP\", \"$DATA_SET\", \"$INPUT_FILE\")" >& log_tr_$INDEX.txt +root -l -q -b "run_digi.C($N_EVENTS, \"$DATA_SET\", 1e7, 1e4, kTRUE)" >& log_digi_$INDEX.txt +root -l -q -b "run_reco_event.C($N_EVENTS, \"$DATA_SET\", \"$SETUP\")" >& log_reco_$INDEX.txt +root -l -q -b "run_treemaker.C(\"\", \"$DATA_SET\")" >& log_$INDEX.txt diff --git a/macro/analysis/common/analysis_tree_converter/batch/run.sh b/macro/analysis/common/analysis_tree_converter/batch/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..f794eec5bed3ea6ed70fd582a8bdb5c5cdb8a2e6 --- /dev/null +++ b/macro/analysis/common/analysis_tree_converter/batch/run.sh @@ -0,0 +1,6 @@ +LOGDIR=/lustre/cbm/users/$USER/log +mkdir -p $LOGDIR +mkdir -p $LOGDIR/out +mkdir -p $LOGDIR/error + +sbatch --partition main -D "$LOGDIR" --export=ALL batch_run.sh diff --git a/macro/analysis/common/analysis_tree_converter/run_treemaker.C b/macro/analysis/common/analysis_tree_converter/run_treemaker.C new file mode 100644 index 0000000000000000000000000000000000000000..328e1dca8a7bb33870d1d7a20d2c2045ebc1aa8f --- /dev/null +++ b/macro/analysis/common/analysis_tree_converter/run_treemaker.C @@ -0,0 +1,171 @@ + +void run_treemaker( + TString inputDir = + "/home/vklochkov/Data/cbm/apr20_fr_18.2.1_fs_jun19p1/cbmroot_new/data/", + TString dataSet = "test", + TString setupName = "sis100_electron") { + const std::string system = "Au+Au"; + const float beam_mom = 12.; + + // --- Logger settings ---------------------------------------------------- + const TString logLevel = "INFO"; + const TString logVerbosity = "LOW"; + // ------------------------------------------------------------------------ + + // ----- Environment -------------------------------------------------- + const TString myName = + "run_treemaker"; // this macro's name for screen output + const TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + // ----- In- and output file names ------------------------------------ + TString traFile = inputDir + dataSet + ".tra.root"; + TString rawFile = inputDir + dataSet + ".event.raw.root"; + TString recFile = inputDir + dataSet + ".rec.root"; + TString geoFile = inputDir + dataSet + ".geo.root"; + TString parFile = inputDir + dataSet + ".par.root"; + const std::string outFile = + dataSet.Data() + std::string(".analysistree.root"); + // ------------------------------------------------------------------------ + + // ----- Remove old CTest runtime dependency file ---------------------- + const TString dataDir = gSystem->DirName(dataSet); + const TString dataName = gSystem->BaseName(dataSet); + const TString testName = ("run_treemaker"); + // ------------------------------------------------------------------------ + + // ----- Load the geometry setup ------------------------------------- + std::cout << std::endl; + const TString setupFile = + srcDir + "/geometry/setup/setup_" + setupName + ".C"; + const TString setupFunct = "setup_" + setupName + "()"; + + std::cout << "-I- " << myName << ": Loading macro " << setupFile << std::endl; + + gROOT->LoadMacro(setupFile); + gROOT->ProcessLine(setupFunct); + CbmSetup* setup = CbmSetup::Instance(); + + // ----- Timer -------------------------------------------------------- + TStopwatch timer; + timer.Start(); + // ------------------------------------------------------------------------ + + TString geoTag; + auto* parFileList = new TList(); + + std::cout << "-I- " << myName << ": Using raw file " << rawFile << std::endl; + std::cout << "-I- " << myName << ": Using parameter file " << parFile + << std::endl; + std::cout << "-I- " << myName << ": Using reco file " << recFile << std::endl; + + // ----- Reconstruction run ------------------------------------------- + auto* run = new FairRunAna(); + auto* inputSource = new FairFileSource(recFile); + inputSource->AddFriend(traFile); + inputSource->AddFriend(rawFile); + run->SetSource(inputSource); + run->SetOutputFile(outFile.c_str()); + run->SetGenerateRunInfo(kTRUE); + // ------------------------------------------------------------------------ + + // ----- Mc Data Manager ------------------------------------------------ + auto* mcManager = new CbmMCDataManager("MCManager", 1); + mcManager->AddFile(traFile); + run->AddTask(mcManager); + // ------------------------------------------------------------------------ + + // --- STS track matching ---------------------------------------------- + auto* matchTask = new CbmMatchRecoToMC(); + run->AddTask(matchTask); + // ------------------------------------------------------------------------ + + auto* KF = new CbmKF(); + run->AddTask(KF); + // needed for tracks extrapolation + auto* l1 = new CbmL1("CbmL1", 1, 3); + if (setup->IsActive(ECbmModuleId::kMvd)) { + setup->GetGeoTag(ECbmModuleId::kMvd, geoTag); + const TString mvdMatBudgetFileName = + srcDir + "/parameters/mvd/mvd_matbudget_" + geoTag + ".root"; + l1->SetMvdMaterialBudgetFileName(mvdMatBudgetFileName.Data()); + } + if (setup->IsActive(ECbmModuleId::kSts)) { + setup->GetGeoTag(ECbmModuleId::kSts, geoTag); + const TString stsMatBudgetFileName = + srcDir + "/parameters/sts/sts_matbudget_" + geoTag + ".root"; + l1->SetStsMaterialBudgetFileName(stsMatBudgetFileName.Data()); + } + run->AddTask(l1); + + // AnalysisTree converter + auto* man = new CbmConverterManager(); + man->SetSystem(system); + man->SetBeamMomentum(beam_mom); + + man->SetOutFileName(outFile); + man->SetOutTreeName("aTree"); + + man->AddTask(new CbmSimEventHeaderConverter("SimEventHeader")); + man->AddTask(new CbmRecEventHeaderConverter("RecEventHeader")); + man->AddTask(new CbmSimTracksConverter("SimParticles")); + + CbmStsTracksConverter* taskCbmStsTracksConverter = + new CbmStsTracksConverter("VtxTracks", "SimParticles"); + taskCbmStsTracksConverter->SetIsWriteKFInfo(); + taskCbmStsTracksConverter->SetIsReproduceCbmKFPF(); + man->AddTask(taskCbmStsTracksConverter); + + man->AddTask(new CbmTofHitsConverter("TofHits", "VtxTracks")); + man->AddTask(new CbmPsdModulesConverter("PsdModules")); + + run->AddTask(man); + + // ----- Parameter database -------------------------------------------- + FairRuntimeDb* rtdb = run->GetRuntimeDb(); + auto* parIo1 = new FairParRootFileIo(); + auto* parIo2 = new FairParAsciiFileIo(); + parIo1->open(parFile.Data()); + parIo2->open(parFileList, "in"); + rtdb->setFirstInput(parIo1); + rtdb->setSecondInput(parIo2); + rtdb->setOutput(parIo1); + rtdb->saveOutput(); + // ------------------------------------------------------------------------ + + // ----- Intialise and run -------------------------------------------- + run->Init(); + + std::cout << "Starting run" << std::endl; + run->Run(); + // ------------------------------------------------------------------------ + + timer.Stop(); + const Double_t rtime = timer.RealTime(); + const Double_t ctime = timer.CpuTime(); + std::cout << "Macro finished succesfully." << std::endl; + std::cout << "Output file is " << outFile << std::endl; + std::cout << "Parameter file is " << parFile << std::endl; + printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime); + + // ----- CTest resource monitoring ------------------------------------ + FairSystemInfo sysInfo; + const Float_t maxMemory = sysInfo.GetMaxMemory(); + std::cout << R"(<DartMeasurement name="MaxMemory" type="numeric/double">)"; + std::cout << maxMemory; + std::cout << "</DartMeasurement>" << std::endl; + std::cout << R"(<DartMeasurement name="WallTime" type="numeric/double">)"; + std::cout << rtime; + std::cout << "</DartMeasurement>" << std::endl; + const Float_t cpuUsage = ctime / rtime; + std::cout << R"(<DartMeasurement name="CpuLoad" type="numeric/double">)"; + std::cout << cpuUsage; + std::cout << "</DartMeasurement>" << std::endl; + // ------------------------------------------------------------------------ + + // ----- Finish ------------------------------------------------------- + std::cout << " Test passed" << std::endl; + std::cout << " All ok " << std::endl; + // Generate_CTest_Dependency_File(depFile); + // ------------------------------------------------------------------------ +}