Belle II Software development
ContinuumSuppressionBuilderModule.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/ContinuumSuppressionBuilder/ContinuumSuppressionBuilderModule.h>
10
11#include <analysis/ContinuumSuppression/Thrust.h>
12#include <analysis/ContinuumSuppression/KsfwMoments.h>
13#include <analysis/ContinuumSuppression/FoxWolfram.h>
14#include <analysis/ContinuumSuppression/CleoCones.h>
15
16#include <analysis/dataobjects/Particle.h>
17#include <analysis/dataobjects/RestOfEvent.h>
18
19#include <analysis/utility/CLHEPToROOT.h>
20#include <analysis/utility/PCmsLabTransform.h>
21#include <analysis/utility/ROOTToCLHEP.h>
22
23#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
24#include <analysis/VertexFitting/KFit/VertexFitKFit.h>
25
26#include <framework/dataobjects/Helix.h>
27
28#include <framework/geometry/BFieldManager.h>
29
30#include <mdst/dataobjects/TrackFitResult.h>
31
32#include <vector>
33#include <Math/Vector3D.h>
34
35using namespace Belle2;
36
37//-----------------------------------------------------------------
38// Register the Module
39//-----------------------------------------------------------------
40REG_MODULE(ContinuumSuppressionBuilder);
41
42//-----------------------------------------------------------------
43// Implementation
44//-----------------------------------------------------------------
45
47{
48 // Set module properties
49 setDescription("Creates for each Particle in the given ParticleLists a ContinuumSuppression dataobject and makes basf2 relation between them.");
51
52 // Parameter definitions
53 addParam("particleList", m_particleListName, "Name of the ParticleList", std::string(""));
54 addParam("ROEMask", m_ROEMask, "ROE mask", std::string(RestOfEvent::c_defaultMaskName));
55 addParam("performIPProfileFit", m_ipProfileFit, "switch to turn on vertex fit of tracks with IP profile constraint", false);
56}
57
59{
60 // Input
61 m_plist.isRequired(m_particleListName);
63
64 if (m_ROEMask.empty()) {
66 }
67
68 if (m_ROEMask == "FS1" or m_ROEMask == "ROE") {
69 B2ERROR("The ROE mask for the continuum suppression must not be called " << m_ROEMask);
70 }
71
72 // Output
73 m_csarray.registerInDataStore(m_ROEMask);
75
76 // set magnetic field
77 m_Bfield = BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
78}
79
81{
82 m_BeamSpotCenter = m_beamSpotDB->getIPPosition();
83 m_beamSpotCov.ResizeTo(3, 3);
84 m_beamSpotCov = m_beamSpotDB->getCovVertex();
85
86 for (unsigned i = 0; i < m_plist->getListSize(); i++) {
87 addContinuumSuppression(m_plist->getParticle(i));
88 }
89}
90
92{
93 // Create ContinuumSuppression object
94 ContinuumSuppression* qqVars = m_csarray.appendNew();
95
96 // Create relation: Particle <-> ContinuumSuppression
97 particle->addRelationTo(qqVars);
98
99 std::vector<ROOT::Math::PxPyPzEVector> p_cms_sigB, p_cms_roe, p_cms_all;
100
101 std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>> p_cms_q_sigA;
102 std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>> p_cms_q_sigB;
103 std::vector<std::pair<ROOT::Math::PxPyPzEVector, int>> p_cms_q_roe;
104
105 std::vector<float> ksfwFS0;
106 std::vector<float> ksfwFS1;
107
108 std::vector<float> cleoConesAll;
109 std::vector<float> cleoConesROE;
110
111 double et[2];
112
113 ROOT::Math::XYZVector thrustB;
114 ROOT::Math::XYZVector thrustO;
115
116 float thrustBm = -1;
117 float thrustOm = -1;
118 float cosTBTO = -1;
119 float cosTBz = -1;
120 float R2 = -1;
121
122
123 // -- B Cand --------------------------------------------------------------------------
125 double BeamEnergy = T.getCMSEnergy() / 2;
126
127 ROOT::Math::PxPyPzEVector p_cms_missA(0, 0, 0, 2 * BeamEnergy);
128 ROOT::Math::PxPyPzEVector p_cms_missB(0, 0, 0, 2 * BeamEnergy);
129 et[0] = et[1] = 0;
130
131 // -- SIG A --- Use B primary daughters - (Belle: use_finalstate_for_sig == 0) --------
132 std::vector<Belle2::Particle*> signalDaughters = particle->getDaughters();
133
134 for (const Belle2::Particle* sigFS0 : signalDaughters) {
135 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS0->get4Vector();
136
137 p_cms_q_sigA.emplace_back(p_cms, sigFS0->getCharge());
138
139 p_cms_missA -= p_cms;
140 et[0] += p_cms.Pt();
141 }
142
143 // -- SIG B --- Use B final-state daughters - (Belle: use_finalstate_for_sig == 1) ----
144 std::vector<const Belle2::Particle*> signalFSParticles = particle->getFinalStateDaughters();
145
146 for (const Belle2::Particle* sigFS1 : signalFSParticles) {
147 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS1->get4Vector();
148
149 p_cms_all.push_back(p_cms);
150 p_cms_sigB.push_back(p_cms);
151
152 p_cms_q_sigB.emplace_back(p_cms, sigFS1->getCharge());
153
154 p_cms_missB -= p_cms;
155 et[1] += p_cms.Pt();
156 }
157
158 // -- ROE -----------------------------------------------------------------------------
159 const RestOfEvent* roe = particle->getRelated<RestOfEvent>();
160
161 if (roe) {
162
163 // Charged tracks
164 //
165 std::vector<const Particle*> chargedROEParticles = roe->getChargedParticles(m_ROEMask);
166
167 for (const Particle* chargedROEParticle : chargedROEParticles) {
168
169 ROOT::Math::PxPyPzEVector p = ipProfileFit(chargedROEParticle);
170
171 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * p;
172
173 p_cms_all.push_back(p_cms);
174 p_cms_roe.push_back(p_cms);
175
176 p_cms_q_roe.emplace_back(p_cms, chargedROEParticle->getCharge());
177
178 p_cms_missA -= p_cms;
179 p_cms_missB -= p_cms;
180 et[0] += p_cms.Pt();
181 et[1] += p_cms.Pt();
182 }
183
184 // ECLCluster
185 //
186 std::vector<const Particle*> roePhotons = roe->getPhotons(m_ROEMask);
187
188 for (const Particle* photon : roePhotons) {
189
190 if (photon->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
191
192 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * photon->get4Vector();
193 p_cms_all.push_back(p_cms);
194 p_cms_roe.push_back(p_cms);
195
196 p_cms_q_roe.emplace_back(p_cms, photon->getCharge());
197
198 p_cms_missA -= p_cms;
199 p_cms_missB -= p_cms;
200 et[0] += p_cms.Pt();
201 et[1] += p_cms.Pt();
202 }
203 }
204
205 // Thrust variables
206 thrustB = Thrust::calculateThrust(p_cms_sigB);
207 thrustO = Thrust::calculateThrust(p_cms_roe);
208 thrustBm = thrustB.R();
209 thrustOm = thrustO.R();
210 cosTBTO = fabs(thrustB.Unit().Dot(thrustO.Unit()));
211 cosTBz = fabs(cos(thrustB.Theta()));
212
213 // Cleo Cones
214 CleoCones cc(p_cms_all, p_cms_roe, thrustB, true, true);
215 cleoConesAll = cc.cleo_cone_with_all();
216 cleoConesROE = cc.cleo_cone_with_roe();
217
218 // Fox-Wolfram Moments: Uses all final-state tracks (= sigB + ROE)
219 FoxWolfram FW(p_cms_all);
221 R2 = FW.getR(2);
222
223 // KSFW moments
224 ROOT::Math::PxPyPzEVector p_cms_B = T.rotateLabToCms() * particle->get4Vector();
225 double Hso0_max(2 * (2 * BeamEnergy - p_cms_B.E()));
226 KsfwMoments KsfwM(Hso0_max,
227 p_cms_q_sigA,
228 p_cms_q_sigB,
229 p_cms_q_roe,
230 p_cms_missA,
231 p_cms_missB,
232 et);
233 // use_finalstate_for_sig == 0
234 KsfwM.usefinal(0);
235 ksfwFS0.push_back(KsfwM.mm2());
236 ksfwFS0.push_back(KsfwM.et());
237 ksfwFS0.push_back(KsfwM.Hso(0, 0));
238 ksfwFS0.push_back(KsfwM.Hso(0, 1));
239 ksfwFS0.push_back(KsfwM.Hso(0, 2));
240 ksfwFS0.push_back(KsfwM.Hso(0, 3));
241 ksfwFS0.push_back(KsfwM.Hso(0, 4));
242 ksfwFS0.push_back(KsfwM.Hso(1, 0));
243 ksfwFS0.push_back(KsfwM.Hso(1, 2));
244 ksfwFS0.push_back(KsfwM.Hso(1, 4));
245 ksfwFS0.push_back(KsfwM.Hso(2, 0));
246 ksfwFS0.push_back(KsfwM.Hso(2, 2));
247 ksfwFS0.push_back(KsfwM.Hso(2, 4));
248 ksfwFS0.push_back(KsfwM.Hoo(0));
249 ksfwFS0.push_back(KsfwM.Hoo(1));
250 ksfwFS0.push_back(KsfwM.Hoo(2));
251 ksfwFS0.push_back(KsfwM.Hoo(3));
252 ksfwFS0.push_back(KsfwM.Hoo(4));
253 // use_finalstate_for_sig == 1
254 KsfwM.usefinal(1);
255 ksfwFS1.push_back(KsfwM.mm2());
256 ksfwFS1.push_back(KsfwM.et());
257 ksfwFS1.push_back(KsfwM.Hso(0, 0));
258 ksfwFS1.push_back(KsfwM.Hso(0, 1));
259 ksfwFS1.push_back(KsfwM.Hso(0, 2));
260 ksfwFS1.push_back(KsfwM.Hso(0, 3));
261 ksfwFS1.push_back(KsfwM.Hso(0, 4));
262 ksfwFS1.push_back(KsfwM.Hso(1, 0));
263 ksfwFS1.push_back(KsfwM.Hso(1, 2));
264 ksfwFS1.push_back(KsfwM.Hso(1, 4));
265 ksfwFS1.push_back(KsfwM.Hso(2, 0));
266 ksfwFS1.push_back(KsfwM.Hso(2, 2));
267 ksfwFS1.push_back(KsfwM.Hso(2, 4));
268 ksfwFS1.push_back(KsfwM.Hoo(0));
269 ksfwFS1.push_back(KsfwM.Hoo(1));
270 ksfwFS1.push_back(KsfwM.Hoo(2));
271 ksfwFS1.push_back(KsfwM.Hoo(3));
272 ksfwFS1.push_back(KsfwM.Hoo(4));
273
274 // TODO: The following is from the original belle ksfwmoments.cc module.
275 // Not sure if necessary here (i.e., will we be using rooksfw in belle II in the same way?).
276 // printf("rooksfw::rooksfw: mm2=%f et=%f hoo2=%f hso02=%f\n",
277 // m_mm2[0], et[0], m_Hoo[0][2], m_Hso[0][0][2]);
278 }
279
280 // Fill ContinuumSuppression with content
281 qqVars->addThrustB(thrustB);
282 qqVars->addThrustO(thrustO);
283 qqVars->addThrustBm(thrustBm);
284 qqVars->addThrustOm(thrustOm);
285 qqVars->addCosTBTO(cosTBTO);
286 qqVars->addCosTBz(cosTBz);
287 qqVars->addR2(R2);
288 qqVars->addKsfwFS0(ksfwFS0);
289 qqVars->addKsfwFS1(ksfwFS1);
290 qqVars->addCleoConesALL(cleoConesAll);
291 qqVars->addCleoConesROE(cleoConesROE);
292}
293
294ROOT::Math::PxPyPzEVector ContinuumSuppressionBuilderModule::ipProfileFit(const Particle* particle)
295{
296 // Return original particle if vertex fit is not requested
297 if (!m_ipProfileFit) return particle->get4Vector();
298
299 // Access track fit result of particle
300 const TrackFitResult* trackFitResult = particle->getTrackFitResult();
301 if (!trackFitResult) return particle->get4Vector();
302
303 // Move helix of TrackFitResult by IP position
304 Helix helix = trackFitResult->getHelix();
305 helix.passiveMoveBy(m_BeamSpotCenter);
306
307 // Reject particle movements that push helix too far away
308 if (std::abs(helix.getD0()) > 5.0 || std::abs(helix.getZ0()) > 10.0) return particle->get4Vector();
309
310 // Create particle with shifted vertex position
311 ROOT::Math::XYZVector p_helix = helix.getMomentum(m_Bfield);
312 double shiftedParticleMass = particle->getPDGMass();
313 double shiftedParticleEnergy = std::sqrt(p_helix.Mag2() + shiftedParticleMass * shiftedParticleMass);
314 ROOT::Math::PxPyPzEVector momentumOfShiftedParticle(p_helix.X(), p_helix.y(), p_helix.Z(), shiftedParticleEnergy);
315 const Particle shiftedParticle(momentumOfShiftedParticle, particle->getPDGCode());
316
317 // Perform vertex fit with IP profile constraint
320 kv.addParticle(&shiftedParticle);
321 kv.setIpProfile(ROOTToCLHEP::getPoint3D(m_BeamSpotCenter), ROOTToCLHEP::getHepSymMatrix(m_beamSpotCov));
322 enum analysis::KFitError::ECode fitError = kv.doFit();
323 if (fitError != analysis::KFitError::kNoError) return particle->get4Vector();
324
325 // Calculate updated particle momentum
330 kmm.setVertex(kv.getVertex());
332 fitError = kmm.doMake();
333 if (fitError != analysis::KFitError::kNoError) return particle->get4Vector();
334 return CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum());
335}
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Helix parameter class.
Definition Helix.h:48
Class to calculate the Cleo clone variables.
Definition CleoCones.h:23
std::vector< float > cleo_cone_with_roe()
Returns calculated Cleo Cones constructed from only ROE tracks.
Definition CleoCones.h:46
std::vector< float > cleo_cone_with_all()
Returns calculated Cleo Cones constructed from all tracks.
Definition CleoCones.h:41
TMatrixDSym m_beamSpotCov
Beam spot covariance matrix.
ROOT::Math::PxPyPzEVector ipProfileFit(const Particle *particle)
get 4Vector for CS calculation with or without IP profile fit
std::string m_particleListName
Name of the ParticleList.
virtual void initialize() override
initialize the module (setup the data store)
bool m_ipProfileFit
Switch to turn on / off vertex fit of tracks with IP profile constraint.
ROOT::Math::XYZVector m_BeamSpotCenter
Beam spot position.
DBObjPtr< BeamSpot > m_beamSpotDB
Beam spot database object.
void addContinuumSuppression(const Particle *particle)
calculate continuum suppression quantities
StoreArray< ContinuumSuppression > m_csarray
StoreArray of ContinuumSuppression.
StoreObjPtr< ParticleList > m_plist
input particle list
This is a class for collecting variables used in continuum suppression.
void addR2(float R2)
Add reduced Fox-Wolfram moment R2.
void addThrustBm(float thrustBm)
Add magnitude of B thrust axis.
void addCleoConesALL(const std::vector< float > &cleoConesALL)
Add vector of Cleo Cones constructed of all final state particles.
void addCleoConesROE(const std::vector< float > &cleoConesROE)
Add vector of Cleo Cones constructed of only ROE particles.
void addThrustB(const ROOT::Math::XYZVector &thrustB)
Add ROE thrust axis.
void addCosTBz(float cosTBz)
Add cosine of the angle between the thrust axis of the B and the z-axis.
void addThrustOm(float thrustOm)
Add magnitude of ROE thrust axis.
void addKsfwFS0(const std::vector< float > &ksfwFS0)
Add vector of KSFW moments, Et, and mm2 for final state = 0.
void addThrustO(const ROOT::Math::XYZVector &thrustO)
Add ROE thrust axis.
void addKsfwFS1(const std::vector< float > &ksfwFS1)
Add vector of KSFW moments, Et, and mm2 for final state = 1.
void addCosTBTO(float cosTBTO)
Add cosine of the angle between the thrust axis of the B and the thrust axis of the ROE.
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
Class to calculate the Fox-Wolfram moments up to order 8.
Definition FoxWolfram.h:28
double getR(int i) const
Returns the i-th moment normalized to the 0th-order moment.
Definition FoxWolfram.h:89
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition FoxWolfram.cc:14
Moment-calculation of the k_sfw improved Super-Fox-Wolfram moments.
Definition KsfwMoments.h:34
int usefinal(int uf)
Sets the flag that specifiies we are using the finalstate for signal.
Definition KsfwMoments.h:77
double et(int uf=-1) const
Returns calculated transverse energy.
Definition KsfwMoments.h:92
double Hso(int i, int j, int uf=-1) const
Returns calculated KSFW Moments.
double mm2(int uf=-1) const
Returns calculated missing mass squared.
Definition KsfwMoments.h:87
double Hoo(int i, int uf=-1) const
Returns calculated KSFW Moments.
Definition KsfwMoments.h:97
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
Module()
Constructor.
Definition Module.cc:30
@ 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.
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
Class to store reconstructed particles.
Definition Particle.h:76
This is a general purpose class for collecting reconstructed MDST data objects that are not used in r...
Definition RestOfEvent.h:55
std::vector< const Particle * > getChargedParticles(const std::string &maskName=c_defaultMaskName, unsigned int pdg=0, bool unpackComposite=true) const
Get charged particles from ROE mask.
static constexpr const char * c_defaultMaskName
Default mask name.
Definition RestOfEvent.h:58
std::vector< const Particle * > getPhotons(const std::string &maskName=c_defaultMaskName, bool unpackComposite=true) const
Get photons from ROE mask.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
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
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::PxPyPzEVector > &momenta)
calculates the thrust axis
Definition Thrust.cc:71
Values of the result of a track fit with a given particle hypothesis.
Helix getHelix() const
Conversion to framework Helix (without covariance).
const CLHEP::HepSymMatrix getTrackError(const int id) const
Get an error matrix of the track.
Definition KFitBase.cc:168
const CLHEP::HepLorentzVector getTrackMomentum(const int id) const
Get a Lorentz vector of the track.
Definition KFitBase.cc:154
const HepPoint3D getTrackPosition(const int id) const
Get a position of the track.
Definition KFitBase.cc:161
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
Definition KFitBase.cc:93
const KFitTrack getTrack(const int id) const
Get a specified track object.
Definition KFitBase.cc:175
enum KFitError::ECode addParticle(const Particle *particle)
Add a particle to the fitter.
Definition KFitBase.cc:59
ECode
ECode is a error code enumerate.
Definition KFitError.h:33
double getCharge(void) const
Get a charge of the track.
Definition KFitTrack.cc:179
MakeMotherKFit is a class to build mother particle from kinematically fitted daughters.
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set a vertex position of the mother particle.
enum KFitError::ECode addTrack(const KFitTrack &kp)
Add a track to the make-mother object.
enum KFitError::ECode doMake(void)
Perform a reconstruction of mother particle.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set a vertex error matrix of the mother particle.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
const CLHEP::HepLorentzVector getMotherMomentum(void) const
Get a Lorentz vector of the mother particle.
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
const CLHEP::HepSymMatrix getVertexError(void) const
Get a fitted vertex error matrix.
enum KFitError::ECode doFit(void)
Perform a vertex-constraint fit.
enum KFitError::ECode setIpProfile(const HepPoint3D &ip, const CLHEP::HepSymMatrix &ipe)
Set an IP-ellipsoid shape for the vertex constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
const CLHEP::HepMatrix getTrackVertexError(const int id) const
Get a vertex error matrix of the track.
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:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.