From 8fe7922f069ebe7b0f5ac2ebbcfcca3921977575 Mon Sep 17 00:00:00 2001
From: Frederic Julian Linz <f.linz@gsi.de>
Date: Mon, 27 Nov 2023 11:22:12 +0000
Subject: [PATCH] Add FSD to AnalysisTreeConverter + Updates to AT binary
---
.../analysis_tree_converter/CMakeLists.txt | 2 +
.../CbmAnalysisTreeInterfaceLinkDef.h | 2 +
.../CbmConverterManager.cxx | 67 ++--
.../CbmFsdHitsConverter.cxx | 306 ++++++++++++++++++
.../CbmFsdHitsConverter.h | 94 ++++++
.../CbmFsdModulesConverter.cxx | 69 ++++
.../CbmFsdModulesConverter.h | 31 ++
.../CbmPsdModulesConverter.cxx | 5 +-
.../config/ATConfig_event_ideal.yaml | 3 +
.../config/ATConfig_timeslice.yaml | 3 +
.../analysis_tree_converter/steer/Config.cxx | 8 +
.../analysis_tree_converter/steer/Config.h | 4 +
.../analysis_tree_converter/steer/Run.cxx | 31 ++
.../analysis_tree_converter/steer/Run.h | 15 +-
.../steer/TaskFactory.cxx | 15 +-
.../run_analysis_tree_maker_json_config.C | 10 +-
16 files changed, 631 insertions(+), 34 deletions(-)
create mode 100644 analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx
create mode 100644 analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h
create mode 100644 analysis/common/analysis_tree_converter/CbmFsdModulesConverter.cxx
create mode 100644 analysis/common/analysis_tree_converter/CbmFsdModulesConverter.h
diff --git a/analysis/common/analysis_tree_converter/CMakeLists.txt b/analysis/common/analysis_tree_converter/CMakeLists.txt
index 273ea100fa..4d92516b76 100644
--- a/analysis/common/analysis_tree_converter/CMakeLists.txt
+++ b/analysis/common/analysis_tree_converter/CMakeLists.txt
@@ -14,6 +14,8 @@ set(SRCS
CbmStsTracksConverter.cxx
CbmSimTracksConverter.cxx
CbmPsdModulesConverter.cxx
+ CbmFsdModulesConverter.cxx
+ CbmFsdHitsConverter.cxx
CbmTofHitsConverter.cxx
CbmTrdTracksConverter.cxx
CbmRichRingsConverter.cxx
diff --git a/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h
index eac02d0865..c497cdba8b 100644
--- a/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h
+++ b/analysis/common/analysis_tree_converter/CbmAnalysisTreeInterfaceLinkDef.h
@@ -16,6 +16,8 @@
#pragma link C++ class CbmStsTracksConverter + ;
#pragma link C++ class CbmTofHitsConverter + ;
#pragma link C++ class CbmPsdModulesConverter + ;
+#pragma link C++ class CbmFsdModulesConverter + ;
+#pragma link C++ class CbmFsdHitsConverter + ;
#pragma link C++ class CbmTrdTracksConverter + ;
#pragma link C++ class CbmRichRingsConverter + ;
//#pragma link C++ class CbmMatchEvents + ;
diff --git a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx
index 6240f8903d..256d5ffde3 100644
--- a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx
+++ b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx
@@ -4,18 +4,18 @@
#include "CbmConverterManager.h"
+#include "AnalysisTree/DataHeader.hpp"
+#include "AnalysisTree/TaskManager.hpp"
#include "CbmConverterTask.h"
+#include "CbmDefs.h"
#include "CbmEvent.h"
-
#include "TClonesArray.h"
#include "TGeoBBox.h"
#include "TGeoManager.h"
+#include "TString.h"
#include <iostream>
-#include "AnalysisTree/DataHeader.hpp"
-#include "AnalysisTree/TaskManager.hpp"
-
ClassImp(CbmConverterManager);
InitStatus CbmConverterManager::Init()
@@ -83,41 +83,56 @@ void CbmConverterManager::FillDataHeader()
data_header->SetBeamMomentum(beam_mom_);
data_header->SetTimeSliceLength(ts_length_);
- auto& psd_mod_pos = data_header->AddDetector();
- const int psd_node_id = 6;
- const char* module_name_prefix = "module";
+ auto& specDet_mod_pos = data_header->AddDetector();
+ TString specDet_names[2] = {(TString) ToString(ECbmModuleId::kPsd), (TString) ToString(ECbmModuleId::kFsd)};
+ const char* module_name = "module"; // PSD and FSD modules always contains "module"
+
+ TVector3 frontFaceGlobal;
+ Int_t nSpecDetModules = 0;
+ TString specDetNameTag = "";
+
+ TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
+ TGeoNode* curNode;
+ geoIterator.Reset(); // safety to reset to "cave" before the loop starts
+ while ((curNode = geoIterator())) {
+ TString nodePath;
+ geoIterator.GetPath(nodePath);
+ if (!(nodePath.Contains(specDet_names[0], TString::kIgnoreCase)
+ || nodePath.Contains(specDet_names[1], TString::kIgnoreCase))) {
+ geoIterator.Skip(); // skip current branch when it is not FSD => should speed up
+ continue; // don't do anything for this branch
+ }
- auto* geoMan = gGeoManager;
- auto* caveNode = geoMan->GetTopNode();
- auto* psdNode = caveNode->GetDaughter(psd_node_id);
- std::cout << "-I- " << psdNode->GetName() << std::endl;
+ TString nodeName(curNode->GetName());
- auto psdGeoMatrix = psdNode->GetMatrix();
- auto psdBox = (TGeoBBox*) psdNode->GetVolume()->GetShape();
- TVector3 frontFaceLocal(0, 0, -psdBox->GetDZ());
+ // spectator detector as a whole
+ if (nodeName.Contains(specDet_names[0], TString::kIgnoreCase)
+ || nodeName.Contains(specDet_names[1], TString::kIgnoreCase)) {
+ specDetNameTag = nodeName;
- TVector3 frontFaceGlobal;
- psdGeoMatrix->LocalToMaster(&frontFaceLocal[0], &frontFaceGlobal[0]);
+ auto specDetGeoMatrix = curNode->GetMatrix();
+ auto specDetBox = (TGeoBBox*) curNode->GetVolume()->GetShape();
+ TVector3 frontFaceLocal(0, 0, -specDetBox->GetDZ());
+ specDetGeoMatrix->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)) {
+ // modules of spectator detector
+ if (nodeName.Contains(module_name)) {
+ nSpecDetModules++;
- auto geoMatrix = daughter->GetMatrix();
+ auto geoMatrix = curNode->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();
+ auto* module = specDet_mod_pos.AddChannel();
module->SetPosition(x, y, frontFaceGlobal[2]);
+ LOG(info) << "Module " << nSpecDetModules << " : " << Form("(%.3f, %.3f)", x, y);
}
}
+ LOG(info) << "Detector " << specDetNameTag << " with " << nSpecDetModules << " modules";
+
task_manager_->SetOutputDataHeader(data_header);
}
CbmConverterManager::~CbmConverterManager() = default;
diff --git a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx
new file mode 100644
index 0000000000..8e60ee5371
--- /dev/null
+++ b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx
@@ -0,0 +1,306 @@
+/* Copyright (C) 2020-2021 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Viktor Klochkov [committer] */
+
+#include "CbmFsdHitsConverter.h"
+
+#include "AnalysisTree/Matching.hpp"
+#include "CbmDefs.h"
+#include "CbmDigiManager.h"
+#include "CbmEvent.h"
+#include "CbmFsdPoint.h"
+#include "CbmGlobalTrack.h"
+#include "CbmKFTrack.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
+#include "CbmMCTrack.h"
+#include "TClonesArray.h"
+#include <CbmFsdDigi.h>
+#include <CbmFsdHit.h>
+#include <CbmGlobalTrack.h>
+#include <CbmStsTrack.h>
+#include <CbmTofHit.h>
+#include <CbmTrackMatchNew.h>
+
+#include <FairMCPoint.h>
+#include <FairRootManager.h>
+#include <Logger.h>
+
+#include <TGeoManager.h> // for TGeoManager, gGeoManager
+#include <TGeoMatrix.h> // for TGeoMatrix
+#include <TGeoNode.h> // for TGeoIterator, TGeoNode
+#include <TGeoShape.h> // for TGeoBBox etc.
+#include <TGeoVolume.h> // for TGeoVolume
+
+#include <AnalysisTree/TaskManager.hpp>
+#include <cassert>
+
+ClassImp(CbmFsdHitsConverter);
+
+void CbmFsdHitsConverter::Init()
+{
+ LOG(debug) << " THIS IS A TEST -> CbmFsdHitsConverter is being initialized LUKAS" << std::endl;
+
+ assert(!out_branch_.empty());
+ auto* ioman = FairRootManager::Instance();
+
+ cbm_fsd_hits_ = (TClonesArray*) ioman->GetObject("FsdHit");
+
+ cbm_global_tracks_ = (TClonesArray*) ioman->GetObject("GlobalTrack");
+ cbm_sts_tracks_ = (TClonesArray*) ioman->GetObject("StsTrack");
+ cbm_tof_hits_ = (TClonesArray*) ioman->GetObject("TofHit");
+ cbm_fsd_hitmatch_ = (TClonesArray*) ioman->GetObject("FsdHitMatch");
+ cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack");
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if (!fDigiMan->IsMatchPresent(ECbmModuleId::kFsd)) {
+ LOG(error) << " No FsdDigiMatch input array present !!";
+ }
+ //cbm_fsd_digis_ = (TClonesArray*) ioman->GetObject("FsdDigi");
+ //cbm_fsd_digimatch_ = (TClonesArray*) ioman->GetObject("FsdDigiMatch");
+
+ cbm_mc_manager_ = dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
+ //cbm_mc_tracks_new_ = cbm_mc_manager_->InitBranch("MCTrack");
+ cbm_fsd_points_new_ = cbm_mc_manager_->InitBranch("FsdPoint");
+
+ AnalysisTree::BranchConfig fsd_branch(out_branch_, AnalysisTree::DetType::kHit);
+
+ fsd_branch.AddField<float>("dEdx", "Energy deposition of given FSD hit [GeV]");
+ fsd_branch.AddField<float>("t", "Reconstructed time of given FSD hit [ps]");
+ fsd_branch.AddField<float>("mass2", "Calculated mass squared from extrapolated global track (by CbmKFTrack) to FSD "
+ "plane and FSD hit time [GeV^2/c^4]");
+ fsd_branch.AddField<float>("l", "Lenght of the extrapolated global track (by CbmKFTrack) to FSD plane [cm]");
+ fsd_branch.AddField<float>(
+ "qp", "charge * momentum of the extrapoleted global track (by CbmKFTrack) to FSD plane [GeV/c]");
+ fsd_branch.AddFields<float>(
+ {"dx", "dy", "dz"},
+ "Component of a 3D distance between FSD hit and extrapolated global track (by CbmKFTrack) [cm]");
+ fsd_branch.AddField<float>(
+ "chi2GtrackHit", "chi2 between extrapolated global track (by CbmKFTrack) to FSD plane (?FIXED Z?) and FSD hit");
+ fsd_branch.AddField<int>(
+ "bestMatchedGtrack2HitId",
+ "Index of best match between extrapolated global track (by CbmKFTrack) and FSD hit based on min chi2GtrackHit");
+ fsd_branch.AddField<int>("multMCtracks", "number of MC particles that cotributed by energy deposition to FSD hit");
+ fsd_branch.AddField<float>("maxWeightMCtrack",
+ "weight of matched link from Hit to Point (?highest energy deposition?)");
+ fsd_branch.AddField<float>("dtHitPoint", "Time difference between FSD hit and matched MC point [ps]");
+ fsd_branch.AddFields<float>({"dxHitPoint", "dyHitPoint", "dzHitPoint"},
+ "Component of a 3D distance between FSD hit and matched MC point [cm]");
+
+ i_edep_ = fsd_branch.GetFieldId("dEdx");
+ i_t_ = fsd_branch.GetFieldId("t");
+ i_mass2_ = fsd_branch.GetFieldId("mass2");
+ i_qp_ = fsd_branch.GetFieldId("qp");
+ i_dx_ = fsd_branch.GetFieldId("dx");
+ i_l_ = fsd_branch.GetFieldId("l");
+ i_dtHP_ = fsd_branch.GetFieldId("dtHitPoint");
+ i_dxHP_ = fsd_branch.GetFieldId("dxHitPoint");
+ i_chi2_ = fsd_branch.GetFieldId("chi2GtrackHit");
+ i_bestMatchedGTrack_ = fsd_branch.GetFieldId("bestMatchedGtrack2HitId");
+ i_multMC_ = fsd_branch.GetFieldId("multMCtracks");
+ i_topW_ = fsd_branch.GetFieldId("maxWeightMCtrack");
+
+ if (fsdgtrack_minChi2_ == 0) fsdgtrack_minChi2_ = -1.;
+ if (fsdgtrack_maxChi2_ == 0) fsdgtrack_maxChi2_ = 10000.;
+
+ auto* man = AnalysisTree::TaskManager::GetInstance();
+ man->AddBranch(fsd_hits_, fsd_branch);
+ man->AddMatching(match_to_, out_branch_, vtx_tracks_2_fsd_);
+ man->AddMatching(out_branch_, mc_tracks_, fsd_hits_2_mc_tracks_);
+}
+
+FairTrackParam CbmFsdHitsConverter::ExtrapolateGtrack(Double_t zpos, FairTrackParam params)
+{
+ // follow inspiration from CbmL1TofMerger
+ FairTrackParam returnParams;
+ CbmKFTrack kfTrack;
+ kfTrack.SetTrackParam(params);
+ kfTrack.Extrapolate(zpos); //this will do straight track extrapolation
+ kfTrack.GetTrackParam(returnParams);
+ return returnParams;
+}
+
+
+Double_t CbmFsdHitsConverter::Chi2FsdhitGtrack(CbmFsdHit* hit, FairTrackParam inputParams)
+{
+ FairTrackParam params = ExtrapolateGtrack(hit->GetZ(), inputParams);
+
+ Float_t dX = params.GetX() - hit->GetX();
+ Float_t dY = params.GetY() - hit->GetY();
+
+ Double_t cxx = params.GetCovariance(0, 0) + hit->GetDx() * hit->GetDx();
+ Double_t cxy = params.GetCovariance(0, 1) + hit->GetDxy() * hit->GetDxy();
+ Double_t cyy = params.GetCovariance(1, 1) + hit->GetDy() * hit->GetDy();
+ Double_t chi2 = 0.5 * (dX * dX * cxx - 2 * dX * dY * cxy + dY * dY * cyy) / (cxx * cyy - cxy * cxy);
+
+ return chi2;
+}
+
+void CbmFsdHitsConverter::ProcessData(CbmEvent* event)
+{
+ assert(cbm_fsd_hits_);
+ fsd_hits_->ClearChannels();
+ vtx_tracks_2_fsd_->Clear();
+ fsd_hits_2_mc_tracks_->Clear();
+
+ auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig();
+ const auto& branch = out_config_->GetBranchConfig(out_branch_);
+
+ auto rec_tracks_map = GetMatchMap(match_to_);
+ auto sim_tracks_map = GetMatchMap(mc_tracks_);
+
+ int file_id{0}, event_id{0};
+ if (event) {
+ auto match = event->GetMatch();
+ if (!match) return;
+ file_id = event->GetMatch()->GetMatchedLink().GetFile();
+ event_id = event->GetMatch()->GetMatchedLink().GetEntry();
+ }
+ else {
+ event_id = FairRootManager::Instance()->GetEntryNr();
+ }
+
+ const int n_fsd_hits = event ? event->GetNofData(ECbmDataType::kFsdHit) : cbm_fsd_hits_->GetEntriesFast();
+ if (n_fsd_hits <= 0) {
+ LOG(warn) << "No FSD hits!";
+ return;
+ }
+
+ const int n_tracks = event ? event->GetNofData(ECbmDataType::kGlobalTrack) : cbm_global_tracks_->GetEntriesFast();
+ if (n_tracks <= 0) {
+ LOG(warn) << "No Global Tracks!";
+ return;
+ }
+
+ fsd_hits_->Reserve(n_fsd_hits);
+
+ // first loop over all FSD hits
+ for (Int_t ifh = 0; ifh < n_fsd_hits; ifh++) {
+ const auto fsdHI = event ? event->GetIndex(ECbmDataType::kFsdHit, ifh) : ifh;
+ auto* fsdHit = static_cast<CbmFsdHit*>(cbm_fsd_hits_->At(fsdHI));
+
+ auto& hit = fsd_hits_->AddChannel(branch);
+
+ const Float_t hitX = fsdHit->GetX();
+ const Float_t hitY = fsdHit->GetY();
+ const Float_t hitZ = fsdHit->GetZ();
+ const Float_t eLoss = fsdHit->GetEdep();
+ const Float_t time = fsdHit->GetTime();
+
+ hit.SetPosition(hitX, hitY, hitZ);
+ hit.SetSignal(time);
+ hit.SetField(eLoss, i_edep_);
+
+ Int_t multMC = 0;
+ Float_t highestWeight = 0.;
+ const auto* fsdHitMatch = dynamic_cast<CbmMatch*>(cbm_fsd_hitmatch_->At(fsdHI));
+ if (fsdHitMatch && fsdHitMatch->GetNofLinks() > 0) {
+ highestWeight = fsdHitMatch->GetMatchedLink().GetWeight();
+
+ for (int32_t ilDigi = 0; ilDigi < fsdHitMatch->GetNofLinks(); ilDigi++) {
+ const auto& digiLink = fsdHitMatch->GetLink(ilDigi);
+ if (digiLink.GetFile() != file_id || digiLink.GetEntry() != event_id) { // match from different event
+ continue;
+ }
+
+ const auto* fsdDigiMatch = fDigiMan->GetMatch(ECbmModuleId::kFsd, digiLink.GetIndex());
+ if (fsdDigiMatch && fsdDigiMatch->GetNofLinks() > 0) {
+ for (int32_t ilPoint = 0; ilPoint < fsdDigiMatch->GetNofLinks(); ilPoint++) {
+ const auto& pointLink = fsdDigiMatch->GetLink(ilPoint);
+ if (pointLink.GetFile() != file_id || pointLink.GetEntry() != event_id) { // match from different event
+ continue;
+ }
+
+ multMC++; // increase counter of MC multiplicity
+
+ // add matching between FSD hit and MC track
+ const auto* fsdPoint = dynamic_cast<FairMCPoint*>(cbm_fsd_points_new_->Get(pointLink));
+ Int_t mc_track_id = fsdPoint->GetTrackID();
+ if (mc_track_id >= 0) {
+ auto it = sim_tracks_map.find(mc_track_id);
+ if (it != sim_tracks_map.end()) { // match is found
+ fsd_hits_2_mc_tracks_->AddMatch(hit.GetId(), it->second);
+ }
+ }
+
+ // for the "matched" links store the difference between Hit and Point
+ if (fsdHitMatch->GetMatchedLink().GetIndex() == ilDigi
+ && fsdDigiMatch->GetMatchedLink().GetIndex() == ilPoint) {
+ hit.SetField(float(fsdPoint->GetTime() - time), i_dtHP_);
+ hit.SetField(float(fsdPoint->GetX() - hitX), i_dxHP_);
+ hit.SetField(float(fsdPoint->GetY() - hitY), i_dxHP_ + 1);
+ hit.SetField(float(fsdPoint->GetZ() - hitZ), i_dxHP_ + 2);
+ }
+
+ } // end of loop over links to points
+ }
+ } // end of loop over links to digis
+ } // end of sanity check of FSD hit matches
+
+ const Int_t hitMCmult = multMC;
+ const Float_t hitTopWeight = highestWeight;
+ hit.SetField(hitMCmult, i_multMC_);
+ hit.SetField(hitTopWeight, i_topW_);
+
+ // now matching with Global (reco) tracks
+ Int_t bestMatchedIndex = -1;
+ Double_t bestChi2 = 0.;
+ for (Int_t igt = 0; igt < n_tracks; igt++) {
+ const auto trackIndex = event ? event->GetIndex(ECbmDataType::kGlobalTrack, igt) : igt;
+ const auto* globalTrack = static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(trackIndex));
+ FairTrackParam param_last = *(globalTrack->GetParamLast());
+
+ Double_t matchingChi2 = Chi2FsdhitGtrack(fsdHit, param_last);
+
+ if (bestChi2 > matchingChi2 || bestMatchedIndex < 0) {
+ bestChi2 = matchingChi2;
+ bestMatchedIndex = trackIndex;
+ }
+ } // end of loop over GlobalTracks
+
+ const Int_t storedIndex = bestMatchedIndex;
+ const Float_t storedChi2 = static_cast<Float_t>(bestChi2);
+ hit.SetField(storedChi2, i_chi2_);
+ hit.SetField(storedIndex, i_bestMatchedGTrack_);
+
+ if ((bestMatchedIndex >= 0) && (bestChi2 > fsdgtrack_minChi2_) && (bestChi2 < fsdgtrack_maxChi2_)) {
+ const auto trackIndex = event ? event->GetIndex(ECbmDataType::kGlobalTrack, bestMatchedIndex) : bestMatchedIndex;
+ const auto* globalTrack = static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(trackIndex));
+ FairTrackParam param_last = *(globalTrack->GetParamLast());
+ FairTrackParam param_fsd = ExtrapolateGtrack(hitZ, param_last);
+
+ TVector3 p_fsd;
+ param_fsd.Momentum(p_fsd);
+
+ const Float_t p = p_fsd.Mag();
+ const Int_t q = param_fsd.GetQp() > 0 ? 1 : -1;
+ const Float_t l = globalTrack->GetLength(); // l is calculated by global tracking
+ const Float_t beta = event ? l / ((time - event->GetTzero()) * 29.9792458) : 0;
+ const Float_t m2 = event ? p * p * (1. / (beta * beta) - 1.) : -1.;
+
+ hit.SetField(m2, i_mass2_);
+ hit.SetField(float(q * p), i_qp_);
+ hit.SetField(float(param_fsd.GetX() - hitX), i_dx_);
+ hit.SetField(float(param_fsd.GetY() - hitY), i_dx_ + 1);
+ hit.SetField(float(param_fsd.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_fsd_->AddMatch(rec_tracks_map.find(stsTrackIndex)->second, hit.GetId());
+ }
+ }
+ } // end of loop over FSD hits
+}
+
+
+CbmFsdHitsConverter::~CbmFsdHitsConverter()
+{
+ delete fsd_hits_;
+ delete vtx_tracks_2_fsd_;
+};
diff --git a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h
new file mode 100644
index 0000000000..2e3ca7abb3
--- /dev/null
+++ b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h
@@ -0,0 +1,94 @@
+/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Lukas Chlad [committer] */
+
+#ifndef ANALYSIS_TREE_FSDHITSCONVERTER_H
+#define ANALYSIS_TREE_FSDHITSCONVERTER_H
+
+#include "AnalysisTree/Detector.hpp"
+#include "CbmConverterTask.h"
+
+class TClonesArray;
+class FairTrackParam;
+class CbmMCDataManager;
+class CbmMCDataArray;
+class CbmDigiManager;
+class CbmFsdHit;
+
+namespace AnalysisTree
+{
+ class Matching;
+}
+
+class CbmFsdHitsConverter final : public CbmConverterTask {
+ public:
+ explicit CbmFsdHitsConverter(std::string out_branch_name, std::string match_to = "")
+ : CbmConverterTask(std::move(out_branch_name), std::move(match_to)){};
+
+ ~CbmFsdHitsConverter() final;
+
+ void Init() final;
+ void ProcessData(CbmEvent* event) final;
+ void Finish() final{};
+
+ Double_t GetMinChi2GtrackHit() { return fsdgtrack_minChi2_; };
+ Double_t GetMaxChi2GtrackHit() { return fsdgtrack_maxChi2_; };
+ void SetMinChi2GtrackHit(Double_t chi2) { fsdgtrack_minChi2_ = chi2; };
+ void SetMaxChi2GtrackHit(Double_t chi2) { fsdgtrack_maxChi2_ = chi2; };
+
+ private:
+ FairTrackParam ExtrapolateGtrack(Double_t zpos, FairTrackParam params);
+ Double_t Chi2FsdhitGtrack(CbmFsdHit* hit, FairTrackParam inputParams);
+
+ const std::map<int, int>& GetMatchMap(const std::string& name) const
+ {
+ const auto& it = indexes_map_->find(name);
+ if (it == indexes_map_->end()) {
+ throw std::runtime_error(name + " is not found to match with FSD hits");
+ }
+ return it->second;
+ }
+ std::string mc_tracks_{"SimParticles"};
+
+ TClonesArray* cbm_global_tracks_{nullptr};
+ TClonesArray* cbm_sts_tracks_{nullptr};
+ TClonesArray* cbm_tof_hits_{nullptr};
+ TClonesArray* cbm_fsd_hits_{nullptr};
+ TClonesArray* cbm_fsd_hitmatch_{nullptr};
+ TClonesArray* cbm_mc_tracks_{nullptr};
+
+ CbmDigiManager* fDigiMan{nullptr};
+ //TClonesArray* cbm_fsd_digimatch_ {nullptr};
+ //TClonesArray* cbm_fsd_digis_ {nullptr};
+
+ CbmMCDataManager* cbm_mc_manager_{nullptr};
+ //CbmMCDataArray* cbm_mc_tracks_new_ {nullptr};
+ CbmMCDataArray* cbm_fsd_points_new_{nullptr};
+
+ AnalysisTree::HitDetector* fsd_hits_{nullptr};
+ AnalysisTree::Matching* vtx_tracks_2_fsd_{
+ nullptr}; // Matching of Reconstructed Global Tracks (extrapolated via CbmKFTrack to fixed Z) to FSD hits, the matching can be configured by cutting on chi2
+ AnalysisTree::Matching* fsd_hits_2_mc_tracks_{
+ nullptr}; // Matching of MC Tracks with highest weight (largest energy deposition) to FSD hits
+
+ int i_mass2_{AnalysisTree::UndefValueInt};
+ int i_qp_{AnalysisTree::UndefValueInt};
+ int i_dx_{AnalysisTree::UndefValueInt};
+ int i_t_{AnalysisTree::UndefValueInt};
+ int i_l_{AnalysisTree::UndefValueInt};
+ int i_edep_{AnalysisTree::UndefValueInt};
+ int i_chi2_{AnalysisTree::UndefValueInt};
+ int i_bestMatchedGTrack_{AnalysisTree::UndefValueInt};
+ int i_multMC_{AnalysisTree::UndefValueInt};
+ int i_topW_{AnalysisTree::UndefValueInt};
+ int i_dxHP_{AnalysisTree::UndefValueInt};
+ int i_dtHP_{AnalysisTree::UndefValueInt};
+
+
+ Double_t fsdgtrack_minChi2_{0.};
+ Double_t fsdgtrack_maxChi2_{0.};
+
+ ClassDef(CbmFsdHitsConverter, 1)
+};
+
+#endif // ANALYSIS_TREE_FSDHITSCONVERTER_H
diff --git a/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.cxx b/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.cxx
new file mode 100644
index 0000000000..1609967b0e
--- /dev/null
+++ b/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.cxx
@@ -0,0 +1,69 @@
+/* Copyright (C) 2020-2021 Physikalisches Institut, Eberhard Karls Universität Tuebingen, Tuebingen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Viktor Klochkov [committer] */
+
+#include "CbmFsdModulesConverter.h"
+
+#include "AnalysisTree/Detector.hpp"
+#include "CbmDefs.h"
+#include "CbmEvent.h"
+#include "CbmFsdHit.h"
+#include "FairRootManager.h"
+#include "TClonesArray.h"
+
+#include <AnalysisTree/TaskManager.hpp>
+#include <cassert>
+#include <vector>
+
+ClassImp(CbmFsdModulesConverter);
+
+void CbmFsdModulesConverter::Init()
+{
+ assert(!out_branch_.empty());
+ auto* ioman = FairRootManager::Instance();
+ assert(ioman != nullptr);
+ cbm_fsd_hits_ = (TClonesArray*) ioman->GetObject("FsdHit");
+
+ AnalysisTree::BranchConfig fsd_branch(out_branch_, AnalysisTree::DetType::kModule);
+
+ auto* man = AnalysisTree::TaskManager::GetInstance();
+ man->AddBranch(fsd_modules_, fsd_branch);
+}
+
+
+void CbmFsdModulesConverter::ProcessData(CbmEvent* event)
+{
+ assert(cbm_fsd_hits_);
+ fsd_modules_->ClearChannels();
+
+ CbmFsdHit* hit{nullptr};
+
+ auto* data_header = AnalysisTree::TaskManager::GetInstance()->GetDataHeader();
+ auto* config = AnalysisTree::TaskManager::GetInstance()->GetConfig();
+ const auto& branch = config->GetBranchConfig(out_branch_);
+
+ const int n_fsd_modules = data_header->GetModulePositions(0).GetNumberOfChannels();
+
+ fsd_modules_->Reserve(n_fsd_modules);
+ for (int i = 0; i < n_fsd_modules; ++i) {
+ auto& module = fsd_modules_->AddChannel(branch);
+ }
+
+ const int nFsdHits = event ? event->GetNofData(ECbmDataType::kFsdHit) : cbm_fsd_hits_->GetEntriesFast();
+ if (nFsdHits <= 0) {
+ LOG(warn) << "No FSD hits!";
+ return;
+ }
+
+ for (int i = 0; i < nFsdHits; ++i) {
+ hit = (CbmFsdHit*) cbm_fsd_hits_->At(i);
+ if (hit == nullptr) continue;
+ auto& module = fsd_modules_->Channel(hit->GetModuleId());
+ module.SetNumber(i);
+ module.SetSignal(hit->GetEdep());
+ }
+}
+
+void CbmFsdModulesConverter::Finish() {}
+
+CbmFsdModulesConverter::~CbmFsdModulesConverter() { delete fsd_modules_; };
diff --git a/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.h b/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.h
new file mode 100644
index 0000000000..65a42d3267
--- /dev/null
+++ b/analysis/common/analysis_tree_converter/CbmFsdModulesConverter.h
@@ -0,0 +1,31 @@
+/* Copyright (C) 2020-2021 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Viktor Klochkov [committer] */
+
+#ifndef ANALYSIS_TREE_FSDMODULESCONVERTER_H_
+#define ANALYSIS_TREE_FSDMODULESCONVERTER_H_
+
+#include "AnalysisTree/Detector.hpp"
+#include "CbmConverterTask.h"
+
+class TClonesArray;
+
+class CbmFsdModulesConverter final : public CbmConverterTask {
+ public:
+ explicit CbmFsdModulesConverter(std::string out_branch_name, std::string match_to = "")
+ : CbmConverterTask(std::move(out_branch_name), std::move(match_to)){};
+
+ ~CbmFsdModulesConverter() final;
+
+ void Init() final;
+ void ProcessData(CbmEvent* event) final;
+ void Finish() final;
+
+ private:
+ AnalysisTree::ModuleDetector* fsd_modules_{nullptr};
+ TClonesArray* cbm_fsd_hits_{nullptr};
+
+ ClassDef(CbmFsdModulesConverter, 1)
+};
+
+#endif // ANALYSIS_TREE_FSDMODULESCONVERTER_H_
diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx
index 14db617c5a..e3121c71cb 100644
--- a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx
+++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx
@@ -50,7 +50,6 @@ void CbmPsdModulesConverter::ProcessData(CbmEvent* event)
psd_modules_->Reserve(n_psd_modules);
for (int i = 0; i < n_psd_modules; ++i) {
auto& module = psd_modules_->AddChannel(branch);
- module.SetSignal(0.f);
}
const int nPsdHits = event ? event->GetNofData(ECbmDataType::kPsdHit) : cbm_psd_hits_->GetEntriesFast();
@@ -62,8 +61,8 @@ void CbmPsdModulesConverter::ProcessData(CbmEvent* event)
for (int i = 0; i < nPsdHits; ++i) {
hit = (CbmPsdHit*) cbm_psd_hits_->At(i);
if (hit == nullptr) continue;
- auto& module = psd_modules_->Channel(hit->GetModuleID() - 1);
- module.SetNumber(i + 1);
+ auto& module = psd_modules_->Channel(i);
+ module.SetNumber(hit->GetModuleID());
module.SetSignal(hit->GetEdep());
}
}
diff --git a/analysis/common/analysis_tree_converter/config/ATConfig_event_ideal.yaml b/analysis/common/analysis_tree_converter/config/ATConfig_event_ideal.yaml
index 16b6c3683d..e1a0779f72 100644
--- a/analysis/common/analysis_tree_converter/config/ATConfig_event_ideal.yaml
+++ b/analysis/common/analysis_tree_converter/config/ATConfig_event_ideal.yaml
@@ -7,3 +7,6 @@ global:
timeslicelength: 1.e5
collisionSystem: "Au+Au"
beamMomentum: 12.
+fsdhits:
+ gtrackmatch_minChi2: -1
+ gtrackmatch_maxChi2: 10000
diff --git a/analysis/common/analysis_tree_converter/config/ATConfig_timeslice.yaml b/analysis/common/analysis_tree_converter/config/ATConfig_timeslice.yaml
index 64e981d91d..a447f5817a 100644
--- a/analysis/common/analysis_tree_converter/config/ATConfig_timeslice.yaml
+++ b/analysis/common/analysis_tree_converter/config/ATConfig_timeslice.yaml
@@ -7,3 +7,6 @@ global:
timeslicelength: 1.e5
collisionSystem: "Au+Au"
beamMomentum: 12.
+fsdhits:
+ gtrackmatch_minChi2: -1
+ gtrackmatch_maxChi2: 10000
diff --git a/analysis/common/analysis_tree_converter/steer/Config.cxx b/analysis/common/analysis_tree_converter/steer/Config.cxx
index f8f12b6664..6105e46600 100644
--- a/analysis/common/analysis_tree_converter/steer/Config.cxx
+++ b/analysis/common/analysis_tree_converter/steer/Config.cxx
@@ -38,6 +38,10 @@ namespace cbm::atconverter
f_glb_tslength = settings["global"]["timeslicelength"].as<float>();
f_glb_system = settings["global"]["collisionSystem"].as<string>();
f_glb_beamMom = settings["global"]["beamMomentum"].as<float>();
+
+ // --- Fsd hits converter settings
+ f_fsd_minChi2match = settings["fsdhits"]["gtrackmatch_minChi2"].as<double>();
+ f_fsd_maxChi2match = settings["fsdhits"]["gtrackmatch_maxChi2"].as<double>();
}
// ----------------------------------------------------------------------------
@@ -83,6 +87,10 @@ namespace cbm::atconverter
settings["global"]["collisionSystem"] = f_glb_system;
settings["global"]["beamMomentum"] = f_glb_beamMom;
+ // --- Fsd hits converter settings
+ settings["fsdhits"]["gtrackmatch_minChi2"] = f_fsd_minChi2match;
+ settings["fsdhits"]["gtrackmatch_maxChi2"] = f_fsd_maxChi2match;
+
return settings;
}
// ----------------------------------------------------------------------------
diff --git a/analysis/common/analysis_tree_converter/steer/Config.h b/analysis/common/analysis_tree_converter/steer/Config.h
index 6406533236..ad9d32c7cd 100644
--- a/analysis/common/analysis_tree_converter/steer/Config.h
+++ b/analysis/common/analysis_tree_converter/steer/Config.h
@@ -94,6 +94,10 @@ namespace cbm::atconverter
// --- Event builder
ECbmEvbuildType f_evbuild_type = ECbmEvbuildType::Undefined;
+
+ // --- Fsd hits converter settings
+ float f_fsd_minChi2match = -1.;
+ float f_fsd_maxChi2match = 10000.;
};
} // namespace cbm::atconverter
diff --git a/analysis/common/analysis_tree_converter/steer/Run.cxx b/analysis/common/analysis_tree_converter/steer/Run.cxx
index 2b300e7f36..70a5aeac6b 100644
--- a/analysis/common/analysis_tree_converter/steer/Run.cxx
+++ b/analysis/common/analysis_tree_converter/steer/Run.cxx
@@ -65,6 +65,33 @@ namespace cbm::atconverter
// --------------------------------------------------------------------------
+ // ----- Check existence of reco branch ---------------------------------
+ void Run::CheckRecoBranch(TTree* tree, ECbmModuleId detector)
+ {
+ TString branchName = ToString(detector);
+ branchName += "Hit";
+ if (tree->GetBranchStatus(branchName.Data())) {
+ LOG(info) << GetName() << ": Found branch " << branchName;
+ fDataPresent.insert(detector);
+ }
+ }
+ // --------------------------------------------------------------------------
+
+
+ // ----- Check which input reco branches are present --------------------
+ void Run::CheckInputBranches(FairFileSource* source)
+ {
+ TFile* inFile = source->GetInFile();
+ if (!inFile) throw std::runtime_error("No input file");
+ TTree* inTree = dynamic_cast<TTree*>(inFile->Get("cbmsim"));
+ if (!inTree) throw std::runtime_error("No input tree");
+
+ for (ECbmModuleId detector = ECbmModuleId::kMvd; detector != ECbmModuleId::kNofSystems; ++detector)
+ CheckRecoBranch(inTree, detector);
+ }
+ // --------------------------------------------------------------------------
+
+
// ----- Create the topology --------------------------------------------
void Run::CreateTopology()
{
@@ -113,6 +140,10 @@ namespace cbm::atconverter
fRun.SetSource(source);
fRun.SetSink(new FairRootFileSink(fOutput));
+ // --- Check presence of input (reco) branches
+ LOG(info) << GetName() << ": Checking reco input...";
+ CheckInputBranches(source);
+
// --- Geometry setup
LOG(info) << GetName() << ": Loading setup " << fSetupTag;
fSetup = CbmSetup::Instance();
diff --git a/analysis/common/analysis_tree_converter/steer/Run.h b/analysis/common/analysis_tree_converter/steer/Run.h
index a9f95a19c4..2c5b10b9df 100644
--- a/analysis/common/analysis_tree_converter/steer/Run.h
+++ b/analysis/common/analysis_tree_converter/steer/Run.h
@@ -11,16 +11,16 @@
#define CBM_ATCONVERTER_STEER_RUN_H 1
#include "CbmDefs.h"
+#include "Config.h"
#include <FairRunAna.h>
#include <TNamed.h>
#include <TString.h>
+#include <TTree.h>
#include <string>
-#include "Config.h"
-
class CbmSetup;
class FairTask;
@@ -128,6 +128,17 @@ namespace cbm::atconverter
bool CheckFile(const char* fileName);
+ /** @brief Check and mark presence of reco branches
+ ** @param tree Pointer to ROOT tree
+ ** @param detector ECbmModuleId
+ **/
+ void CheckRecoBranch(TTree* tree, ECbmModuleId detector);
+
+
+ /** @brief Check the presence of reco input branches **/
+ void CheckInputBranches(FairFileSource* source);
+
+
/** @brief Create the reconstruction task topology (chain) **/
void CreateTopology();
diff --git a/analysis/common/analysis_tree_converter/steer/TaskFactory.cxx b/analysis/common/analysis_tree_converter/steer/TaskFactory.cxx
index 1ccc65e2fe..519a8f933b 100644
--- a/analysis/common/analysis_tree_converter/steer/TaskFactory.cxx
+++ b/analysis/common/analysis_tree_converter/steer/TaskFactory.cxx
@@ -10,6 +10,8 @@
#include "TaskFactory.h"
#include "CbmConverterManager.h"
+#include "CbmFsdHitsConverter.h"
+#include "CbmFsdModulesConverter.h"
#include "CbmKF.h"
#include "CbmL1.h"
#include "CbmMCDataManager.h"
@@ -83,10 +85,19 @@ namespace cbm::atconverter
ststracksconverter->SetIsReproduceCbmKFPF();
man->AddTask(ststracksconverter);
- man->AddTask(new CbmRichRingsConverter("RichRings", "VtxTracks"));
+ if (fRun->IsDataPresent(ECbmModuleId::kRich)) man->AddTask(new CbmRichRingsConverter("RichRings", "VtxTracks"));
man->AddTask(new CbmTofHitsConverter("TofHits", "VtxTracks"));
man->AddTask(new CbmTrdTracksConverter("TrdTracks", "VtxTracks"));
- man->AddTask(new CbmPsdModulesConverter("PsdModules"));
+ if (fRun->IsDataPresent(ECbmModuleId::kPsd)) man->AddTask(new CbmPsdModulesConverter("PsdModules"));
+
+ if (fRun->IsDataPresent(ECbmModuleId::kFsd)) {
+ man->AddTask(new CbmFsdModulesConverter("FsdModules"));
+
+ CbmFsdHitsConverter* fsdHitsConverter = new CbmFsdHitsConverter("FsdHits", "VtxTracks");
+ fsdHitsConverter->SetMinChi2GtrackHit(fRun->fConfig.f_fsd_minChi2match);
+ fsdHitsConverter->SetMaxChi2GtrackHit(fRun->fConfig.f_fsd_maxChi2match);
+ man->AddTask(fsdHitsConverter);
+ }
fRun->AddTask(man);
}
diff --git a/macro/PWG/common/production/run_analysis_tree_maker_json_config.C b/macro/PWG/common/production/run_analysis_tree_maker_json_config.C
index 2b4dfa2985..4707b43c1a 100644
--- a/macro/PWG/common/production/run_analysis_tree_maker_json_config.C
+++ b/macro/PWG/common/production/run_analysis_tree_maker_json_config.C
@@ -158,7 +158,15 @@ void run_analysis_tree_maker_json_config(TString traPath = "test", TString rawPa
man->AddTask(new CbmRichRingsConverter("RichRings", "VtxTracks"));
man->AddTask(new CbmTofHitsConverter("TofHits", "VtxTracks"));
man->AddTask(new CbmTrdTracksConverter("TrdTracks", "VtxTracks"));
- man->AddTask(new CbmPsdModulesConverter("PsdModules"));
+ if (setup->IsActive(ECbmModuleId::kPsd)) man->AddTask(new CbmPsdModulesConverter("PsdModules"));
+ if (setup->IsActive(ECbmModuleId::kFsd)) {
+ man->AddTask(new CbmFsdModulesConverter("FsdModules"));
+
+ CbmFsdHitsConverter* taskCbmFsdHitsConverter = new CbmFsdHitsConverter("FsdHits", "VtxTracks");
+ taskCbmFsdHitsConverter->SetMinChi2GtrackHit(-1.);
+ taskCbmFsdHitsConverter->SetMaxChi2GtrackHit(10000.);
+ man->AddTask(taskCbmFsdHitsConverter);
+ }
run->AddTask(man);
--
GitLab