Belle II Software light-2406-ragdoll
TreeFitterModule.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// Implementation of Decay Tree Fitter based on arXiv:physics/0503191
10// Main module implementation
11
12#include <analysis/modules/TreeFitter/TreeFitterModule.h>
13#include <analysis/VertexFitting/TreeFitter/FitManager.h>
14
15#include <framework/datastore/StoreArray.h>
16#include <framework/particledb/EvtGenDatabasePDG.h>
17#include <framework/database/DBObjPtr.h>
18
19#include <analysis/utility/ParticleCopy.h>
20#include <analysis/utility/PCmsLabTransform.h>
21
22#include <analysis/VertexFitting/TreeFitter/FitParameterDimensionException.h>
23
24using namespace Belle2;
25
26REG_MODULE(TreeFitter);
27
28TreeFitterModule::TreeFitterModule() : Module(), m_nCandidatesBeforeFit(-1), m_nCandidatesAfter(-1)
29{
30 setDescription("Tree Fitter module. Performs simultaneous fit of all vertices in a decay chain. Can also be used to just fit a single vertex.");
32 //
33 addParam("particleList", m_particleList,
34 "Type::[string]. Input mother of the decay tree to fit. For example 'B0:myB0particleList'.");
35 addParam("confidenceLevel", m_confidenceLevel,
36 "Type::[double]. Confidence level to accept fitted decay tree. Candidates with < confidenceLevel will be removed from the particle list! Typical Values: -1: keep all particle candidates, 0: remove all that fail the fit, 0.001: standard cut, 0.1: (too) tight cut. Optimise using a figure of merit (for example S/(sqrt{S+B}) ) for your analysis.",
37 0.0);
38 addParam("convergencePrecision", m_precision,
39 "Type::[double]. Fractional upper limit for chi2 fluctuations to accept result. Larger value = less signal rejection but also less background rejection. Optimized for FOM on different topologies - don't touch unless you REALLY want this.",
40 0.01);
41 addParam("massConstraintList", m_massConstraintList,
42 "Type::[int]. List of particles to mass constrain with int = pdg code. Note that the variables 'M': fit result for the particle and 'InvM': calculated from the daughter momenta, will look different (especially if you don't update the daughters!).", {});
43 addParam("massConstraintListParticlename", m_massConstraintListParticlename,
44 "Type::[string]. List of particles to mass constrain with string = particle name.", {});
45
46
47 addParam("geoConstraintList", m_geoConstraintListPDG,
48 "Type::[int], if 'autoSetGeoConstraintAndMergeVertices==False' you can manually set the particles that will be geometrically constrained here.", {});
49 addParam("sharedVertexList", m_fixedToMotherVertexListPDG,
50 "Type::[int], if 'autoSetGeoConstraintAndMergeVertices==False' you can manually set the particles that share the vertex with their mother here.", {});
51 addParam("autoSetGeoConstraintAndMergeVertices", m_automatic_vertex_constraining,
52 "Type::bool, shall vertices of strong resonance be merged with their mothers? Can the particles vertex be constraint geometrically?",
53 true);
54
55 addParam("customOriginVertex", m_customOriginVertex,
56 "Type::[double]. List of vertex coordinates to be used in the custom origin constraint.", {0.001, 0, 0.0116});
57 addParam("customOriginCovariance", m_customOriginCovariance,
58 "Type::[double]. List vertex covariance elements used in the custom origin constraint (as a vector). Default is meant for B0 decays and is taken from 100k generated B0 to mumu events.",
59 {
60 0.0048, 0, 0,
61 0, 0.003567, 0,
62 0, 0, 0.0400
63 }
64 );
65 addParam("customOriginConstraint", m_customOrigin,
66 "Type::[bool]. Use a custom vertex as the production point of the highest hierarchy particle (register this as the mother of the list you specify). Like the beam constraint but you can specify the position its covariance yourself. ",
67 false);
68 addParam("ipConstraint", m_ipConstraint,
69 "Type::[bool]. Use the IP as the origin of the tree. This registers an internal IP particle as the mother of the list you give. Or in other words forces the PRODUCTION vertex of your particle to be the IP and its covariance as specified in the database.",
70 false);
71 addParam("originDimension", m_originDimension,
72 "Type int, default 3. If origin or ip constraint used, specify the dimension of the constraint 3->x,y,z; 2->x,y. This also changes the dimension of the geometric constraints! So you might want to turn them off for some particles. (That means turn auto off and manually on for the ones you want to constrain)",
73 3);
74 addParam("updateAllDaughters", m_updateDaughters,
75 "Type::[bool]. Update all daughters (vertex position and momenta) in the tree. If not set only the 4-momenta for the head of the tree will be updated. We also update the vertex position of the daughters regardless of what you put here, because otherwise the default when the particle list is created is {0,0,0}.",
76 false);
77 addParam("expertBeamConstraintPDG", m_beamConstraintPDG,
78 "Type int, default 0. The 4-momentum of particles with the given PDG will be constrained to the 4-momentum of the initial e+e- system.",
79 0);
80 addParam("expertMassConstraintType", m_massConstraintType,
81 "Type::[int]. False(0): use particles parameters in mass constraint; True: use sum of daughter parameters for mass constraint. WAARNING not even guaranteed that it works.",
82 0);
83 addParam("expertRemoveConstraintList", m_removeConstraintList,
84 "Type::[string]. List of constraints that you do not want to be used in the fit. WARNING don't use if you don't know exactly what it does.", {});
85 addParam("expertUseReferencing", m_useReferencing,
86 "Type::[bool]. Use the Extended Kalman Filter. This implementation linearises around the previous state vector which gives smoother convergence.",
87 true);
88 addParam("inflationFactorCovZ", m_inflationFactorCovZ,
89 "Inflate the covariance of the beamspot by this number so that the 3d beam constraint becomes weaker in Z.And: thisnumber->infinity : dim(beamspot constr) 3d->2d.",
90 1);
91 addParam("treatAsInvisible", m_treatAsInvisible,
92 "Type::[string]. Decay string to select one particle that will be ignored in the fit.", {});
93
94 addParam("ignoreFromVertexFit", m_ignoreFromVertexFit,
95 "Type::[string]. Decay string to select particles that will be ignored to determine the vertex position while kept for kinematics determination.", {});
96}
97
99{
100 m_plist.isRequired(m_particleList);
104
105 if ((m_massConstraintList.size()) == 0 && (m_massConstraintListParticlename.size()) > 0) {
106 for (auto& containedParticle : m_massConstraintListParticlename) {
107 TParticlePDG* particletemp = TDatabasePDG::Instance()->GetParticle((containedParticle).c_str());
108 m_massConstraintList.push_back(particletemp->PdgCode());
109 }
110 }
111
112 if (!m_treatAsInvisible.empty()) {
114 if (!valid)
115 B2ERROR("TreeFitterModule::initialize Invalid Decay Descriptor: " << m_treatAsInvisible);
116 else if (m_pDDescriptorInvisibles.getSelectionPDGCodes().size() != 1)
117 B2ERROR("TreeFitterModule::initialize Please select exactly one particle to ignore: " << m_treatAsInvisible);
118 }
119
120 if (!m_ignoreFromVertexFit.empty()) {
122 if (!valid)
123 B2ERROR("TreeFitterModule::initialize Invalid Decay Descriptor: " << m_ignoreFromVertexFit);
124 }
125
126}
127
129{
131 const ROOT::Math::PxPyPzEVector cms = T.getBeamFourMomentum();
132
133 m_beamMomE(0) = cms.X();
134 m_beamMomE(1) = cms.Y();
135 m_beamMomE(2) = cms.Z();
136 m_beamMomE(3) = cms.E();
137
138 m_beamCovariance = Eigen::Matrix4d::Zero();
139 const double covE = 3.19575e-05; // TODO Avoid using the hard coded value.
140 // It is taked from the BeamParameter which was used previously. The uncertainty from the CollisionInvariantMass is a possibility.
141
142 for (size_t i = 0; i < 4; ++i) {
143 m_beamCovariance(i, i) = covE;
144 // TODO Currently, we do not get a full covariance matrix from beamparams, and the py value is zero, which means there is no constraint on py. Therefore, we approximate it by a diagonal matrix using the energy value for all components. This is based on the assumption that the components of the beam four-momentum are independent and of comparable size.
145 }
146}
147
149{
150 if (!m_plist) {
151 B2ERROR("ParticleList " << m_particleList << " not found");
152 return;
153 }
154
155 std::vector<unsigned int> toRemove;
156 const unsigned int nParticles = m_plist->getListSize();
157 m_nCandidatesBeforeFit += nParticles;
158
159 TMatrixFSym dummyCovMatrix(7);
160 for (int row = 0; row < 7; ++row) { //diag
161 dummyCovMatrix(row, row) = 10000;
162 }
163
164 TMatrixFSym dummyCovMatrix_smallMomError(dummyCovMatrix);
165 for (int row = 0; row < 4; ++row) {
166 dummyCovMatrix_smallMomError(row, row) = 1e-10;
167 }
168
169 for (unsigned iPart = 0; iPart < nParticles; iPart++) {
170 Particle* particle = m_plist->getParticle(iPart);
171
172 if (m_updateDaughters == true) {
174 }
175
176 if (!m_treatAsInvisible.empty()) {
177 std::vector<const Particle*> selParticlesTarget = m_pDDescriptorInvisibles.getSelectionParticles(particle);
178 Particle* targetD = m_particles[selParticlesTarget[0]->getArrayIndex()];
179 Particle* daughterCopy = ParticleCopy::copyParticle(targetD);
180 daughterCopy->writeExtraInfo("treeFitterTreatMeAsInvisible", 1);
181 daughterCopy->setMomentumVertexErrorMatrix(dummyCovMatrix);
182
183 bool isReplaced = particle->replaceDaughterRecursively(targetD, daughterCopy);
184 if (!isReplaced)
185 B2ERROR("TreeFitterModule::event No target particle found for " << m_treatAsInvisible);
186 }
187
188 if (!m_ignoreFromVertexFit.empty()) {
189 std::vector<const Particle*> selParticlesTarget = m_pDDescriptorForIgnoring.getSelectionParticles(particle);
190
191 for (auto part : selParticlesTarget) {
192 Particle* targetD = m_particles[part->getArrayIndex()];
193 Particle* daughterCopy = ParticleCopy::copyParticle(targetD);
194 daughterCopy->writeExtraInfo("treeFitterTreatMeAsInvisible", 1);
195 daughterCopy->setMomentumVertexErrorMatrix(dummyCovMatrix_smallMomError);
196
197 bool isReplaced = particle->replaceDaughterRecursively(targetD, daughterCopy);
198 if (!isReplaced)
199 B2ERROR("TreeFitterModule::event No target particle found for " << m_ignoreFromVertexFit);
200 }
201 }
202
203 try {
204 const bool ok = fitTree(particle);
205 if (!ok) { particle->setPValue(-1); }
207 B2ERROR(e.what());
208 }
209
210 if (particle->getPValue() < m_confidenceLevel) {
211 toRemove.push_back(particle->getArrayIndex());
212 }
213
214 }
215 m_plist->removeParticles(toRemove);
216 m_nCandidatesAfter += m_plist->getListSize();
217}
218
219
221{
222 if (m_nCandidatesAfter > 0) {
224 } else {
225 B2WARNING("Not a single candidate survived the fit. Candidates before fit: " << m_nCandidatesBeforeFit << " after: " <<
227 }
228}
229
231{
232 const TreeFitter::ConstraintConfiguration constrConfig(
248 );
249
250 std::unique_ptr<TreeFitter::FitManager> TreeFitter(
252 head,
253 constrConfig,
257 )
258 );
259 bool rc = TreeFitter->fit();
260 return rc;
261}
262
264{
265 B2INFO("\033[1;35m================================================================================\033[0m");
266 B2INFO("\033[40;97m ,., \033[0m");
267 B2INFO("\033[40;97m ;%&M%;_ ,.., \033[0m");
268 B2INFO("\033[40;97m \"_ \"__\" % M % M %; , ..., , \033[0m");
269 B2INFO("\033[40;97m , ..., __.\" --\" , ., _ - \" % &W % WM %; \033[0m");
270 B2INFO("\033[40;97m ; % M&$M % \"___ \"_._ %M%\"_.\"\" _ \"\"\"\"\"\" \033[0m");
271 B2INFO("\033[40;97m \"\"\"\"\" \"\" , \\_. \"_. .\" \033[0m");
272 B2INFO("\033[40;97m , ., _\"__ \\__./ .\" \033[0m");
273 B2INFO("\033[40;97m ___ __ | y , .., \033[97;40mThank you for using TreeFitter. \033[0m");
274 B2INFO("\033[40;97m /)'\\ ''''| u \\ %W%W%%; \033[0m");
275 B2INFO("\033[40;97m ___)/ \"---\\_ \\ |____\" \033[97;40mPlease cite both: \033[0m");
276 B2INFO("\033[40;97m ;&&%%; (|__.|)./ ,.., \033[97;40m10.1016/j.nima.2020.164269 \033[0m");
277 B2INFO("\033[40;97m ,.., ___\\ |/ &&\" \033[97;40m 10.1016/j.nima.2005.06.078 \033[0m");
278 B2INFO("\033[40;97m &&%%& (| Uo / '\" \033[97;40mEmail: \033[0m");
279 B2INFO("\033[40;97m '''' \\ 7 \\ \033[97;40mfrancesco.tenchini@desy.de \033[0m");
280 B2INFO("\033[40;97m ._______________.-'____\"\"-.____. \033[97;40mjo-frederik.krohn@desy.de \033[0m");
281 B2INFO("\033[40;97m \\ / \033[0m");
282 B2INFO("\033[40;97m \\ \033[0m\033[32;40mTREEFITTER\033[0m\033[40;97m / \033[0m");
283 B2INFO("\033[40;97m \\_______________________/ \033[0m");
284 B2INFO("\033[40;97m (_) (_) \033[0m");
285 B2INFO("\033[40;97m \033[0m");
286 B2INFO("\033[1;35m============= TREEFIT STATISTICS ===============================================\033[0m");
287 B2INFO("\033[1;39mTarget particle list: " << m_particleList << "\033[0m");
288 B2INFO("\033[1;39mCandidates before fit: " << m_nCandidatesBeforeFit << "\033[0m");
289 B2INFO("\033[1;39mCandidates after fit: " << m_nCandidatesAfter << "\033[0m");
290 B2INFO("\033[1;39mA total of " << m_nCandidatesBeforeFit - m_nCandidatesAfter <<
291 " candidates were removed during the fit.\033[0m");
292 B2INFO("\033[1;39m" << (double)m_nCandidatesAfter / (double)m_nCandidatesBeforeFit * 100.0 <<
293 "% of candidates survived the fit.\033[0m");
294 B2INFO("\033[1;39m" << 100. - (double)m_nCandidatesAfter / (double)m_nCandidatesBeforeFit * 100.0 <<
295 "% of candidates did not.\033[0m");
296 B2INFO("\033[1;39mYou chose to drop all candidates with pValue < " << m_confidenceLevel << ".\033[0m");
297 B2INFO("\033[1;35m================================================================================\033[0m");
298}
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
std::vector< int > getSelectionPDGCodes()
Return list of PDG codes of selected particles.
std::vector< const Particle * > getSelectionParticles(const Particle *particle)
Get a vector of pointers with selected daughters in the decay tree.
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.
ROOT::Math::PxPyPzEVector getBeamFourMomentum() const
Returns LAB four-momentum of e+e-, i.e.
Class to store reconstructed particles.
Definition: Particle.h:75
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1308
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:667
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:393
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:366
bool replaceDaughterRecursively(const Particle *oldDaughter, Particle *newDaughter)
Apply replaceDaughter to all Particles in the decay tree by looping recursively through it,...
Definition: Particle.cc:723
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
DecayDescriptor m_pDDescriptorForIgnoring
Decay descriptor of the ignored particles.
void plotFancyASCII()
plot ascii art and statistics
std::string m_ignoreFromVertexFit
decay string to select one particle that will be ignored to determine the vertex position
std::vector< double > m_customOriginVertex
vertex coordinates of the custom origin
virtual void initialize() override
initialize
virtual void event() override
performed for each event
Eigen::Matrix< double, 4, 1 > m_beamMomE
beam four-momentum
bool m_updateDaughters
flag if you want to update all particle momenta in the decay tree.
DecayDescriptor m_pDDescriptorInvisibles
Decay descriptor of the invisible particles.
std::vector< std::string > m_massConstraintListParticlename
vector carrying the names of the particles to be mass constraint
unsigned int m_nCandidatesAfter
after the fit
std::vector< int > m_massConstraintList
vector carrying the PDG codes of the particles to be mass constraint
std::vector< int > m_fixedToMotherVertexListPDG
list of pdg codes of particles where we use the same vertex for production and decay which is the ver...
StoreArray< Particle > m_particles
StoreArray of Particles.
virtual void terminate() override
stuff at the end
bool m_automatic_vertex_constraining
should the vertex be joined with the mother and should it be geometrically constrained?...
virtual void beginRun() override
performed at the start of run
bool m_ipConstraint
Use x-y-z beamspot constraint.
int m_originDimension
dimension to use for beam/origin constraint
int m_inflationFactorCovZ
inflate beamspot covariance of z by this number
std::vector< std::string > m_removeConstraintList
list of constraints not to be applied in tree fit WARNING only use if you know what you are doing
int m_massConstraintType
type of the mass constraint false: use normal one.
std::vector< int > m_geoConstraintListPDG
list of pdg codes of particles to use a geo constraint for
std::vector< double > m_customOriginCovariance
covariance of the custom origin
std::string m_particleList
name of the particle list fed to the fitter
double m_confidenceLevel
minimum confidence level to accept fit calculated as f(chiSquared, NDF) -2: accept all 0: only accept...
bool m_customOrigin
use a custom vertex as the production vertex of the highest hierarchy particle
std::string m_treatAsInvisible
decay string to select one particle that will be treated as invisible
Eigen::Matrix< double, 4, 4 > m_beamCovariance
beam covariance matrix
int m_beamConstraintPDG
PDG code of particle to be constrained to the beam 4-momentum.
bool m_useReferencing
linearise around a previous state of the Kalman Filter
bool fitTree(Particle *head)
this fits all particle candidates contained in the m_particleList
double m_precision
convergence precision for the newton method When the delta chiSquared between 2 iterations divided by...
StoreObjPtr< ParticleList > m_plist
input particle list
unsigned int m_nCandidatesBeforeFit
before the fit
constraint configuration class
exception template, runtime_error implements what()
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
void copyDaughters(Particle *mother)
Function copies all (grand-)^n-daughter particles of the argument mother Particle.
Definition: ParticleCopy.cc:56
Particle * copyParticle(const Particle *original)
Function takes argument Particle and creates a copy of it and copies of all its (grand-)^n-daughters.
Definition: ParticleCopy.cc:18
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24