/
TrackFastSimEval.cc
103 lines (75 loc) · 3.29 KB
/
TrackFastSimEval.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#include <iostream>
#include <TVector3.h>
#include <TH1D.h>
#include <trackbase_historic/SvtxTrackMap.h>
#include <trackbase_historic/SvtxTrack_FastSim.h>
#include <g4main/PHG4Particle.h>
#include <g4main/PHG4TruthInfoContainer.h>
#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/PHTFileServer.h>
#include <fun4all/SubsysReco.h>
#include <phool/getClass.h>
#include "TrackFastSimEval.h"
using namespace std;
// ---------------------------------------------------------------------------------------
TrackFastSimEval::TrackFastSimEval(const string &name, const string &filename,
const string &trackmapname)
: SubsysReco(name)
, _outfile_name(filename)
, _trackmapname(trackmapname)
, _event(0)
, _h1d_Delta_mom(nullptr)
{
} // TrackFastSimEval::TrackFastSimEval()
// ---------------------------------------------------------------------------------------
int TrackFastSimEval::Init(PHCompositeNode *topNode)
{
cout << PHWHERE << " Opening file " << _outfile_name << endl;
PHTFileServer::get().open(_outfile_name, "RECREATE");
_h1d_Delta_mom = new TH1D("_h1d_Delta_mom", "#frac{#Delta p}{truth p}", 100, -0.2, 0.2);
return Fun4AllReturnCodes::EVENT_OK;
} // TrackFastSimEval::Init()
// ---------------------------------------------------------------------------------------
int TrackFastSimEval::process_event(PHCompositeNode *topNode)
{
if (++_event % 100 == 0) cout << PHWHERE << "Events processed: " << _event << endl;
auto _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
auto _trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname);
if (!_truth_container || !_trackmap) {
cout << PHWHERE << " Either PHG4TruthInfoContainer or SvtxTrackMap node not found on node tree"
<< endl;
return Fun4AllReturnCodes::ABORTEVENT;
} //if
{
auto range = _truth_container->GetPrimaryParticleRange();
// Loop through all the truth particles;
for (auto truth_itr = range.first; truth_itr != range.second; ++truth_itr) {
auto g4particle = truth_itr->second;
if (!g4particle) continue;
// Loop through all the reconstructed particles and find a match; how about std::map?;
for (auto track_itr = _trackmap->begin(); track_itr != _trackmap->end(); track_itr++) {
auto track = dynamic_cast<SvtxTrack_FastSim *>(track_itr->second);
if (!track) {
std::cout << "ERROR CASTING PARTICLE!" << std::endl;
continue;
} //if
// Matching reconstructed particle found, use it and break;
if ((track->get_truth_track_id() - g4particle->get_track_id()) == 0) {
TVector3 truth_mom(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz());
TVector3 reco_mom ( track->get_px(), track->get_py(), track->get_pz());
_h1d_Delta_mom->Fill((reco_mom.Mag() - truth_mom.Mag()) / truth_mom.Mag());
break;
} //if
} //for
}
}
return Fun4AllReturnCodes::EVENT_OK;
} // TrackFastSimEval::process_event()
// ---------------------------------------------------------------------------------------
int TrackFastSimEval::End(PHCompositeNode *topNode)
{
PHTFileServer::get().cd(_outfile_name);
_h1d_Delta_mom->Write();
return Fun4AllReturnCodes::EVENT_OK;
} // TrackFastSimEval::End()
// ---------------------------------------------------------------------------------------