Belle II Software development
KlongValidationModule.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 <reconstruction/modules/KlId/KlongValidation/KlongValidationModule.h>
10#include <mdst/dataobjects/KlId.h>
11
12#include <framework/datastore/StoreArray.h>
13#include <framework/logging/Logger.h>
14#include <framework/gearbox/Const.h>
15
16#include <mdst/dataobjects/MCParticle.h>
17#include <mdst/dataobjects/KLMCluster.h>
18
19// here's where the functions are hidden
20#include <reconstruction/modules/KlId/KLMExpert/KlId.h>
21#include <tracking/dataobjects/TrackClusterSeparation.h>
22
23#include <TFile.h>
24#include <TGraph.h>
25#include <cstring>
26
27using namespace std;
28using namespace Belle2;
29using namespace Belle2::KlongId;
30
31REG_MODULE(KlongValidation);
32
34{
35 setDescription("Module used by the validation server to generate root files for the K0L validation. It calculates fake rates and efficiencies.");
36 addParam("outputName", m_outputName,
37 "Name of the output root file (must end with .root)",
39 addParam("KlIdCut", m_KlIdCut,
40 "If cut < 0, then only the neutral KLMClusters will be used. Otherwise, only the KLMClusters that satisfies the selection will be used.",
41 m_KlIdCut);
42 addParam("contact", m_contact, "Email address of contact for validation plots.", m_contact);
43}
44
48
50{
51 // require existence of necessary datastore obj
52 m_klmClusters.isRequired();
53 m_mcParticles.isRequired();
54
55 // initialize root tree to write stuff into
56 m_f = new TFile(m_outputName.c_str(), "recreate");
58 m_Phi_Pass = new TH1F("Phi Passed", "Passed #Phi;#Phi [rad];Count", 32, -3.2, 3.2);
59 m_Phi_all = new TH1F("Phi All", "All #Phi;#Phi [rad];Count", 32, -3.2, 3.2);
60 m_effPhi = new TH1F("Phi Efficiency", "Efficiency #Phi;#Phi [rad];Efficiency", 32, -3.2, 3.2);
61
62 m_Theta_Pass = new TH1F("Theta Passed", "Passed #Theta;#Theta [rad];Count", 32, 0, 3.2);
63 m_Theta_all = new TH1F("Theta All", "All #Theta;#Theta [rad];Count", 32, 0, 3.2);
64 m_effTheta = new TH1F("Theta Efficiency", "Efficiency #Theta;#Theta [rad];Efficiency", 32, 0, 3.2);
65
66 m_Mom_Pass = new TH1F("Momentum Efficiency", "Efficiency Momentum;Momentum [GeV];Efficiency", 25, 0, 5);
67 m_Mom_all = new TH1F("Momentum Efficiency normalise", "Efficiency Momentum;Momentum [GeV];Count", 25, 0, 5);
68 m_Mom_all_plot = new TH1F("Momentum Efficiency all, obtained from clusters", "Efficiency Momentum;Momentum [GeV];Count", 25, 0,
69 5);
70 m_effMom = new TH1F("Momentum Efficiency obtained from cluster", "Efficiency Momentum;Momentum [GeV];Efficiency", 25, 0, 5);
71
72 m_fakePhi_Pass = new TH1F("Phi Fake Passed", "Fake Passed #Phi;#Phi [rad];Count", 32, -3.2, 3.2);
73 m_fakePhi = new TH1F("Phi Fake Rate", "Fake Rate #Phi;#Phi [rad];Fake Rate", 32, -3.2, 3.2);
74
75 m_fakeTheta_Pass = new TH1F("Theta Fake Passed", "Fake Passed #Theta;#Theta [rad];Count", 32, 0, 3.2);
76 m_fakeTheta = new TH1F("Theta Fake Rate", "Fake Rate #Theta;#Theta [rad];Fake Rate", 32, 0, 3.2);
77
78 m_fakeMom_Pass = new TH1F("Momentum Fake Passed", "Momentum Fake Passed;Momentum [GeV];Count", 25, 0, 5);
79 m_fakeMom = new TH1F("Momentum Fake Rate", "Momentum Fake Rate;Momentum [GeV];Fake Rate", 25, 0, 5);
80
81 m_time = new TH1F("KLM Cluster Time", "Cluster Timing;Cluster time [ns];Count", 100, -30, 70);
82 m_trackSep = new TH1F("KLM trackSeperation Distance", "KLM trackSeperation Distance;Distance [mm];Count", 100, 0, 4000);
83 m_nLayer = new TH1F("KLM N-Layer", "N-layer;N-layer;count", 20, 0, 20);
84
85 m_bkgPhi = new TH1F("Phi Beam BKG", "BeamBKG #Phi;#Phi [rad];Background count", 32, -3.2, 3.2);
86 m_bkgTheta = new TH1F("Theta Beam BKG", "BeamBKG #Theta;#Theta [rad];Background count", 32, 0, 3.2);
87 m_bkgMom = new TH1F("Momentum Beam BKG", "BeamBKG #Momentum;#Momentum [GeV];Background count", 25, 0, 5);
88
89 m_innermostLayer = new TH1F("Innermost layer", "Innermost layer;Innermost layer; Count", 20, 0, 20);
90 m_trackFlag = new TH1F("Track flag", "TrackFlag;Flag; Count", 2, 0, 1);
91 m_ECLFlag = new TH1F("ECL flag", "ECLFlag;Flag; Count", 2, 0, 1);
92
93 m_klidFake = new TH1F("Klid Fake", "klid Fake;klid; Count", m_xbins.size() - 1, m_xbins.data());
94 m_klidTrue = new TH1F("Klid True", "Klid True;klid; Count", m_xbins.size() - 1, m_xbins.data());
95 m_klidAll = new TH1F("Klid all", "klid all clusters;klid; Count", m_xbins.size() - 1, m_xbins.data());
96}
97
101
105
107{
108 // the performance of the classifier depends on the cut on the classifier
109 // output. This is generally arbitrary as it depends on what the user wants.
110 // Finding K_L or rejecting?
111 // Therefore an arbitrary cut is chosen so that something happens and we can judge
112 // if the behavior changes.
113 for (const KLMCluster& cluster : m_klmClusters) {
114
115 MCParticle* mcpart = NULL;
116
117 auto klidObj = cluster.getRelatedTo<KlId>();
118
119 if (m_KlIdCut >= 0) {
120 if (klidObj) {
121 if (klidObj->getKlId() > m_KlIdCut) {
122 m_reconstructedAsKl = true;
123 } else {
124 m_reconstructedAsKl = false;
125 }
126 }
127 } else {
128 if (!cluster.getAssociatedTrackFlag()) {
129 m_reconstructedAsKl = true;
130 }
131 }
132 mcpart = cluster.getRelatedTo<MCParticle>();
133 // find mc truth
134 // go thru all mcparts mothers up to highest particle and check if its a Klong
135 m_isKl = 0;
136 if (mcpart == nullptr) {
137 m_isKl = 0; // this is the case for beambkg
138 m_isBeamBKG = 1;
139 } else {
140 m_isBeamBKG = 0;
141 while (!(mcpart -> getMother() == nullptr)) {
142 if (mcpart -> getPDG() == Const::Klong.getPDGCode()) {
143 m_isKl = 1;
144 break;
145 }
146 mcpart = mcpart -> getMother();
147 }// while
148 }
149 // second param is cut on
150 m_isKl = isKLMClusterSignal(cluster, 0);
151
152 double time = cluster.getTime();
153 double nLayer = cluster.getLayers();
154 double trackSep = 1.0e20;
155 auto trackSeperations = cluster.getRelationsTo<TrackClusterSeparation>();
156 for (auto trackSeperation : trackSeperations) {
157 double dist = trackSeperation.getDistance();
158 if (dist < trackSep) {
159 trackSep = dist;
160 }
161 }
162 // only fill if cluster is correctly reconstructed Kl
163 if (m_reconstructedAsKl && m_isKl) {m_passed = true;}
164 else {m_passed = false;}
165
166 // used for fake rate
167 if (m_reconstructedAsKl && (!m_isKl)) {m_faked = true;}
168 else {m_faked = false;}
169
170 m_phi = cluster.getMomentum().Phi();
171 m_theta = cluster.getMomentum().Theta();
172 m_momentum = std::abs(cluster.getMomentumMag());
173
174 if (m_passed) {
175 m_Phi_Pass -> Fill(m_phi);
176 m_Theta_Pass -> Fill(m_theta);
177 m_Mom_Pass -> Fill(m_momentum);
178 }
179
180 if (m_faked) {
181 m_fakePhi_Pass -> Fill(m_phi);
182 m_fakeTheta_Pass -> Fill(m_theta);
183 m_fakeMom_Pass -> Fill(m_momentum);
184 }
185
186 if (m_isBeamBKG) {
187 m_bkgPhi-> Fill(m_phi);
188 m_bkgTheta-> Fill(m_theta);
189 m_bkgMom-> Fill(m_momentum);
190 }
191
192 m_klidAll->SetMinimum(0.);
193 m_Mom_all_plot->SetMinimum(0.);
194 m_effPhi->SetMinimum(0.);
195 m_effTheta->SetMinimum(0.);
196 m_effMom->SetMinimum(0.);
197 m_fakePhi->SetMinimum(0.);
198 m_fakeTheta->SetMinimum(0.);
199 m_fakeMom->SetMinimum(0.);
200 m_time->SetMinimum(0.);
201 m_trackSep->SetMinimum(0.);
202 m_nLayer->SetMinimum(0.);
203 m_bkgPhi->SetMinimum(0.);
204 m_bkgTheta->SetMinimum(0.);
205 m_bkgMom->SetMinimum(0.);
206 m_innermostLayer->SetMinimum(0.);
207 m_trackFlag->SetMinimum(0.);
208 m_ECLFlag->SetMinimum(0.);
209
210 //fill all to normalise later
211 m_Phi_all -> Fill(m_phi);
212 m_Theta_all -> Fill(m_theta);
213 m_Mom_all -> Fill(m_momentum);
214 m_Mom_all_plot -> Fill(m_momentum);
215
216 m_time ->Fill(time);
217 m_trackSep->Fill(trackSep);
218 m_nLayer->Fill(nLayer);
219 m_innermostLayer->Fill(cluster.getInnermostLayer());
220 m_trackFlag->Fill(cluster.getAssociatedTrackFlag());
221 m_ECLFlag->Fill(cluster.getAssociatedEclClusterFlag());
222
223 if (klidObj) {
224 if (!m_isKl) {
225 m_klidFake ->Fill(klidObj->getKlId());
226 }
227 if (m_isKl) {
228 m_klidTrue ->Fill(klidObj->getKlId());
229 }
230 m_klidAll ->Fill(klidObj->getKlId());
231 }
232 }// for klm clusters
233
234}
235
237{
238 // write TEff to root file ,
239 m_f -> cd();
240 // TH1F is not compatible with the validation server
241 m_effPhi->Divide(m_Phi_Pass, m_Phi_all);
243 m_effMom->Divide(m_Mom_Pass, m_Mom_all);
244
248
249 std::vector<double> efficiency(m_xbins.size());
250 std::vector<double> purity(m_xbins.size());
251 std::vector<double> backrej(m_xbins.size());
252 double NallKlong = m_klidTrue->GetEntries();
253 double NallFakes = m_klidFake->GetEntries();
254 for (unsigned int i = 0; i < m_xbins.size() - 1; ++i) {
255 double NtruePositive = m_klidTrue->Integral(i, m_xbins.size() - 1);
256 double NfalsePositive = m_klidFake->Integral(i, m_xbins.size() - 1);
257 if ((NtruePositive + NfalsePositive) <= 0) {
258 purity[i] = 0;
259 } else {
260 purity[i] = NtruePositive / (NtruePositive + NfalsePositive);
261 }
262 if (NallFakes) {
263 backrej[i] = 1 - NfalsePositive / NallFakes;
264 } else {
265 backrej[i] = 0;
266 }
267 if (NallKlong <= 0) {
268 B2WARNING("No Klongs found in ROC calculation.");
269 } else {
270 efficiency[i] = NtruePositive / NallKlong;
271 }
272 }
273
274 m_ROC = new TGraph(m_xbins.size(), efficiency.data(), purity.data());
275 m_backRej = new TGraph(m_xbins.size(), efficiency.data(), backrej.data());
276
277 m_ROC->SetMarkerStyle(3);
278 m_backRej->SetMarkerStyle(3);
279 m_ROC->SetLineColorAlpha(kBlue, 0.0);
280 m_backRej->SetLineColorAlpha(kBlue, 0.0);
281 m_ROC->GetXaxis()->SetTitle("Efficiency");
282 m_backRej->GetXaxis()->SetTitle("Efficiency");
283 m_ROC->GetYaxis()->SetTitle("Purity");
284 m_backRej->GetYaxis()->SetTitle("Background rejection");
285 m_ROC-> SetTitle("Klong Purity vs Efficiency");
286 m_backRej-> SetTitle("Klong background rejection vs Efficiency");
287 m_ROC-> GetListOfFunctions() -> Add(new TNamed("Description",
288 "Purity vs Efficiency each point represents a cut on the klong ID."));
289 m_ROC -> GetListOfFunctions() -> Add(new TNamed("Check", "Should be as high as possible"));
290 m_backRej -> GetListOfFunctions() -> Add(new TNamed("Check", "Should be as high as possible"));
291 m_ROC -> GetListOfFunctions() -> Add(new TNamed("Contact", m_contact));
292 m_backRej -> GetListOfFunctions() -> Add(new TNamed("Contact", m_contact));
293
294 // tuple: pointer to the plot, name of the plot, true for shifter plots
295 std::vector<std::tuple<TH1F*, std::string, std::string, bool>> histograms;
296 std::string defaultCheck{"Nightly result should not differ significantly from the reference"};
297 histograms.push_back(make_tuple(m_klidAll, "KlId distribution", defaultCheck, true));
298 histograms.push_back(make_tuple(m_Mom_all_plot, "All Momenta generated", defaultCheck, false));
299 histograms.push_back(make_tuple(m_effPhi, "KlId efficiency in Phi", defaultCheck, false));
300 histograms.push_back(make_tuple(m_effTheta, "KlId efficiency in Theta",
301 defaultCheck + "; efficiency must be 0 below 0.3 rad and above 2.7 rad", true));
302 histograms.push_back(make_tuple(m_effMom, "KlId efficiency in Momentum", defaultCheck, true));
303 histograms.push_back(make_tuple(m_fakePhi, "KlId fake rate in Phi", defaultCheck, false));
304 histograms.push_back(make_tuple(m_fakeTheta, "KlId fake rate in Theta",
305 defaultCheck + "; fake rate must be 0 below 0.3 rad and above 2.7 rad", true));
306 histograms.push_back(make_tuple(m_fakeMom, "KlId fake rate in Momentum", defaultCheck, true));
307 histograms.push_back(make_tuple(m_time, "KLMClusters time", defaultCheck + "; no secondary peaks", true));
308 histograms.push_back(make_tuple(m_trackSep, "Distance of KLMCluster to closest Track", defaultCheck, true));
309 histograms.push_back(make_tuple(m_nLayer, "KLMClusters nLayer", defaultCheck, false));
310 histograms.push_back(make_tuple(m_innermostLayer, "KLMClusters innermostLayer", defaultCheck, false));
311 histograms.push_back(make_tuple(m_bkgPhi, "Beam bkg in Phi", defaultCheck, false));
312 histograms.push_back(make_tuple(m_bkgTheta, "Beam bkg in Theta",
313 defaultCheck + "; background must be 0 below 0.3 rad and above 2.7 rad", false));
314 histograms.push_back(make_tuple(m_bkgMom, "Beam bkg in Momentum", defaultCheck, false));
315 histograms.push_back(make_tuple(m_trackFlag, "Track flag", defaultCheck, false));
316 histograms.push_back(make_tuple(m_ECLFlag, "ECLCluster flag", defaultCheck, false));
317
318 for (auto hist : histograms) {
319 std::get<0>(hist) -> SetTitle(std::get<1>(hist).c_str());
320 std::get<0>(hist) -> GetListOfFunctions() -> Add(new TNamed("Description", std::get<1>(hist)));
321 std::get<0>(hist) -> GetListOfFunctions() -> Add(new TNamed("Check", std::get<2>(hist).c_str()));
322 std::get<0>(hist) -> GetListOfFunctions() -> Add(new TNamed("Contact", m_contact));
323 if (std::get<3>(hist))
324 std::get<0>(hist) -> GetListOfFunctions() -> Add(new TNamed("MetaOptions", "shifter,pvalue-warn=0.99,pvalue-error=0.02"));
325 else
326 std::get<0>(hist) -> GetListOfFunctions() -> Add(new TNamed("MetaOptions", "pvalue-warn=0.99,pvalue-error=0.1"));
327 std::get<0>(hist) -> Write();
328 }
329
330 //this guy is a tgraph
331 m_ROC -> Write();
332 m_backRej -> Write();
333 m_f -> Close();
334}
static const ParticleType Klong
K^0_L particle.
Definition Const.h:678
KLM cluster data.
Definition KLMCluster.h:29
Klong identification (KlId) datastore object to store results from KlId calculations.
Definition KlId.h:23
TH1F * m_Mom_Pass
momentum efficiency
TH1F * m_effTheta
efficiency in angle to z
bool m_passed
did cluster pass selection of algorythm?
TH1F * m_trackSep
track separation distance plot
TH1F * m_Mom_all_plot
momentum efficiency
TH1F * m_Theta_Pass
efficiency in angle to z
TH1F * m_fakeMom_Pass
fake momentum plot
StoreArray< KLMCluster > m_klmClusters
storearrays
std::string m_contact
contact email address
virtual void initialize() override
initialize
TGraph * m_backRej
background rejection
TH1F * m_fakePhi_Pass
fake phi, angle in x-y
virtual void event() override
process event
const std::vector< double > m_xbins
bins used for the ROC plots
virtual void endRun() override
end run
virtual void terminate() override
terminate
TH1F * m_fakeTheta_Pass
fake theta, angle to z
TH1F * m_Phi_Pass
efficiency in x-y plane
TH1F * m_Phi_all
efficiency in x-y plane
virtual void beginRun() override
beginn run
std::string m_outputName
output file name
TH1F * m_time
cluster timing plot
float m_KlIdCut
of > 0 use Klid else use trackflag as reconstruction criterion
virtual ~KlongValidationModule()
Destructor.
StoreArray< MCParticle > m_mcParticles
storearrays
bool m_faked
cluster wrongly reconstructed as K0L?
bool m_reconstructedAsKl
cluster reconstructed as K0L?
TH1F * m_fakeTheta
fake theta, angle to z
TH1F * m_effMom
momentum efficiency
TH1F * m_effPhi
efficiency in x-y plane
TH1F * m_Mom_all
momentum efficiency
TH1F * m_fakeMom
fake momentum plot
TH1F * m_fakePhi
fake phi, angle in x-y
TH1F * m_Theta_all
efficiency in angle to z
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
Store one Track-KLMCluster separation as a ROOT object.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Helper functions for all klid modules to improve readability of the code.
Definition KlId.h:27
bool isKLMClusterSignal(const Belle2::KLMCluster &cluster, float mcWeigthCut=0.66)
checks if a cluster is signal under the mcWeightcondition (mcWeight = energy deposition)
Definition KlId.h:100
Abstract base class for different kinds of events.
STL namespace.