Belle II Software light-2406-ragdoll
ParticleStatsModule.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 <analysis/modules/ParticleStats/ParticleStatsModule.h>
10#include <analysis/dataobjects/ParticleList.h>
11#include <framework/core/Environment.h>
12#include <framework/datastore/StoreObjPtr.h>
13#include <nlohmann/json.hpp>
14#include <fstream>
15
16
17using namespace Belle2;
18
19// Register module in the framework
20REG_MODULE(ParticleStats);
21
23{
24 // Set module properties and parameters
25 setDescription("Make a summary of specific ParticleLists.");
26 addParam("outputFile", m_outputFile, "Name of output file", std::string());
27 addParam("printPassMatrix", m_printPassMatrix, "Should we also calculate and print the pass matrix?", true);
28 addParam("particleLists", m_strParticleLists, "List of ParticleLists", std::vector<std::string>());
29}
30
32{
35
36 // check all particle lists are valid
37 for (unsigned i = 0; i < m_nLists; ++i) {
39 if (!valid)
40 B2ERROR("Invalid input list name: " << m_strParticleLists[i]);
41
43 B2ERROR("ParticleStatsModule::initialize Invalid input DecayString " << m_strParticleLists[i]
44 << ". DecayString should not contain any daughters, only the mother particle.");
45 }
46 B2INFO("Number of ParticleLists to be studied: " << m_nLists << " ");
47
48 m_PassMatrix = new TMatrix(m_nLists, m_nLists + 1);
49 m_MultiplicityMatrix = new TMatrix(m_nLists, 4); // 0 All particles; 1 Negative; 2 Positive; 3 SelfConjugated
50}
51
53{
54 bool unique = true;
55 bool pass = false;
56
57 // first loop over all particle lists we are interested in
58 for (unsigned iList = 0; iList < m_nLists; ++iList) {
59
60 // check that the list exists and is not empty
62 if (!particlelist) {
63 B2ERROR("ParticleListi " << m_strParticleLists[iList] << " not found");
64 continue;
65 }
66 if (!particlelist->getListSize()) continue;
67
68 pass = true;
69
70 /**************************************************************************
71 * For the multiplicity table, count the number of candidates in each
72 * list. Add to the total for all events processed so far.
73 */
74
75 // all particles and anti-particles, or self-conjugate
76 (*m_MultiplicityMatrix)(iList, 0) = (*m_MultiplicityMatrix)(iList, 0)
77 + particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle)
78 + particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle, true)
79 + particlelist->getNParticlesOfType(ParticleList::c_SelfConjugatedParticle);
80
81 // particles alone
82 if (particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle))
83 (*m_MultiplicityMatrix)(iList, 1) = (*m_MultiplicityMatrix)(iList, 1)
84 + particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle);
85
86 // anti-particles alone
87 if (particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle, true))
88 (*m_MultiplicityMatrix)(iList, 2) = (*m_MultiplicityMatrix)(iList, 2)
89 + particlelist->getNParticlesOfType(ParticleList::c_FlavorSpecificParticle, true);
90
91 // self-conjugate particles alone
92 if (particlelist->getNParticlesOfType(ParticleList::c_SelfConjugatedParticle))
93 (*m_MultiplicityMatrix)(iList, 3) = (*m_MultiplicityMatrix)(iList, 3)
94 + particlelist->getNParticlesOfType(ParticleList::c_SelfConjugatedParticle);
95
96 /**************************************************************************
97 * Now, loop again over all particle lists to fill the pass matrix for this
98 * event. Add to the totals for all events processed so far.
99 */
100 for (unsigned jList = 0; jList < m_nLists; ++jList) {
101
102 // check that the list exists and is not empty
103 StoreObjPtr<ParticleList> particlelistj(m_strParticleLists[jList]);
104 if (!particlelistj) {
105 B2INFO("ParticleListj " << m_strParticleLists[jList] << " not found");
106 continue;
107 }
108 if (!particlelistj->getListSize()) continue;
109
110 // if the event passed both the i and j list selection then increment
111 (*m_PassMatrix)(iList, jList) = (*m_PassMatrix)(iList, jList) + 1.;
112 if (iList != jList) unique = false;
113
114 } // end inner loop over all particle lists
115
116 if (unique)(*m_PassMatrix)(iList, m_nLists) = (*m_PassMatrix)(iList, m_nLists) + 1.;
117
118 } // end outer loop over all particle lists
119
120 // finally count the particles created in this event, and the number of events
121 // passing any of the list selections that we were interested in
123 if (pass) m_nPass++;
124}
125
126// finish up the calculation and print a table to B2INFO
128{
129 auto nEvents = (float)Environment::Instance().getNumberOfEvents();
130 B2INFO("ParticleStats Summary:");
131 std::ostringstream stream;
132
133 nlohmann::json json;
134
135 // divide every entry of the pass matrix by the number of events processed
136 // (entrywise) to get the fraction of events
137 for (unsigned iList = 0; iList < m_nLists; ++iList)
138 for (unsigned jList = 0; jList < m_nLists + 1; ++jList)
139 (*m_PassMatrix)(iList, jList) = (*m_PassMatrix)(iList, jList) / nEvents;
140
141 for (unsigned iList = 0; iList < m_nLists; ++iList)
142 for (unsigned jList = 0; jList < m_nLists + 1; ++jList)
143 if (iList != jList && (*m_PassMatrix)(iList, iList) > 0.)
144 (*m_PassMatrix)(iList, jList) = (*m_PassMatrix)(iList, jList) / (*m_PassMatrix)(iList, iList) ;
145
146 if (m_printPassMatrix && m_nLists > 1) {
147
148 // format the pass matrix table
149 stream << "=======================================================================\n";
150 stream << " Pass matrix (for i, j: fraction passing selection i also passes selection j)\n\n";
151 stream << " ";
152
153 // table headings
154 for (unsigned iList = 0; iList < m_nLists; ++iList)
155 stream << "|" << Form(" %5d", iList);
156 stream << "| Unique|\n";
157
158 // table data
159 for (unsigned iList = 0; iList < m_nLists; ++iList) {
160 // first column is particle list name
161 stream << Form("%14s(%2d)", m_strParticleLists[iList].c_str(), iList) << "|";
162 // next N columns are particle list survival fractions (need to use tabs here)
163 for (unsigned jList = 0; jList < m_nLists + 1; ++jList) { // matrix
164 if (iList != jList) {
165 stream << Form(" %6.4f|", (*m_PassMatrix)(iList, jList));
166 std::string jName = (jList < m_nLists ? m_strParticleLists[jList].c_str() : "Unique");
167 json["Pass matrix"][m_strParticleLists[iList].c_str()][jName.c_str()] = (*m_PassMatrix)(iList, jList);
168 }
169 if (iList == jList) {
170 stream << Form(" %6.4f|", 1.0);
171 json["Pass matrix"][m_strParticleLists[iList].c_str()][m_strParticleLists[jList].c_str()] = 1.0;
172 }
173
174 }
175 stream << "\n";
176 }
177 B2INFO(stream.str()); // print the pass matrix table
178
179
180 } // end if
181
182 // clear the stream and format the multiplicity table
183 stream.str("");
184 stream << "\n======================================================================\n";
185 stream << " Average Candidate Multiplicity (ACM) and ACM for Passed Events (ACMPE) \n\n";
186
187 // table headings
188 stream << " | All Particles | Particles | Anti-particles | Self-conjugates |\n";
189 stream << " | Retention| ACM| ACMPE| ACM| ACMPE| ACM| ACMPE| ACM| ACMPE|\n";
190
191 // table data
192 for (unsigned iList = 0; iList < m_nLists; ++iList) {
193 // first column is particle list name
194 stream << Form("%14s(%2d)", m_strParticleLists[iList].c_str(), iList) << "|";
195 // second column is retention (no need to use tabs here because it can't be
196 // more than 1.{followed by 4 digits})
197 stream << Form(" %6.4f|", (*m_PassMatrix)(iList, iList));
198
199 std::string pName = m_strParticleLists[iList].c_str();
200 float retRate = (*m_PassMatrix)(iList, iList);
201 json["Retention"][pName]["Retention"] = retRate;
202 std::string flavs[4] = {"All Particles", "Particles", "Anti Particles", "Self-conjugates"};
203 // now the ACM and ACMPE
204 for (int iFlav = 0; iFlav < 4; ++iFlav) {
205 stream << Form(" %8.4f|", (*m_MultiplicityMatrix)(iList, iFlav) / nEvents);
206 json["Retention"][pName][Form("%s ACM", flavs[iFlav].c_str())] = (*m_MultiplicityMatrix)(iList, iFlav) / nEvents;
207
208 stream << Form(" %8.4f|", (*m_MultiplicityMatrix)(iList, iFlav) / nEvents / (*m_PassMatrix)(iList, iList));
209 json["Retention"][pName][Form("%s ACPME", flavs[iFlav].c_str())] = (*m_MultiplicityMatrix)(iList,
210 iFlav) / nEvents / (*m_PassMatrix)(iList,
211 iList);
212 }
213 stream << "\n";
214 }
215 B2INFO(stream.str()); // print the multiplicity table
216
217 // now print some global information
218 stream.str("");
219 stream << "\n======================================================================\n";
220 stream << "Total Retention: " << m_nPass << " events passing / " << nEvents << " events processed = " << Form("%6.4f\n",
221 (float)m_nPass / (float)nEvents);
222 stream << "Total Number of Particles created in the DataStore: " << m_nParticles;
223 stream << "\n======================================================================\n";
224 json["Total retention"] = m_nPass / nEvents;
225 json["Events passing"] = m_nPass;
226 json["Events processed"] = nEvents;
227 json["Total particles number"] = m_nParticles;
228 B2INFO(stream.str());
229 if (m_outputFile != "") {
230 std::ofstream jsonFile(m_outputFile);
231 jsonFile << json.dump(2) << std::endl;
232 }
233
234 delete m_PassMatrix;
236}
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
int getNDaughters() const
return number of direct daughters.
unsigned int getNumberOfEvents() const
Return the number of events, from either input or EventInfoSetter, or -n command line override (if le...
Definition: Environment.cc:39
static Environment & Instance()
Static method to get a reference to the Environment instance.
Definition: Environment.cc:28
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
TMatrix * m_PassMatrix
Pass matrix for the particle lists.
TMatrix * m_MultiplicityMatrix
Particle multiplicities for the particle lists.
bool m_printPassMatrix
Shall we also calculate and print out the pass matrix?
virtual void initialize() override
Initialises the module.
virtual void event() override
Fill the particle stats matrices for each event.
std::string m_outputFile
Name of output file.
StoreArray< Particle > m_particles
StoreArray of Particles.
virtual void terminate() override
Print out the arrays as a table in B2INFO.
std::vector< std::string > m_strParticleLists
The particle lists to produce statistics about.
int m_nParticles
Count the total number of particles.
unsigned m_nLists
The number of particle lists.
int m_nPass
Number of events with Particle candidates.
DecayDescriptor m_decaydescriptor
Decay descriptor of the particle being selected.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24