Belle II Software light-2406-ragdoll
HelixErrorScalerModule.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/modules/HelixErrorScaler/HelixErrorScalerModule.h>
10
11#include <analysis/DecayDescriptor/ParticleListName.h>
12#include <analysis/utility/ParticleCopy.h>
13#include <mdst/dataobjects/HitPatternCDC.h>
14#include <mdst/dataobjects/HitPatternVXD.h>
15#include <framework/datastore/RelationArray.h>
16#include <framework/gearbox/Const.h>
17#include <framework/geometry/B2Vector3.h>
18
19#include <Math/Vector3D.h>
20#include <vector>
21
22using namespace Belle2;
23
24//-----------------------------------------------------------------
25// Register the Module
26//-----------------------------------------------------------------
27REG_MODULE(HelixErrorScaler);
28
29//-----------------------------------------------------------------
30// Implementation
31//-----------------------------------------------------------------
32
33HelixErrorScalerModule::HelixErrorScalerModule() : Module(), m_pdgCode(0), m_scaleKshort(false)
34{
35 // Set module properties
36 setDescription(R"DOC(scale the error of helix parameters
37
38Creates a new charged particle list whose helix errors are scaled by constant factors.
39Different sets of scale factors are defined for tracks with/without a PXD hit.
40For tracks with a PXD hit, in order to avoid severe underestimation of d0 and z0 errors,
41lower limits (best resolution) can be set in a momentum-dependent form.
42The module also accepts a V0 Kshort particle list as input and applies the error correction to its daughters.
43Note the difference in impact parameter resolution between V0 daughters and tracks from IP,
44as V0 daughters are free from multiple scattering through the beam pipe.
45 )DOC");
46
47 // Parameter definitions
48 addParam("inputListName", m_inputListName, "The name of input particle list (charged stable or V0 Kshort)", std::string(""));
49 addParam("outputListName", m_outputListName, "The name of output charged particle list", std::string(""));
50 addParam("scaleFactors_PXD", m_scaleFactors_PXD,
51 "vector of five scale factors for helix parameter errors (for tracks with a PXD hit)", {1.0, 1.0, 1.0, 1.0, 1.0});
52 addParam("scaleFactors_noPXD", m_scaleFactors_noPXD,
53 "vector of five scale factors for helix parameter errors (for tracks without a PXD hit)", {1.0, 1.0, 1.0, 1.0, 1.0});
54 addParam("d0ResolutionParameters", m_d0ResolPars, "d0 best resolution parameters", {0.0, 0.0});
55 addParam("z0ResolutionParameters", m_z0ResolPars, "z0 best resolution parameters", {0.0, 0.0});
56 addParam("d0MomentumThreshold", m_d0MomThr, "d0 best resolution is kept constant below this momentum.", 0.0);
57 addParam("z0MomentumThreshold", m_z0MomThr, "z0 best resolution is kept constant below this momentum.", 0.0);
58
59}
60
62{
63 // check the validity of output ParticleList name
65 if (!valid)
66 B2ERROR("Invalid output ParticleList name: " << m_outputListName);
67
68 // output particle
70 m_pdgCode = mother->getPDGCode();
71 if (m_pdgCode == Const::Kshort.getPDGCode()) {
72 m_scaleKshort = true;
74 B2ERROR("Invalid input ParticleList PDG code (must be ChargedStable): " << m_pdgCode);
75 }
76
77 // get existing particle lists
79 B2ERROR("Input and output particle list names are the same: " << m_inputListName);
81 B2ERROR("Invalid input particle list name: " << m_inputListName);
82 } else {
84 }
85
86 // make output list
87 m_outputparticleList.registerInDataStore(m_outputListName);
88 if (! m_scaleKshort) {
91 }
92
95}
96
98{
99 RelationArray particlesToMCParticles(m_particles, m_mcparticles);
100
101 // new output particle list
102 if (! m_outputparticleList.isValid()) {
103 m_outputparticleList.create();
105 }
106
107 if (! m_scaleKshort) {
110 m_outputAntiparticleList->bindAntiParticleList(*(m_outputparticleList));
111 }
112
113 if (m_scaleKshort) { // scale V0 Kshort
114
115 // loop over Kshort
116 const unsigned int nPar = m_inputparticleList->getListSize();
117 for (unsigned i = 0; i < nPar; i++) {
118 const Particle* particle = m_inputparticleList->getParticle(i);
119 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_V0) {
120 B2WARNING(" Input ParticleList " << m_inputListName <<
121 " contains a particle which is not from V0. It will not be copied to output list.");
122 continue;
123 }
124
125 if (particle->getNDaughters() != 2)
126 B2ERROR("V0 particle should have exactly two daughters");
127
128 Particle* newV0 = ParticleCopy::copyParticle(particle);
129
130 const Particle* dauP = newV0->getDaughter(0);
131 const Particle* dauM = newV0->getDaughter(1);
132
133 Particle* newDauP = getChargedWithScaledError(dauP);
134 Particle* newDauM = getChargedWithScaledError(dauM);
135
136 ROOT::Math::PxPyPzEVector v0Momentum = newDauP->get4Vector() + newDauM->get4Vector();
137 newV0->set4VectorDividingByMomentumScaling(v0Momentum);
138
139 newV0->replaceDaughter(dauP, newDauP);
140 newV0->replaceDaughter(dauM, newDauM);
141
142 m_outputparticleList->addParticle(newV0);
143
144 } // loop over Kshort
145
146 } else { // scale charged particles
147
148 // loop over charged particles
149 const unsigned int nPar = m_inputparticleList->getListSize();
150 for (unsigned i = 0; i < nPar; i++) {
151 const Particle* charged = m_inputparticleList->getParticle(i);
152 Particle* newCharged = getChargedWithScaledError(charged);
153 m_outputparticleList->addParticle(newCharged);
154 }
155 }
156
157}
158
160{
161 const TrackFitResult* new_trkfit = getTrackFitResultWithScaledError(particle);
162
163 Const::ChargedStable chargedtype(abs(particle->getPDGCode()));
164 Particle new_charged(particle->getMdstArrayIndex(), new_trkfit, chargedtype);
165 Particle* newCharged = m_particles.appendNew(new_charged);
166
167 // add relation to PID, MCParticle (if any) and TrackFitResult
168 const PIDLikelihood* pid = particle->getPIDLikelihood();
169 const MCParticle* mcCharged = particle->getRelated<MCParticle>();
170 newCharged->addRelationTo(new_trkfit);
171 if (pid) newCharged->addRelationTo(pid);
172 if (mcCharged != nullptr) newCharged->addRelationTo(mcCharged);
173
174 return newCharged;
175}
176
178{
179 const TrackFitResult* trkfit = particle->getTrackFitResult();
180
181 const std::vector<float> helix = trkfit->getTau();
182 const std::vector<float> cov5 = trkfit->getCov();
183 const Const::ParticleType ptype = trkfit->getParticleType();
184 const double pvalue = trkfit->getPValue();
185 const ULong64_t hitCDC = trkfit->getHitPatternCDC().getInteger();
186 const ULong64_t hitVXD = trkfit->getHitPatternVXD().getInteger();
187 const int ndf = trkfit->getNDF();
188
189 std::vector<double> scaleFactors = getScaleFactors(particle, trkfit);
190 std::vector<float> cov5_scaled;
191 unsigned int counter = 0;
192 for (unsigned int j = 0; j < 5; j++) {
193 for (unsigned int k = j; k < 5; k++) {
194 cov5_scaled.push_back(cov5[counter++] * scaleFactors[j] * scaleFactors[k]);
195 }
196 }
197 TrackFitResult* new_trkfit = m_trackfitresults.appendNew(helix, cov5_scaled, ptype, pvalue, hitCDC, hitVXD, ndf);
198 return new_trkfit;
199}
200
201std::vector<double> HelixErrorScalerModule::getScaleFactors(const Particle* particle, const TrackFitResult* trkfit)
202{
203 if (trkfit->getHitPatternVXD().getNPXDHits() > 0) {
204
205 // d0, z0 resolution = a (+) b / pseudo-momentum
206 B2Vector3D p = particle->getMomentum();
207 double sinTheta = TMath::Sin(p.Theta());
208 double pD0 = p.Mag2() / (particle->getEnergy()) * TMath::Power(sinTheta, 1.5); // p*beta*sinTheta**1.5
209 double pZ0 = pD0 * sinTheta; // p*beta*sinTheta**2.5
210
211 pD0 = TMath::Max(pD0, m_d0MomThr); // if pD0 is smaller than the threshold, *overwrite* it with the threshold.
212 pZ0 = TMath::Max(pZ0, m_z0MomThr); // if pZ0 is smaller than the threshold, *overwrite* it with the threshold.
213 double d0Resol = TMath::Sqrt(TMath::Power(m_d0ResolPars[0], 2) + TMath::Power(m_d0ResolPars[1] / pD0, 2));
214 double z0Resol = TMath::Sqrt(TMath::Power(m_z0ResolPars[0], 2) + TMath::Power(m_z0ResolPars[1] / pZ0, 2));
215 double d0Err = TMath::Sqrt(trkfit->getCovariance5()[0][0]);
216 double z0Err = TMath::Sqrt(trkfit->getCovariance5()[3][3]);
217
218 std::vector<double> scaleFactors = { TMath::Max(d0Resol / d0Err, m_scaleFactors_PXD[0]),
221 TMath::Max(z0Resol / z0Err, m_scaleFactors_PXD[3]),
223 };
224 return scaleFactors;
225 } else {
227 }
228}
229
230
231
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:589
The ParticleType class for identifying different particle types.
Definition: Const.h:408
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
static const ParticleType invalidParticle
Invalid particle, used internally.
Definition: Const.h:681
static const ParticleType Kshort
K^0_S particle.
Definition: Const.h:677
Represents a particle in the DecayDescriptor.
int getPDGCode() const
Return PDG code.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
bool m_scaleKshort
Whether the input particle list is Kshort or not.
StoreArray< TrackFitResult > m_trackfitresults
StoreArray of TrackFitResult objects.
std::vector< double > m_d0ResolPars
parameters (a,b) to define d0 best resolution = a (+) b / (p*beta*sinTheta**1.5)
StoreObjPtr< ParticleList > m_outputparticleList
StoreObjptr for output particlelist.
virtual void initialize() override
Register input and output data.
virtual void event() override
loop over the input charged particles
StoreArray< MCParticle > m_mcparticles
StoreArray of MCParticle objects.
double m_d0MomThr
d0 best resolution is kept constant below this momentum.
std::vector< double > m_scaleFactors_PXD
vector of five scale factors for helix parameter errors (for tracks with a PXD hit)
StoreArray< Particle > m_particles
StoreArray of Particle objects.
std::vector< double > getScaleFactors(const Particle *particle, const TrackFitResult *trkfit)
get scale factors
StoreArray< PIDLikelihood > m_pidlikelihoods
StoreArray of PIDLikelihood objects.
HelixErrorScalerModule()
Constructor: Sets the description, the properties and the parameters of the module.
const TrackFitResult * getTrackFitResultWithScaledError(const Particle *particle)
create a TrackFitResult with scaled errors
std::vector< double > m_z0ResolPars
parameters (a,b) to define z0 best resolution = a (+) b / (p*beta*sinTheta**2.5)
StoreObjPtr< ParticleList > m_inputparticleList
StoreObjptr for input charged particlelist.
Particle * getChargedWithScaledError(const Particle *particle)
create a Particle with scaled errors
StoreObjPtr< ParticleList > m_outputAntiparticleList
StoreObjptr for output antiparticlelist.
DecayDescriptor m_decaydescriptor
Decay descriptor of the charged particle.
double m_z0MomThr
z0 best resolution is kept constant below this momentum.
int m_pdgCode
PDG code of the charged particle to be scaled.
std::string m_outputAntiListName
output anti-particle list name
std::vector< double > m_scaleFactors_noPXD
vector of five scale factors for helix parameter errors (for tracks without a PXD hit)
std::string m_inputListName
The name of input charged particle list.
std::string m_outputListName
The name of output charged particle list.
ULong64_t getInteger() const
Getter for underlying integer type.
unsigned int getInteger() const
Getter for the underlying integer.
Definition: HitPatternVXD.h:58
unsigned short getNPXDHits() const
Get total number of hits in the PXD.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:29
Class to store reconstructed particles.
Definition: Particle.h:75
bool replaceDaughter(const Particle *oldDaughter, Particle *newDaughter)
Replace index of given daughter with new daughter, return true if a replacement is made.
Definition: Particle.cc:705
void set4VectorDividingByMomentumScaling(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector dividing by the momentum scaling factor.
Definition: Particle.h:283
double getEnergy() const
Returns total energy.
Definition: Particle.h:535
unsigned getMdstArrayIndex(void) const
Returns 0-based index of MDST store array object (0 for composite particles)
Definition: Particle.h:487
const PIDLikelihood * getPIDLikelihood() const
Returns the pointer to the PIDLikelihood object that is related to the Track, which was used to creat...
Definition: Particle.cc:871
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition: Particle.h:727
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
EParticleSourceObject getParticleSource() const
Returns particle source as defined with enum EParticleSourceObject.
Definition: Particle.h:478
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition: Particle.cc:854
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
Values of the result of a track fit with a given particle hypothesis.
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
float getNDF() const
Getter for number of degrees of freedom of the track fit.
double getPValue() const
Getter for Chi2 Probability of the track fit.
std::vector< float > getTau() const
Getter for all perigee parameters.
std::vector< float > getCov() const
Getter for all covariance matrix elements of perigee parameters.
Const::ParticleType getParticleType() const
Getter for ParticleType of the mass hypothesis of the track fit.
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
HitPatternVXD getHitPatternVXD() const
Getter for the hit pattern in the VXD;.
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
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
std::string antiParticleListName(const std::string &listName)
Returns name of anti-particle-list corresponding to listName.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24