12 #include <analysis/modules/BremsCorrection/BremsFinderModule.h>
15 #include <framework/gearbox/Const.h>
16 #include <framework/logging/Logger.h>
17 #include <framework/datastore/RelationArray.h>
18 #include <framework/datastore/RelationVector.h>
19 #include <framework/datastore/StoreArray.h>
22 #include <analysis/dataobjects/Particle.h>
23 #include <mdst/dataobjects/MCParticle.h>
24 #include <mdst/dataobjects/PIDLikelihood.h>
25 #include <mdst/dataobjects/Track.h>
26 #include <mdst/dataobjects/ECLCluster.h>
29 #include <analysis/DecayDescriptor/ParticleListName.h>
30 #include <analysis/DecayDescriptor/DecayDescriptor.h>
33 #include <analysis/variables/ECLVariables.h>
37 #include <TMatrixFSym.h>
61 This module copies each particle in the ``inputList`` to the ``outputList`` and uses
62 the results of the **eclTrackBremFinder** module to look for possible bremsstrahlung photons; if these
63 photons exists, it adds their four momentum to the particle in the ``outputList``.
64 It also adds the original particle and these photons as daughters of the new, corrected particle.
65 Track and PID information of the original particle are copied onto the new one to facilitate their access
66 in the analysis scripts.
68 The **eclTrackBremFinder** module uses the lepton track PXD and SVD hits and extrapolates them to the ECL;
69 then looks for ECL clusters with energies between 0.2 and 1 times the track energy and without associated
70 tracks, and checks if the distance between each of these clusters
71 and the extrapolated hit is smaller than 0.5 mm. If it is, a *Bremsstrahlung* weighted relation
72 between said cluster and the track is established. The weight is determined as
76 \text{max}\left(\frac{\left|\phi_{\text{cluster}}-\phi_{\text{hit}}\right|}{\Delta\phi_{\text{cluster}}+\Delta\phi_{\text{hit}}}, \, \frac{\left|\theta_{\text{cluster}}-\theta_{\text{hit}}\right|}{\Delta\theta_{\text{cluster}}+\Delta\theta_{\text{hit}}}\right)
78 where :math:`\phi_i` and :math:`\theta_i` are the azimutal and polar angles of the ECL cluster and the
79 extrapolated hit, and :math:`\Delta x` represents the uncertainty of the value :math:`x`. The details of
80 the calculation of these quantities are `here`_. By default, only relations with a weight
81 smaller than 3.0 are stored. The user can further determine the maximum value of this weight required in order
82 to perform the bremsstrahlung correction.
84 This module looks for photons in the ``gammaList`` whose clusters have a *Bremsstrahlung* relation with the track
85 of one of the particles in the ``inputList``, and adds their 4-momentum to the particle's one. It also stores the value
86 of each relation weight as ``extraInfo`` of the corrected particle, under the name ``"bremsWeightWithPhotonN"``, where
87 N is the index of the photon as daughter of the corrected particle; thus ``"bremsWeightWithPhoton0"`` gives the weight
88 of the Bremsstrahlung relation between the new, corrected particle, and the first photon daughter.
91 Even in the event of no bremsstrahlung photons found, a new particle is still created, and the original one is still
92 added as its daughter.
95 `eclTrackBremFinder module`_
97 .. _eclTrackBremFinder module: https://stash.desy.de/projects/B2/repos/software/browse/ecl/modules/eclTrackBremFinder
98 .. _here: https://stash.desy.de/projects/B2/repos/software/browse/ecl/modules/eclTrackBremFinder/src/BremFindingMatchCompute.cc)DOC");
99 setPropertyFlags(c_ParallelProcessingCertified);
102 addParam(
"outputList", m_outputListName,
"The output particle list name.");
103 addParam(
"inputList", m_inputListName,
104 R
"DOC(The initial particle list name containing the particles to correct. *It should already exist* and *the particles must have an associated track.*)DOC");
105 addParam("gammaList", m_gammaListName,
106 R
"DOC(The photon list containing the preselected bremsstrahlung candidates. *It should already exist* and *the particles in the list must be photons*)DOC");
107 addParam("maximumAcceptance", m_maximumAcceptance,
108 "The maximum value of the relation weight between a bremsstrahlung cluster and a particle track",
109 m_maximumAcceptance);
110 addParam(
"multiplePhotons", m_addMultiplePhotons,
"If true, use all possible photons to correct the particle's 4-momentum",
111 m_addMultiplePhotons);
112 addParam(
"usePhotonOnlyOnce", m_usePhotonOnlyOnce,
113 R
"DOC(If true, each brems candidate is used to correct maximum 1 particle (the one with the lowest relation weight among all in the ``inputList``).)DOC",
114 m_usePhotonOnlyOnce);
115 addParam("writeOut", m_writeOut,
116 R
"DOC(If true, the output ``ParticleList`` will be saved by `RootOutput`. If false, it will be ignored when writing the file.)DOC",
120 void BremsFinderModule::initialize()
123 m_inputList.isRequired(m_inputListName);
126 decDes.
init(m_inputListName);
130 if (Const::chargedStableSet.find(abs(m_pdgCode)) == Const::invalidParticle)
131 B2ERROR(
"[BremsFinderModule] Invalid particle list. the particles in " << m_inputListName <<
" should have an associated track.");
134 m_gammaList.isRequired(m_gammaListName);
136 decDes.
init(m_gammaListName);
140 if (temp_pdg != Const::photon.getPDGCode())
141 B2ERROR(
"[BremsFinderModule] Invalid particle list. the particles in " << m_gammaListName <<
" should be photons!");
143 decDes.
init(m_outputListName);
146 if (m_inputListName == m_outputListName)
147 B2ERROR(
"[BremsFinderModule] Input and output particle list names are the same: " << m_inputListName);
148 else if (temp_pdg != m_pdgCode) {
149 B2ERROR(
"[BremsFinderModule] The input and output particle list correspond to different particles: " << m_inputListName <<
" --> "
150 << m_outputListName);
154 m_outputAntiListName = ParticleListName::antiParticleListName(m_outputListName);
158 m_outputList.registerInDataStore(m_outputListName, flags);
159 m_outputAntiList.registerInDataStore(m_outputAntiListName, flags);
163 particles.registerRelationTo(pidlikelihoods);
166 void BremsFinderModule::event()
171 RelationArray particlesToMCParticles(particles, mcParticles);
174 m_outputList.create();
175 m_outputList->initialize(m_pdgCode, m_outputListName);
177 m_outputAntiList.create();
178 m_outputAntiList->initialize(-1 * m_pdgCode, m_outputAntiListName);
179 m_outputAntiList->bindAntiParticleList(*(m_outputList));
182 const unsigned int nGamma = m_gammaList->getListSize();
185 const unsigned int nLep = m_inputList->getListSize();
187 const std::string relationName =
"Bremsstrahlung";
190 if (m_usePhotonOnlyOnce) {
191 for (
unsigned n = 0; n < nGamma; n++) {
192 Particle* gamma = m_gammaList->getParticle(n);
195 if (gamma->hasExtraInfo(
"bestMatchIndex"))
continue;
197 auto cluster = gamma->getECLCluster();
200 double bestWeight = m_maximumAcceptance;
201 unsigned bestMatchIndex = 0;
202 unsigned trkIndex = 0;
204 for (
auto trk = relatedTracks.
begin(); trk != relatedTracks.
end(); trk++, trkIndex++) {
206 for (
unsigned i = 0; i < nLep; i++) {
207 const Particle* lepton = m_inputList->getParticle(i);
208 auto leptonTrack = lepton->
getTrack();
210 if (leptonTrack->getArrayIndex() == trk->getArrayIndex()) {
211 double weight = relatedTracks.
weight(trkIndex);
212 if (weight < bestWeight) {
215 bestMatchIndex = trk->getArrayIndex();
222 gamma->addExtraInfo(
"bestMatchIndex", bestMatchIndex);
228 for (
unsigned i = 0; i < nLep; i++) {
229 const Particle* lepton = m_inputList->getParticle(i);
237 std::vector<std::pair <double, Particle*> > selectedGammas;
240 for (
auto bremCluster = bremClusters.
begin(); bremCluster != bremClusters.
end(); bremCluster++, j++) {
241 double weight = bremClusters.
weight(j);
243 if (weight > m_maximumAcceptance)
continue;
245 for (
unsigned k = 0; k < nGamma; k++) {
247 Particle* gamma = m_gammaList->getParticle(k);
248 auto cluster = gamma->getECLCluster();
250 if (bremCluster->getClusterId() == cluster->getClusterId()) {
251 if (m_usePhotonOnlyOnce) {
252 if (track->getArrayIndex() ==
static_cast<int>
253 (gamma->getExtraInfo(
"bestMatchIndex")))
254 selectedGammas.push_back(std::make_pair(weight, gamma));
256 selectedGammas.push_back(std::make_pair(weight, gamma));
265 TLorentzVector new4Vec = lepton->
get4Vector();
268 std::sort(selectedGammas.begin(), selectedGammas.end());
271 for (
auto const& bremsPair : selectedGammas) {
273 new4Vec += g->get4Vector();
274 if (! m_addMultiplePhotons)
break;
278 Particle correctedLepton(new4Vec, lepton->
getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
279 track->getArrayIndex());
285 TMatrixFSym corLepMatrix(lepErrorMatrix);
287 double bremsGammaEnergySum = 0.0;
290 for (
auto const& bremsPair : selectedGammas) {
292 Particle* bremsGamma = bremsPair.second;
293 std::string extraInfoName =
"bremsWeightWithPhoton" + std::to_string(photonIndex);
294 correctedLepton.
addExtraInfo(extraInfoName, bremsPair.first);
296 bremsGammaEnergySum += Variable::eclClusterE(bremsGamma);
299 for (
int irow = 0; irow <= 3; irow++) {
300 for (
int icol = irow; icol <= 3; icol++) corLepMatrix(irow, icol) += gammaErrorMatrix(irow, icol);
303 B2DEBUG(10,
"[BremsFinderModule] Found a bremsstrahlung gamma and added its 4-vector to the charged particle");
304 if (! m_addMultiplePhotons)
break;
312 correctedLepton.
addExtraInfo(
"bremsCorrected",
float(selectedGammas.size() > 0));
313 correctedLepton.
addExtraInfo(
"bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
316 Particle* newLepton = particles.appendNew(correctedLepton);
324 m_outputList->addParticle(newLepton);