Belle II Software  light-2212-foldex
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 
16 #include <framework/logging/Logger.h>
17 
18 #include <iostream>
19 
20 #include <analysis/ContinuumSuppression/Thrust.h>
21 #include <analysis/ContinuumSuppression/HarmonicMoments.h>
22 #include <analysis/ContinuumSuppression/CleoCones.h>
23 #include <analysis/ContinuumSuppression/FoxWolfram.h>
24 #include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
25 
26 #include <Math/Vector4D.h>
27 
28 
29 using namespace std;
30 using namespace Belle2;
31 
32 //-----------------------------------------------------------------
33 // Register the Module
34 //-----------------------------------------------------------------
35 REG_MODULE(EventShapeCalculator);
36 
37 //-----------------------------------------------------------------
38 // Implementation
39 //-----------------------------------------------------------------
40 
41 EventShapeCalculatorModule::EventShapeCalculatorModule() : Module()
42 {
43  // Set module properties
44  setDescription("Module to compute event shape attributes starting from particlelists. The core algorithms are not implemented in this module, but in dedicated basf2 classes.");
46  // Parameter definitions
47  addParam("particleListNames", m_particleListNames, "List of the ParticleLists to be used for the calculation of the EventShapes.",
48  vector<string>());
49  addParam("enableThrust", m_enableThrust, "Enables the calculation of thust-related quantities.", true);
50  addParam("enableCollisionAxis", m_enableCollisionAxis, "Enables the calculation of the quantities related to the collision axis.",
51  true);
52  addParam("enableFoxWolfram", m_enableFW, "Enables the calculation of the Fox-Wolfram moments.", true);
53  addParam("enableHarmonicMoments", m_enableHarmonicMoments, "Enables the calculation of the Harmonic moments.", true);
54  addParam("enableJets", m_enableJets, "Enables the calculation of jet-related quantities.", true);
55  addParam("enableSphericity", m_enableSphericity, "Enables the calculation of the sphericity-related quantities.", true);
56  addParam("enableCleoCones", m_enableCleoCones, "Enables the calculation of the CLEO cones.", true);
57  addParam("enableAllMoments", m_enableAllMoments, "Enables the calculation of FW and harmonic moments from 5 to 8", false);
58  addParam("checkForDuplicates", m_checkForDuplicates,
59  "Enables the check for duplicates in the input list. If a duplicate entry is found, the first one is kept.", false);
60 }
61 
62 
64 {
65  m_eventShapeContainer.registerInDataStore();
66 }
67 
68 
70 {
71 
73  double sqrtS = T.getCMSEnergy();
74 
76 
78 
79 
80  // --------------------
81  // Calculates the FW moments
82  // --------------------
83  if (m_enableFW) {
84  FoxWolfram fw(m_p3List);
85  if (m_enableAllMoments) {
87  for (short i = 0; i < 9; i++) {
88  m_eventShapeContainer->setFWMoment(i, fw.getH(i));
89  }
90  } else {
92  for (short i = 0; i < 5; i++) {
93  m_eventShapeContainer->setFWMoment(i, fw.getH(i));
94  }
95  }
96  }
97 
98  // --------------------
99  // Calculates the sphericity quantities
100  // --------------------
101  if (m_enableSphericity) {
103  Sph.calculateEigenvalues();
104  if (Sph.getEigenvalue(0) < Sph.getEigenvalue(1) || Sph.getEigenvalue(0) < Sph.getEigenvalue(2)
105  || Sph.getEigenvalue(1) < Sph.getEigenvalue(2))
106  B2WARNING("Eigenvalues not ordered!!!!!!!!!!");
107 
108  for (short i = 0; i < 3; i++) {
109  m_eventShapeContainer->setSphericityEigenvalue(i, Sph.getEigenvalue(i));
110  m_eventShapeContainer->setSphericityEigenvector(i, Sph.getEigenvector(i));
111  }
112  }
113 
114 
115  // --------------------
116  // Calculates thrust and thrust-related quantities
117  // --------------------
118  if (m_enableThrust) {
119  ROOT::Math::XYZVector thrust = Thrust::calculateThrust(m_p3List);
120  float thrustVal = thrust.R();
121  thrust = thrust.Unit();
122  m_eventShapeContainer->setThrustAxis(thrust);
123  m_eventShapeContainer->setThrust(thrustVal);
124 
125  // --- If required, calculates the HarmonicMoments ---
127  HarmonicMoments MM(m_p3List, thrust);
128  if (m_enableAllMoments) {
129  MM.calculateAllMoments();
130  for (short i = 0; i < 9; i++) {
131  auto moment = MM.getMoment(i, sqrtS);
132  m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
133  }
134  } else {
136  for (short i = 0; i < 5; i++) {
137  auto moment = MM.getMoment(i, sqrtS);
138  m_eventShapeContainer->setHarmonicMomentThrust(i, moment);
139  }
140  }
141  }
142 
143  // --- If required, calculates the cleo cones w/ respect to the thrust axis ---
144  if (m_enableCleoCones) {
145  // Cleo cone class constructor. Unfortunately this class is designed
146  // to use the ROE, so the constructor takes two std::vector of momenta ("all" and "ROE"),
147  // then a vector to be used as axis, and finally two flags that determine if the cleo cones
148  // are calculated using the ROE, all the particles or both. Here we use the m_p3List as dummy
149  // list of the ROE momenta, that is however not used at all since the calculate only the
150  // cones with all the particles. This whole class would need some heavy restructuring...
151  CleoCones cleoCones(m_p3List, m_p3List, thrust, true, false);
152  std::vector<float> cones;
153  cones = cleoCones.cleo_cone_with_all();
154  for (short i = 0; i < 10; i++) {
155  m_eventShapeContainer->setCleoConeThrust(i, cones[i]);
156  }
157  } // end of if m_enableCleoCones
158 
159 
160  // --- If required, calculates the jet 4-momentum using the thrust axis ---
161  if (m_enableJets) {
162  ROOT::Math::PxPyPzEVector p4FWD;
163  ROOT::Math::PxPyPzEVector p4BKW;
164  for (const auto& p4 : m_p4List) {
165  if (p4.Vect().Dot(thrust) > 0)
166  p4FWD += p4;
167  else
168  p4BKW += p4;
169  }
170  m_eventShapeContainer->setForwardHemisphere4Momentum(p4FWD);
171  m_eventShapeContainer->setBackwardHemisphere4Momentum(p4BKW);
172  } // end of if m_enableJets
173  }// end of if m_enableThrust
174 
175 
176 
177  // --------------------
178  // Calculates the collision axis quantities
179  // --------------------
180  if (m_enableCollisionAxis) {
181 
182  ROOT::Math::XYZVector collisionAxis(0., 0., 1.);
183 
184  // --- If required, calculates the cleo cones w/ respect to the collision axis ---
185  if (m_enableCleoCones) {
186  CleoCones cleoCones(m_p3List, m_p3List, collisionAxis, true, false);
187  std::vector<float> cones;
188  cones = cleoCones.cleo_cone_with_all();
189  for (short i = 0; i < 10; i++) {
190  m_eventShapeContainer->setCleoConeCollision(i, cones[i]);
191  }
192  }
193 
194  // --- If required, calculates the HarmonicMoments ---
196  HarmonicMoments MM(m_p3List, collisionAxis);
197  if (m_enableAllMoments) {
198  MM.calculateAllMoments();
199  for (short i = 0; i < 9; i++) {
200  auto moment = MM.getMoment(i, sqrtS);
201  m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
202  }
203  } else {
205  for (short i = 0; i < 5; i++) {
206  auto moment = MM.getMoment(i, sqrtS);
207  m_eventShapeContainer->setHarmonicMomentCollision(i, moment);
208  }
209  }
210  } // end of m_enableHarmonicMoments
211  } // end of m_enableCollisionAxis
212 } // end of event()
213 
214 
215 
216 int EventShapeCalculatorModule::parseParticleLists(vector<string> particleListNames)
217 {
218  int nPart = 0; // number of particles
219 
221  m_p4List.clear();
222  m_p3List.clear();
223 
224  unsigned short nParticleLists = particleListNames.size();
225  if (nParticleLists == 0)
226  B2ERROR("No particle lists found. EventShape calculation not performed.");
227 
228  // This vector temporary stores the 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<Particle> tmpParticles;
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  for (unsigned int iPart = 0; iPart < particleList->getListSize(); iPart++) {
241  const Particle* part = particleList->getParticle(iPart);
242 
243  // Flag to check for duplicates across the lists.
244  // It can be true only if m_checkForDuplicates is enabled
245  bool isDuplicate = false;
246 
247  if (m_checkForDuplicates) {
248  // loops over all the particles loaded so far
249  for (const auto& testPart : tmpParticles) {
250  if (testPart.isCopyOf(part)) {
251  B2WARNING("Duplicate particle found. The new one won't be used for the calculation of the event shape variables. Please, double check your input lists and try to make them mutually exclusive.");
252  isDuplicate = true;
253  break;
254  }
255  }
256  tmpParticles.push_back(*part);
257  }
258 
259  if (!isDuplicate) {
260  ROOT::Math::PxPyPzEVector p4CMS = T.rotateLabToCms() * part->get4Vector();
261  // it need to fill an std::vector of XYZVector to use the current FW routines.
262  // It will hopefully change in release 3
263  m_p4List.push_back(p4CMS);
264  m_p3List.push_back(p4CMS.Vect());
265  }
266  }
267  }
268  return nPart;
269 }
270 
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.
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.
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:74
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:532
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:95
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
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:23