Belle II Software  release-05-02-19
ParticleListManipulatorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Anze Zupanc *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <analysis/modules/ParticleListManipulator/ParticleListManipulatorModule.h>
13 
14 
15 // framework - DataStore
16 #include <framework/datastore/DataStore.h>
17 #include <framework/datastore/StoreArray.h>
18 #include <framework/datastore/StoreObjPtr.h>
19 
20 // framework aux
21 #include <framework/logging/Logger.h>
22 
23 // dataobjects
24 #include <analysis/dataobjects/Particle.h>
25 #include <analysis/dataobjects/ParticleList.h>
26 
27 // utilities
28 #include <analysis/DecayDescriptor/ParticleListName.h>
29 #include <analysis/utility/ValueIndexPairSorting.h>
30 
31 using namespace std;
32 
33 namespace Belle2 {
39  //-----------------------------------------------------------------
40  // Register module
41  //-----------------------------------------------------------------
42 
43  REG_MODULE(ParticleListManipulator)
44 
45  //-----------------------------------------------------------------
46  // Implementation
47  //-----------------------------------------------------------------
48 
50  m_variable(nullptr)
51 
52  {
53  setDescription("Manipulates ParticleLists: copies/merges/performs particle selection");
54  setPropertyFlags(c_ParallelProcessingCertified);
55 
56  // Add parameters
57  addParam("outputListName", m_outputListName, "Output ParticleList name");
58 
59  vector<string> defaultList;
60  addParam("inputListNames", m_inputListNames,
61  "list of input ParticleList names", defaultList);
62 
63  addParam("cut", m_cutParameter, "Selection criteria to be applied", std::string(""));
64 
65  addParam("variable", m_variableName, "Variable which defines the best duplicate (see ``selectLowest`` for ordering)",
66  std::string("mdstIndex"));
67  addParam("preferLowest", m_preferLowest,
68  "If true, duplicate with lowest value of ``variable`` is accepted, otherwise higher one.", true);
69 
70  addParam("writeOut", m_writeOut,
71  "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
72 
73  // initializing the rest of private members
74  m_pdgCode = 0;
75  m_isSelfConjugatedParticle = false;
76  }
77 
78  void ParticleListManipulatorModule::initialize()
79  {
80  m_pdgCode = 0;
81 
82  // check the validity of output ParticleList name
83  bool valid = m_decaydescriptor.init(m_outputListName);
84  if (!valid)
85  B2ERROR("ParticleListManipulatorModule::initialize Invalid output ParticleList name: " << m_outputListName);
86 
87  // Output particle
88  const DecayDescriptorParticle* mother = m_decaydescriptor.getMother();
89 
90  m_pdgCode = mother->getPDGCode();
91 
92  m_outputAntiListName = ParticleListName::antiParticleListName(m_outputListName);
93  m_isSelfConjugatedParticle = (m_outputListName == m_outputAntiListName);
94 
95  // Input lists
96  for (const std::string& listName : m_inputListNames) {
97  if (listName == m_outputListName) {
98  B2ERROR("ParticleListManipulatorModule: cannot copy Particles from " << listName <<
99  " to itself! Use applyCuts() (ParticleSelector module) instead.");
100  } else if (!m_decaydescriptor.init(listName)) {
101  B2ERROR("Invalid input ParticleList name: " << listName);
102  } else {
103  StoreObjPtr<ParticleList>().isRequired(listName);
104  }
105  }
106 
107  StoreObjPtr<ParticleList> particleList(m_outputListName);
108  DataStore::EStoreFlags flags = m_writeOut ? DataStore::c_WriteOut : DataStore::c_DontWriteOut;
109  particleList.registerInDataStore(flags);
110  if (!m_isSelfConjugatedParticle) {
111  StoreObjPtr<ParticleList> antiParticleList(m_outputAntiListName);
112  antiParticleList.registerInDataStore(flags);
113  }
114 
115  m_variable = Variable::Manager::Instance().getVariable(m_variableName);
116  if (!m_variable) {
117  B2ERROR("Variable '" << m_variableName << "' is not available in Variable::Manager!");
118  }
119  m_cut = Variable::Cut::compile(m_cutParameter);
120  }
121 
122  void ParticleListManipulatorModule::event()
123  {
124  // clear the list
125  m_particlesInTheList.clear();
126 
127  const StoreArray<Particle> particles;
128  StoreObjPtr<ParticleList> plist(m_outputListName);
129  bool existingList = plist.isValid();
130 
131  if (!existingList) {
132  // new particle list: create it
133  plist.create();
134  plist->initialize(m_pdgCode, m_outputListName);
135 
136  if (!m_isSelfConjugatedParticle) {
137  StoreObjPtr<ParticleList> antiPlist(m_outputAntiListName);
138  antiPlist.create();
139  antiPlist->initialize(-1 * m_pdgCode, m_outputAntiListName);
140 
141  antiPlist->bindAntiParticleList(*(plist));
142  }
143  } else {
144  // output list already contains Particles
145  // fill m_particlesInTheList with unique
146  // identifiers of particles already in
147  for (unsigned i = 0; i < plist->getListSize(); i++) {
148  const Particle* particle = plist->getParticle(i);
149 
150  std::vector<int> idSeq;
151  fillUniqueIdentifier(particle, idSeq);
152  m_particlesInTheList.push_back(idSeq);
153  }
154  }
155 
156  // create list of candidate indices and corresponding sorting values
157  typedef std::pair<double, unsigned int> ValueIndexPair;
158  std::vector<ValueIndexPair> valueToIndex;
159 
160  // fill all particles from input lists that pass selection criteria into comparison list
161  for (const auto& inputListName : m_inputListNames) {
162  const StoreObjPtr<ParticleList> inPList(inputListName);
163 
164  std::vector<int> fsParticles = inPList->getList(ParticleList::EParticleType::c_FlavorSpecificParticle, false);
165  const std::vector<int>& scParticles = inPList->getList(ParticleList::EParticleType::c_SelfConjugatedParticle, false);
166  const std::vector<int>& fsAntiParticles = inPList->getList(ParticleList::EParticleType::c_FlavorSpecificParticle, true);
167 
168  fsParticles.insert(fsParticles.end(), scParticles.begin(), scParticles.end());
169  fsParticles.insert(fsParticles.end(), fsAntiParticles.begin(), fsAntiParticles.end());
170 
171  for (int fsParticle : fsParticles) {
172  const Particle* part = particles[fsParticle];
173 
174  if (m_cut->check(part)) {
175  valueToIndex.emplace_back(m_variable->function(part), part->getArrayIndex());
176  }
177  }
178  }
179 
180  // use stable sort to make sure we keep the relative order of elements with
181  // same value as it was before
182  if (m_preferLowest) {
183  std::stable_sort(valueToIndex.begin(), valueToIndex.end(), ValueIndexPairSorting::lowerPair<ValueIndexPair>);
184  } else {
185  std::stable_sort(valueToIndex.begin(), valueToIndex.end(), ValueIndexPairSorting::higherPair<ValueIndexPair>);
186  }
187 
188  // starting from the best candidate add all particles to output list that are not already in it
189  for (const auto& candidate : valueToIndex) {
190  const Particle* part = particles[candidate.second];
191 
192  std::vector<int> idSeq;
193  fillUniqueIdentifier(part, idSeq);
194  bool uniqueSeq = isUnique(idSeq);
195 
196  if (uniqueSeq) {
197  plist->addParticle(part);
198  m_particlesInTheList.push_back(idSeq);
199  }
200  }
201  }
202 
203  void ParticleListManipulatorModule::fillUniqueIdentifier(const Particle* p, std::vector<int>& idSequence)
204  {
205  idSequence.push_back(p->getPDGCode());
206 
207  if (p->getNDaughters() == 0) {
208  idSequence.push_back(p->getMdstArrayIndex());
209  } else {
210  idSequence.push_back(p->getNDaughters());
211  auto daughters = p->getDaughters();
212  // sorting the daughters by their pdgCode to identify decay chains only differing by the order of their daughters
213  sort(daughters.begin(), daughters.end(), [](const auto a, const auto b) {
214  return a->getPDGCode() > b->getPDGCode();
215  });
216  // this is not FSP (go one level down)
217  for (const auto& daughter : daughters)
218  fillUniqueIdentifier(daughter, idSequence);
219  }
220  }
221 
222  bool ParticleListManipulatorModule::isUnique(const std::vector<int>& idSeqOUT)
223  {
224  for (const auto& idSeqIN : m_particlesInTheList) {
225  bool sameSeq = (idSeqIN == idSeqOUT);
226  if (sameSeq)
227  return false;
228  }
229 
230  return true;
231  }
233 } // end Belle2 namespace
Belle2::ParticleListManipulatorModule
Module for copying Particles (actually their indices) from two or more ParticleLists(s) to another Pa...
Definition: ParticleListManipulatorModule.h:47
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
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
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::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::StoreArray< Particle >
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120