Belle II Software light-2406-ragdoll
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 <iostream>
20
21#include <analysis/ContinuumSuppression/Thrust.h>
22#include <analysis/ContinuumSuppression/HarmonicMoments.h>
23#include <analysis/ContinuumSuppression/CleoCones.h>
24#include <analysis/ContinuumSuppression/FoxWolfram.h>
25#include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
26
27#include <Math/Vector4D.h>
28
29
30using namespace std;
31using namespace Belle2;
32
33//-----------------------------------------------------------------
34// Register the Module
35//-----------------------------------------------------------------
36REG_MODULE(EventShapeCalculator);
37
38//-----------------------------------------------------------------
39// Implementation
40//-----------------------------------------------------------------
41
43{
44 // Set module properties
45 setDescription("Module to compute event shape attributes starting from particlelists. The core algorithms are not implemented in this module, but in dedicated basf2 classes.");
47 // Parameter definitions
48 addParam("particleListNames", m_particleListNames, "List of the ParticleLists to be used for the calculation of the EventShapes.",
49 vector<string>());
50 addParam("enableThrust", m_enableThrust, "Enables the calculation of thust-related quantities.", true);
51 addParam("enableCollisionAxis", m_enableCollisionAxis, "Enables the calculation of the quantities related to the collision axis.",
52 true);
53 addParam("enableFoxWolfram", m_enableFW, "Enables the calculation of the Fox-Wolfram moments.", true);
54 addParam("enableHarmonicMoments", m_enableHarmonicMoments, "Enables the calculation of the Harmonic moments.", true);
55 addParam("enableJets", m_enableJets, "Enables the calculation of jet-related quantities.", true);
56 addParam("enableSphericity", m_enableSphericity, "Enables the calculation of the sphericity-related quantities.", true);
57 addParam("enableCleoCones", m_enableCleoCones, "Enables the calculation of the CLEO cones.", true);
58 addParam("enableAllMoments", m_enableAllMoments, "Enables the calculation of FW and harmonic moments from 5 to 8", false);
59 addParam("checkForDuplicates", m_checkForDuplicates,
60 "Enables the check for duplicates in the input list. If a duplicate entry is found, the first one is kept.", false);
61}
62
63
65{
66 m_eventShapeContainer.registerInDataStore();
67}
68
69
71{
72
74 double sqrtS = T.getCMSEnergy();
75
77
78 const int nPart = parseParticleLists(m_particleListNames);
79 if (nPart == 0) return;
80
81 // --------------------
82 // Calculates the FW moments
83 // --------------------
84 if (m_enableFW) {
88 for (short i = 0; i < 9; i++) {
89 m_eventShapeContainer->setFWMoment(i, fw.getH(i));
90 }
91 } else {
93 for (short i = 0; i < 5; i++) {
94 m_eventShapeContainer->setFWMoment(i, fw.getH(i));
95 }
96 }
97 }
98
99 // --------------------
100 // Calculates the sphericity quantities
101 // --------------------
102 if (m_enableSphericity) {
105 if (Sph.getEigenvalue(0) < Sph.getEigenvalue(1) || Sph.getEigenvalue(0) < Sph.getEigenvalue(2)
106 || Sph.getEigenvalue(1) < Sph.getEigenvalue(2))
107 B2WARNING("Eigenvalues not ordered!!!!!!!!!!");
108
109 for (short i = 0; i < 3; i++) {
110 m_eventShapeContainer->setSphericityEigenvalue(i, Sph.getEigenvalue(i));
111 m_eventShapeContainer->setSphericityEigenvector(i, Sph.getEigenvector(i));
112 }
113 }
114
115
116 // --------------------
117 // Calculates thrust and thrust-related quantities
118 // --------------------
119 if (m_enableThrust) {
120 ROOT::Math::XYZVector thrust = Thrust::calculateThrust(m_p3List);
121 float thrustVal = thrust.R();
122 thrust = thrust.Unit();
123 m_eventShapeContainer->setThrustAxis(thrust);
124 m_eventShapeContainer->setThrust(thrustVal);
125
126 // --- If required, calculates the HarmonicMoments ---
128 HarmonicMoments MM(m_p3List, thrust);
129 if (m_enableAllMoments) {
131 for (short i = 0; i < 9; i++) {
132 auto moment = MM.getMoment(i, sqrtS);
133 m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
134 }
135 } else {
137 for (short i = 0; i < 5; i++) {
138 auto moment = MM.getMoment(i, sqrtS);
139 m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
140 }
141 }
142 }
143
144 // --- If required, calculates the cleo cones w/ respect to the thrust axis ---
145 if (m_enableCleoCones) {
146 // Cleo cone class constructor. Unfortunately this class is designed
147 // to use the ROE, so the constructor takes two std::vector of momenta ("all" and "ROE"),
148 // then a vector to be used as axis, and finally two flags that determine if the cleo cones
149 // are calculated using the ROE, all the particles or both. Here we use the m_p3List as dummy
150 // list of the ROE momenta, that is however not used at all since the calculate only the
151 // cones with all the particles. This whole class would need some heavy restructuring...
152 CleoCones cleoCones(m_p3List, m_p3List, thrust, true, false);
153 std::vector<float> cones;
154 cones = cleoCones.cleo_cone_with_all();
155 for (short i = 0; i < 10; i++) {
156 m_eventShapeContainer->setCleoConeThrust(i, cones[i]);
157 }
158 } // end of if m_enableCleoCones
159
160
161 // --- If required, calculates the jet 4-momentum using the thrust axis ---
162 if (m_enableJets) {
163 ROOT::Math::PxPyPzEVector p4FWD;
164 ROOT::Math::PxPyPzEVector p4BKW;
165 for (const auto& p4 : m_p4List) {
166 if (p4.Vect().Dot(thrust) > 0)
167 p4FWD += p4;
168 else
169 p4BKW += p4;
170 }
171 m_eventShapeContainer->setForwardHemisphere4Momentum(p4FWD);
172 m_eventShapeContainer->setBackwardHemisphere4Momentum(p4BKW);
173 } // end of if m_enableJets
174 }// end of if m_enableThrust
175
176
177
178 // --------------------
179 // Calculates the collision axis quantities
180 // --------------------
182
183 ROOT::Math::XYZVector collisionAxis(0., 0., 1.);
184
185 // --- If required, calculates the cleo cones w/ respect to the collision axis ---
186 if (m_enableCleoCones) {
187 CleoCones cleoCones(m_p3List, m_p3List, collisionAxis, true, false);
188 std::vector<float> cones;
189 cones = cleoCones.cleo_cone_with_all();
190 for (short i = 0; i < 10; i++) {
191 m_eventShapeContainer->setCleoConeCollision(i, cones[i]);
192 }
193 }
194
195 // --- If required, calculates the HarmonicMoments ---
197 HarmonicMoments MM(m_p3List, collisionAxis);
198 if (m_enableAllMoments) {
200 for (short i = 0; i < 9; i++) {
201 auto moment = MM.getMoment(i, sqrtS);
202 m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
203 }
204 } else {
206 for (short i = 0; i < 5; i++) {
207 auto moment = MM.getMoment(i, sqrtS);
208 m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
209 }
210 }
211 } // end of m_enableHarmonicMoments
212 } // end of m_enableCollisionAxis
213} // end of event()
214
215
216
217int EventShapeCalculatorModule::parseParticleLists(vector<string> particleListNames)
218{
220 m_p4List.clear();
221 m_p3List.clear();
222
223 unsigned int nParticlesInAllLists = 0;
224 unsigned short nParticleLists = particleListNames.size();
225 if (nParticleLists == 0)
226 B2WARNING("No particle lists found. EventShape calculation not performed.");
227
228 // This vector temporary stores the mdstSource of particle objects
229 // that have been processed so far (not only the momenta)
230 // in order to check for duplicates before pushing the 3- and 4- vectors
231 // in the corresponding lists
232 std::vector<int> usedMdstSources;
233
234 // Loops over the number of particle lists
235 for (unsigned short iList = 0; iList < nParticleLists; iList++) {
236 string particleListName = particleListNames[iList];
237 StoreObjPtr<ParticleList> particleList(particleListName);
238
239 // Loops over the number of particles in the list
240 nParticlesInAllLists += particleList->getListSize();
241
242 for (unsigned int iPart = 0; iPart < particleList->getListSize(); iPart++) {
243 const Particle* part = particleList->getParticle(iPart);
244 const MCParticle* mcParticle = part->getMCParticle();
245 if (mcParticle and mcParticle->isInitial()) continue;
246
247 // Flag to check for duplicates across the lists.
248 // It can be true only if m_checkForDuplicates is enabled
249 bool isDuplicate = false;
250
252 int mdstSource = part->getMdstSource();
253
254 auto result = std::find(usedMdstSources.begin(), usedMdstSources.end(), mdstSource);
255 if (result == usedMdstSources.end()) {
256 usedMdstSources.push_back(mdstSource);
257 } else {
258 B2WARNING("Duplicate particle found. The new one won't be used for the calculation of the event shape variables. "
259 "Please, double check your input lists and try to make them mutually exclusive.");
260 isDuplicate = true;
261 }
262 }
263
264 if (!isDuplicate) {
265 ROOT::Math::PxPyPzEVector p4CMS = T.rotateLabToCms() * part->get4Vector();
266 // it need to fill an std::vector of XYZVector to use the current FW routines.
267 // It will hopefully change in release 3
268 m_p4List.push_back(p4CMS);
269 m_p3List.push_back(p4CMS.Vect());
270 }
271 }
272 }
273
274 return nParticlesInAllLists;
275}
276
Class to calculate the Cleo clone variables.
Definition: CleoCones.h:22
int parseParticleLists(std::vector< std::string >)
Turns the ParticleLists provided as inputs in std::vectors of PxPyPzEVector and B2Vector3D,...
bool m_enableCollisionAxis
Enables the calculation of the quantities related to the collision axis.
virtual void initialize() override
Define the physical parameters.
std::vector< ROOT::Math::XYZVector > m_p3List
vector containing all the 3-momenta of the particles contained in the input lists.
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 thust-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
Base class for Modules.
Definition: Module.h:72
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
@ 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.
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
const ROOT::Math::LorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Class to store reconstructed particles.
Definition: Particle.h:75
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:942
int getMdstSource() const
Returns unique identifier of final state particle (needed in particle combiner)
Definition: Particle.cc:369
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
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::XYZVector > &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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
bool isInitial() const
Check if particle is an initial particle such as ISR.
Definition: MCParticle.h:590
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.