Belle II Software development
EventShapeCalculatorModule.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/utility/PCmsLabTransform.h>
10
11#include <analysis/modules/EventShapeCalculator/EventShapeCalculatorModule.h>
12
13#include <analysis/dataobjects/ParticleList.h>
14#include <analysis/dataobjects/Particle.h>
15#include <mdst/dataobjects/MCParticle.h>
16
17#include <framework/logging/Logger.h>
18
19#include <analysis/ContinuumSuppression/Thrust.h>
20#include <analysis/ContinuumSuppression/HarmonicMoments.h>
21#include <analysis/ContinuumSuppression/CleoCones.h>
22#include <analysis/ContinuumSuppression/FoxWolfram.h>
23#include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
24
25using namespace std;
26using namespace Belle2;
27
28//-----------------------------------------------------------------
29// Register the Module
30//-----------------------------------------------------------------
31REG_MODULE(EventShapeCalculator);
32
33//-----------------------------------------------------------------
34// Implementation
35//-----------------------------------------------------------------
36
38{
39 // Set module properties
40 setDescription("Module to compute event shape attributes starting from particlelists. The core algorithms are not implemented in this module, but in dedicated basf2 classes.");
42 // Parameter definitions
43 addParam("particleListNames", m_particleListNames, "List of the ParticleLists to be used for the calculation of the EventShapes.",
44 vector<string>());
45 addParam("enableThrust", m_enableThrust, "Enables the calculation of thrust-related quantities.", true);
46 addParam("enableCollisionAxis", m_enableCollisionAxis, "Enables the calculation of the quantities related to the collision axis.",
47 true);
48 addParam("enableFoxWolfram", m_enableFW, "Enables the calculation of the Fox-Wolfram moments.", true);
49 addParam("enableHarmonicMoments", m_enableHarmonicMoments, "Enables the calculation of the Harmonic moments.", true);
50 addParam("enableJets", m_enableJets, "Enables the calculation of jet-related quantities.", true);
51 addParam("enableSphericity", m_enableSphericity, "Enables the calculation of the sphericity-related quantities.", true);
52 addParam("enableCleoCones", m_enableCleoCones, "Enables the calculation of the CLEO cones.", true);
53 addParam("enableAllMoments", m_enableAllMoments, "Enables the calculation of FW and harmonic moments from 5 to 8", false);
54 addParam("checkForDuplicates", m_checkForDuplicates,
55 "Enables the check for duplicates in the input list. If a duplicate entry is found, the first one is kept.", false);
56}
57
58
60{
61 m_eventShapeContainer.registerInDataStore();
62 if (m_enableJets and not m_enableThrust) {
63 B2WARNING("The jet-related quantities can only be calculated if the thrust calculation is activated as well.");
64 m_enableThrust = true;
65 }
67 B2WARNING("The CLEO cones can only be calculated if either the thrust or the collision axis calculation are activated as well.");
69 B2WARNING("The harmonic moments can only be calculated if either the thrust or the collision axis calculation are activated as well.");
70}
71
72
74{
75
77 double sqrtS = T.getCMSEnergy();
78
80
81 const int nPart = parseParticleLists(m_particleListNames);
82 if (nPart == 0) return;
83
84 // --------------------
85 // Calculates the FW moments
86 // --------------------
87 if (m_enableFW) {
91 for (short i = 0; i < 9; i++) {
92 m_eventShapeContainer->setFWMoment(i, fw.getH(i));
93 }
94 } else {
96 for (short i = 0; i < 5; i++) {
97 m_eventShapeContainer->setFWMoment(i, fw.getH(i));
98 }
99 }
100 }
101
102 // --------------------
103 // Calculates the sphericity quantities
104 // --------------------
105 if (m_enableSphericity) {
108 if (Sph.getEigenvalue(0) < Sph.getEigenvalue(1) || Sph.getEigenvalue(0) < Sph.getEigenvalue(2)
109 || Sph.getEigenvalue(1) < Sph.getEigenvalue(2))
110 B2WARNING("Eigenvalues not ordered!!!!!!!!!!");
111
112 for (short i = 0; i < 3; i++) {
113 m_eventShapeContainer->setSphericityEigenvalue(i, Sph.getEigenvalue(i));
114 m_eventShapeContainer->setSphericityEigenvector(i, Sph.getEigenvector(i));
115 }
116 }
117
118
119 // --------------------
120 // Calculates thrust and thrust-related quantities
121 // --------------------
122 if (m_enableThrust) {
123 ROOT::Math::XYZVector thrust = Thrust::calculateThrust(m_p4List);
124 float thrustVal = thrust.R();
125 thrust = thrust.Unit();
126 m_eventShapeContainer->setThrustAxis(thrust);
127 m_eventShapeContainer->setThrust(thrustVal);
128
129 // --- If required, calculates the HarmonicMoments ---
131 HarmonicMoments MM(m_p4List, thrust);
132 if (m_enableAllMoments) {
134 for (short i = 0; i < 9; i++) {
135 auto moment = MM.getMoment(i, sqrtS);
136 m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
137 }
138 } else {
140 for (short i = 0; i < 5; i++) {
141 auto moment = MM.getMoment(i, sqrtS);
142 m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
143 }
144 }
145 }
146
147 // --- If required, calculates the cleo cones w/ respect to the thrust axis ---
148 if (m_enableCleoCones) {
149 // Cleo cone class constructor. Unfortunately this class is designed
150 // to use the ROE, so the constructor takes two std::vector of momenta ("all" and "ROE"),
151 // then a vector to be used as axis, and finally two flags that determine if the cleo cones
152 // are calculated using the ROE, all the particles or both. Here we use the m_p4List as dummy
153 // list of the ROE momenta, that is however not used at all since the calculate only the
154 // cones with all the particles. This whole class would need some heavy restructuring...
155 CleoCones cleoCones(m_p4List, m_p4List, thrust, true, false);
156 std::vector<float> cones;
157 cones = cleoCones.cleo_cone_with_all();
158 for (short i = 0; i < 10; i++) {
159 m_eventShapeContainer->setCleoConeThrust(i, cones[i]);
160 }
161 } // end of if m_enableCleoCones
162
163
164 // --- If required, calculates the jet 4-momentum using the thrust axis ---
165 if (m_enableJets) {
166 ROOT::Math::PxPyPzEVector p4FWD;
167 ROOT::Math::PxPyPzEVector p4BKW;
168 for (const auto& p4 : m_p4List) {
169 if (p4.Vect().Dot(thrust) > 0)
170 p4FWD += p4;
171 else
172 p4BKW += p4;
173 }
174 m_eventShapeContainer->setForwardHemisphere4Momentum(p4FWD);
175 m_eventShapeContainer->setBackwardHemisphere4Momentum(p4BKW);
176 } // end of if m_enableJets
177 }// end of if m_enableThrust
178
179
180
181 // --------------------
182 // Calculates the collision axis quantities
183 // --------------------
185
186 ROOT::Math::XYZVector collisionAxis(0., 0., 1.);
187
188 // --- If required, calculates the cleo cones w/ respect to the collision axis ---
189 if (m_enableCleoCones) {
190 CleoCones cleoCones(m_p4List, m_p4List, collisionAxis, true, false);
191 std::vector<float> cones;
192 cones = cleoCones.cleo_cone_with_all();
193 for (short i = 0; i < 10; i++) {
194 m_eventShapeContainer->setCleoConeCollision(i, cones[i]);
195 }
196 }
197
198 // --- If required, calculates the HarmonicMoments ---
200 HarmonicMoments MM(m_p4List, collisionAxis);
201 if (m_enableAllMoments) {
203 for (short i = 0; i < 9; i++) {
204 auto moment = MM.getMoment(i, sqrtS);
205 m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
206 }
207 } else {
209 for (short i = 0; i < 5; i++) {
210 auto moment = MM.getMoment(i, sqrtS);
211 m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
212 }
213 }
214 } // end of m_enableHarmonicMoments
215 } // end of m_enableCollisionAxis
216} // end of event()
217
218
219
220int EventShapeCalculatorModule::parseParticleLists(vector<string> particleListNames)
221{
223 m_p4List.clear();
224
225 unsigned int nParticlesInAllLists = 0;
226 unsigned short nParticleLists = particleListNames.size();
227 if (nParticleLists == 0)
228 B2WARNING("No particle lists found. EventShape calculation not performed.");
229
230 // This vector temporarily stores the mdstSource of particle objects
231 // that have been processed so far (not only the momenta)
232 // in order to check for duplicates before pushing the 3- and 4- vectors
233 // in the corresponding lists
234 std::vector<int> usedMdstSources;
235
236 // Loops over the number of particle lists
237 for (unsigned short iList = 0; iList < nParticleLists; iList++) {
238 string particleListName = particleListNames[iList];
239 StoreObjPtr<ParticleList> particleList(particleListName);
240
241 // Loops over the number of particles in the list
242 nParticlesInAllLists += particleList->getListSize();
243
244 for (unsigned int iPart = 0; iPart < particleList->getListSize(); iPart++) {
245 const Particle* part = particleList->getParticle(iPart);
246 const MCParticle* mcParticle = part->getMCParticle();
247 if (mcParticle and mcParticle->isInitial()) continue;
248
250
251 std::vector<const Particle*> finalStateDaughters = part->getFinalStateDaughters();
252
253 for (const auto fsp : finalStateDaughters) {
254 int mdstSource = fsp->getMdstSource();
255 auto result = std::find(usedMdstSources.begin(), usedMdstSources.end(), mdstSource);
256 if (result == usedMdstSources.end()) {
257 usedMdstSources.push_back(mdstSource);
258 ROOT::Math::PxPyPzEVector p4CMS = T.rotateLabToCms() * fsp->get4Vector();
259 m_p4List.push_back(p4CMS);
260 B2DEBUG(19, "non-duplicate has pdgCode " << fsp->getPDGCode() << " and mdstSource " << mdstSource);
261 } else {
262 B2DEBUG(19, "duplicate has pdgCode " << fsp->getPDGCode() << " and mdstSource " << mdstSource);
263 B2DEBUG(19, "Duplicate particle found. The new one won't be used for the calculation of the event shape variables. "
264 "Please, double check your input lists and try to make them mutually exclusive.");
265 }
266 }
267 } else {
268 ROOT::Math::PxPyPzEVector p4CMS = T.rotateLabToCms() * part->get4Vector();
269 m_p4List.push_back(p4CMS);
270 }
271 }
272 }
273
274 return nParticlesInAllLists;
275}
276
Class to calculate the Cleo clone variables.
Definition CleoCones.h:23
int parseParticleLists(std::vector< std::string >)
Turns the ParticleLists provided as inputs into std::vector of PxPyPzEVector, after boosting them int...
bool m_enableCollisionAxis
Enables the calculation of the quantities related to the collision axis.
virtual void initialize() override
Define the physical parameters.
virtual void event() override
Main method, called for each events.
bool m_enableHarmonicMoments
Enables the calculation of the Harmonic moments.
std::vector< std::string > m_particleListNames
Names of the ParticleLists (inputs).
bool m_enableAllMoments
Enables the calculation of the FW moments from 5 to 8.
bool m_enableThrust
Enables the calculation of thrust-related quantities.
bool m_enableCleoCones
Enables the calculation of the Cleo Cones.
bool m_enableJets
Enables the calculation of the Jet-related quantities.
std::vector< ROOT::Math::PxPyPzEVector > m_p4List
vector containing all the 4-momenta of the particles contained in the input lists.
bool m_enableFW
Enables the calculation of the FW moments.
bool m_checkForDuplicates
Enables the check for the duplicates in the input list.
EventShapeCalculatorModule()
Constructor: Sets the description, the properties and the parameters of the module.
StoreObjPtr< EventShapeContainer > m_eventShapeContainer
event shape container object pointer
bool m_enableSphericity
Enables the calculation of the Sphericity matrix.
Class to calculate the Fox-Wolfram moments up to order 8.
Definition FoxWolfram.h:28
void calculateAllMoments()
Method to perform the calculation of the moments up to order 8.
Definition FoxWolfram.cc:58
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition FoxWolfram.cc:14
double getH(int i) const
Returns the i-th moment.
Definition FoxWolfram.h:83
Class to calculate the Harmonic moments up to order 8 with respect to a given axis.
void calculateAllMoments()
Calculates the moments up to order 8.
double getMoment(short i, double sqrts) const
Returns the moment of order i.
void calculateBasicMoments()
Calculates the moments up to order 4.
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
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
Class to hold Lorentz transformations from/to CMS and boost vector.
Class to store reconstructed particles.
Definition Particle.h:76
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition Particle.cc:1027
std::vector< const Particle * > getFinalStateDaughters() const
Returns a vector of pointers to Final State daughter particles.
Definition Particle.cc:680
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition Particle.h:567
Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-mom...
void calculateEigenvalues()
Calculates eigenvalues and eigenvectors.
ROOT::Math::XYZVector getEigenvector(short i) const
Returns the i-th Eigenvector.
double getEigenvalue(short i) const
Returns the i-th Eigenvalue.
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::PxPyPzEVector > &momenta)
calculates the thrust axis
Definition Thrust.cc:71
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
bool isInitial() const
Check if particle is an initial particle such as ISR.
Definition MCParticle.h:580
Abstract base class for different kinds of events.
STL namespace.