Belle II Software  release-05-02-19
ParticleCombinerFromMCModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Yo Sato *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <analysis/modules/ParticleCombinerFromMC/ParticleCombinerFromMCModule.h>
13 
14 // framework aux
15 #include <framework/logging/Logger.h>
16 
17 // dataobjects
18 #include <analysis/dataobjects/Particle.h>
19 
20 // decay descriptor
21 #include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
22 
23 // utilities
24 #include <analysis/DecayDescriptor/ParticleListName.h>
25 #include <analysis/utility/PCmsLabTransform.h>
26 #include <analysis/utility/MCMatching.h>
27 
28 #include <memory>
29 
30 using namespace std;
31 
32 namespace Belle2 {
38 //-----------------------------------------------------------------
39 // Register module
40 //-----------------------------------------------------------------
41 
42  REG_MODULE(ParticleCombinerFromMC)
43 
44 //-----------------------------------------------------------------
45 // Implementation
46 //-----------------------------------------------------------------
47 
49  Module()
50 
51  {
52  // set module description (e.g. insert text)
53  setDescription("Makes particle combinations");
54  setPropertyFlags(c_ParallelProcessingCertified);
55 
56  // Add parameters
57  addParam("decayString", m_decayString,
58  "Input DecayDescriptor string (see :ref:`DecayString`).");
59  addParam("cut", m_cutParameter, "Selection criteria to be applied", std::string(""));
60  addParam("decayMode", m_decayModeID, "User-specified decay mode identifier (saved in 'decayModeID' extra-info for each Particle)",
61  0);
62 
63  addParam("writeOut", m_writeOut,
64  "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
65  addParam("chargeConjugation", m_chargeConjugation,
66  "If true, the charge-conjugated mode will be reconstructed as well", true);
67 
68  // initializing the rest of private members
69  m_pdgCode = 0;
70  m_isSelfConjugatedParticle = false;
71  m_generator = nullptr;
72  m_cut = nullptr;
73  }
74 
75  void ParticleCombinerFromMCModule::initialize()
76  {
77  // clear everything
78  m_pdgCode = 0;
79  m_listName = "";
80 
81  // obtain the input and output particle lists from the decay string
82  bool valid = m_decaydescriptor.init(m_decayString);
83  if (!valid)
84  B2ERROR("Invalid input DecayString: " << m_decayString);
85 
86  // Mother particle
87  const DecayDescriptorParticle* mother = m_decaydescriptor.getMother();
88 
89  m_pdgCode = mother->getPDGCode();
90  m_listName = mother->getFullName();
91 
92  m_antiListName = ParticleListName::antiParticleListName(m_listName);
93  m_isSelfConjugatedParticle = (m_listName == m_antiListName);
94 
95  m_cut = Variable::Cut::compile(m_cutParameter);
96 
97  // register particles which have (sub-)decay recursively
98  registerParticleRecursively(m_decaydescriptor);
99 
100  }
101 
102  void ParticleCombinerFromMCModule::event()
103  {
104  B2DEBUG(10, "event() started !!!");
105 
106  // combine particles recursively
107  combineRecursively(m_decaydescriptor);
108 
109  StoreObjPtr<ParticleList> plist(m_listName);
110  bool existingList = plist.isValid();
111 
112  if (!existingList) {
113  B2WARNING("Output list " << m_listName << " was not created");
114  return;
115  }
116 
117  // loop over list only if cuts should be applied
118  if (!m_cutParameter.empty()) {
119  std::vector<unsigned int> toRemove;
120  unsigned int n = plist->getListSize();
121  for (unsigned i = 0; i < n; i++) {
122  const Particle* part = plist->getParticle(i);
123 
124  MCMatching::setMCTruth(part);
125  MCMatching::getMCErrors(part);
126 
127  if (!m_cut->check(part)) toRemove.push_back(part->getArrayIndex());
128  }
129  plist->removeParticles(toRemove);
130  }
131  }
132 
133 
134  void ParticleCombinerFromMCModule::registerParticleRecursively(const DecayDescriptor& decaydescriptor)
135  {
136  int nProducts = decaydescriptor.getNDaughters();
137 
138  for (int i = 0; i < nProducts; ++i) {
139  if (decaydescriptor.getDaughter(i)->getNDaughters() == 0) {
140  // if daughter does not have daughters, it must be already registered
141  const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
142  StoreObjPtr<ParticleList>().isRequired(daughter->getFullName());
143  } else {
144  // if daughter has daughters, call the function recursively
145  registerParticleRecursively(*(decaydescriptor.getDaughter(i)));
146  }
147  }
148 
149  // Mother particle
150  const DecayDescriptorParticle* mother = decaydescriptor.getMother();
151  std::string listName = mother->getFullName();
152  StoreObjPtr<ParticleList> particleList(listName);
153 
154  // if particleList already exists
155  auto existList = std::find(m_vector_listName.begin(), m_vector_listName.end(), listName);
156  if (existList != m_vector_listName.end()) {
157  B2ERROR(listName << " already exist ! You cannot write same sub-decay twice !!!");
158  return;
159  }
160 
161  std::string antiListName = ParticleListName::antiParticleListName(listName);
162  bool isSelfConjugatedParticle = (listName == antiListName);
163 
164  DataStore::EStoreFlags flags = m_writeOut ? DataStore::c_WriteOut : DataStore::c_DontWriteOut;
165  particleList.registerInDataStore(flags);
166  if (!isSelfConjugatedParticle && m_chargeConjugation) {
167  StoreObjPtr<ParticleList> antiParticleList(antiListName);
168  antiParticleList.registerInDataStore(flags);
169  }
170 
171  m_vector_listName.push_back(listName);
172 
173  }
174 
175  void ParticleCombinerFromMCModule::combineRecursively(const DecayDescriptor& decaydescriptor)
176  {
177  // Mother particle
178  const DecayDescriptorParticle* mother = decaydescriptor.getMother();
179  int pdgCode = mother->getPDGCode();
180  std::string listName = mother->getFullName();
181  std::string antiListName = ParticleListName::antiParticleListName(listName);
182  bool isSelfConjugatedParticle = (listName == antiListName);
183 
184  StoreArray<Particle> particles;
185 
186  StoreObjPtr<ParticleList> outputList(listName);
187  outputList.create();
188  outputList->initialize(pdgCode, listName);
189 
190  if (!isSelfConjugatedParticle && m_chargeConjugation) {
191  StoreObjPtr<ParticleList> outputAntiList(antiListName);
192  outputAntiList.create();
193  outputAntiList->initialize(-1 * pdgCode, antiListName);
194 
195  outputList->bindAntiParticleList(*(outputAntiList));
196  }
197 
198  unsigned int numberOfLists = decaydescriptor.getNDaughters();
199 
200  for (unsigned int i = 0; i < numberOfLists; ++i) {
201  const DecayDescriptor* dDaughter = decaydescriptor.getDaughter(i);
202 
203  if (dDaughter->getNDaughters() == 0) {
204  // if daughter does not have daughters, check if it is not reconstructed particle.
205  const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
206  StoreObjPtr<ParticleList> plist(daughter->getFullName());
207  // if daughter is not created, returns error
208  if (!plist.isValid())
209  B2ERROR(daughter->getFullName() << " is not created");
210  // if daughter contains reconstructed particles, returns error.
211  unsigned nPart = plist->getListSize();
212  for (unsigned iPart = 0; iPart < nPart; iPart++) {
213  const Particle* part = plist->getParticle(iPart);
214  Particle::EParticleSourceObject particleType = part->getParticleSource();
215  if (particleType == Particle::c_Track or
216  particleType == Particle::c_ECLCluster or
217  particleType == Particle::c_KLMCluster)
218  B2ERROR(daughter->getFullName() << " contains a reconstructed particle! It is not accepted in ParticleCombinerFromMCModule!");
219  }
220  // if not, do nothing.
221  } else {
222  // if daughter has daughter, call the function recursively
223  combineRecursively(*dDaughter);
224  }
225  }
226 
227  // initialize the generator
228  m_generator = std::make_unique<ParticleGenerator>(decaydescriptor, "");
229  m_generator->init();
230 
231  while (m_generator->loadNext(m_chargeConjugation)) {
232  Particle&& particle = m_generator->getCurrentParticle();
233 
234  Particle* newParticle = particles.appendNew(particle);
235  // append to the created particle the user specified decay mode ID
236  newParticle->addExtraInfo("decayModeID", m_decayModeID);
237 
238  int iparticle = particles.getEntries() - 1;
239  outputList->addParticle(iparticle, particle.getPDGCode(), particle.getFlavorType());
240  }
241 
242  // select only signal particles
243  std::unique_ptr<Variable::Cut> cutIsSignal = Variable::Cut::compile("isSignal");
244 
245  std::vector<unsigned int> toRemove;
246  unsigned int n = outputList->getListSize();
247  for (unsigned i = 0; i < n; i++) {
248  const Particle* part = outputList->getParticle(i);
249 
250  MCMatching::setMCTruth(part);
251  MCMatching::getMCErrors(part);
252 
253  if (!cutIsSignal->check(part)) toRemove.push_back(part->getArrayIndex());
254  }
255  outputList->removeParticles(toRemove);
256 
257  }
258 
260 } // end Belle2 namespace
261 
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::Particle::addExtraInfo
void addExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1224
Belle2::DataStore::EStoreFlags
EStoreFlags
Flags describing behaviours of objects etc.
Definition: DataStore.h:71
Belle2::DecayDescriptorParticle
Represents a particle in the DecayDescriptor.
Definition: DecayDescriptorParticle.h:37
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::DecayDescriptor::getDaughter
const DecayDescriptor * getDaughter(int i) const
return i-th daughter (0 based index).
Definition: DecayDescriptor.h:146
Belle2::DecayDescriptor::getNDaughters
int getNDaughters() const
return number of direct daughters.
Definition: DecayDescriptor.h:141
Belle2::Particle::EParticleSourceObject
EParticleSourceObject
particle source enumerators
Definition: Particle.h:84
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::DecayDescriptor
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
Definition: DecayDescriptor.h:43
Belle2::StoreArray< Particle >
Belle2::DecayDescriptor::getMother
const DecayDescriptorParticle * getMother() const
return mother.
Definition: DecayDescriptor.h:136
Belle2::ParticleCombinerFromMCModule
particle combiner module
Definition: ParticleCombinerFromMCModule.h:41
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120