CbmRoot
CbmStsTracksConverter.cxx
Go to the documentation of this file.
1 
2 #include <cassert>
3 #include <cmath>
4 
5 #include "TClonesArray.h"
6 
7 #include "FairRootManager.h"
8 
9 #include "AnalysisTree/Matching.hpp"
10 
11 #include "CbmMCTrack.h"
12 #include "CbmStsTrack.h"
13 #include "CbmTrackMatchNew.h"
14 #include "CbmVertex.h"
15 #include "Interface/CbmKFVertex.h"
16 #include "L1Field.h"
17 #include "ParticleFinder/CbmL1PFFitter.h"
18 
19 #include "CbmStsTracksConverter.h"
20 
21 ClassImp(CbmStsTracksConverter)
22 
23  void CbmStsTracksConverter::Exec() {
24  assert(cbm_sts_tracks_ != nullptr && cbm_mc_tracks_ != nullptr);
25 
26  vtx_tracks_2_sim_->Clear();
27  out_indexes_map_.clear();
28 
29  ReadVertexTracks();
30  MapTracks();
31 }
32 
33 CbmStsTracksConverter::~CbmStsTracksConverter() {
34  delete vtx_tracks_;
35  delete vtx_tracks_2_sim_;
36 };
37 
38 // TODO misleading name, move field filling somewhere else?
39 float CbmStsTracksConverter::ExtrapolateToVertex(CbmStsTrack* sts_track,
40  AnalysisTree::Track* track,
41  int pdg) {
42 
43  vector<CbmStsTrack> tracks = {*sts_track};
44  CbmL1PFFitter fitter;
45  vector<float> chi2_to_vtx;
46  vector<L1FieldRegion> field;
47  CbmKFVertex kfVertex = CbmKFVertex(*cbm_prim_vertex_);
48  if (is_reproduce_cbmkfpf_) {
49  std::vector<int> pdgVector = {pdg};
50  fitter.Fit(tracks, pdgVector);
51  }
52  fitter.GetChiToVertex(tracks, field, chi2_to_vtx, kfVertex, 3.);
53  *sts_track = tracks[0];
54 
55  if (is_write_kfinfo_) {
56  track->SetField(float(field.at(0).cx0[0]), imf_);
57  track->SetField(float(field.at(0).cx1[0]), imf_ + 1);
58  track->SetField(float(field.at(0).cx2[0]), imf_ + 2);
59  track->SetField(float(field.at(0).cy0[0]), imf_ + 3);
60  track->SetField(float(field.at(0).cy1[0]), imf_ + 4);
61  track->SetField(float(field.at(0).cy2[0]), imf_ + 5);
62  track->SetField(float(field.at(0).cz0[0]), imf_ + 6);
63  track->SetField(float(field.at(0).cz1[0]), imf_ + 7);
64  track->SetField(float(field.at(0).cz2[0]), imf_ + 8);
65  track->SetField(float(field.at(0).z0[0]), imf_ + 9);
66  }
67 
68  return chi2_to_vtx[0];
69 }
70 
71 void CbmStsTracksConverter::ReadVertexTracks() {
72  assert(cbm_prim_vertex_ && cbm_sts_tracks_);
73 
74  vtx_tracks_->ClearChannels();
75 
76  const int iq =
77  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("q");
78  const int indf =
79  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("ndf");
80  const int ichi2 =
81  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("chi2");
82  const int inhits =
83  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("nhits");
84  const int inhits_mvd =
85  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("nhits_mvd");
86  const int idcax =
87  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("dcax");
88  const int ivtx_chi2 =
89  out_config_->GetBranchConfig(vtx_tracks_->GetId()).GetFieldId("vtx_chi2");
90 
91  const int n_sts_tracks = cbm_sts_tracks_->GetEntries();
92  vtx_tracks_->Reserve(n_sts_tracks);
93 
94  for (short i_track = 0; i_track < n_sts_tracks; ++i_track) {
95  auto* sts_track = (CbmStsTrack*) cbm_sts_tracks_->At(i_track);
96  if (!sts_track) { throw std::runtime_error("empty track!"); }
97 
98  auto* track = vtx_tracks_->AddChannel();
99  track->Init(out_config_->GetBranchConfig(vtx_tracks_->GetId()));
100 
101  int pdg = GetMcPid((CbmTrackMatchNew*) cbm_sts_match_->At(i_track), track);
102 
103  bool is_good_track = IsGoodCovMatrix(sts_track);
104 
105  float chi2_vertex = ExtrapolateToVertex(sts_track, track, pdg);
106  const FairTrackParam* trackParamFirst = sts_track->GetParamFirst();
107  TVector3 momRec;
108  trackParamFirst->Momentum(momRec);
109  const Int_t q = trackParamFirst->GetQp() > 0 ? 1 : -1;
110 
111  track->SetMomentum3(momRec);
112  track->SetField(int(q), iq);
113  track->SetField(int(sts_track->GetNDF()), indf);
114  track->SetField(float(sts_track->GetChiSq()), ichi2);
115  track->SetField(int(sts_track->GetNofHits()), inhits);
116  track->SetField(float(trackParamFirst->GetX() - cbm_prim_vertex_->GetX()),
117  idcax);
118  track->SetField(float(trackParamFirst->GetY() - cbm_prim_vertex_->GetY()),
119  idcax + 1);
120  track->SetField(float(trackParamFirst->GetZ() - cbm_prim_vertex_->GetZ()),
121  idcax + 2);
122  track->SetField(int(sts_track->GetNofMvdHits()), inhits_mvd);
123  track->SetField(float(chi2_vertex), ivtx_chi2);
124 
125  out_indexes_map_.insert(std::make_pair(i_track, track->GetId()));
126 
127  if (is_write_kfinfo_) { WriteKFInfo(track, sts_track, is_good_track); }
128  }
129 }
130 
131 void CbmStsTracksConverter::WriteKFInfo(AnalysisTree::Track* track,
132  const CbmStsTrack* sts_track,
133  bool is_good_track) const {
134  assert(track && sts_track);
135  const FairTrackParam* trackParamFirst = sts_track->GetParamFirst();
136 
137  track->SetField(float(trackParamFirst->GetX()), ipar_);
138  track->SetField(float(trackParamFirst->GetY()), ipar_ + 1);
139  track->SetField(float(trackParamFirst->GetZ()), ipar_ + 2);
140  track->SetField(float(trackParamFirst->GetTx()), ipar_ + 3);
141  track->SetField(float(trackParamFirst->GetTy()), ipar_ + 4);
142  track->SetField(float(trackParamFirst->GetQp()), ipar_ + 5);
143 
144  for (Int_t i = 0, iCov = 0; i < 5; i++) {
145  for (Int_t j = 0; j <= i; j++, iCov++) {
146  track->SetField(float(trackParamFirst->GetCovariance(i, j)),
147  icov_ + iCov);
148  }
149  }
150  track->SetField(is_good_track, ipasscuts_);
151 }
152 
153 bool CbmStsTracksConverter::IsGoodCovMatrix(
154  const CbmStsTrack* sts_track) const {
155 
156  if (!is_reproduce_cbmkfpf_) return true;
157 
158  assert(sts_track);
159  const FairTrackParam* trackParamFirst = sts_track->GetParamFirst();
160 
161  Double_t cov_matrix[15] = {0.f};
162  for (Int_t i = 0, iCov = 0; i < 5; i++) {
163  for (Int_t j = 0; j <= i; j++, iCov++) {
164  cov_matrix[iCov] = trackParamFirst->GetCovariance(i, j);
165  }
166  }
167  // Cuts, coded in MZ's CbmKFParticleFinder.cxx
168  bool ok = true;
169  ok = ok && isfinite(sts_track->GetParamFirst()->GetX());
170  ok = ok && isfinite(sts_track->GetParamFirst()->GetY());
171  ok = ok && isfinite(sts_track->GetParamFirst()->GetZ());
172  ok = ok && isfinite(sts_track->GetParamFirst()->GetTx());
173  ok = ok && isfinite(sts_track->GetParamFirst()->GetTy());
174  ok = ok && isfinite(sts_track->GetParamFirst()->GetQp());
175 
176  for (auto element : cov_matrix) {
177  ok = ok && isfinite(element);
178  }
179  ok = ok && (cov_matrix[0] < 1. && cov_matrix[0] > 0.)
180  && (cov_matrix[2] < 1. && cov_matrix[2] > 0.)
181  && (cov_matrix[5] < 1. && cov_matrix[5] > 0.)
182  && (cov_matrix[9] < 1. && cov_matrix[9] > 0.)
183  && (cov_matrix[14] < 1. && cov_matrix[14] > 0.);
184 
185  ok = ok && sts_track->GetChiSq() < 10. * sts_track->GetNDF();
186 
187  return ok;
188 }
189 
190 int CbmStsTracksConverter::GetMcPid(const CbmTrackMatchNew* match,
191  AnalysisTree::Track* track) const {
192 
193  if (!is_write_kfinfo_) { return -2; }
194  //-----------PID as in MZ's
195  //CbmKFPF----------------------------------------------------------
196  Int_t nMCTracks = cbm_mc_tracks_->GetEntriesFast();
197 
198  int pdgCode = -2;
199  if (match->GetNofLinks() > 0) { // there is at least one matched MC-track
200  Float_t bestWeight = 0.f;
201  Float_t totalWeight = 0.f;
202  Int_t mcTrackId = -1;
203 
204  for (int iLink = 0; iLink < match->GetNofLinks(); iLink++) {
205  totalWeight += match->GetLink(iLink).GetWeight();
206  if (match->GetLink(iLink).GetWeight() > bestWeight) {
207  bestWeight = match->GetLink(iLink).GetWeight();
208  mcTrackId = match->GetLink(iLink).GetIndex();
209  }
210  }
211 
212  if (!((bestWeight / totalWeight < 0.7)
213  || (mcTrackId >= nMCTracks || mcTrackId < 0))) {
214  auto* mctrack = static_cast<CbmMCTrack*>(cbm_mc_tracks_->At(mcTrackId));
215 
216  if (!(TMath::Abs(mctrack->GetPdgCode()) == 11
217  || TMath::Abs(mctrack->GetPdgCode()) == 13
218  || TMath::Abs(mctrack->GetPdgCode()) == 211
219  || TMath::Abs(mctrack->GetPdgCode()) == 321
220  || TMath::Abs(mctrack->GetPdgCode()) == 2212
221  || TMath::Abs(mctrack->GetPdgCode()) == 1000010020
222  || TMath::Abs(mctrack->GetPdgCode()) == 1000010030
223  || TMath::Abs(mctrack->GetPdgCode()) == 1000020030
224  || TMath::Abs(mctrack->GetPdgCode()) == 1000020040)) {
225  pdgCode = -1;
226  } else {
227  pdgCode = mctrack->GetPdgCode();
228  }
229  if (mctrack->GetMotherId() > -1) {
230  track->SetField(
231  int(((CbmMCTrack*) cbm_mc_tracks_->At(mctrack->GetMotherId()))
232  ->GetPdgCode()),
233  imother_pdg_);
234  }
235  }
236  }
237  track->SetField(pdgCode, imc_pdg_);
238 
239  return pdgCode;
240 }
241 
242 void CbmStsTracksConverter::InitInput() {
243  auto* ioman = FairRootManager::Instance();
244 
245  cbm_prim_vertex_ =
246  (CbmVertex*) ioman->GetObject(in_branches_.at(ePrimiryVertex).c_str());
247  cbm_sts_tracks_ =
248  (TClonesArray*) ioman->GetObject(in_branches_.at(eStsTracks).c_str());
249  cbm_mc_tracks_ =
250  (TClonesArray*) ioman->GetObject(in_branches_.at(eSimTracks).c_str());
251  cbm_sts_match_ = (TClonesArray*) ioman->GetObject("StsTrackMatch");
252 }
253 
254 void CbmStsTracksConverter::Init(map<std::string, void*>& branches) {
255  InitInput();
256 
257  AnalysisTree::BranchConfig vtx_tracks_config(out_branch_,
258  AnalysisTree::DetType::kTrack);
259  vtx_tracks_config.AddField<float>("chi2");
260  vtx_tracks_config.AddField<float>("vtx_chi2");
261  vtx_tracks_config.AddFields<float>({"dcax", "dcay", "dcaz"});
262  vtx_tracks_config.AddField<int>("ndf");
263  vtx_tracks_config.AddField<int>("q");
264  vtx_tracks_config.AddField<int>("nhits");
265  vtx_tracks_config.AddField<int>("nhits_mvd");
266 
267  if (is_write_kfinfo_) {
268  vtx_tracks_config.AddFields<float>({"x", "y", "z", "tx", "ty", "qp"});
269  vtx_tracks_config.AddFields<float>(
270  {"cx0", "cx1", "cx2", "cy0", "cy1", "cy2", "cz0", "cz1", "cz2", "z0"});
271  vtx_tracks_config.AddFields<float>({"cov1",
272  "cov2",
273  "cov3",
274  "cov4",
275  "cov5",
276  "cov6",
277  "cov7",
278  "cov8",
279  "cov9",
280  "cov10",
281  "cov11",
282  "cov12",
283  "cov13",
284  "cov14",
285  "cov15"});
286 
287  vtx_tracks_config.AddField<int>("mother_pdg");
288  vtx_tracks_config.AddField<int>("mc_pdg");
289  vtx_tracks_config.AddField<bool>("pass_cuts");
290 
291  ipar_ = vtx_tracks_config.GetFieldId("x");
292  imf_ = vtx_tracks_config.GetFieldId("cx0");
293  icov_ = vtx_tracks_config.GetFieldId("cov1");
294  imc_pdg_ = vtx_tracks_config.GetFieldId("mc_pdg");
295  imother_pdg_ = vtx_tracks_config.GetFieldId("mother_pdg");
296  ipasscuts_ = vtx_tracks_config.GetFieldId("pass_cuts");
297  }
298 
299  out_config_->AddBranchConfig(std::move(vtx_tracks_config));
300  vtx_tracks_ = new AnalysisTree::TrackDetector(out_config_->GetLastId());
301  vtx_tracks_2_sim_ = new AnalysisTree::Matching(
302  out_config_->GetLastId(), out_config_->GetBranchConfig(match_to_).GetId());
303 
304  out_config_->AddMatch(vtx_tracks_2_sim_);
305 
306  out_tree_->Branch(
307  out_branch_.c_str(), "AnalysisTree::TrackDetector", &vtx_tracks_);
308  out_tree_->Branch((out_branch_ + "2" + match_to_).c_str(),
309  "AnalysisTree::Matching",
310  &vtx_tracks_2_sim_);
311 
312  branches.emplace(out_branch_, vtx_tracks_);
313 }
314 
315 void CbmStsTracksConverter::MapTracks() {
316  const auto it = indexes_map_->find(match_to_);
317  if (it == indexes_map_->end()) {
318  throw std::runtime_error(match_to_
319  + " is not found to match with vertex tracks");
320  }
321  auto sim_tracks_map = it->second;
322 
323  if (out_indexes_map_.empty() || sim_tracks_map.empty()) return;
324  CbmTrackMatchNew* match {nullptr};
325 
326  for (const auto& track_id : out_indexes_map_) {
327  const int cbm_id = track_id.first;
328  const int out_id = track_id.second;
329  match = (CbmTrackMatchNew*) cbm_sts_match_->At(cbm_id);
330  if (match->GetNofLinks() > 0) { // there is at least one matched MC-track
331  const int mc_track_id = match->GetMatchedLink().GetIndex();
332 
333  if (match->GetTrueOverAllHitsRatio()
334  < 0.5) // not enough hits to be a match
335  continue;
336 
337  auto p = sim_tracks_map.find(mc_track_id);
338  if (p == sim_tracks_map.end()) // match is not found
339  continue;
340 
341  vtx_tracks_2_sim_->AddMatch(out_id, p->second);
342  }
343  }
344 }
CbmTrack::GetChiSq
Double_t GetChiSq() const
Definition: CbmTrack.h:58
CbmVertex.h
CbmStsTracksConverter::is_write_kfinfo_
bool is_write_kfinfo_
Definition: CbmStsTracksConverter.h:68
CbmStsTracksConverter.h
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmConverterTask::match_to_
std::string match_to_
AT branch to match.
Definition: CbmConverterTask.h:36
CbmStsTracksConverter::icov_
int icov_
Definition: CbmStsTracksConverter.h:73
L1Field.h
CbmL1PFFitter
Definition: CbmL1PFFitter.h:31
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, std::vector< int > &pidHypo)
Definition: CbmL1PFFitter.cxx:81
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmStsTracksConverter::IsGoodCovMatrix
bool IsGoodCovMatrix(const CbmStsTrack *sts_track) const
Definition: CbmStsTracksConverter.cxx:153
CbmStsTracksConverter::ePrimiryVertex
@ ePrimiryVertex
Definition: CbmStsTracksConverter.h:22
CbmStsTracksConverter::imc_pdg_
int imc_pdg_
Definition: CbmStsTracksConverter.h:74
CbmStsTracksConverter::~CbmStsTracksConverter
~CbmStsTracksConverter() final
Definition: CbmStsTracksConverter.cxx:33
CbmVertex::GetX
Double_t GetX() const
Definition: CbmVertex.h:68
CbmConverterTask::indexes_map_
std::map< std::string, std::map< int, int > > * indexes_map_
Definition: CbmConverterTask.h:34
CbmStsTrack.h
Data class for STS tracks.
CbmStsTracksConverter::vtx_tracks_
AnalysisTree::TrackDetector * vtx_tracks_
raw pointers are needed for TTree::Branch
Definition: CbmStsTracksConverter.h:58
CbmStsTracksConverter::ExtrapolateToVertex
float ExtrapolateToVertex(CbmStsTrack *sts_track, AnalysisTree::Track *track, int pdg)
Definition: CbmStsTracksConverter.cxx:39
ClassImp
ClassImp(CbmStsTracksConverter) void CbmStsTracksConverter
Definition: CbmStsTracksConverter.cxx:21
CbmStsTracksConverter
Definition: CbmStsTracksConverter.h:17
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition: CbmL1PFFitter.cxx:403
tracks
TClonesArray * tracks
Definition: Analyze_matching.h:17
CbmStsTracksConverter::cbm_sts_match_
TClonesArray * cbm_sts_match_
non-owning pointer
Definition: CbmStsTracksConverter.h:66
CbmTrackMatchNew.h
CbmVertex
Definition: CbmVertex.h:26
CbmL1PFFitter.h
nMCTracks
Int_t nMCTracks
Definition: CbmHadronAnalysis.cxx:63
CbmStsTracksConverter::eStsTracks
@ eStsTracks
Definition: CbmStsTracksConverter.h:21
CbmTrack::GetNDF
Int_t GetNDF() const
Definition: CbmTrack.h:59
CbmVertex::GetZ
Double_t GetZ() const
Definition: CbmVertex.h:70
CbmStsTracksConverter::cbm_mc_tracks_
TClonesArray * cbm_mc_tracks_
non-owning pointer
Definition: CbmStsTracksConverter.h:64
CbmStsTracksConverter::is_reproduce_cbmkfpf_
bool is_reproduce_cbmkfpf_
Definition: CbmStsTracksConverter.h:69
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmConverterTask::out_indexes_map_
std::map< int, int > out_indexes_map_
Definition: CbmConverterTask.h:31
CbmStsTracksConverter::imf_
int imf_
Definition: CbmStsTracksConverter.h:72
CbmStsTracksConverter::cbm_prim_vertex_
CbmVertex * cbm_prim_vertex_
non-owning pointer
Definition: CbmStsTracksConverter.h:63
CbmMCTrack.h
CbmVertex::GetY
Double_t GetY() const
Definition: CbmVertex.h:69
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmStsTracksConverter::ReadVertexTracks
void ReadVertexTracks()
Definition: CbmStsTracksConverter.cxx:71
CbmStsTracksConverter::vtx_tracks_2_sim_
AnalysisTree::Matching * vtx_tracks_2_sim_
raw pointers are needed for TTree::Branch
Definition: CbmStsTracksConverter.h:60
CbmStsTracksConverter::Exec
void Exec() final
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmStsTracksConverter::MapTracks
void MapTracks()
Definition: CbmStsTracksConverter.cxx:315
CbmStsTracksConverter::ipar_
int ipar_
Definition: CbmStsTracksConverter.h:71
CbmStsTracksConverter::WriteKFInfo
void WriteKFInfo(AnalysisTree::Track *track, const CbmStsTrack *sts_track, bool is_good_track) const
Definition: CbmStsTracksConverter.cxx:131
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmKFVertex.h
CbmStsTracksConverter::imother_pdg_
int imother_pdg_
Definition: CbmStsTracksConverter.h:75
CbmStsTracksConverter::InitInput
void InitInput()
Definition: CbmStsTracksConverter.cxx:242
CbmStsTracksConverter::ipasscuts_
int ipasscuts_
Definition: CbmStsTracksConverter.h:76
CbmStsTracksConverter::GetMcPid
int GetMcPid(const CbmTrackMatchNew *match, AnalysisTree::Track *track) const
Definition: CbmStsTracksConverter.cxx:190
CbmKFVertex
Definition: CbmKFVertex.h:6
CbmStsTracksConverter::eSimTracks
@ eSimTracks
Definition: CbmStsTracksConverter.h:23
CbmStsTracksConverter::Init
void Init(std::map< std::string, void * > &) final
Definition: CbmStsTracksConverter.cxx:254
CbmStsTracksConverter::cbm_sts_tracks_
TClonesArray * cbm_sts_tracks_
non-owning pointer
Definition: CbmStsTracksConverter.h:65