Belle II Software light-2406-ragdoll
ContinuumSuppression.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/ContinuumSuppression/ContinuumSuppression.h>
10#include <analysis/ContinuumSuppression/Thrust.h>
11#include <analysis/ContinuumSuppression/KsfwMoments.h>
12#include <analysis/ContinuumSuppression/FoxWolfram.h>
13#include <analysis/ContinuumSuppression/CleoCones.h>
14#include <analysis/dataobjects/RestOfEvent.h>
15#include <analysis/dataobjects/ContinuumSuppression.h>
16#include <analysis/utility/PCmsLabTransform.h>
17#include <framework/datastore/StoreArray.h>
18#include <framework/logging/Logger.h>
19
20#include <vector>
21#include <Math/Vector3D.h>
22
23namespace Belle2 {
29 void addContinuumSuppression(const Particle* particle, const std::string& maskName)
30 {
31 // Output
32 StoreArray<ContinuumSuppression> qqArray(maskName);
33 // Create ContinuumSuppression object
34 ContinuumSuppression* qqVars = qqArray.appendNew();
35
36 // Create relation: Particle <-> ContinuumSuppression
37 particle->addRelationTo(qqVars);
38
39 std::vector<ROOT::Math::XYZVector> p3_cms_sigB, p3_cms_roe, p3_cms_all;
40
41 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigA;
42 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigB;
43 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_roe;
44
45 std::vector<float> ksfwFS0;
46 std::vector<float> ksfwFS1;
47
48 std::vector<float> cleoConesAll;
49 std::vector<float> cleoConesROE;
50
51 double et[2];
52
53 ROOT::Math::XYZVector thrustB;
54 ROOT::Math::XYZVector thrustO;
55
56 float thrustBm = -1;
57 float thrustOm = -1;
58 float cosTBTO = -1;
59 float cosTBz = -1;
60 float R2 = -1;
61
62
63 // -- B Cand --------------------------------------------------------------------------
65 double BeamEnergy = T.getCMSEnergy() / 2;
66
67 ROOT::Math::PxPyPzEVector p_cms_missA(0, 0, 0, 2 * BeamEnergy);
68 ROOT::Math::PxPyPzEVector p_cms_missB(0, 0, 0, 2 * BeamEnergy);
69 et[0] = et[1] = 0;
70
71 // -- SIG A --- Use B primary daughters - (Belle: use_finalstate_for_sig == 0) --------
72 std::vector<Belle2::Particle*> signalDaughters = particle->getDaughters();
73
74 for (const Belle2::Particle* sigFS0 : signalDaughters) {
75 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS0->get4Vector();
76
77 p3_cms_q_sigA.emplace_back(p_cms.Vect(), sigFS0->getCharge());
78
79 p_cms_missA -= p_cms;
80 et[0] += p_cms.Pt();
81 }
82
83 // -- SIG B --- Use B final-state daughters - (Belle: use_finalstate_for_sig == 1) ----
84 std::vector<const Belle2::Particle*> signalFSParticles = particle->getFinalStateDaughters();
85
86 for (const Belle2::Particle* sigFS1 : signalFSParticles) {
87 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS1->get4Vector();
88
89 p3_cms_all.push_back(p_cms.Vect());
90 p3_cms_sigB.push_back(p_cms.Vect());
91
92 p3_cms_q_sigB.emplace_back(p_cms.Vect(), sigFS1->getCharge());
93
94 p_cms_missB -= p_cms;
95 et[1] += p_cms.Pt();
96 }
97
98 // -- ROE -----------------------------------------------------------------------------
99 const RestOfEvent* roe = particle->getRelated<RestOfEvent>();
100
101 if (roe) {
102
103 // Charged tracks
104 //
105 std::vector<const Particle*> chargedROEParticles = roe->getChargedParticles(maskName);
106
107 for (const Particle* chargedROEParticle : chargedROEParticles) {
108
109 // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
110 // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
111
112 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * chargedROEParticle->get4Vector();
113
114 p3_cms_all.push_back(p_cms.Vect());
115 p3_cms_roe.push_back(p_cms.Vect());
116
117 p3_cms_q_roe.emplace_back(p_cms.Vect(), chargedROEParticle->getCharge());
118
119 p_cms_missA -= p_cms;
120 p_cms_missB -= p_cms;
121 et[0] += p_cms.Pt();
122 et[1] += p_cms.Pt();
123 }
124
125 // ECLCluster
126 //
127 std::vector<const Particle*> roePhotons = roe->getPhotons(maskName);
128
129 for (const Particle* photon : roePhotons) {
130
131 if (photon->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
132
133 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * photon->get4Vector();
134 p3_cms_all.push_back(p_cms.Vect());
135 p3_cms_roe.push_back(p_cms.Vect());
136
137 p3_cms_q_roe.emplace_back(p_cms.Vect(), photon->getCharge());
138
139 p_cms_missA -= p_cms;
140 p_cms_missB -= p_cms;
141 et[0] += p_cms.Pt();
142 et[1] += p_cms.Pt();
143 }
144 }
145
146 // Thrust variables
147 thrustB = Thrust::calculateThrust(p3_cms_sigB);
148 thrustO = Thrust::calculateThrust(p3_cms_roe);
149 thrustBm = thrustB.R();
150 thrustOm = thrustO.R();
151 cosTBTO = fabs(thrustB.Unit().Dot(thrustO.Unit()));
152 cosTBz = fabs(cos(thrustB.Theta()));
153
154 // Cleo Cones
155 CleoCones cc(p3_cms_all, p3_cms_roe, thrustB, true, true);
156 cleoConesAll = cc.cleo_cone_with_all();
157 cleoConesROE = cc.cleo_cone_with_roe();
158
159 // Fox-Wolfram Moments: Uses all final-state tracks (= sigB + ROE)
160 FoxWolfram FW(p3_cms_all);
162 R2 = FW.getR(2);
163
164 // KSFW moments
165 ROOT::Math::PxPyPzEVector p_cms_B = T.rotateLabToCms() * particle->get4Vector();
166 double Hso0_max(2 * (2 * BeamEnergy - p_cms_B.E()));
167 KsfwMoments KsfwM(Hso0_max,
168 p3_cms_q_sigA,
169 p3_cms_q_sigB,
170 p3_cms_q_roe,
171 p_cms_missA,
172 p_cms_missB,
173 et);
174 // use_finalstate_for_sig == 0
175 KsfwM.usefinal(0);
176 ksfwFS0.push_back(KsfwM.mm2());
177 ksfwFS0.push_back(KsfwM.et());
178 ksfwFS0.push_back(KsfwM.Hso(0, 0));
179 ksfwFS0.push_back(KsfwM.Hso(0, 1));
180 ksfwFS0.push_back(KsfwM.Hso(0, 2));
181 ksfwFS0.push_back(KsfwM.Hso(0, 3));
182 ksfwFS0.push_back(KsfwM.Hso(0, 4));
183 ksfwFS0.push_back(KsfwM.Hso(1, 0));
184 ksfwFS0.push_back(KsfwM.Hso(1, 2));
185 ksfwFS0.push_back(KsfwM.Hso(1, 4));
186 ksfwFS0.push_back(KsfwM.Hso(2, 0));
187 ksfwFS0.push_back(KsfwM.Hso(2, 2));
188 ksfwFS0.push_back(KsfwM.Hso(2, 4));
189 ksfwFS0.push_back(KsfwM.Hoo(0));
190 ksfwFS0.push_back(KsfwM.Hoo(1));
191 ksfwFS0.push_back(KsfwM.Hoo(2));
192 ksfwFS0.push_back(KsfwM.Hoo(3));
193 ksfwFS0.push_back(KsfwM.Hoo(4));
194 // use_finalstate_for_sig == 1
195 KsfwM.usefinal(1);
196 ksfwFS1.push_back(KsfwM.mm2());
197 ksfwFS1.push_back(KsfwM.et());
198 ksfwFS1.push_back(KsfwM.Hso(0, 0));
199 ksfwFS1.push_back(KsfwM.Hso(0, 1));
200 ksfwFS1.push_back(KsfwM.Hso(0, 2));
201 ksfwFS1.push_back(KsfwM.Hso(0, 3));
202 ksfwFS1.push_back(KsfwM.Hso(0, 4));
203 ksfwFS1.push_back(KsfwM.Hso(1, 0));
204 ksfwFS1.push_back(KsfwM.Hso(1, 2));
205 ksfwFS1.push_back(KsfwM.Hso(1, 4));
206 ksfwFS1.push_back(KsfwM.Hso(2, 0));
207 ksfwFS1.push_back(KsfwM.Hso(2, 2));
208 ksfwFS1.push_back(KsfwM.Hso(2, 4));
209 ksfwFS1.push_back(KsfwM.Hoo(0));
210 ksfwFS1.push_back(KsfwM.Hoo(1));
211 ksfwFS1.push_back(KsfwM.Hoo(2));
212 ksfwFS1.push_back(KsfwM.Hoo(3));
213 ksfwFS1.push_back(KsfwM.Hoo(4));
214
215 // TODO: The following is from the original belle ksfwmoments.cc module.
216 // Not sure if necessary here (i.e., will we be using rooksfw in belle II in the same way?).
217 // printf("rooksfw::rooksfw: mm2=%f et=%f hoo2=%f hso02=%f\n",
218 // m_mm2[0], et[0], m_Hoo[0][2], m_Hso[0][0][2]);
219 }
220
221 // Fill ContinuumSuppression with content
222 qqVars->addThrustB(thrustB);
223 qqVars->addThrustO(thrustO);
224 qqVars->addThrustBm(thrustBm);
225 qqVars->addThrustOm(thrustOm);
226 qqVars->addCosTBTO(cosTBTO);
227 qqVars->addCosTBz(cosTBz);
228 qqVars->addR2(R2);
229 qqVars->addKsfwFS0(ksfwFS0);
230 qqVars->addKsfwFS1(ksfwFS1);
231 qqVars->addCleoConesALL(cleoConesAll);
232 qqVars->addCleoConesROE(cleoConesROE);
233 }
235}
Class to calculate the Cleo clone variables.
Definition: CleoCones.h:22
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)
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:35
int usefinal(int uf)
Sets the flag that specifiies we are using the finalstate for signal.
Definition: KsfwMoments.h:78
double et(int uf=-1) const
Returns calculated transverse energy.
Definition: KsfwMoments.h:93
double Hso(int i, int j, int uf=-1) const
Returns calculated KSFW Moments.
Definition: KsfwMoments.h:103
double mm2(int uf=-1) const
Returns calculated missing mass squared.
Definition: KsfwMoments.h:88
double Hoo(int i, int uf=-1) const
Returns calculated KSFW Moments.
Definition: KsfwMoments.h:98
Class to hold Lorentz transformations from/to CMS and boost vector.
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
const ROOT::Math::LorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Class to store reconstructed particles.
Definition: Particle.h:75
std::vector< const Belle2::Particle * > getFinalStateDaughters() const
Returns a vector of pointers to Final State daughter particles.
Definition: Particle.cc:649
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:637
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
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.
This is a general purpose class for collecting reconstructed MDST data objects that are not used in r...
Definition: RestOfEvent.h:57
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.
Definition: RestOfEvent.cc:108
std::vector< const Particle * > getPhotons(const std::string &maskName=c_defaultMaskName, bool unpackComposite=true) const
Get photons from ROE mask.
Definition: RestOfEvent.cc:84
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::XYZVector > &momenta)
calculates the thrust axis
Definition: Thrust.cc:71
void addContinuumSuppression(const Particle *particle, const std::string &maskName)
Adds continuum suppression variables.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24