Belle II Software development
ParticleCombinerModule.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// Own header.
10#include <analysis/modules/ParticleCombiner/ParticleCombinerModule.h>
11
12// framework aux
13#include <framework/logging/Logger.h>
14
15// decay descriptor
16#include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
17#include <analysis/DecayDescriptor/ParticleListName.h>
18
19#include <analysis/ParticleCombiner/ParticleCombiner.h>
20
21// utilities
22#include <analysis/utility/EvtPDLUtil.h>
23#include <analysis/utility/PCmsLabTransform.h>
24
25using namespace std;
26using namespace Belle2;
27
28//-----------------------------------------------------------------
29// Register module
30//-----------------------------------------------------------------
31
32REG_MODULE(ParticleCombiner);
33
34//-----------------------------------------------------------------
35// Implementation
36//-----------------------------------------------------------------
37
39 Module()
40
41{
42 // set module description (e.g. insert text)
43 setDescription("Makes particle combinations");
45
46 // Add parameters
47 addParam("decayString", m_decayString,
48 "Input DecayDescriptor string (see :ref:`DecayString`).");
49 addParam("cut", m_cutParameter, "Selection criteria to be applied", std::string(""));
50 addParam("maximumNumberOfCandidates", m_maximumNumberOfCandidates,
51 "Max. number of candidates reconstructed. By default, if the limit is reached no candidates will be produced.\n"
52 "This behaviour can be changed by \'ignoreIfTooManyCandidates\' flag.", 10000);
53
54 addParam("ignoreIfTooManyCandidates", m_ignoreIfTooManyCandidates,
55 "Don't reconstruct channel if more candidates than given by \'maximumNumberOfCandidates\' are produced.", true);
56
57 addParam("decayMode", m_decayModeID, "User-specified decay mode identifier (saved in 'decayModeID' extra-info for each Particle)",
58 0);
59 addParam("writeOut", m_writeOut,
60 "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
61 addParam("recoilParticleType", m_recoilParticleType, R"DOC(If not equal 0,
62the mother Particle is reconstructed in the recoil against the daughter particles.
63In the case of the following decay chain M -> D1 D2 ... Dn and
64
65- recoilParticleType = 1:
66 - the mother momentum is given by: p(M) = p(e+e-) - p(D1) - p(D2) - ... - p(DN)
67 - D1, D2, ..., DN are attached as daughters of M
68- recoilParticleType = 2:
69 - the mother momentum is given by: p(M) = p(D1) - p(D2) - ... - p(DN)
70 - D1, D2, ..., DN are attached as daughters of M
71 )DOC", 0);
72 addParam("chargeConjugation", m_chargeConjugation,
73 "If true, the charge-conjugated mode will be reconstructed as well", true);
74 addParam("allowChargeViolation", m_allowChargeViolation,
75 "If true the decay string does not have to conserve electric charge", false);
76
77 // initializing the rest of private members
78 m_pdgCode = 0;
80 m_generator = nullptr;
81}
82
84{
85 // clear everything
86 m_pdgCode = 0;
87 m_listName = "";
88
89 // obtain the input and output particle lists from the decay string
90 bool valid = m_decaydescriptor.init(m_decayString);
91 if (!valid)
92 B2ERROR("Invalid input DecayString: " << m_decayString);
93
94 // Mother particle
95 const DecayDescriptorParticle* mother = m_decaydescriptor.getMother();
96
97 m_pdgCode = mother->getPDGCode();
98 m_listName = mother->getFullName();
99
102
103 // Daughters
104 int nProducts = m_decaydescriptor.getNDaughters();
105 int daughtersNetCharge = 0;
106 for (int i = 0; i < nProducts; ++i) {
107 const DecayDescriptorParticle* daughter =
108 m_decaydescriptor.getDaughter(i)->getMother();
109 StoreObjPtr<ParticleList>().isRequired(daughter->getFullName());
110 int daughterPDGCode = daughter->getPDGCode();
111 if (m_recoilParticleType == 2 && i == 0) {
112 daughtersNetCharge -= EvtPDLUtil::charge(daughterPDGCode);
113 } else {
114 daughtersNetCharge += EvtPDLUtil::charge(daughterPDGCode);
115 }
116 }
117
118 // The electric charge check makes no sense for an inclusive decay
119 if (!m_decaydescriptor.isIgnoreMassive() and daughtersNetCharge != EvtPDLUtil::charge(m_pdgCode)) {
121 B2FATAL("Your decay string " << m_decayString << " violates electric charge conservation!\n"
122 "If you want to allow this you can set the argument 'allowChargeViolation' to True. Something like:\n"
123 "modularAnalysis.reconstructDecay(" << m_decayString << ", your_cuts, allowChargeViolation=True, path=mypath)");
124 }
125 B2WARNING("Your decay string " << m_decayString << " violates electric charge conservation!\n"
126 "Processing is continued assuming that you allowed this deliberately, e.g. for systematic studies etc.");
127 }
128
129 if (m_recoilParticleType == 0) {
130 m_generator = std::make_unique<ParticleGenerator>(m_decayString, m_cutParameter);
131 } else {
132 string noCutInParticleCombiner = "";
133 m_generator = std::make_unique<ParticleGenerator>(m_decayString, noCutInParticleCombiner);
135 }
136
138 m_outputList.registerInDataStore(m_listName, flags);
140 m_outputAntiList.registerInDataStore(m_antiListName, flags);
141 }
142
144 B2FATAL("Invalid recoil particle type = " << m_recoilParticleType <<
145 "! Valid values are 0 (not a recoil), 1 (recoiling against e+e- and daughters), 2 (daughter of a recoil)");
146}
147
149{
150 m_outputList.create();
151 m_outputList->initialize(m_pdgCode, m_listName);
152
154 m_outputAntiList.create();
155 m_outputAntiList->initialize(-1 * m_pdgCode, m_antiListName);
156
157 m_outputList->bindAntiParticleList(*(m_outputAntiList));
158 }
159
160 m_generator->init();
161
162 int numberOfCandidates = 0;
163 while (m_generator->loadNext(m_chargeConjugation)) {
164
165 Particle&& particle = m_generator->getCurrentParticle();
166
167 // if particle is reconstructed in the recoil,
168 // its 4-momentum vector needs to be fixed
169 // at this stage its 4 momentum is:
170 // p(mother) = Sum_i p(daughter_i)
171 // but it needs to be
172 // a) in the case of recoilParticleType == 1
173 // - p(mother) = p(e-) + p(e+) - Sum_i p(daughter_i)
174 // b) in the case of recoilParticleType == 2
175 // - p(mother) = p(daughter_0) - Sum_i p(daughter_i) (where i > 0)
176 if (m_recoilParticleType == 1) {
178 ROOT::Math::PxPyPzEVector recoilMomentum = T.getBeamFourMomentum() - particle.get4Vector();
179 particle.set4Vector(recoilMomentum);
180 } else if (m_recoilParticleType == 2) {
181 const std::vector<Particle*> daughters = particle.getDaughters();
182
183 if (daughters.size() < 2)
184 B2FATAL("Reconstructing particle as a daughter of a recoil with less then 2 daughters!");
185
186 ROOT::Math::PxPyPzEVector pDaughters;
187 for (unsigned i = 1; i < daughters.size(); i++) {
188 pDaughters += daughters[i]->get4Vector();
189 }
190
191 ROOT::Math::PxPyPzEVector mom = daughters[0]->get4Vector() - pDaughters;
192 particle.set4Vector(mom);
193 }
194 if (m_recoilParticleType > 0) {
195 if (!m_cut->check(&particle)) {
196 continue;
197 }
198 }
199
200 numberOfCandidates++;
201
202 if (m_maximumNumberOfCandidates > 0 and numberOfCandidates > m_maximumNumberOfCandidates) {
204 B2WARNING("Maximum number of " << m_maximumNumberOfCandidates << " candidates reached, skipping event");
205 m_outputList->clear();
206 } else {
207 B2WARNING("Maximum number of " << m_maximumNumberOfCandidates << " candidates reached. Ignoring others");
208 }
209 break;
210 }
211
212 Particle* newParticle = m_particles.appendNew(particle);
213
214 m_outputList->addParticle(newParticle);
215
216 // append to the created particle the user specified decay mode ID
217 newParticle->addExtraInfo("decayModeID", m_decayModeID);
218 }
219}
EStoreFlags
Flags describing behaviours of objects etc.
Definition DataStore.h:69
@ c_WriteOut
Object/array should be saved by output modules.
Definition DataStore.h:70
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition DataStore.h:71
Represents a particle in the DecayDescriptor.
int getPDGCode() const
Return PDG code.
std::string getFullName() const
returns the full name of the particle full_name = name:label
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Definition GeneralCut.h:84
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.
int m_recoilParticleType
type of recoil particle: 0 - not recoil (normal reconstruction); 1 - recoil against e+e- and all daug...
bool m_isSelfConjugatedParticle
flag that indicates whether an anti-particle mother does not exist and should not be reconstructed as...
int m_maximumNumberOfCandidates
maximum number of reconstructed candidates
std::string m_antiListName
output anti-particle list name
virtual void initialize() override
Initialize the Module.
virtual void event() override
Event processor.
std::string m_decayString
Input DecayString specifying the decay being reconstructed.
std::string m_listName
output particle list name
StoreObjPtr< ParticleList > m_outputAntiList
output anti-particle list
StoreArray< Particle > m_particles
StoreArray of Particles.
StoreObjPtr< ParticleList > m_outputList
output particle list
std::unique_ptr< ParticleGenerator > m_generator
Generates the combinations.
std::unique_ptr< Variable::Cut > m_cut
cut object which performs the cuts
int m_decayModeID
user specified decay mode identifier
bool m_allowChargeViolation
switch to turn on and off the requirement of electric charge conservation
std::string m_cutParameter
selection criteria
DecayDescriptor m_decaydescriptor
Decay descriptor of the decay being reconstructed.
bool m_writeOut
toggle output particle list btw.
bool m_ignoreIfTooManyCandidates
drop all candidates if max.
int m_pdgCode
PDG code of the combined mother particle.
bool m_chargeConjugation
boolean to control whether charge conjugated decay should be reconstructed as well
Class to store reconstructed particles.
Definition Particle.h:76
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition Particle.cc:1421
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
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:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition EvtPDLUtil.cc:44
std::string antiParticleListName(const std::string &listName)
Returns name of anti-particle-list corresponding to listName.
Abstract base class for different kinds of events.
STL namespace.