diff --git a/analysis/common/analysis_tree_converter/CMakeLists.txt b/analysis/common/analysis_tree_converter/CMakeLists.txt index 6053191de13420942de89d5351e78368de424643..f918cbeec0d0616ad22a007520bf102f51d49de1 100644 --- a/analysis/common/analysis_tree_converter/CMakeLists.txt +++ b/analysis/common/analysis_tree_converter/CMakeLists.txt @@ -9,6 +9,8 @@ set(SRCS CbmSimTracksConverter.cxx CbmPsdModulesConverter.cxx CbmTofHitsConverter.cxx + CbmTrdTracksConverter.cxx + CbmRichRingsConverter.cxx ) @@ -23,6 +25,7 @@ Set(INCLUDE_DIRECTORIES ${CBMROOT_SOURCE_DIR}/core/data/sts ${CBMROOT_SOURCE_DIR}/core/data/tof ${CBMROOT_SOURCE_DIR}/core/data/psd + ${CBMROOT_SOURCE_DIR}/core/data/trd ${CBMROOT_SOURCE_DIR}/sim/transport/generators ${CBMROOT_SOURCE_DIR}/sim/transport/steer diff --git a/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h index 0c8265b512cac8d7a2faf7f9f3d05fe5b6756c7f..f7d0fbb956b198abed1411bcc81e0f19279d7ac8 100644 --- a/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h +++ b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h @@ -12,9 +12,9 @@ #pragma link C++ class CbmStsTracksConverter + ; #pragma link C++ class CbmTofHitsConverter + ; #pragma link C++ class CbmPsdModulesConverter + ; +#pragma link C++ class CbmTrdTracksConverter + ; +#pragma link C++ class CbmRichRingsConverter + ; //#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 index 067a372d12c2a4c345b5e1be479f911c0d6ca12e..0a3b792cda899bbe1c10e58a9fdc143443014ace 100644 --- a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx +++ b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx @@ -12,69 +12,46 @@ 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"); - + task_manager_->Init(); 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(); + task_manager_->FillOutput(); } void CbmConverterManager::Finish() { - TDirectory* old_dir = gDirectory; + TDirectory* curr = gDirectory; // TODO check why this is needed + TFile* currentFile = gFile; - out_file_->cd(); - for (auto* task : tasks_) { - task->Finish(); - } - out_tree_->Write(); - out_file_->Close(); + task_manager_->GetChain()->Write(); + + task_manager_->Finish(); - old_dir->cd(); + gFile = currentFile; + gDirectory = curr; } void CbmConverterManager::FillDataHeader() { + // Force user to write data info //TODO is there a way to read it from a file automatically? + assert(!system_.empty() && beam_mom_); - 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; + auto* data_header = new AnalysisTree::DataHeader(); std::cout << "ReadDataHeader" << std::endl; - data_header_->SetSystem(system_); - data_header_->SetBeamMomentum(beam_mom_); + data_header->SetSystem(system_); + data_header->SetBeamMomentum(beam_mom_); - auto& psd_mod_pos = data_header_->AddDetector(); + auto& psd_mod_pos = data_header->AddDetector(); const int psd_node_id = 6; const char* module_name_prefix = "module"; @@ -110,18 +87,6 @@ void CbmConverterManager::FillDataHeader() { } } - 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_; + task_manager_->SetOutputDataHeader(data_header); } +CbmConverterManager::~CbmConverterManager() = default; diff --git a/analysis/common/analysis_tree_converter/CbmConverterManager.h b/analysis/common/analysis_tree_converter/CbmConverterManager.h index 0eb217ff98dffb9076a8a27a73f84e87057739d6..98e15edc0d4cd3468efbd34688f006de06a62203 100644 --- a/analysis/common/analysis_tree_converter/CbmConverterManager.h +++ b/analysis/common/analysis_tree_converter/CbmConverterManager.h @@ -3,6 +3,8 @@ #include "FairTask.h" +#include "AnalysisTree/TaskManager.hpp" + namespace AnalysisTree { class Configuration; class DataHeader; @@ -20,26 +22,27 @@ public: 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 AddTask(CbmConverterTask* task) { + tasks_.emplace_back(task); + task_manager_->AddTask(reinterpret_cast<AnalysisTree::Task*>(task)); + } void SetSystem(const std::string& system) { system_ = system; } void SetBeamMomentum(float beam_mom) { beam_mom_ = beam_mom; } + void SetOutputName(std::string file, std::string tree = "aTree") { + task_manager_->SetOutputName(file, tree); + } + private: void FillDataHeader(); - TFile* out_file_ {nullptr}; - TTree* out_tree_ {nullptr}; - std::string out_file_name_ {""}; - std::string out_tree_name_ {""}; + AnalysisTree::TaskManager* task_manager_ { + AnalysisTree::TaskManager::GetInstance()}; - std::string system_ {""}; + 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>> diff --git a/analysis/common/analysis_tree_converter/CbmConverterTask.h b/analysis/common/analysis_tree_converter/CbmConverterTask.h index e37d0cf74282056953452ddde288f56dbeb54b66..d88e6302759236d35698ef71b09a6fd198ddfbde 100644 --- a/analysis/common/analysis_tree_converter/CbmConverterTask.h +++ b/analysis/common/analysis_tree_converter/CbmConverterTask.h @@ -1,11 +1,11 @@ #ifndef ANALYSIS_TREE_CONVERTERTASK_H_ #define ANALYSIS_TREE_CONVERTERTASK_H_ -#include "AnalysisTree/FillTask.hpp" +#include "AnalysisTree/Task.hpp" class FairRootManager; -class CbmConverterTask : public AnalysisTree::FillTask { +class CbmConverterTask : public AnalysisTree::Task { public: CbmConverterTask() = delete; @@ -17,7 +17,6 @@ public: ~CbmConverterTask() override = default; - // void SetIoMan(FairRootManager* man) { ioman_ = man; } const std::map<int, int>& GetOutIndexesMap() const { return out_indexes_map_; } @@ -26,13 +25,15 @@ public: indexes_map_ = indexes_map; } + const std::string& GetOutputBranchName() const { return out_branch_; } + protected: - // FairRootManager* ioman_{nullptr}; std::map<int, int> out_indexes_map_ {}; ///< CbmRoot to AnalysisTree indexes ///< map for output branch + std::string out_branch_ {}; std::map<std::string, std::map<int, int>>* indexes_map_ {}; ///< CbmRoot to AnalysisTree indexes map for branches - ///< from other tasks + ///< from other tasks std::string match_to_ {}; ///< AT branch to match }; diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx index 02dbf1f8eeb67434f2766fc06ec762a8711c165f..0b1e02539394f73fdc49c7ab19c169e049f9a4a5 100644 --- a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx @@ -1,3 +1,4 @@ +#include <AnalysisTree/TaskManager.hpp> #include <cassert> #include <vector> @@ -13,21 +14,20 @@ ClassImp(CbmPsdModulesConverter) - void CbmPsdModulesConverter::Init(std::map<std::string, void*>&) { - assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ - && out_tree_); + void CbmPsdModulesConverter::Init() { + assert(!out_branch_.empty()); 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_); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, psd_modules_, psd_branch); } + void CbmPsdModulesConverter::Exec() { assert(cbm_psd_hits_); psd_modules_->ClearChannels(); @@ -35,22 +35,23 @@ void CbmPsdModulesConverter::Exec() { 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; + auto* data_header = AnalysisTree::TaskManager::GetInstance()->GetDataHeader(); + auto* config = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = config->GetBranchConfig(out_branch_); + const int n_psd_modules = - data_header_->GetModulePositions(0).GetNumberOfChannels(); + 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); + auto& module = psd_modules_->AddChannel(branch); + 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); + auto& module = psd_modules_->Channel(hit->GetModuleID() - 1); module.SetNumber(i + 1); module.SetSignal(hit->GetEdep()); psd_energy += hit->GetEdep(); @@ -60,7 +61,6 @@ void CbmPsdModulesConverter::Exec() { } void CbmPsdModulesConverter::Finish() { - delete psd_modules_; - delete cbm_psd_hits_; + // delete psd_modules_; } CbmPsdModulesConverter::~CbmPsdModulesConverter() { delete psd_modules_; }; diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h index a96028c6756d0830e3e07737bf64f9fd557a0b05..b57e72c9c593b3a0b674ade4e5d9ccb8802b8527 100644 --- a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h +++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.h @@ -12,12 +12,11 @@ public: 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 Init() final; void Exec() final; void Finish() final; diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx index e5bdcc73a94347829be6252acc894602f7123db4..0d744b92e6594c9411d2e8b42e1f34d6718b159d 100644 --- a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx @@ -3,51 +3,57 @@ #include "FairRootManager.h" #include "TClonesArray.h" #include "cassert" +#include "iostream" +#include <AnalysisTree/TaskManager.hpp> #include "CbmRecEventHeaderConverter.h" ClassImp(CbmRecEventHeaderConverter) - void CbmRecEventHeaderConverter::Init(std::map<std::string, void*>&) { - assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ - && out_tree_); + void CbmRecEventHeaderConverter::Init() { + assert(!out_branch_.empty()); 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()); + cbm_header_ = (FairMCEventHeader*) ioman->GetObject("MCEventHeader."); + cbm_prim_vertex_ = (CbmVertex*) ioman->GetObject("PrimaryVertex."); // ***** 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()); + RecEventHeaderBranch.AddField<float>("vtx_chi2", "primiry vertex fit chi2"); + RecEventHeaderBranch.AddField<float>("Epsd", + "GeV, full energy deposit in PSD"); + RecEventHeaderBranch.AddField<int>("M", "total multiplicity in STS"); + RecEventHeaderBranch.AddField<int>("evt_id", "event identifier"); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, rec_event_header_, RecEventHeaderBranch); 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()); + + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); + 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")); + rec_event_header_->SetField(float(-999.), branch.GetFieldId("vtx_chi2")); return; } + + // std::cout << cbm_prim_vertex_->GetX() << " " << cbm_prim_vertex_->GetY() << " " << cbm_prim_vertex_->GetZ() << std::endl; + 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")); + branch.GetFieldId("evt_id")); rec_event_header_->SetField( float(cbm_prim_vertex_->GetChi2() / cbm_prim_vertex_->GetNDF()), - conf.GetFieldId("vtx_chi2")); + branch.GetFieldId("vtx_chi2")); } diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h index 1208bf7361f0f2123f47fb00daf223967da0dd06..ef6c163930f4944ae68e11e6a4a2a6be4923a5df 100644 --- a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h @@ -11,13 +11,10 @@ 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."); - }; + : CbmConverterTask(std::move(out_branch_name)) {}; ~CbmRecEventHeaderConverter() final = default; - void Init(std::map<std::string, void*>&) final; + void Init() final; void Exec() final; void Finish() final { delete rec_event_header_; }; diff --git a/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cc7a411ffca24f84cffd9f3814b30b71c1e5f352 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx @@ -0,0 +1,114 @@ +#include <cassert> +#include <vector> + +#include "TClonesArray.h" + +#include <FairMCPoint.h> +#include <FairRootManager.h> + +#include "AnalysisTree/Matching.hpp" + +#include <AnalysisTree/TaskManager.hpp> +#include <CbmGlobalTrack.h> +#include <CbmTrackMatchNew.h> +#include <rich/CbmRichRing.h> + +#include "CbmRichRingsConverter.h" + +void CbmRichRingsConverter::Init() { + + assert(!out_branch_.empty()); + auto* ioman = FairRootManager::Instance(); + + cbm_rich_rings_ = (TClonesArray*) ioman->GetObject("RichRing"); + cbm_global_tracks_ = (TClonesArray*) ioman->GetObject("GlobalTrack"); + + AnalysisTree::BranchConfig rich_branch(out_branch_, + AnalysisTree::DetType::kHit); + rich_branch.AddField<float>("radius"); + rich_branch.AddFields<int>({"n_hits", "n_hits_on_ring"}); + rich_branch.AddFields<float>({"axis_a", "axis_b"}); + rich_branch.AddFields<float>({"center_x", "center_y"}); + rich_branch.AddField<float>("chi2_ov_ndf", "chi2/ndf ring fit"); + rich_branch.AddField<float>("phi_ellipse", "phi rotation angle of ellipse"); + rich_branch.AddField<float>("radial_pos", "sqrt(x**2+abs(y-110)**2)"); + rich_branch.AddField<float>("radial_angle", + "(0||1||2)*pi +- atan( abs((+-100-y)/-x) )"); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, rich_rings_, rich_branch); + man->AddMatching(match_to_, out_branch_, vtx_tracks_2_rich_); +} + +void CbmRichRingsConverter::FillRichRings() { + + assert(cbm_rich_rings_); + rich_rings_->ClearChannels(); + vtx_tracks_2_rich_->Clear(); + + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); + + const int i_r = branch.GetFieldId("radius"); + const int i_n_hits = branch.GetFieldId("n_hits"); + const int i_axis = branch.GetFieldId("axis_a"); + const int i_center = branch.GetFieldId("center_x"); + const int i_chi2 = branch.GetFieldId("chi2_ov_ndf"); + const int i_radial_angle = branch.GetFieldId("radial_angle"); + const int i_radial_pos = branch.GetFieldId("radial_pos"); + const int i_phi_ellipse = branch.GetFieldId("phi_ellipse"); + + 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; + + rich_rings_->Reserve(cbm_global_tracks_->GetEntries()); + + for (Int_t igt = 0; igt < cbm_global_tracks_->GetEntries(); igt++) { + const auto* global_track = + static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(igt)); + + Int_t i_rich = global_track->GetRichRingIndex(); + if (i_rich < 0) continue; + auto rich_ring = static_cast<CbmRichRing*>(cbm_rich_rings_->At(i_rich)); + // auto rich_match = static_cast<CbmTrackMatchNew*>( cbm_rich_rings_->At(i_rich)); + // Int_t itrdMC = (rich_match ? rich_match->GetMatchedLink().GetIndex() : -1); + // richProj = static_cast<FairTrackParam*>(fRichProjection->At(i)); + + auto& ring = rich_rings_->AddChannel(branch); + + ring.SetPosition(rich_ring->GetCenterX(), rich_ring->GetCenterY(), 0.f); + + ring.SetField(int(rich_ring->GetNofHits()), i_n_hits); + ring.SetField(int(rich_ring->GetNofHitsOnRing()), i_n_hits + 1); + ring.SetField(float(rich_ring->GetAaxis()), i_axis); + ring.SetField(float(rich_ring->GetBaxis()), i_axis + 1); + ring.SetField(float(rich_ring->GetCenterX()), i_center); + ring.SetField(float(rich_ring->GetCenterY()), i_center + 1); + ring.SetField(float(rich_ring->GetRadius()), i_r); + ring.SetField(float(rich_ring->GetNDF() > 0. + ? rich_ring->GetChi2() / rich_ring->GetNDF() + : -999.), + i_chi2); + ring.SetField(float(rich_ring->GetRadialAngle()), i_radial_angle); + ring.SetField(float(rich_ring->GetRadialPosition()), i_radial_pos); + ring.SetField(float(rich_ring->GetPhi()), i_phi_ellipse); + + if (rec_tracks_map.empty()) { continue; } + const Int_t stsTrackIndex = global_track->GetStsTrackIndex(); + if (rec_tracks_map.find(stsTrackIndex) != rec_tracks_map.end()) { + vtx_tracks_2_rich_->AddMatch(rec_tracks_map.find(stsTrackIndex)->second, + ring.GetId()); + } + } +} + +void CbmRichRingsConverter::Exec() { FillRichRings(); } + +CbmRichRingsConverter::~CbmRichRingsConverter() { + delete rich_rings_; + delete vtx_tracks_2_rich_; +}; diff --git a/analysis/common/analysis_tree_converter/CbmRichRingsConverter.h b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..5f8a9b9e844112d241980278ee2a82b707df4220 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.h @@ -0,0 +1,39 @@ +#ifndef ANALYSIS_TREE_RICHRINGSCONVERTER_H_ +#define ANALYSIS_TREE_RICHRINGSCONVERTER_H_ + +#include "CbmConverterTask.h" + +#include "AnalysisTree/Detector.hpp" + +class TClonesArray; + +namespace AnalysisTree { + class Matching; +} + +class CbmRichRingsConverter final : public CbmConverterTask { +public: + explicit CbmRichRingsConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) {}; + + ~CbmRichRingsConverter() final; + + void Init() final; + void Exec() final; + void Finish() final {}; + +private: + void FillRichRings(); + + TClonesArray* cbm_global_tracks_ {nullptr}; + TClonesArray* cbm_rich_rings_ {nullptr}; + + AnalysisTree::HitDetector* rich_rings_ {nullptr}; + AnalysisTree::Matching* vtx_tracks_2_rich_ {nullptr}; + + ClassDef(CbmRichRingsConverter, 1) +}; + + +#endif //ANALYSIS_TREE_RICHRINGSCONVERTER_H_ diff --git a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx index 77f7bdd263550976143d8d17667f4702021b531b..46f1bb964f17c16fa60e9d2a38b6c6fb855ce42c 100644 --- a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.cxx @@ -2,29 +2,31 @@ #include "FairRootManager.h" #include "TClonesArray.h" #include "cassert" +#include "iostream" +#include <AnalysisTree/TaskManager.hpp> #include "CbmSimEventHeaderConverter.h" ClassImp(CbmSimEventHeaderConverter) - void CbmSimEventHeaderConverter::Init(std::map<std::string, void*>&) { - assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ - && out_tree_); + void CbmSimEventHeaderConverter::Init() { + assert(!out_branch_.empty()); auto* ioman = FairRootManager::Instance(); assert(ioman != nullptr); - cbm_header_ = (FairMCEventHeader*) ioman->GetObject(in_branches_[0].c_str()); + cbm_header_ = (FairMCEventHeader*) ioman->GetObject("MCEventHeader."); // ***** 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); + SimEventHeaderBranch.AddField<float>("psi_RP", "reaction plane orientation"); + SimEventHeaderBranch.AddField<float>("b", "impact parameter"); + SimEventHeaderBranch.AddField<float>("time", "ns, event time"); + SimEventHeaderBranch.AddField<int>("run_id", "run identifier"); + SimEventHeaderBranch.AddField<int>("event_id", "event identifier"); - out_tree_->Branch( - out_branch_.c_str(), "AnalysisTree::EventHeader", &sim_event_header_); + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, sim_event_header_, SimEventHeaderBranch); + sim_event_header_->Init(SimEventHeaderBranch); } void CbmSimEventHeaderConverter::Exec() { @@ -32,15 +34,20 @@ void CbmSimEventHeaderConverter::Exec() { throw std::runtime_error( "CbmSimEventHeaderConverter::Exec - ERROR! No fHeader!"); } + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); - TVector3 pos = { - cbm_header_->GetX(), cbm_header_->GetY(), cbm_header_->GetZ()}; + 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")); + + sim_event_header_->SetField(float(cbm_header_->GetRotZ()), + branch.GetFieldId("psi_RP")); + sim_event_header_->SetField(float(cbm_header_->GetB()), + branch.GetFieldId("b")); + sim_event_header_->SetField(int(cbm_header_->GetEventID()), + branch.GetFieldId("event_id")); + sim_event_header_->SetField(int(cbm_header_->GetRunID()), + branch.GetFieldId("run_id")); + sim_event_header_->SetField(float(cbm_header_->GetT()), + branch.GetFieldId("time")); } diff --git a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h index db81811f9f81bae9c2c6ec7f3219748eab1dcc3c..e99870f42a3257a3d2b7628dd04c5aba2503f400 100644 --- a/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h +++ b/analysis/common/analysis_tree_converter/CbmSimEventHeaderConverter.h @@ -11,12 +11,10 @@ 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."); - }; + : CbmConverterTask(std::move(out_branch_name)) {}; ~CbmSimEventHeaderConverter() final = default; - void Init(std::map<std::string, void*>&) final; + void Init() final; void Exec() final; void Finish() final { delete sim_event_header_; }; diff --git a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx index 20122644f0e6bf25c1487140b3e8cd5490d3c239..10c7e56bd8a80488447af7103a5b0fb37288837f 100644 --- a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.cxx @@ -1,6 +1,8 @@ #include <cassert> #include <vector> +#include <AnalysisTree/TaskManager.hpp> + #include "TClonesArray.h" #include "FairRootManager.h" @@ -10,58 +12,75 @@ ClassImp(CbmSimTracksConverter) - void CbmSimTracksConverter::Init(std::map<std::string, void*>&) { + void CbmSimTracksConverter::Init() { - assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ - && out_tree_); + assert(!out_branch_.empty()); 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_); + cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); + + AnalysisTree::BranchConfig sim_particles_branch( + out_branch_, AnalysisTree::DetType::kParticle); + sim_particles_branch.AddField<int>("mother_id", + "id of mother particle, -1 for primaries"); + sim_particles_branch.AddField<int>("cbmroot_id", + "track id in CbmRoot transport file"); + sim_particles_branch.AddField<int>("geant_process_id", ""); + sim_particles_branch.AddFields<int>( + {"n_hits_mvd", "n_hits_sts", "n_hits_trd"}, + "Number of hits in the detector"); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, sim_tracks_, sim_particles_branch); } void CbmSimTracksConverter::Exec() { assert(cbm_mc_tracks_); out_indexes_map_.clear(); - std::cout << "ReadMcTracks" << std::endl; + // std::cout << "ReadMcTracks" << std::endl; sim_tracks_->ClearChannels(); + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); + + const int nMcTracks = cbm_mc_tracks_->GetEntries(); + const int imother_id = branch.GetFieldId("mother_id"); + const int igeant_id = branch.GetFieldId("geant_process_id"); + const int in_hits = branch.GetFieldId("n_hits_mvd"); + const int icbm_id = branch.GetFieldId("cbmroot_id"); + + // std::cout << "cbm_mc_tracks_->GetEntries() = " << nMcTracks << std::endl; + // std::cout << imother_id << " " << igeant_id << " " << in_hits << std::endl; - 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())); + if (mctrack->GetPdgCode() == 50000050) { //Cherenkov + continue; + } + auto& track = sim_tracks_->AddChannel(branch); - out_indexes_map_.insert(std::make_pair(iMcTrack, Track->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())); + track.SetMomentum(mctrack->GetPx(), mctrack->GetPy(), mctrack->GetPz()); + track.SetMass(float(mctrack->GetMass())); + track.SetPid(int(mctrack->GetPdgCode())); + track.SetField(int(mctrack->GetGeantProcessId()), igeant_id); + track.SetField(int(mctrack->GetNPoints(ECbmModuleId::kMvd)), in_hits); + track.SetField(int(mctrack->GetNPoints(ECbmModuleId::kSts)), in_hits + 1); + track.SetField(int(mctrack->GetNPoints(ECbmModuleId::kTrd)), in_hits + 2); + track.SetField(int(mctrack->GetUniqueID()), icbm_id); if (mctrack->GetMotherId() == -1) { // mother id should < than track id, so we can do it here - Track->SetField(int(-1), imother_id); + 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); + track.SetField(int(-999), imother_id); + else { + track.SetField(int(p->second), imother_id); + } } } } diff --git a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h index dd3c92a4bd601db7308d93e02b1b9b56f15b3a01..a1e0676bdf73f1daaa6b38d297d000df541515b6 100644 --- a/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h +++ b/analysis/common/analysis_tree_converter/CbmSimTracksConverter.h @@ -11,13 +11,11 @@ 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"); - }; + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) {}; ~CbmSimTracksConverter() final; - void Init(std::map<std::string, void*>&) final; + void Init() final; void Exec() final; void Finish() final {}; diff --git a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx index 4d8f645197d39f866d403a8169197d4f3dd3052d..6aef8251ae07e08ae53f848fa27338ddf384a245 100644 --- a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx @@ -1,12 +1,4 @@ - -#include <cassert> -#include <cmath> - -#include "TClonesArray.h" - -#include "FairRootManager.h" - -#include "AnalysisTree/Matching.hpp" +#include "CbmStsTracksConverter.h" #include "CbmMCTrack.h" #include "CbmStsTrack.h" @@ -16,7 +8,16 @@ #include "L1Field.h" #include "ParticleFinder/CbmL1PFFitter.h" -#include "CbmStsTracksConverter.h" +#include "FairRootManager.h" + +#include "AnalysisTree/Matching.hpp" +#include <AnalysisTree/TaskManager.hpp> + +#include "TClonesArray.h" + +#include <cassert> +#include <cmath> + ClassImp(CbmStsTracksConverter) @@ -35,9 +36,74 @@ CbmStsTracksConverter::~CbmStsTracksConverter() { delete vtx_tracks_2_sim_; }; +void CbmStsTracksConverter::InitInput() { + auto* ioman = FairRootManager::Instance(); + + cbm_prim_vertex_ = (CbmVertex*) ioman->GetObject("PrimaryVertex."); + cbm_sts_tracks_ = (TClonesArray*) ioman->GetObject("StsTrack"); + cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); + cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch"); +} + +void CbmStsTracksConverter::Init() { + InitInput(); + + AnalysisTree::BranchConfig vtx_tracks_config(out_branch_, + AnalysisTree::DetType::kTrack); + vtx_tracks_config.AddField<float>("chi2", "chi2 of the track fit"); + vtx_tracks_config.AddField<float>("vtx_chi2", + "chi2 to to the primary vertex"); + vtx_tracks_config.AddFields<float>( + {"dcax", "dcay", "dcaz"}, + "not actuall Distance of Closest Approach, but extrapolated to z=z_vtx"); + vtx_tracks_config.AddField<int>("ndf", "number degrees of freedom"); + vtx_tracks_config.AddField<int>("q", "charge"); + vtx_tracks_config.AddField<int>("nhits", "number of hits (total MVD+STS)"); + vtx_tracks_config.AddField<int>("nhits_mvd", "number of hits in MVD"); + + if (is_write_kfinfo_) { + vtx_tracks_config.AddFields<float>({"x", "y", "z", "tx", "ty", "qp"}, + "track parameters"); + vtx_tracks_config.AddFields<float>( + {"cx0", "cx1", "cx2", "cy0", "cy1", "cy2", "cz0", "cz1", "cz2", "z0"}, + "magnetic field approximation"); + vtx_tracks_config.AddFields<float>({"cov1", + "cov2", + "cov3", + "cov4", + "cov5", + "cov6", + "cov7", + "cov8", + "cov9", + "cov10", + "cov11", + "cov12", + "cov13", + "cov14", + "cov15"}, + "covarience matrix"); + + vtx_tracks_config.AddField<int>("mother_pdg", + "PDG code of mother particle"); + vtx_tracks_config.AddField<int>("mc_pdg", "MC-true PDG code"); + vtx_tracks_config.AddField<bool>("pass_cuts", "ask Oleksii"); + + 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"); + } + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, vtx_tracks_, vtx_tracks_config); + man->AddMatching(out_branch_, match_to_, vtx_tracks_2_sim_); +} + // TODO misleading name, move field filling somewhere else? float CbmStsTracksConverter::ExtrapolateToVertex(CbmStsTrack* sts_track, - AnalysisTree::Track* track, + AnalysisTree::Track& track, int pdg) { vector<CbmStsTrack> tracks = {*sts_track}; @@ -53,16 +119,16 @@ float CbmStsTracksConverter::ExtrapolateToVertex(CbmStsTrack* sts_track, *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); + 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]; @@ -72,21 +138,16 @@ void CbmStsTracksConverter::ReadVertexTracks() { assert(cbm_prim_vertex_ && cbm_sts_tracks_); vtx_tracks_->ClearChannels(); + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); - 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 iq = branch.GetFieldId("q"); + const int indf = branch.GetFieldId("ndf"); + const int ichi2 = branch.GetFieldId("chi2"); + const int inhits = branch.GetFieldId("nhits"); + const int inhits_mvd = branch.GetFieldId("nhits_mvd"); + const int idcax = branch.GetFieldId("dcax"); + const int ivtx_chi2 = branch.GetFieldId("vtx_chi2"); const int n_sts_tracks = cbm_sts_tracks_->GetEntries(); vtx_tracks_->Reserve(n_sts_tracks); @@ -95,11 +156,9 @@ void CbmStsTracksConverter::ReadVertexTracks() { 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())); + auto& track = vtx_tracks_->AddChannel(branch); int pdg = GetMcPid((CbmTrackMatchNew*) cbm_sts_match_->At(i_track), track); - bool is_good_track = IsGoodCovMatrix(sts_track); float chi2_vertex = ExtrapolateToVertex(sts_track, track, pdg); @@ -108,46 +167,45 @@ void CbmStsTracksConverter::ReadVertexTracks() { 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())); + 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, is_good_track); } } } -void CbmStsTracksConverter::WriteKFInfo(AnalysisTree::Track* track, +void CbmStsTracksConverter::WriteKFInfo(AnalysisTree::Track& track, const CbmStsTrack* sts_track, bool is_good_track) const { - assert(track && sts_track); + assert(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); + 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(float(trackParamFirst->GetCovariance(i, j)), icov_ + iCov); } } - track->SetField(is_good_track, ipasscuts_); + track.SetField(is_good_track, ipasscuts_); } bool CbmStsTracksConverter::IsGoodCovMatrix( @@ -188,7 +246,7 @@ bool CbmStsTracksConverter::IsGoodCovMatrix( } int CbmStsTracksConverter::GetMcPid(const CbmTrackMatchNew* match, - AnalysisTree::Track* track) const { + AnalysisTree::Track& track) const { if (!is_write_kfinfo_) { return -2; } //-----------PID as in MZ's @@ -227,91 +285,18 @@ int CbmStsTracksConverter::GetMcPid(const CbmTrackMatchNew* match, pdgCode = mctrack->GetPdgCode(); } if (mctrack->GetMotherId() > -1) { - track->SetField( + track.SetField( int(((CbmMCTrack*) cbm_mc_tracks_->At(mctrack->GetMotherId())) ->GetPdgCode()), imother_pdg_); } } } - track->SetField(pdgCode, imc_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()) { diff --git a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h index e9cd9db1e5e2458307fee3f3026e2531bc714b04..4eef6df46589edc796cf568601141b5b1eaa0a81 100644 --- a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h +++ b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.h @@ -17,25 +17,13 @@ namespace AnalysisTree { 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"; - } + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) {} ~CbmStsTracksConverter() final; - void Init(std::map<std::string, void*>&) final; + void Init() final; void Exec() final; void Finish() final {} @@ -47,13 +35,14 @@ private: void MapTracks(); void InitInput(); float ExtrapolateToVertex(CbmStsTrack* sts_track, - AnalysisTree::Track* track, + AnalysisTree::Track& track, int pdg); - void WriteKFInfo(AnalysisTree::Track* track, + + void WriteKFInfo(AnalysisTree::Track& track, const CbmStsTrack* sts_track, bool is_good_track) const; bool IsGoodCovMatrix(const CbmStsTrack* sts_track) const; - int GetMcPid(const CbmTrackMatchNew* match, AnalysisTree::Track* track) const; + int GetMcPid(const CbmTrackMatchNew* match, AnalysisTree::Track& track) const; AnalysisTree::TrackDetector* vtx_tracks_ { nullptr}; ///< raw pointers are needed for TTree::Branch diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx index 9819cea643e1292cbcb8dda5b0741d04de4291d5..dd8035246ec00b61c1e45ca5fd1bfacff27f61a8 100644 --- a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx @@ -8,6 +8,7 @@ #include "AnalysisTree/Matching.hpp" +#include <AnalysisTree/TaskManager.hpp> #include <CbmGlobalTrack.h> #include <CbmMCTrack.h> #include <CbmTofHit.h> @@ -17,10 +18,9 @@ ClassImp(CbmTofHitsConverter) - void CbmTofHitsConverter::Init(std::map<std::string, void*>&) { + void CbmTofHitsConverter::Init() { - assert(!in_branches_.empty() && !out_branch_.empty() && out_config_ - && out_tree_); + assert(!out_branch_.empty()); auto* ioman = FairRootManager::Instance(); cbm_tof_hits_ = (TClonesArray*) ioman->GetObject("TofHit"); @@ -30,7 +30,8 @@ ClassImp(CbmTofHitsConverter) // cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); // cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch"); - AnalysisTree::BranchConfig tof_branch("TofHits", AnalysisTree::DetType::kHit); + AnalysisTree::BranchConfig tof_branch(out_branch_, + AnalysisTree::DetType::kHit); tof_branch.AddField<float>("mass2"); tof_branch.AddField<float>("l"); tof_branch.AddField<float>("t"); @@ -38,18 +39,10 @@ ClassImp(CbmTofHitsConverter) 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_); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, tof_hits_, tof_branch); + man->AddMatching(match_to_, out_branch_, vtx_tracks_2_tof_); } void CbmTofHitsConverter::ExtrapolateStraightLine(FairTrackParam* params, @@ -71,20 +64,14 @@ void CbmTofHitsConverter::FillTofHits() { 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"); + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); + + const int i_mass2 = branch.GetFieldId("mass2"); + const int i_qp = branch.GetFieldId("qp_tof"); + const int i_dx = branch.GetFieldId("dx"); + const int i_t = branch.GetFieldId("t"); + const int i_l = branch.GetFieldId("l"); const auto it = indexes_map_->find(match_to_); if (it == indexes_map_->end()) { @@ -122,23 +109,22 @@ void CbmTofHitsConverter::FillTofHits() { 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; } + auto& hit = tof_hits_->AddChannel(branch); + 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()) { continue; } 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()); + hit.GetId()); } // const auto* tofMatch = @@ -152,7 +138,7 @@ void CbmTofHitsConverter::FillTofHits() { // if(itofMC >= 0){ // const auto* mc_track = dynamic_cast<const // CbmMCTrack*>(cbm_mc_tracks_->At(itofMC)); - // Hit->SetField(mc_track->GetPdgCode(), i_pdg); + // hit.SetField(mc_track->GetPdgCode(), i_pdg); // // const Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); // if(stsTrackIndex<0) return; @@ -160,7 +146,7 @@ void CbmTofHitsConverter::FillTofHits() { // 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); + // hit.SetField( bool(mc_id_sts == itofMC), i_is_correct); // } // } } diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h index 990102074b51d67f84ede97113b9ea58dbdb6f57..de8eeabe7ff2cdaa3ab7ade74ee0e48cb6fcdaa2 100644 --- a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h @@ -16,14 +16,12 @@ 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"); - }; + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) {}; ~CbmTofHitsConverter() final; ; - void Init(std::map<std::string, void*>&) final; + void Init() final; void Exec() final; void Finish() final {}; diff --git a/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a350bbee95489b0b9c465bca62c8bfa93c1b2bc8 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx @@ -0,0 +1,135 @@ +#include "CbmTrdTracksConverter.h" + +#include <cassert> +#include <vector> + +#include "TClonesArray.h" + +#include <FairMCPoint.h> +#include <FairRootManager.h> + +#include "AnalysisTree/Matching.hpp" +#include <AnalysisTree/TaskManager.hpp> + +#include <CbmGlobalTrack.h> +#include <CbmTrackMatchNew.h> +#include <CbmTrdHit.h> +#include <CbmTrdTrack.h> + +#include "CbmTrdTracksConverter.h" + +ClassImp(CbmTrdTracksConverter) + + void CbmTrdTracksConverter::Init() { + + assert(!out_branch_.empty()); + auto* ioman = FairRootManager::Instance(); + + cbm_trd_tracks_ = (TClonesArray*) ioman->GetObject("TrdTrack"); + cbm_global_tracks_ = (TClonesArray*) ioman->GetObject("GlobalTrack"); + cbm_trd_hits_ = (TClonesArray*) ioman->GetObject("TrdHit"); + + AnalysisTree::BranchConfig trd_branch(out_branch_, + AnalysisTree::DetType::kTrack); + trd_branch.AddFields<float>( + {"energy_loss_0", "energy_loss_1", "energy_loss_2", "energy_loss_3"}); + trd_branch.AddFields<float>( + {"pid_like_e", "pid_like_pi", "pid_like_k", "pid_like_p", "pid_like_mu"}); + trd_branch.AddField<float>("chi2_ov_ndf"); + trd_branch.AddFields<float>({"pT_out", "p_out"}); + trd_branch.AddField<int>("n_hits"); + + auto* man = AnalysisTree::TaskManager::GetInstance(); + man->AddBranch(out_branch_, trd_tracks_, trd_branch); + man->AddMatching(match_to_, out_branch_, vtx_tracks_2_trd_); +} + +void CbmTrdTracksConverter::FillTrdTracks() { + + assert(cbm_trd_tracks_); + trd_tracks_->ClearChannels(); + vtx_tracks_2_trd_->Clear(); + + auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); + const auto& branch = out_config_->GetBranchConfig(out_branch_); + + const int i_e_loss_i = branch.GetFieldId("energy_loss_0"); + // const int i_ann = branch.GetFieldId("pid_ann"); + const int i_pid_like = branch.GetFieldId("pid_like_e"); + const int i_chi2_ov_ndf = branch.GetFieldId("chi2_ov_ndf"); + const int i_pT_out = branch.GetFieldId("pT_out"); + const int i_n_hits = branch.GetFieldId("n_hits"); + + 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; + + trd_tracks_->Reserve(cbm_global_tracks_->GetEntries()); + + for (Int_t igt = 0; igt < cbm_global_tracks_->GetEntries(); igt++) { + const auto* global_track = + static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(igt)); + + Int_t itrd = global_track->GetTrdTrackIndex(); + if (itrd < 0) continue; + + auto trd_track = static_cast<CbmTrdTrack*>(cbm_trd_tracks_->At(itrd)); + auto trd_match = static_cast<CbmTrackMatchNew*>(cbm_trd_tracks_->At(itrd)); + Int_t itrdMC = (trd_match ? trd_match->GetMatchedLink().GetIndex() : -1); + + auto& track = trd_tracks_->AddChannel(branch); + TVector3 mom, mom_last; + trd_track->GetParamFirst()->Momentum(mom); + trd_track->GetParamLast()->Momentum(mom_last); + + track.SetMomentum3(mom); + track.SetField(int(trd_track->GetNofHits()), i_n_hits); + // track.SetField(float(trd_track->GetPidANN()), i_ann); + // track.SetField(float(trd_track->GetPidWkn()), i_ann+1); + + track.SetField(float(trd_track->GetPidLikeEL()), i_pid_like); + track.SetField(float(trd_track->GetPidLikePI()), i_pid_like + 1); + track.SetField(float(trd_track->GetPidLikeKA()), i_pid_like + 2); + track.SetField(float(trd_track->GetPidLikePR()), i_pid_like + 3); + track.SetField(float(trd_track->GetPidLikeMU()), i_pid_like + 4); + + track.SetField(float(trd_track->GetNDF() > 0. + ? trd_track->GetChiSq() / trd_track->GetNDF() + : -999.), + i_chi2_ov_ndf); + + track.SetField(float(mom_last.Pt()), i_pT_out); + track.SetField(float(mom_last.Mag()), i_pT_out + 1); + + for (int i = 0; i < 4; ++i) { + track.SetField(0.f, i_e_loss_i + i); + } + + for (Int_t ihit = 0; ihit < trd_track->GetNofHits(); ihit++) { + Int_t idx = trd_track->GetHitIndex(ihit); + auto* hit = (CbmTrdHit*) cbm_trd_hits_->At(idx); + if (hit) { + // std::cout << hit->GetELoss()*1e6 << " " << hit->GetPlaneId() << std::endl; + track.SetField(float(hit->GetELoss() * 1e6), + i_e_loss_i + hit->GetPlaneId()); + } + } + + if (rec_tracks_map.empty()) { continue; } + const Int_t stsTrackIndex = global_track->GetStsTrackIndex(); + if (rec_tracks_map.find(stsTrackIndex) != rec_tracks_map.end()) { + vtx_tracks_2_trd_->AddMatch(rec_tracks_map.find(stsTrackIndex)->second, + track.GetId()); + } + } +} + +void CbmTrdTracksConverter::Exec() { FillTrdTracks(); } + +CbmTrdTracksConverter::~CbmTrdTracksConverter() { + delete trd_tracks_; + delete vtx_tracks_2_trd_; +}; diff --git a/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.h b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..ff194782af298b22c635c4d1506cf1e2efdecb43 --- /dev/null +++ b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.h @@ -0,0 +1,43 @@ +#ifndef ANALYSIS_TREE_TRDTRACKSCONVERTER_H_ +#define ANALYSIS_TREE_TRDTRACKSCONVERTER_H_ + +#include "CbmConverterTask.h" + +#include "AnalysisTree/Detector.hpp" + +class TClonesArray; + +namespace AnalysisTree { + class Matching; +} + +class CbmTrdTracksConverter final : public CbmConverterTask { + + enum kInBranches { kTrdTrack = 0, kGlobalTrack, kTrdHit }; + +public: + explicit CbmTrdTracksConverter(std::string out_branch_name, + std::string match_to = "") + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) {} + + ~CbmTrdTracksConverter() final; + + void Init() final; + void Exec() final; + void Finish() final {}; + +private: + void FillTrdTracks(); + + TClonesArray* cbm_global_tracks_ {nullptr}; + TClonesArray* cbm_trd_tracks_ {nullptr}; + TClonesArray* cbm_trd_hits_ {nullptr}; + + AnalysisTree::TrackDetector* trd_tracks_ {nullptr}; + AnalysisTree::Matching* vtx_tracks_2_trd_ {nullptr}; + + ClassDef(CbmTrdTracksConverter, 1) +}; + + +#endif //ANALYSIS_TREE_TRDTRACKSCONVERTER_H_ diff --git a/external/CMakeLists.txt b/external/CMakeLists.txt index 8c6f00abf406db5917913466836a3baa636d4e9e..056adc3e9b9986aa1b751607da64c28136cdc19a 100644 --- a/external/CMakeLists.txt +++ b/external/CMakeLists.txt @@ -25,7 +25,7 @@ if(DOWNLOAD_EXTERNALS) Include(InstallKFParticle.cmake) Include(InstallNicaFemto.cmake) Include(InstallAnalysisTree.cmake) -# Include(InstallAnalysisTreeQA.cmake) + Include(InstallAnalysisTreeQA.cmake) Include(InstallParameter.cmake) Include(InstallInput.cmake) @@ -34,6 +34,7 @@ else() # Define targets which are needed by CbmRoot but are not available # whithout the external packages add_library(ANALYSISTREE SHARED IMPORTED GLOBAL) +# add_library(ANALYSISTREEQA SHARED IMPORTED GLOBAL) add_library(NICAFEMTO SHARED IMPORTED GLOBAL) add_library(KFPARTICLE SHARED IMPORTED GLOBAL) endif() diff --git a/external/InstallAnalysisTree.cmake b/external/InstallAnalysisTree.cmake index 1a04182f8b65733a00b062e41da435e94239b3ad..48c178c2f7a89dcc1dd2c0ebfb8bab4aaced5854 100644 --- a/external/InstallAnalysisTree.cmake +++ b/external/InstallAnalysisTree.cmake @@ -1,12 +1,10 @@ -set(ANALYSISTREE_VERSION 78086bc3b1ea62054fcb2d56cc2b2f0002f14b2f) +set(ANALYSISTREE_VERSION 9f49db8c1cd0cd1079a3f9267358a9a7362bde41) 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} @@ -60,11 +58,6 @@ foreach(LIB_NAME ${ANALYSISTREE_LIBNAME}) 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 -# ) Install(DIRECTORY ${CMAKE_BINARY_DIR}/include/AnalysisTree DESTINATION include - ) \ No newline at end of file + ) diff --git a/external/InstallAnalysisTreeQA.cmake b/external/InstallAnalysisTreeQA.cmake index 797633641a3b9f9a5fe7ba6b4047fd189f41b438..2e01864d2cf241d5daddd6a0e04247e3e028ab8a 100644 --- a/external/InstallAnalysisTreeQA.cmake +++ b/external/InstallAnalysisTreeQA.cmake @@ -1,4 +1,4 @@ -set(ANALYSISTREEQA_VERSION 634a0200ba4110eb31328c327a462d11f96f7e14) +set(ANALYSISTREEQA_VERSION be56549d0ee90b628e1ec49131924a03982365d6) set(ANALYSISTREEQA_SRC_URL "https://github.com/HeavyIonAnalysis/AnalysisTreeQA.git") set(ANALYSISTREEQA_DESTDIR "${CMAKE_BINARY_DIR}/external/ANALYSISTREEQA-prefix") @@ -17,7 +17,7 @@ If(ProjectUpdated) EndIf() ExternalProject_Add(ANALYSISTREEQA - DEPENDS AnalysisTreeBase + DEPENDS AnalysisTreeInfra BUILD_IN_SOURCE 0 SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/AnalysisTreeQA BUILD_BYPRODUCTS ${ANALYSISTREEQA_LIBRARY} @@ -31,9 +31,9 @@ ExternalProject_Add(ANALYSISTREEQA -DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR} -DCMAKE_CXX_STANDARD=11 -DROOTSYS=${SIMPATH} - -DUSEBOOST=TRUE -DBOOST_ROOT=${SIMPATH} -DBoost_NO_BOOST_CMAKE=ON + -DAnalysisTreeQA_BUNDLED_AT=OFF -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON INSTALL_COMMAND ${CMAKE_COMMAND} --build . --target install ) @@ -47,8 +47,16 @@ 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 +Install( FILES + ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${ANALYSISTREEQA_LIBNAME}${CMAKE_SHARED_LIBRARY_SUFFIX} + ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${ANALYSISTREEQA_LIBNAME}_rdict.pcm + ${CMAKE_BINARY_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${ANALYSISTREEQA_LIBNAME}.rootmap + DESTINATION lib) + +Install(DIRECTORY ${CMAKE_BINARY_DIR}/include/AnalysisTreeQA + DESTINATION include + ) + +Install(FILES {CMAKE_SOURCE_DIR}/macro/analysis/common/qa/run_analysistree_qa.C + DESTINATION share/cbmroot/macro ) diff --git a/macro/C2F/CMakeLists.txt b/macro/C2F/CMakeLists.txt index 39445e65a5c88e6f9c308ebc2f0c07df9d8306e1..0fa85f4e7afca4590b8a9ca3fda99db719627d51 100644 --- a/macro/C2F/CMakeLists.txt +++ b/macro/C2F/CMakeLists.txt @@ -5,6 +5,7 @@ GENERATE_CBM_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_digi.C ${MACRO_DIR} GENERATE_CBM_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_reco_event.C ${MACRO_DIR}) GENERATE_CBM_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/analysis/common/analysis_tree_converter/run_analysis_tree_maker.C ${MACRO_DIR}) GENERATE_CBM_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/analysis/common/at_kfpf_interface/run_at_kfpf.C ${MACRO_DIR}) +GENERATE_CBM_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/analysis/common/qa/run_analysistree_qa.C ${MACRO_DIR}) # Put the .rootrc file into the directory from which root is executed. # Otherwise the rootalias file is not loaded @@ -30,47 +31,63 @@ EndIf() # --- Transport run, using run_transport.C Set(testname c2f_transport_${setup}) Add_Test(${testname} ${MACRO_DIR}/c2f_transport.sh ${NumEvents} \"${setup}\" \"data/${setup}_test\" \"\" \"TGeant4\" false) - Set_Tests_Properties(${testname} PROPERTIES TIMEOUT "900") - Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Test Passed;All ok") - Set(fixture_c2f_digi fixture_c2f_digi_${testname}) - set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_c2f_digi}) + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "900" + PASS_REGULAR_EXPRESSION "Test Passed;All ok" + FIXTURES_SETUP fixture_c2f_digi_${setup} + ) # --- Digitisation run, event-by-event Set(testname c2f_digi_event_${setup}) Add_Test(${testname} ${MACRO_DIR}/run_digi.sh \"data/${setup}_test\" ${NumEvents} \"data/${setup}_test.event\" -1.) - Set_Tests_Properties(${testname} PROPERTIES TIMEOUT "300") - Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Macro finished successfully") - set_tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_c2f_digi}) - Set(fixture_c2f_reco_event fixture_c2f_reco_event_${testname}) - set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_c2f_reco_event}) + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "300" + PASS_REGULAR_EXPRESSION "Macro finished successfully" + FIXTURES_REQUIRED fixture_c2f_digi_${setup} + FIXTURES_SETUP fixture_c2f_reco_event_${setup} + ) # --- Reconstruction run, event-by-event Set(testname c2f_reco_event_${setup}) Add_Test(${testname} ${MACRO_DIR}/run_reco_event.sh ${NumEvents} \"data/${setup}_test\" \"${setup}\") - Set_Tests_Properties(${testname} PROPERTIES TIMEOUT "1300") - Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Test Passed;All ok") - set_tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_c2f_reco_event}) - Set(fixture_c2f_ana fixture_c2f_ana_${testname}) - set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_c2f_ana}) + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "1300" + PASS_REGULAR_EXPRESSION "Test Passed;All ok" + FIXTURES_REQUIRED fixture_c2f_reco_event_${setup} + FIXTURES_SETUP fixture_c2f_ana_${setup} + ) # --- AnalysisTree Set(testname analysis_tree_maker_${setup}) - Add_Test(${testname} ${MACRO_DIR}/run_analysis_tree_maker.sh ${NumEvents} \"data/${setup}_test\" \"${setup}\") - Set_Tests_Properties(${testname} PROPERTIES TIMEOUT "300") - Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Test Passed;All ok") - set_tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_c2f_ana}) - Set(fixture_c2f_analysistree fixture_c2f_analysistree_${testname}) - set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_c2f_analysistree}) - + Add_Test(${testname} ${MACRO_DIR}/run_analysis_tree_maker.sh \"data/${setup}_test\" \"${setup}\") + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "300" + PASS_REGULAR_EXPRESSION "Test Passed;All ok" + FIXTURES_REQUIRED fixture_c2f_ana_${setup} + FIXTURES_SETUP fixture_c2f_analysistree_${setup} + ) + # --- AnalysisTree to KFPF interface Set(testname at_kfpf_${setup}) Add_Test(${testname} ${MACRO_DIR}/run_at_kfpf.sh ${NumEvents} \"data/${setup}_test\") - Set_Tests_Properties(${testname} PROPERTIES TIMEOUT "300") - Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Test Passed;All ok") - set_tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_c2f_analysistree}) - Set(fixture_c2f_at_kfpf_interface fixture_c2f_at_kfpf_interface_${testname}) - set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_c2f_at_kfpf_interface}) + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "300" + PASS_REGULAR_EXPRESSION "Test Passed;All ok" + FIXTURES_REQUIRED fixture_c2f_analysistree_${setup} + FIXTURES_SETUP fixture_c2f_at_kfpf_interface_${setup} + ) + + # --- AnalysisTreeQA + Set(testname analysistree_qa_${setup}) + Add_Test(${testname} ${MACRO_DIR}/run_analysistree_qa.sh \"data/${setup}_test.analysistree.root\" true) + Set_Tests_Properties(${testname} PROPERTIES + TIMEOUT "300" + PASS_REGULAR_EXPRESSION "Test Passed;All ok" + FIXTURES_REQUIRED fixture_c2f_at_kfpf_interface_${setup} +# FIXTURES_SETUP fixture_c2f_at_qa_${setup} + ) + EndForEach(setup IN LISTS cbm_setup) # end of test CBM setups from geometry/setup @@ -79,3 +96,5 @@ Install(FILES .rootrc c2f_transport.C DESTINATION share/cbmroot/macro/c2f ) Install(CODE "FILE(MAKE_DIRECTORY \${CMAKE_INSTALL_PREFIX}/share/cbmroot/macro/run/data)") + + diff --git a/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh b/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh index 733836b97d2b7fbb541a6cf4ad926053862548e9..305b6354b25e2bf78bba62276f3b8d377e580a5d 100755 --- a/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh +++ b/macro/analysis/common/analysis_tree_converter/batch/batch_run.sh @@ -2,10 +2,10 @@ #SBATCH -J CbmAnalysisTree #SBATCH -o out/%j.out.log #SBATCH -e error/%j.err.log -#SBATCH --time=2:00:00 -#SBATCH --array=1-50 +#SBATCH --time=8:00:00 +#SBATCH --array=1-100 -N_EVENTS=100 +N_EVENTS=500 SETUP="sis100_electron" CBMROOT_DIR=/lustre/cbm/users/klochkov/soft/cbmroot/c2f_fork/ @@ -15,7 +15,7 @@ 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 +OUTPUT_DIR=/lustre/nyx/cbm/users/klochkov/cbm/test_prod_AT/urqmd/auau/12agev/mbias/TGeant3/$INDEX INPUT_FILE=/lustre/cbm/users/ogolosov/mc/generators/urqmd/v3.4/auau/pbeam12agev_eos0/mbias/root/urqmd_$INDEX.root DATA_SET=$INDEX @@ -31,9 +31,9 @@ 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 ./ +cp $MACRO_DIR/run_analysis_tree_maker.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_digi.C(\"$DATA_SET\", $N_EVENTS, \"$DATA_SET\", -1)" >& 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 +root -l -q -b "run_analysis_tree_maker.C(\"$DATA_SET\", \"$SETUP\")" >& log_$INDEX.txt diff --git a/macro/analysis/common/analysis_tree_converter/run_analysis_tree_maker.C b/macro/analysis/common/analysis_tree_converter/run_analysis_tree_maker.C index 8e9c9f94bdc9dbd0ed2a969fd24d7aa2e3372910..a5d35d3126ed692e4c96f02c6f4f15560cd73074 100644 --- a/macro/analysis/common/analysis_tree_converter/run_analysis_tree_maker.C +++ b/macro/analysis/common/analysis_tree_converter/run_analysis_tree_maker.C @@ -1,6 +1,5 @@ -void run_analysis_tree_maker(Int_t nEvents = 2, - TString dataSet = "test", +void run_analysis_tree_maker(TString dataSet = "../../../run/test", TString setupName = "sis100_electron") { const std::string system = "Au+Au"; // TODO can we read it automatically? const float beam_mom = 12.; @@ -95,13 +94,25 @@ void run_analysis_tree_maker(Int_t nEvents = 2, } run->AddTask(l1); + // --- TRD pid tasks + if (setup->IsActive(ECbmModuleId::kTrd)) { + + CbmTrdSetTracksPidLike* trdLI = + new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood"); + trdLI->SetUseMCInfo(kTRUE); + trdLI->SetUseMomDependence(kTRUE); + run->AddTask(trdLI); + std::cout << "-I- : Added task " << trdLI->GetName() << std::endl; + // ------------------------------------------------------------------------ + } + // AnalysisTree converter auto* man = new CbmConverterManager(); man->SetSystem(system); man->SetBeamMomentum(beam_mom); - man->SetOutFileName(outFile); - man->SetOutTreeName("aTree"); + man->SetOutputName(outFile, "rTree"); + // man->SetOutTreeName("aTree"); man->AddTask(new CbmSimEventHeaderConverter("SimEventHeader")); man->AddTask(new CbmRecEventHeaderConverter("RecEventHeader")); @@ -113,7 +124,9 @@ void run_analysis_tree_maker(Int_t nEvents = 2, taskCbmStsTracksConverter->SetIsReproduceCbmKFPF(); man->AddTask(taskCbmStsTracksConverter); + man->AddTask(new CbmRichRingsConverter("RichRings", "VtxTracks")); man->AddTask(new CbmTofHitsConverter("TofHits", "VtxTracks")); + man->AddTask(new CbmTrdTracksConverter("TrdTracks", "VtxTracks")); man->AddTask(new CbmPsdModulesConverter("PsdModules")); run->AddTask(man); @@ -134,7 +147,7 @@ void run_analysis_tree_maker(Int_t nEvents = 2, run->Init(); std::cout << "Starting run" << std::endl; - run->Run(0, nEvents); + run->Run(0); // ------------------------------------------------------------------------ timer.Stop(); @@ -143,6 +156,7 @@ void run_analysis_tree_maker(Int_t nEvents = 2, 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 ------------------------------------ @@ -166,5 +180,5 @@ void run_analysis_tree_maker(Int_t nEvents = 2, // Generate_CTest_Dependency_File(depFile); // ------------------------------------------------------------------------ - RemoveGeoManager(); + // RemoveGeoManager(); } diff --git a/macro/analysis/common/qa/run_analysistree_qa.C b/macro/analysis/common/qa/run_analysistree_qa.C new file mode 100644 index 0000000000000000000000000000000000000000..ec6d9f8b1aaacecd8a1e0f3a192ac83575521568 --- /dev/null +++ b/macro/analysis/common/qa/run_analysistree_qa.C @@ -0,0 +1,311 @@ +//#include <AnalysisTree/TaskManager.hpp> +//#include "AnalysisTree/Cuts.hpp" +// +//#include "src/EntryConfig.hpp" +//#include "src/Task.hpp" +//#include "src/Utils.hpp" + +typedef AnalysisTree::QA::EntryConfig::PlotType PlotType; +using AnalysisTree::QA::gNbins; + +using namespace AnalysisTree; + +void KFPFTracksQA(QA::Task& task); +void VertexTracksQA(QA::Task& task); +void TofHitsQA(QA::Task& task); +void SimParticlesQA(QA::Task& task); +void SimEventHeaderQA(QA::Task& task); +void RecEventHeaderQA(QA::Task& task); +void EfficiencyMaps(QA::Task& task); +void TrdTracksQA(QA::Task& task); +void RichRingsQA(QA::Task& task); + +void run_analysistree_qa(std::string filelist, bool is_single_file = false); + +const std::string sim_event_header = "SimEventHeader"; +const std::string rec_event_header = "RecEventHeader"; +const std::string sim_particles = "SimParticles"; +const std::string rec_tracks = "VtxTracks"; +const std::string tof_hits = "TofHits"; +const std::string trd_tracks = "TrdTracks"; +const std::string rich_rings = "RichRings"; + +int main(int argc, char** argv) { + if (argc <= 1) { + std::cout << "Not enough arguments! Please use:" << std::endl; + std::cout << " ./cbm_qa filelist" << std::endl; + return -1; + } + + const std::string filelist = argv[1]; + run_analysistree_qa(filelist); + return 0; +} + +void run_analysistree_qa(std::string filelist, bool is_single_file) { + + if (is_single_file) { + std::ofstream fl("fl_temp.txt"); + fl << filelist << "\n"; + fl.close(); + filelist = "fl_temp.txt"; + } + + TaskManager* man = TaskManager::GetInstance(); + + auto* task = new QA::Task; + task->SetOutputFileName("cbm_qa.root"); + + RichRingsQA(*task); + TrdTracksQA(*task); + // KFPFTracksQA(*task); + VertexTracksQA(*task); + TofHitsQA(*task); + SimParticlesQA(*task); + SimEventHeaderQA(*task); + RecEventHeaderQA(*task); + // EfficiencyMaps(*task); + // AddParticlesFlowQA(task, sim_particles, {sim_event_header, "psi_RP"}, {2212, 211, -211}); + + man->AddTask(task); + + man->Init({filelist}, {"rTree"}); + man->Run(-1); + man->Finish(); + + if (is_single_file) { + // ----- Finish ------------------------------------------------------- + std::cout << " Test passed" << std::endl; + std::cout << " All ok " << std::endl; + } +} + +void TrdTracksQA(QA::Task& task) { + AddTrackQA(&task, trd_tracks); + AddTracksMatchQA(&task, trd_tracks, rec_tracks); + + task.AddH1({"TRD total energy loss, keV", + {trd_tracks, "energy_loss_total"}, + {gNbins, 0, 1}}); + + task.AddH1({"TRD energy loss in 1st station, keV", + {trd_tracks, "energy_loss_0"}, + {gNbins, 0, 1}}); + task.AddH1({"TRD energy loss in 2nd station", + {trd_tracks, "energy_loss_1"}, + {gNbins, 0, 1}}); + task.AddH1({"TRD energy loss in 3rd station", + {trd_tracks, "energy_loss_2"}, + {gNbins, 0, 1}}); + task.AddH1({"TRD energy loss in 4th station", + {trd_tracks, "energy_loss_3"}, + {gNbins, 0, 1}}); + + task.AddH1({"Number of hits in TRD", {trd_tracks, "n_hits"}, {6, 0, 6}}); + + task.AddH1({"PID ANN", {trd_tracks, "pid_ann"}, {gNbins, -1, 1}}); + // task.AddH1({"PID WKN", {trd_tracks, "pid_wkn"}, {gNbins, -1, 1}}); + task.AddH1( + {"PID Likelihood, e^{-}", {trd_tracks, "pid_like_e"}, {gNbins, -1, 1}}); + task.AddH1( + {"PID Likelihood, #pi", {trd_tracks, "pid_like_pi"}, {gNbins, -1, 1}}); + task.AddH1( + {"PID Likelihood, K", {trd_tracks, "pid_like_k"}, {gNbins, -1, 1}}); + task.AddH1( + {"PID Likelihood, p", {trd_tracks, "pid_like_p"}, {gNbins, -1, 1}}); + task.AddH1( + {"PID Likelihood, #mu", {trd_tracks, "pid_like_mu"}, {gNbins, -1, 1}}); + + task.AddH1({"#chi^{2}/NDF", {trd_tracks, "chi2_ov_ndf"}, {gNbins, 0, 30}}); + task.AddH1({"p_{T}^{last} (GeV/c)", {trd_tracks, "pT_out"}, {gNbins, 0, 10}}); + task.AddH1({"p^{last} (GeV/c)", {trd_tracks, "p_out"}, {gNbins, 0, 10}}); +} + +void RichRingsQA(QA::Task& task) { + task.AddH1( + {"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}}); + task.AddH1( + {"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}}); + task.AddH2( + {"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}}, + {"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}}); + + task.AddH1({"Ring radius", {rich_rings, "radius"}, {gNbins, 0, 10}}); + + task.AddH1({"n_hits", {rich_rings, "n_hits"}, {gNbins, 0, 60}}); + task.AddH1( + {"n_hits_on_ring", {rich_rings, "n_hits_on_ring"}, {gNbins, 0, 100}}); + task.AddH1({"axis_a", {rich_rings, "axis_a"}, {gNbins, 0, 10}}); + task.AddH1({"axis_b", {rich_rings, "axis_b"}, {gNbins, 0, 10}}); + + task.AddH1({"chi2_ov_ndf", {rich_rings, "chi2_ov_ndf"}, {gNbins, 0, 1}}); + task.AddH1({"phi_ellipse", {rich_rings, "phi_ellipse"}, {gNbins, -3.2, 3.2}}); +} + +void VertexTracksQA(QA::Task& task) { + AddTrackQA(&task, rec_tracks); + if (!sim_particles.empty()) { + AddTracksMatchQA(&task, rec_tracks, sim_particles); + } + + Variable chi2_over_ndf( + "chi2_ndf", + {{rec_tracks, "chi2"}, {rec_tracks, "ndf"}}, + [](std::vector<double>& var) { return var.at(0) / var.at(1); }); + + task.AddH1({"DCA_{x}, cm", {rec_tracks, "dcax"}, {gNbins, -1, 1}}); + task.AddH1({"DCA_{y}, cm", {rec_tracks, "dcay"}, {gNbins, -1, 1}}); + task.AddH1({"DCA_{z}, cm", {rec_tracks, "dcaz"}, {gNbins, -1, 1}}); + task.AddH1({"NDF", {rec_tracks, "ndf"}, {30, 0, 30}}); + task.AddH1({"Number of hits (STS+MVD)", {rec_tracks, "nhits"}, {15, 0, 15}}); + task.AddH1({"Number of hits (MVD)", {rec_tracks, "nhits_mvd"}, {5, 0, 5}}); + task.AddH1({"#chi^{2}_{vertex}", {rec_tracks, "vtx_chi2"}, {500, 0, 100}}); + task.AddH1({"#chi^{2}/NDF", chi2_over_ndf, {gNbins, 0, 10}}); + task.AddH2({"DCA_{x}, cm", {rec_tracks, "dcax"}, {gNbins, -1, 1}}, + {"DCA_{y}, cm", {rec_tracks, "dcay"}, {gNbins, -1, 1}}); +} + +void TofHitsQA(QA::Task& task) { + task.AddH1({"TOF hit x-position (cm)", {tof_hits, "x"}, {gNbins, -600, 600}}); + task.AddH1({"TOF hit y-position (cm)", {tof_hits, "y"}, {gNbins, -400, 400}}); + task.AddH1({"TOF hit z-position (cm)", {tof_hits, "z"}, {gNbins, 600, 800}}); + + task.AddH2({"TOF hit x-position (cm)", {tof_hits, "x"}, {gNbins, -600, 600}}, + {"TOF hit y-position (cm)", {tof_hits, "y"}, {gNbins, -600, 600}}); + + Variable qp_sts("qp_reco", + {{rec_tracks, "q"}, {rec_tracks, "p"}}, + [](std::vector<double>& qp) { return qp.at(0) * qp.at(1); }); + task.AddH2({"sign(q)*p, GeV/c", qp_sts, {gNbins, -10, 10}}, + {"m^{2}, GeV^{2}/c^{2}", {tof_hits, "mass2"}, {gNbins, -1, 2}}); + task.AddH2({"sign(q)*p STS, GeV/c", qp_sts, {gNbins, -10, 10}}, + {"sign(q)*p TOF, GeV/c", {tof_hits, "qp_tof"}, {gNbins, -10, 10}}); +} + +void SimParticlesQA(QA::Task& task) { + AddParticleQA(&task, sim_particles); + + // task.AddH1({"PDG code", {sim_particles, "pid"}, {8001, -3999.5, 4000.5}}); +} + +void SimEventHeaderQA(QA::Task& task) { + task.AddH1( + {"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}}); + task.AddH1( + {"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}}); + task.AddH1( + {"z_{vertex}^{MC} (cm)", {sim_event_header, "vtx_z"}, {gNbins, -1, 1}}); + task.AddH1({"b (fm)", {sim_event_header, "b"}, {gNbins, 0, 20}}); + task.AddH1({"#Psi_{RP}", {sim_event_header, "psi_RP"}, {gNbins, 0, 6.5}}); + + task.AddH2( + {"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}}, + {"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}}); +} + +void RecEventHeaderQA(QA::Task& task) { + task.AddH1({"x_{vertex} (cm)", {rec_event_header, "vtx_x"}, {gNbins, -1, 1}}); + task.AddH1({"y_{vertex} (cm)", {rec_event_header, "vtx_y"}, {gNbins, -1, 1}}); + task.AddH1({"z_{vertex} (cm)", {rec_event_header, "vtx_z"}, {gNbins, -1, 1}}); + task.AddH1( + {"#chi^{2}_{vertex fit}", {rec_event_header, "vtx_chi2"}, {gNbins, 0, 5}}); + + task.AddH1({"E_{PSD} (GeV)", {rec_event_header, "Epsd"}, {gNbins, 0, 60}}); + task.AddH1({"M_{tracks}", {rec_event_header, "M"}, {800, 0, 800}}); + task.AddH1({"Event ID", {rec_event_header, "evt_id"}, {gNbins, 0, 2000}}); + + task.AddH2({"x_{vertex} (cm)", {rec_event_header, "vtx_x"}, {gNbins, -1, 1}}, + {"y_{vertex} (cm)", {rec_event_header, "vtx_y"}, {gNbins, -1, 1}}); + + task.AddH2({"x_{vertex} (cm)", {rec_event_header, "vtx_x"}, {gNbins, -1, 1}}, + {"x_{vertex} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}}); + task.AddH2({"y_{vertex} (cm)", {rec_event_header, "vtx_y"}, {gNbins, -1, 1}}, + {"y_{vertex} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}}); + task.AddH2({"z_{vertex} (cm)", {rec_event_header, "vtx_z"}, {gNbins, -1, 1}}, + {"z_{vertex} (cm)", {sim_event_header, "vtx_z"}, {gNbins, -1, 1}}); +} + +void EfficiencyMaps(QA::Task& task) { + + const float y_beam = 1.62179f; // TODO from DataHeader + const float p_mass = 0.938; + const float pi_mass = 0.14; + + Variable proton_y("y-y_{beam}", + {{rec_tracks, "p"}, {rec_tracks, "pz"}}, + [y_beam, p_mass](std::vector<double>& var) { + const float e = sqrt(p_mass * p_mass + var[0] * var[0]); + return 0.5 * log((e + var[1]) / (e - var[1])); + }); + + Variable pion_y("y-y_{beam}", + {{rec_tracks, "p"}, {rec_tracks, "pz"}}, + [y_beam, pi_mass](std::vector<double>& var) { + const float e = sqrt(pi_mass * pi_mass + var[0] * var[0]); + return 0.5 * log((e + var[1]) / (e - var[1])); + }); + + Cuts* mc_protons = + new Cuts("McProtons", {EqualsCut({sim_particles + ".pid"}, 2212)}); + Cuts* mc_pions_neg = + new Cuts("McPionsNeg", {EqualsCut({sim_particles + ".pid"}, -211)}); + Cuts* mc_pions_pos = + new Cuts("McPionsPos", {EqualsCut({sim_particles + ".pid"}, -211)}); + + task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, + mc_protons); + task.AddH2({"#it{y}_{Lab}", proton_y, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}}, + mc_protons); + + task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, + mc_pions_neg); + task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}}, + mc_pions_neg); + + task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, + mc_pions_pos); + task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}}, + {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}}, + mc_pions_pos); +} + +void KFPFTracksQA(QA::Task& task) { + const std::vector<std::string> field_par_names_ { + "cx0", "cx1", "cx2", "cy1", "cy2", "cz0", "cz1", "cz2"}; + const std::vector<std::string> cov_names_ {"cov1", + "cov2", + "cov3", + "cov4", + "cov5", + "cov6", + "cov7", + "cov8", + "cov9", + "cov10", + "cov11", + "cov12", + "cov13", + "cov14", + "cov15"}; + // KFPF tracks + task.AddH1({"x, cm", {rec_tracks, "x"}, {gNbins, -10, 10}}); + task.AddH1({"y, cm", {rec_tracks, "y"}, {gNbins, -10, 10}}); + task.AddH1({"z, cm", {rec_tracks, "z"}, {gNbins, 0, 40}}); + task.AddH1({"t_{x}", {rec_tracks, "tx"}, {gNbins, -2, 2}}); + task.AddH1({"t_{y}", {rec_tracks, "ty"}, {gNbins, -2, 2}}); + task.AddH1({"q*p, GeV/c", {rec_tracks, "qp"}, {gNbins, -10, 10}}); + + for (const auto& field_par_name : field_par_names_) + task.AddH1({field_par_name, {rec_tracks, field_par_name}, {gNbins, -1, 1}}); + + for (const auto& cov_name : cov_names_) + task.AddH1({cov_name, {rec_tracks, cov_name}, {1000, -.1, .1}}); + + task.AddH1({"cy0", {rec_tracks, "cy0"}, {gNbins, -12, -8}}); + task.AddH1({"z0", {rec_tracks, "z0"}, {gNbins, 0, 40}}); +}