diff --git a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx index 49ce439ec8211e62e0c0c90bb993a9dd115ee730..3b62e5cd8083836e1682345dc0e24315d4fd3e31 100644 --- a/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmPsdModulesConverter.cxx @@ -56,8 +56,6 @@ void CbmPsdModulesConverter::Exec() module.SetSignal(hit->GetEdep()); psd_energy += hit->GetEdep(); } - // rec_event_header_ -> SetField(psd_energy, - // config_.GetBranchConfig(rec_event_header_->GetId() ).GetFieldId("Epsd")); } void CbmPsdModulesConverter::Finish() {} diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx index 5e01f623c34ad3e09f10c85033d07210a245aca1..ccfbb646d9cdaa7f6286ff6785a534c1cadcc9ad 100644 --- a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.cxx @@ -1,5 +1,6 @@ #include "CbmRecEventHeaderConverter.h" +#include "CbmPsdHit.h" #include "CbmVertex.h" #include "FairMCEventHeader.h" @@ -21,12 +22,14 @@ ClassImp(CbmRecEventHeaderConverter) assert(ioman != nullptr); cbm_header_ = (FairMCEventHeader*) ioman->GetObject("MCEventHeader."); cbm_prim_vertex_ = (CbmVertex*) ioman->GetObject("PrimaryVertex."); + cbm_sts_tracks_ = (TClonesArray*) ioman->GetObject("StsTrack"); + cbm_psd_hits_ = (TClonesArray*) ioman->GetObject("PsdHit"); // ***** RecEventHeader ******* AnalysisTree::BranchConfig RecEventHeaderBranch("RecEventHeader", AnalysisTree::DetType::kEventHeader); - RecEventHeaderBranch.AddField<float>("vtx_chi2", "primiry vertex fit chi2"); + RecEventHeaderBranch.AddField<float>("vtx_chi2", "primiry vertex fit chi^2/NDF"); RecEventHeaderBranch.AddField<float>("Epsd", "GeV, full energy deposit in PSD"); - RecEventHeaderBranch.AddField<int>("M", "total multiplicity in STS"); + RecEventHeaderBranch.AddField<int>("M", "total multiplicity in STS(+MVD)"); RecEventHeaderBranch.AddField<int>("evt_id", "event identifier"); auto* man = AnalysisTree::TaskManager::GetInstance(); @@ -40,18 +43,25 @@ void CbmRecEventHeaderConverter::Exec() 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.), branch.GetFieldId("vtx_chi2")); - return; - } + if (!cbm_prim_vertex_) { throw std::runtime_error("No fPrimVtx"); } 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_sts_tracks_->GetEntries()), branch.GetFieldId("M")); + rec_event_header_->SetField(GetPsdEnergy(), branch.GetFieldId("Epsd")); rec_event_header_->SetField(int(cbm_header_->GetEventID()), branch.GetFieldId("evt_id")); rec_event_header_->SetField(float(cbm_prim_vertex_->GetChi2() / cbm_prim_vertex_->GetNDF()), branch.GetFieldId("vtx_chi2")); } + +float CbmRecEventHeaderConverter::GetPsdEnergy() +{ + //TODO avoid duplicating the code with PsdModulesConverter + float psd_energy {0.f}; + const int nPsdHits = cbm_psd_hits_->GetEntriesFast(); + for (int i = 0; i < nPsdHits; ++i) { + auto* hit = (CbmPsdHit*) cbm_psd_hits_->At(i); + if (hit == nullptr) continue; + psd_energy += hit->GetEdep(); + } + return psd_energy; +} \ No newline at end of file diff --git a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h index 51e56b33e72c7be8c061818ce64d603e8317d9e4..53ef4c641a7af955332883dfeb081ad45081cb14 100644 --- a/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h +++ b/analysis/common/analysis_tree_converter/CbmRecEventHeaderConverter.h @@ -7,6 +7,7 @@ class FairMCEventHeader; class CbmVertex; +class TClonesArray; class CbmRecEventHeaderConverter final : public CbmConverterTask { public: @@ -18,10 +19,14 @@ public: void Finish() final { delete rec_event_header_; }; private: + float GetPsdEnergy(); + AnalysisTree::EventHeader* rec_event_header_ {nullptr}; - FairMCEventHeader* cbm_header_ {nullptr}; - CbmVertex* cbm_prim_vertex_ {nullptr}; + TClonesArray* cbm_psd_hits_ {nullptr}; + TClonesArray* cbm_sts_tracks_ {nullptr}; ///< non-owning pointer + FairMCEventHeader* cbm_header_ {nullptr}; ///< non-owning pointer + CbmVertex* cbm_prim_vertex_ {nullptr}; ///< non-owning pointer ClassDef(CbmRecEventHeaderConverter, 1) }; diff --git a/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx index 974236dfa0707da80db4e9855657abfeeb17bc09..893f2d6e1cfa84d6ddf6584716f39e3e5c05987c 100644 --- a/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmRichRingsConverter.cxx @@ -1,6 +1,7 @@ #include "CbmRichRingsConverter.h" #include <CbmGlobalTrack.h> +#include <rich/CbmRichRing.h> #include <FairMCPoint.h> #include <FairRootManager.h> @@ -11,8 +12,6 @@ #include <cassert> #include <vector> -#include <rich/CbmRichRing.h> - #include "AnalysisTree/Matching.hpp" void CbmRichRingsConverter::Init() diff --git a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx index ec2f2d6959d6b2c08bdb29c5ca4869eb38d5af45..4a3eb5a86745f3aaed13c2cc9e83d8a7b31d7a3b 100644 --- a/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmStsTracksConverter.cxx @@ -4,6 +4,8 @@ #include "CbmStsTrack.h" #include "CbmTrackMatchNew.h" #include "CbmVertex.h" +#include "Interface/CbmKFVertex.h" +#include "ParticleFinder/CbmL1PFFitter.h" #include "FairRootManager.h" @@ -15,9 +17,7 @@ #include <cmath> #include "AnalysisTree/Matching.hpp" -#include "Interface/CbmKFVertex.h" #include "L1Field.h" -#include "ParticleFinder/CbmL1PFFitter.h" ClassImp(CbmStsTracksConverter) @@ -90,7 +90,6 @@ void CbmStsTracksConverter::Init() // TODO misleading name, move field filling somewhere else? float CbmStsTracksConverter::ExtrapolateToVertex(CbmStsTrack* sts_track, AnalysisTree::Track& track, int pdg) { - vector<CbmStsTrack> tracks = {*sts_track}; CbmL1PFFitter fitter; vector<float> chi2_to_vtx; diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx index 341bb7b3b4eb23de8dccff885d38163afb77985b..7ffe8a942c87be5e3170e9ca62a5862e58a717ea 100644 --- a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.cxx @@ -1,6 +1,7 @@ #include "CbmTofHitsConverter.h" #include <CbmGlobalTrack.h> +#include <CbmMCTrack.h> #include <CbmTofHit.h> #include <CbmTrackMatchNew.h> @@ -15,20 +16,18 @@ #include "AnalysisTree/Matching.hpp" -ClassImp(CbmTofHitsConverter) +ClassImp(CbmTofHitsConverter); - void CbmTofHitsConverter::Init() +void CbmTofHitsConverter::Init() { - assert(!out_branch_.empty()); auto* ioman = FairRootManager::Instance(); cbm_tof_hits_ = (TClonesArray*) ioman->GetObject("TofHit"); cbm_global_tracks_ = (TClonesArray*) ioman->GetObject("GlobalTrack"); - // cbm_tof_match_ = (TClonesArray*) ioman->GetObject("TofHitMatch"); - // cbm_tof_points_ = (TClonesArray*) ioman->GetObject("TofPoint"); - // cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); - // cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch"); + cbm_tof_match_ = (TClonesArray*) ioman->GetObject("TofHitMatch"); + cbm_tof_points_ = (TClonesArray*) ioman->GetObject("TofPoint"); + cbm_mc_tracks_ = (TClonesArray*) ioman->GetObject("MCTrack"); AnalysisTree::BranchConfig tof_branch(out_branch_, AnalysisTree::DetType::kHit); tof_branch.AddField<float>("mass2", "Mass squared"); @@ -37,16 +36,15 @@ ClassImp(CbmTofHitsConverter) tof_branch.AddField<float>("qp_tof", "charge * momentum extrapoleted to TOF"); tof_branch.AddFields<float>({"dx", "dy", "dz"}, "Distance between TOF hit and extrapolated global track, cm"); tof_branch.AddField<int>("mc_pdg", "MC-true PDG code of particle with highest contribution to TOF hit"); - tof_branch.AddField<bool>("is_correct_match", "is the matched track corresponds to MC-true track from TOF hit"); auto* man = AnalysisTree::TaskManager::GetInstance(); man->AddBranch(out_branch_, tof_hits_, tof_branch); man->AddMatching(match_to_, out_branch_, vtx_tracks_2_tof_); + man->AddMatching(out_branch_, mc_tracks_, tof_hits_2_mc_tracks_); } void CbmTofHitsConverter::ExtrapolateStraightLine(FairTrackParam* params, float z) { - const Float_t Tx = params->GetTx(); const Float_t Ty = params->GetTy(); const Float_t old_z = params->GetZ(); @@ -58,11 +56,13 @@ void CbmTofHitsConverter::ExtrapolateStraightLine(FairTrackParam* params, float params->SetPosition({x, y, z}); } + void CbmTofHitsConverter::FillTofHits() { assert(cbm_tof_hits_); tof_hits_->ClearChannels(); vtx_tracks_2_tof_->Clear(); + tof_hits_2_mc_tracks_->Clear(); auto* out_config_ = AnalysisTree::TaskManager::GetInstance()->GetConfig(); const auto& branch = out_config_->GetBranchConfig(out_branch_); @@ -73,9 +73,8 @@ void CbmTofHitsConverter::FillTofHits() 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()) { throw std::runtime_error(match_to_ + " is not found to match with TOF hits"); } - auto rec_tracks_map = it->second; + auto rec_tracks_map = GetMatchMap(match_to_); + auto sim_tracks_map = GetMatchMap(mc_tracks_); tof_hits_->Reserve(cbm_global_tracks_->GetEntriesFast()); @@ -120,28 +119,19 @@ void CbmTofHitsConverter::FillTofHits() vtx_tracks_2_tof_->AddMatch(rec_tracks_map.find(stsTrackIndex)->second, hit.GetId()); } - // const auto* tofMatch = - // dynamic_cast<CbmMatch*>(cbm_tof_match_->At(tofHitIndex)); if(tofMatch - // != nullptr && tofMatch->GetNofLinks()>0) { - // - // const auto* tofPoint = dynamic_cast<FairMCPoint*>( - // cbm_tof_points_->At(tofMatch->GetMatchedLink().GetIndex()) ); - // - // Int_t itofMC = (tofPoint ? tofPoint->GetTrackID() : -1 ); - // if(itofMC >= 0){ - // const auto* mc_track = dynamic_cast<const - // CbmMCTrack*>(cbm_mc_tracks_->At(itofMC)); - // hit.SetField(mc_track->GetPdgCode(), i_pdg); - // - // const Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); - // if(stsTrackIndex<0) return; - // - // auto* match = (CbmTrackMatchNew*) cbm_sts_match_->At(stsTrackIndex); - // if(match == nullptr || match->GetNofLinks() == 0) continue; - // const int mc_id_sts = match->GetMatchedLink().GetIndex(); - // hit.SetField( bool(mc_id_sts == itofMC), i_is_correct); - // } - // } + const auto* tofMatch = dynamic_cast<CbmMatch*>(cbm_tof_match_->At(tofHitIndex)); + if (tofMatch && tofMatch->GetNofLinks() > 0) { + const auto* tofPoint = dynamic_cast<FairMCPoint*>(cbm_tof_points_->At(tofMatch->GetMatchedLink().GetIndex())); + if (!tofPoint) { throw std::runtime_error("no TOF point"); } + + Int_t mc_track_id = tofPoint->GetTrackID(); + if (mc_track_id >= 0) { + auto it = sim_tracks_map.find(mc_track_id); + if (it != sim_tracks_map.end()) { // match is found + tof_hits_2_mc_tracks_->AddMatch(hit.GetId(), it->second); + } + } + } } } diff --git a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h index 796e679ef43fc535f6dfafbeefc8d2c2e390fed6..8b2d3c5065ebc4068b55a3ee623a62c3abb66acd 100644 --- a/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h +++ b/analysis/common/analysis_tree_converter/CbmTofHitsConverter.h @@ -16,28 +16,37 @@ namespace AnalysisTree 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)) {}; + : CbmConverterTask(std::move(out_branch_name), std::move(match_to)) + { + } ~CbmTofHitsConverter() final; - ; void Init() final; void Exec() final; - void Finish() final {}; + void Finish() final {} private: void FillTofHits(); static void ExtrapolateStraightLine(FairTrackParam* params, float z); + const std::map<int, int>& GetMatchMap(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 TOF hits"); } + return it->second; + } + std::string mc_tracks_ {"SimParticles"}; + TClonesArray* cbm_global_tracks_ {nullptr}; TClonesArray* cbm_tof_hits_ {nullptr}; - // TClonesArray *cbm_tof_points_{nullptr}; - // TClonesArray *cbm_tof_match_{nullptr}; - // TClonesArray* cbm_mc_tracks_{nullptr}; - // TClonesArray* cbm_sts_match_{nullptr}; + TClonesArray* cbm_tof_points_ {nullptr}; + TClonesArray* cbm_tof_match_ {nullptr}; + TClonesArray* cbm_mc_tracks_ {nullptr}; AnalysisTree::HitDetector* tof_hits_ {nullptr}; AnalysisTree::Matching* vtx_tracks_2_tof_ {nullptr}; + AnalysisTree::Matching* tof_hits_2_mc_tracks_ {nullptr}; ClassDef(CbmTofHitsConverter, 1) }; diff --git a/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx index 98125e44dbc4865bc9ac3bd414a78df0f89b11ba..2d393653fb987647ded86599e57d60409ea1e105 100644 --- a/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmTrdTracksConverter.cxx @@ -52,8 +52,7 @@ void CbmTrdTracksConverter::FillTrdTracks() 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_e_loss_i = branch.GetFieldId("energy_loss_0"); 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"); @@ -72,8 +71,6 @@ void CbmTrdTracksConverter::FillTrdTracks() 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; @@ -82,14 +79,11 @@ void CbmTrdTracksConverter::FillTrdTracks() 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);