Belle II Software development
ECLTrackClusterMatchingPerformanceModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <ecl/modules/eclTrackClusterMatching/ECLTrackClusterMatchingPerformanceModule.h>
10
11#include <framework/datastore/RelationVector.h>
12#include <mdst/dataobjects/HitPatternCDC.h>
13#include <mdst/dataobjects/HitPatternVXD.h>
14#include <framework/gearbox/Const.h>
15
16#include <root/TFile.h>
17#include <root/TTree.h>
18
19using namespace Belle2;
20
21//-----------------------------------------------------------------
22// Register the Module
23//-----------------------------------------------------------------
24REG_MODULE(ECLTrackClusterMatchingPerformance);
25
26ECLTrackClusterMatchingPerformanceModule::ECLTrackClusterMatchingPerformanceModule() :
27 Module(), m_trackProperties()
28{
29 setDescription("Module to test the track cluster matching efficiency. Writes information about the tracks and MCParticles in a ROOT file.");
30 setPropertyFlags(c_ParallelProcessingCertified);
31
32 addParam("outputFileName", m_outputFileName, "Name of output root file.",
33 std::string("ECLTrackClusterMatchingPerformance.root"));
34 addParam("minClusterEnergy", m_minClusterEnergy, "Minimal cluster energy in units of particle's true energy for MC match",
35 0.5);
36 addParam("minWeight", m_minWeight, "Fraction of cluster energy required to originate from particle to be matched to cluster.",
37 0.5);
38 addParam("trackClusterRelationName", m_trackClusterRelationName, "Name of relation array between tracks and ECL clusters",
39 std::string(""));
40}
41
43{
44 // MCParticles and Tracks needed for this module
45 m_mcParticles.isRequired();
46 m_tracks.isRequired();
47 m_trackFitResults.isRequired();
48 m_eclClusters.isRequired();
49
50 m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE");
51 TDirectory* oldDir = gDirectory;
52 m_outputFile->cd();
53 m_tracksTree = new TTree("tracks", "tracks");
54 m_clusterTree = new TTree("cluster", "cluster");
55 oldDir->cd();
56
57 setupTree();
58}
59
61{
62 m_iEvent = m_EventMetaData->getEvent();
63 m_iRun = m_EventMetaData->getRun();
64 m_iExperiment = m_EventMetaData->getExperiment();
65
66 for (const ECLCluster& eclCluster : m_eclClusters) {
68 bool found_photon = false, found_mcmatch = false, found_charged_stable = false;
69 // find all MCParticles matched to the ECLCluster
70 const auto& relatedMCParticles = eclCluster.getRelationsWith<MCParticle>();
71 for (unsigned int index = 0; index < relatedMCParticles.size() && (!found_photon || !found_charged_stable); ++index) {
72 const MCParticle* relatedMCParticle = relatedMCParticles.object(index);
73 // check if matched MCParticle is primary photon
74 if (!(isPrimaryMcParticle(*relatedMCParticle))) continue;
75 found_mcmatch = true;
76 // get total energy of MCParticle deposited in this ECLCluster
77 const double weight = relatedMCParticles.weight(index);
78 // check that at least 50% of the generated energy of the particle is contained in this ECLCluster
79 // and check that the total cluster energy is greater than 50% of the energy coming from the particle
80
82 if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron)
83 and not eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
85 }
86
87 if (eclCluster.getEnergy(hypo) >= 0.5 * relatedMCParticle->getEnergy()) {
88 if (isChargedStable(*relatedMCParticle) && weight >= 0.5 * relatedMCParticle->getEnergy()) {
89 found_charged_stable = true;
90 } else if (relatedMCParticle->getPDG() == Const::photon.getPDGCode() && weight >= 0.5 * eclCluster.getEnergy(hypo)
91 && !found_photon) {
92 found_photon = true;
93 m_photonEnergy = relatedMCParticle->getEnergy();
94 }
95 }
96 }
97 if (found_mcmatch) {
98 if (eclCluster.isTrack()) {
100 }
101 m_clusterIsPhoton = int(found_photon);
102 m_clusterIsChargedStable = int(found_charged_stable);
103 m_clusterPhi = eclCluster.getPhi();
104 m_clusterTheta = eclCluster.getTheta();
105 if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
106 if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron)) {
108 } else {
110 }
112 } else if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron)) {
115 }
116
117 m_clusterErrorTiming = eclCluster.getDeltaTime99();
118 m_clusterE1E9 = eclCluster.getE1oE9();
119 m_clusterDetectorRegion = eclCluster.getDetectorRegion();
120 m_clusterMinTrkDistance = eclCluster.getMinTrkDistance();
121 m_clusterDeltaL = eclCluster.getDeltaL();
122 m_clusterTree->Fill();
123 }
124 }
125
126 for (const MCParticle& mcParticle : m_mcParticles) {
127 // check status of mcParticle
128 if (isPrimaryMcParticle(mcParticle) && isChargedStable(mcParticle) && mcParticle.hasStatus(MCParticle::c_StableInGenerator)) {
130
131 int pdgCode = mcParticle.getPDG();
132 B2DEBUG(29, "Primary MCParticle has PDG code " << pdgCode);
133 m_trackProperties.pdg_gen = pdgCode;
134
135 m_trackProperties.cosTheta_gen = cos(mcParticle.getMomentum().Theta());
136 m_trackProperties.phi_gen = mcParticle.getMomentum().Phi();
137 m_trackProperties.ptot_gen = mcParticle.getMomentum().R();
138 m_trackProperties.pt_gen = mcParticle.getMomentum().Rho();
139 m_trackProperties.px_gen = mcParticle.getMomentum().X();
140 m_trackProperties.py_gen = mcParticle.getMomentum().Y();
141 m_trackProperties.pz_gen = mcParticle.getMomentum().Z();
142 m_trackProperties.x_gen = mcParticle.getVertex().X();
143 m_trackProperties.y_gen = mcParticle.getVertex().Y();
144 m_trackProperties.z_gen = mcParticle.getVertex().Z();
145
146 const Track* b2Track = nullptr;
147 double maximumWeight = -2;
148 // find highest rated Track
149 const auto& relatedTracks = mcParticle.getRelationsWith<Track>();
150 for (unsigned int index = 0; index < relatedTracks.size(); ++index) {
151 const Track* relatedTrack = relatedTracks.object(index);
152 const double weight = relatedTracks.weight(index);
153
154 if (weight > maximumWeight) {
155 maximumWeight = weight;
156 b2Track = relatedTrack;
157 }
158 }
159
160 if (b2Track) {
161 const TrackFitResult* fitResult = b2Track->getTrackFitResultWithClosestMass(Const::ChargedStable(std::abs(pdgCode)));
162 B2ASSERT("Related Belle2 Track has no related track fit result!", fitResult);
163
164 // write some data to the root tree
165 ROOT::Math::XYZVector mom = fitResult->getMomentum();
166 m_trackProperties.cosTheta = cos(mom.Theta());
167 m_trackProperties.phi = mom.Phi();
168 m_trackProperties.ptot = mom.R();
169 m_trackProperties.pt = mom.Rho();
170 m_trackProperties.px = mom.X();
171 m_trackProperties.py = mom.Y();
172 m_trackProperties.pz = mom.Z();
173 m_trackProperties.x = fitResult->getPosition().X();
174 m_trackProperties.y = fitResult->getPosition().Y();
175 m_trackProperties.z = fitResult->getPosition().Z();
176
177 m_pValue = fitResult->getPValue();
178 m_charge = (int)fitResult->getChargeSign();
179 m_d0 = fitResult->getD0();
180 m_z0 = fitResult->getZ0();
181
182 // Count hits
183 m_trackProperties.nPXDhits = fitResult->getHitPatternVXD().getNPXDHits();
184 m_trackProperties.nSVDhits = fitResult->getHitPatternVXD().getNSVDHits();
185 m_trackProperties.nCDChits = fitResult->getHitPatternCDC().getNHits();
186
188
189 const ECLCluster* eclCluster_PhotonHypothesis_track_related = nullptr;
190 const ECLCluster* eclCluster_HadronHypothesis_track_related = nullptr;
191 for (const auto& eclCluster : b2Track->getRelationsTo<ECLCluster>("", m_trackClusterRelationName)) {
192 if (!(eclCluster.isTrack())) continue;
193 if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
195 m_matchedPhotonHypothesisClusterDetectorRegion = eclCluster.getDetectorRegion();
196 m_matchedPhotonHypothesisClusterTheta = eclCluster.getTheta();
197 m_matchedPhotonHypothesisClusterPhi = eclCluster.getPhi();
198 m_matchedPhotonHypothesisClusterMinTrkDistance = eclCluster.getMinTrkDistance();
199 m_matchedPhotonHypothesisClusterDeltaL = eclCluster.getDeltaL();
200 eclCluster_PhotonHypothesis_track_related = &eclCluster;
201 }
202 if (eclCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron)) {
204 m_matchedHadronHypothesisClusterDetectorRegion = eclCluster.getDetectorRegion();
205 m_matchedHadronHypothesisClusterTheta = eclCluster.getTheta();
206 m_matchedHadronHypothesisClusterPhi = eclCluster.getPhi();
207 m_matchedHadronHypothesisClusterMinTrkDistance = eclCluster.getMinTrkDistance();
208 m_matchedHadronHypothesisClusterDeltaL = eclCluster.getDeltaL();
209 eclCluster_HadronHypothesis_track_related = &eclCluster;
210 }
211 }
212 // find ECLCluster in which the particle has deposited the majority of its energy
213 maximumWeight = -2.;
214 const ECLCluster* eclCluster_matchedBestToMCParticle = nullptr;
215 const auto& relatedECLClusters = mcParticle.getRelationsFrom<ECLCluster>();
216 for (unsigned int index = 0; index < relatedECLClusters.size(); ++index) {
217 const ECLCluster* relatedECLCluster = relatedECLClusters.object(index);
218 if (!relatedECLCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) continue;
219 const double weight = relatedECLClusters.weight(index);
220 if (weight > maximumWeight) {
221 eclCluster_matchedBestToMCParticle = relatedECLCluster;
222 maximumWeight = weight;
223 }
224 }
225 if (eclCluster_matchedBestToMCParticle != nullptr) {
227 m_mcparticle_cluster_energy = maximumWeight;
228 m_mcparticle_cluster_theta = eclCluster_matchedBestToMCParticle->getTheta();
229 m_mcparticle_cluster_phi = eclCluster_matchedBestToMCParticle->getPhi();
230 m_mcparticle_cluster_detectorregion = eclCluster_matchedBestToMCParticle->getDetectorRegion();
231 if ((eclCluster_PhotonHypothesis_track_related == eclCluster_matchedBestToMCParticle
232 && eclCluster_matchedBestToMCParticle->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
233 || (eclCluster_HadronHypothesis_track_related == eclCluster_matchedBestToMCParticle
234 && eclCluster_matchedBestToMCParticle->hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron))) {
235 m_sameclusters = 1;
236 }
237 }
238 }
239 m_tracksTree->Fill(); // write data to tree
240 }
241 }
242}
243
244
249
254
256{
257 return Const::chargedStableSet.find(abs(mcParticle.getPDG()))
259}
260
262{
263 if (m_tracksTree == nullptr || m_clusterTree == nullptr) {
264 B2FATAL("Data tree or event tree was not created.");
265 }
269
271
272 addVariableToTree("cosTheta", m_trackProperties.cosTheta, m_tracksTree);
273 addVariableToTree("cosTheta_gen", m_trackProperties.cosTheta_gen, m_tracksTree);
274
277
280
283
286
289
292
295
298
300 addVariableToTree("ptot_gen", m_trackProperties.ptot_gen, m_tracksTree);
301
303
305
308
309 addVariableToTree("nPXDhits", m_trackProperties.nPXDhits, m_tracksTree);
310 addVariableToTree("nSVDhits", m_trackProperties.nSVDhits, m_tracksTree);
311 addVariableToTree("nCDChits", m_trackProperties.nCDChits, m_tracksTree);
312
314
318
320
326
332
337
338
342
353
357}
358
360{
361 if (m_tracksTree != nullptr && m_clusterTree != nullptr) {
362 TDirectory* oldDir = gDirectory;
363 if (m_outputFile)
364 m_outputFile->cd();
365 m_tracksTree->Write();
366 delete m_tracksTree;
367 m_clusterTree->Write();
368 delete m_clusterTree;
369 oldDir->cd();
370 }
371 if (m_outputFile != nullptr) {
372 m_outputFile->Close();
373 delete m_outputFile;
374 }
375}
376
425
427{
428 m_clusterPhi = -999;
429
430 m_clusterTheta = -999;
431
433
434 m_clusterIsPhoton = -999;
435
437
439
440 m_clusterEnergy = -999;
441
442 m_photonEnergy = -999;
443
444 m_clusterE1E9 = -999;
445
447
449
451
452 m_clusterDeltaL = -999;
453}
454
455void ECLTrackClusterMatchingPerformanceModule::addVariableToTree(const std::string& varName, double& varReference, TTree* tree)
456{
457 std::stringstream leaf;
458 leaf << varName << "/D";
459 tree->Branch(varName.c_str(), &varReference, leaf.str().c_str());
460}
461
462void ECLTrackClusterMatchingPerformanceModule::addVariableToTree(const std::string& varName, int& varReference, TTree* tree)
463{
464 std::stringstream leaf;
465 leaf << varName << "/I";
466 tree->Branch(varName.c_str(), &varReference, leaf.str().c_str());
467}
Provides a type-safe way to pass members of the chargedStableSet set.
Definition Const.h:589
static const ParticleSet chargedStableSet
set of charged stable particles
Definition Const.h:618
static const ParticleType invalidParticle
Invalid particle, used internally.
Definition Const.h:681
static const ParticleType photon
photon particle
Definition Const.h:673
ECL cluster data.
Definition ECLCluster.h:27
bool hasHypothesis(EHypothesisBit bitmask) const
Return if specific hypothesis bit is set.
Definition ECLCluster.h:351
double getPhi() const
Return Corrected Phi of Shower (radian).
Definition ECLCluster.h:304
int getDetectorRegion() const
Return detector region: 0: below acceptance, 1: FWD, 2: BRL, 3: BWD, 11: FWDGAP, 13: BWDGAP.
Definition ECLCluster.cc:70
double getTheta() const
Return Corrected Theta of Shower (radian).
Definition ECLCluster.h:307
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition ECLCluster.h:31
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
@ c_neutralHadron
CR is reconstructed as a neutral hadron (N2)
Definition ECLCluster.h:43
int m_sameclusters
boolean whether matched to ECL cluster with highest weight
double m_mcparticle_cluster_theta
theta of cluster matched to MCParticle
double m_matchedPhotonHypothesisClusterDeltaL
delta l of cluster with photon hypothesis
double m_d0
signed distance of the track to the IP in the r-phi plane
std::string m_trackClusterRelationName
name of relation array between tracks and ECL clusters
int m_mcparticle_cluster_match
boolean for match between MCParticle and ECL cluster
double m_matchedPhotonHypothesisClusterMinTrkDistance
minimal distance between cluster with photon hypothesis and track (not necessarily the matched one)
TTree * m_tracksTree
MCParticle based root tree with all output data.
void initialize() override
Register the needed StoreArrays and open th output TFile.
int m_matchedHadronHypothesisClusterDetectorRegion
detector region of cluster with hadron hypothesis matched to track
void addVariableToTree(const std::string &varName, double &varReference, TTree *tree)
add a variable with double format
void setClusterVariablesToDefaultValue()
sets cluster related variables to default values
double m_matchedHadronHypothesisClusterPhi
phi of cluster with hadron hypothesis matched to track
double m_matchedHadronHypothesisClusterMinTrkDistance
minimal distance between cluster with hadron hypothesis and track (not necessarily the matched one)
StoreArray< TrackFitResult > m_trackFitResults
Required input array of TrackFitResults.
void terminate() override
Write the tree into the opened root file.
double m_clusterErrorTiming
cluster's timing uncertainty containing 99% of true photons
ParticleProperties m_trackProperties
properties of a reconstructed track
bool isChargedStable(const MCParticle &mcParticle)
Tests if MCParticle is a charged stable particle.
int m_mcparticle_cluster_detectorregion
detector region of cluster matched to MCParticle
void setVariablesToDefaultValue()
Sets all variables to the default value, here -999.
int m_clusterIsPhoton
cluster fulfills requirements for being product of a photon
int m_matchedToPhotonHypothesisECLCluster
boolean for match between track and ECL cluster with photon hypothesis
double m_clusterE1E9
energy sum of central crystal over 3x3 array around central crystal
int m_matchedToHadronHypothesisECLCluster
boolean for match between track and ECL cluster with hadron hypothesis
StoreArray< Track > m_tracks
Required input array of Tracks.
bool isPrimaryMcParticle(const MCParticle &mcParticle)
Tests if MCParticle is a primary one.
StoreArray< ECLCluster > m_eclClusters
Required input array of ECLClusters.
double m_mcparticle_cluster_energy
amount of particle energy contained in cluster matched to MCParticle
double m_matchedHadronHypothesisClusterDeltaL
delta l of cluster with hadron hypothesis
TTree * m_clusterTree
root tree containing information on all truth-matched photon clusters.
double m_matchedPhotonHypothesisClusterPhi
phi of cluster with photon hypothesis matched to track
StoreArray< MCParticle > m_mcParticles
Required input array of MCParticles.
double m_matchedPhotonHypothesisClusterTheta
theta of cluster with photon hypothesis matched to track
double m_z0
distance of the track to the IP along the beam axis
void writeData()
write root tree to output file and close the file
double m_matchedHadronHypothesisClusterTheta
theta of cluster with hadron hypothesis matched to track
int m_matchedPhotonHypothesisClusterDetectorRegion
detector region of cluster with photon hypothesis matched to track
int m_clusterIsChargedStable
cluster has related MCParticle which is charged and stable
short getLastLayer() const
Returns the index of the last layer with a hit.
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
unsigned short getNPXDHits() const
Get total number of hits in the PXD.
unsigned short getNSVDHits() const
Get total number of hits in the SVD.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
float getEnergy() const
Return particle energy in GeV.
Definition MCParticle.h:136
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition MCParticle.h:49
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
Definition MCParticle.h:118
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
Base class for Modules.
Definition Module.h:72
RelationVector< FROM > getRelationsFrom(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from another store array to this object.
RelationVector< T > getRelationsWith(const std::string &name="", const std::string &namedRelation="") const
Get the relations between this object and another store array.
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
Values of the result of a track fit with a given particle hypothesis.
short getChargeSign() const
Return track charge (1 or -1).
double getPValue() const
Getter for Chi2 Probability of the track fit.
double getD0() const
Getter for d0.
double getZ0() const
Getter for z0.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
ROOT::Math::XYZVector getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
HitPatternVXD getHitPatternVXD() const
Getter for the hit pattern in the VXD;.
Class that bundles various TrackFitResults.
Definition Track.h:25
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for the fit hypothesis with the closest mass.
Definition Track.cc:104
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.