Belle II Software  release-05-01-25
ContinuumSuppression.cc
1 /* BASF2 (Belle Analysis Framework 2) *
2  * Copyright(C) 2013 - Belle II Collaboration *
3  * *
4  * Author: The Belle II Collaboration *
5  * Contributors: Pablo Goldenzweig, *
6  * *
7  * This software is provided "as is" without any warranty. *
8  **************************************************************************/
9 
10 #include <analysis/ContinuumSuppression/ContinuumSuppression.h>
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 #include <analysis/dataobjects/RestOfEvent.h>
16 #include <analysis/dataobjects/ContinuumSuppression.h>
17 #include <analysis/utility/PCmsLabTransform.h>
18 #include <mdst/dataobjects/Track.h>
19 #include <mdst/dataobjects/ECLCluster.h>
20 #include <mdst/dataobjects/KLMCluster.h>
21 #include <mdst/dataobjects/PIDLikelihood.h>
22 #include <framework/datastore/StoreArray.h>
23 
24 #include <vector>
25 
26 
27 namespace Belle2 {
33  void addContinuumSuppression(const Particle* particle, const std::string& maskName)
34  {
35  // Output
36  StoreArray<ContinuumSuppression> qqArray;
37  // Create ContinuumSuppression object
38  ContinuumSuppression* qqVars = qqArray.appendNew();
39 
40  // Create relation: Particle <-> ContinuumSuppression
41  particle->addRelationTo(qqVars);
42 
43  std::vector<TVector3> p3_cms_sigB, p3_cms_roe, p3_cms_all;
44 
45  std::vector<std::pair<TVector3, int>> p3_cms_q_sigA;
46  std::vector<std::pair<TVector3, int>> p3_cms_q_sigB;
47  std::vector<std::pair<TVector3, int>> p3_cms_q_roe;
48 
49  std::vector<float> ksfwFS0;
50  std::vector<float> ksfwFS1;
51 
52  std::vector<float> cleoConesAll;
53  std::vector<float> cleoConesROE;
54 
55  double et[2];
56 
57  TVector3 thrustB;
58  TVector3 thrustO;
59 
60  float thrustBm = -1;
61  float thrustOm = -1;
62  float cosTBTO = -1;
63  float cosTBz = -1;
64  float R2 = -1;
65 
66 
67  // -- B Cand --------------------------------------------------------------------------
68  PCmsLabTransform T;
69  double BeamEnergy = T.getCMSEnergy() / 2;
70 
71  TLorentzVector p_cms_missA(0, 0, 0, 2 * BeamEnergy);
72  TLorentzVector p_cms_missB(0, 0, 0, 2 * BeamEnergy);
73  et[0] = et[1] = 0;
74 
75  // -- SIG A --- Use B primary daughters - (Belle: use_finalstate_for_sig == 0) --------
76  std::vector<Belle2::Particle*> signalDaughters = particle->getDaughters();
77 
78  for (const Belle2::Particle* sigFS0 : signalDaughters) {
79  TLorentzVector p_cms = T.rotateLabToCms() * sigFS0->get4Vector();
80 
81  p3_cms_q_sigA.emplace_back(p_cms.Vect(), sigFS0->getCharge());
82 
83  p_cms_missA -= p_cms;
84  et[0] += p_cms.Perp();
85  }
86 
87  // -- SIG B --- Use B final-state daughters - (Belle: use_finalstate_for_sig == 1) ----
88  std::vector<const Belle2::Particle*> signalFSParticles = particle->getFinalStateDaughters();
89 
90  for (const Belle2::Particle* sigFS1 : signalFSParticles) {
91  TLorentzVector p_cms = T.rotateLabToCms() * sigFS1->get4Vector();
92 
93  p3_cms_all.push_back(p_cms.Vect());
94  p3_cms_sigB.push_back(p_cms.Vect());
95 
96  p3_cms_q_sigB.emplace_back(p_cms.Vect(), sigFS1->getCharge());
97 
98  p_cms_missB -= p_cms;
99  et[1] += p_cms.Perp();
100  }
101 
102  // -- ROE -----------------------------------------------------------------------------
103  const RestOfEvent* roe = particle->getRelated<RestOfEvent>();
104 
105  if (roe) {
106 
107  // Charged tracks -> Pion
108  //
109  std::vector<const Track*> roeTracks = roe->getTracks(maskName);
110 
111  for (const Track* track : roeTracks) {
112 
113  // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
114  // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
115 
116  // Create particle from track with most probable hypothesis
117  const PIDLikelihood* iPidLikelihood = track->getRelated<PIDLikelihood>();
118  const Const::ChargedStable charged = iPidLikelihood ? iPidLikelihood->getMostLikely() : Const::pion;
119  // Here we skip tracks with 0 charge
120  if (track->getTrackFitResultWithClosestMass(charged)->getChargeSign() == 0) {
121  B2WARNING("Track with charge = 0 skipped! From the ContinuumSuppression");
122  continue;
123  }
124  Particle charged_particle(track, charged);
125  if (charged_particle.getParticleSource() == Particle::c_Track) {
126  TLorentzVector p_cms = T.rotateLabToCms() * charged_particle.get4Vector();
127 
128  p3_cms_all.push_back(p_cms.Vect());
129  p3_cms_roe.push_back(p_cms.Vect());
130 
131  p3_cms_q_roe.emplace_back(p_cms.Vect(), charged_particle.getCharge());
132 
133  p_cms_missA -= p_cms;
134  p_cms_missB -= p_cms;
135  et[0] += p_cms.Perp();
136  et[1] += p_cms.Perp();
137  }
138  }
139 
140  // ECLCluster -> Gamma
141  //
142  std::vector<const ECLCluster*> roeECLClusters = roe->getECLClusters(maskName);
143 
144  for (const ECLCluster* cluster : roeECLClusters) {
145 
146  if (cluster->isNeutral() and cluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
147  // Create particle from ECLCluster with gamma hypothesis
148  Particle gamma_particle(cluster, Const::photon);
149 
150  TLorentzVector p_cms = T.rotateLabToCms() * gamma_particle.get4Vector();
151  p3_cms_all.push_back(p_cms.Vect());
152  p3_cms_roe.push_back(p_cms.Vect());
153 
154  p3_cms_q_roe.emplace_back(p_cms.Vect(), gamma_particle.getCharge());
155 
156  p_cms_missA -= p_cms;
157  p_cms_missB -= p_cms;
158  et[0] += p_cms.Perp();
159  et[1] += p_cms.Perp();
160  }
161  }
162 
163  // Thrust variables
164  thrustB = Thrust::calculateThrust(p3_cms_sigB);
165  thrustO = Thrust::calculateThrust(p3_cms_roe);
166  thrustBm = thrustB.Mag();
167  thrustOm = thrustO.Mag();
168  cosTBTO = fabs(cos(thrustB.Angle(thrustO)));
169  cosTBz = fabs(thrustB.CosTheta());
170 
171  // Cleo Cones
172  CleoCones cc(p3_cms_all, p3_cms_roe, thrustB, true, true);
173  cleoConesAll = cc.cleo_cone_with_all();
174  cleoConesROE = cc.cleo_cone_with_roe();
175 
176  // Fox-Wolfram Moments: Uses all final-state tracks (= sigB + ROE)
177  FoxWolfram FW(p3_cms_all);
178  FW.calculateBasicMoments();
179  R2 = FW.getR(2);
180 
181  // KSFW moments
182  TLorentzVector p_cms_B = T.rotateLabToCms() * particle->get4Vector();
183  double Hso0_max(2 * (2 * BeamEnergy - p_cms_B.E()));
184  KsfwMoments KsfwM(Hso0_max,
185  p3_cms_q_sigA,
186  p3_cms_q_sigB,
187  p3_cms_q_roe,
188  p_cms_missA,
189  p_cms_missB,
190  et);
191  // use_finalstate_for_sig == 0
192  KsfwM.usefinal(0);
193  ksfwFS0.push_back(KsfwM.mm2());
194  ksfwFS0.push_back(KsfwM.et());
195  ksfwFS0.push_back(KsfwM.Hso(0, 0));
196  ksfwFS0.push_back(KsfwM.Hso(0, 1));
197  ksfwFS0.push_back(KsfwM.Hso(0, 2));
198  ksfwFS0.push_back(KsfwM.Hso(0, 3));
199  ksfwFS0.push_back(KsfwM.Hso(0, 4));
200  ksfwFS0.push_back(KsfwM.Hso(1, 0));
201  ksfwFS0.push_back(KsfwM.Hso(1, 2));
202  ksfwFS0.push_back(KsfwM.Hso(1, 4));
203  ksfwFS0.push_back(KsfwM.Hso(2, 0));
204  ksfwFS0.push_back(KsfwM.Hso(2, 2));
205  ksfwFS0.push_back(KsfwM.Hso(2, 4));
206  ksfwFS0.push_back(KsfwM.Hoo(0));
207  ksfwFS0.push_back(KsfwM.Hoo(1));
208  ksfwFS0.push_back(KsfwM.Hoo(2));
209  ksfwFS0.push_back(KsfwM.Hoo(3));
210  ksfwFS0.push_back(KsfwM.Hoo(4));
211  // use_finalstate_for_sig == 1
212  KsfwM.usefinal(1);
213  ksfwFS1.push_back(KsfwM.mm2());
214  ksfwFS1.push_back(KsfwM.et());
215  ksfwFS1.push_back(KsfwM.Hso(0, 0));
216  ksfwFS1.push_back(KsfwM.Hso(0, 1));
217  ksfwFS1.push_back(KsfwM.Hso(0, 2));
218  ksfwFS1.push_back(KsfwM.Hso(0, 3));
219  ksfwFS1.push_back(KsfwM.Hso(0, 4));
220  ksfwFS1.push_back(KsfwM.Hso(1, 0));
221  ksfwFS1.push_back(KsfwM.Hso(1, 2));
222  ksfwFS1.push_back(KsfwM.Hso(1, 4));
223  ksfwFS1.push_back(KsfwM.Hso(2, 0));
224  ksfwFS1.push_back(KsfwM.Hso(2, 2));
225  ksfwFS1.push_back(KsfwM.Hso(2, 4));
226  ksfwFS1.push_back(KsfwM.Hoo(0));
227  ksfwFS1.push_back(KsfwM.Hoo(1));
228  ksfwFS1.push_back(KsfwM.Hoo(2));
229  ksfwFS1.push_back(KsfwM.Hoo(3));
230  ksfwFS1.push_back(KsfwM.Hoo(4));
231 
232  // TODO: The following is from the original belle ksfwmoments.cc module.
233  // Not sure if necessary here (i.e., will we be using rooksfw in belle II in the same way?).
234  // printf("rooksfw::rooksfw: mm2=%f et=%f hoo2=%f hso02=%f\n",
235  // m_mm2[0], et[0], m_Hoo[0][2], m_Hso[0][0][2]);
236  }
237 
238  // Fill ContinuumSuppression with content
239  qqVars->addThrustB(thrustB);
240  qqVars->addThrustO(thrustO);
241  qqVars->addThrustBm(thrustBm);
242  qqVars->addThrustOm(thrustOm);
243  qqVars->addCosTBTO(cosTBTO);
244  qqVars->addCosTBz(cosTBz);
245  qqVars->addR2(R2);
246  qqVars->addKsfwFS0(ksfwFS0);
247  qqVars->addKsfwFS1(ksfwFS1);
248  qqVars->addCleoConesALL(cleoConesAll);
249  qqVars->addCleoConesROE(cleoConesROE);
250  }
252 }
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::Const::photon
static const ParticleType photon
photon particle
Definition: Const.h:547
Belle2::ContinuumSuppression::addCosTBTO
void addCosTBTO(float cosTBTO)
Add cosine of the angle between the thrust axis of the B and the thrust axis of the ROE.
Definition: ContinuumSuppression.cc:35
Belle2::ContinuumSuppression::addR2
void addR2(float R2)
Add reduced Fox-Wolfram moment R2.
Definition: ContinuumSuppression.cc:45
Belle2::PIDLikelihood::getMostLikely
Const::ChargedStable getMostLikely(const double *fractions=0, Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Return most likely particle among chargedStableSet; if prior fractions not given equal prior probabil...
Definition: PIDLikelihood.cc:108
Belle2::ContinuumSuppression::addCleoConesROE
void addCleoConesROE(const std::vector< float > &cleoConesROE)
Add vector of Cleo Cones constructed of only ROE particles.
Definition: ContinuumSuppression.cc:65
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::Thrust::calculateThrust
static TVector3 calculateThrust(const std::vector< TVector3 > &momenta)
calculates the thrust axis
Definition: Thrust.cc:5
Belle2::ContinuumSuppression::addCosTBz
void addCosTBz(float cosTBz)
Add cosine of the angle between the thrust axis of the B and the z-axis.
Definition: ContinuumSuppression.cc:40
Belle2::PCmsLabTransform::getCMSEnergy
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
Definition: PCmsLabTransform.h:57
Belle2::KsfwMoments::Hso
double Hso(int i, int j, int uf=-1) const
Returns calculated KSFW Moments.
Definition: KsfwMoments.h:117
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::FoxWolfram::calculateBasicMoments
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition: FoxWolfram.cc:20
Belle2::RestOfEvent
This is a general purpose class for collecting reconstructed MDST data objects that are not used in r...
Definition: RestOfEvent.h:59
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::KsfwMoments::usefinal
int usefinal(int uf)
Sets the flag that specifiies we are using the finalstate for signal.
Definition: KsfwMoments.h:92
Belle2::FoxWolfram::getR
double getR(int i) const
Returns the i-th moment normalized to the 0th-order moment.
Definition: FoxWolfram.h:109
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::Particle::getCharge
float getCharge(void) const
Returns particle charge.
Definition: Particle.cc:586
Belle2::ContinuumSuppression::addKsfwFS1
void addKsfwFS1(const std::vector< float > &ksfwFS1)
Add vector of KSFW moments, Et, and mm2 for final state = 1.
Definition: ContinuumSuppression.cc:55
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ContinuumSuppression::addThrustO
void addThrustO(const TVector3 &thrustO)
Add ROE thrust axis.
Definition: ContinuumSuppression.cc:20
Belle2::ContinuumSuppression::addCleoConesALL
void addCleoConesALL(const std::vector< float > &cleoConesALL)
Add vector of Cleo Cones constructed of all final state particles.
Definition: ContinuumSuppression.cc:60
Belle2::ContinuumSuppression::addThrustBm
void addThrustBm(float thrustBm)
Add magnitude of B thrust axis.
Definition: ContinuumSuppression.cc:25
Belle2::KsfwMoments::et
double et(int uf=-1) const
Returns calculated transverse energy.
Definition: KsfwMoments.h:107
Belle2::ContinuumSuppression
This is a class for collecting variables used in continuum suppression.
Definition: ContinuumSuppression.h:52
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::RestOfEvent::getTracks
std::vector< const Track * > getTracks(const std::string &maskName="") const
Get vector of all (no mask) or a subset (use mask) of all Tracks in ROE.
Definition: RestOfEvent.cc:308
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:37
Belle2::PCmsLabTransform::rotateLabToCms
const TLorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:74
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::KsfwMoments::Hoo
double Hoo(int i, int uf=-1) const
Returns calculated KSFW Moments.
Definition: KsfwMoments.h:112
Belle2::RestOfEvent::getECLClusters
std::vector< const ECLCluster * > getECLClusters(const std::string &maskName="") const
Get vector of all (no mask) or a subset (use mask) of all ECLClusters in ROE.
Definition: RestOfEvent.cc:319
Belle2::KsfwMoments::mm2
double mm2(int uf=-1) const
Returns calculated missing mass squared.
Definition: KsfwMoments.h:102
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::CleoCones
Class to calculate the Cleo clone variables.
Definition: CleoCones.h:32
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::Particle::getParticleSource
EParticleSourceObject getParticleSource() const
Returns particle source as defined with enum EParticleSourceObject.
Definition: Particle.h:404
Belle2::KsfwMoments
Moment-calculation of the k_sfw improved Super-Fox-Wolfram moments.
Definition: KsfwMoments.h:49
Belle2::ContinuumSuppression::addThrustOm
void addThrustOm(float thrustOm)
Add magnitude of ROE thrust axis.
Definition: ContinuumSuppression.cc:30
Belle2::ContinuumSuppression::addKsfwFS0
void addKsfwFS0(const std::vector< float > &ksfwFS0)
Add vector of KSFW moments, Et, and mm2 for final state = 0.
Definition: ContinuumSuppression.cc:50
Belle2::ContinuumSuppression::addThrustB
void addThrustB(const TVector3 &thrustB)
Add ROE thrust axis.
Definition: ContinuumSuppression.cc:15
Belle2::FoxWolfram
Class to calculate the Fox-Wolfram moments up to order 8.
Definition: FoxWolfram.h:48
Belle2::addContinuumSuppression
void addContinuumSuppression(const Particle *particle, const std::string &maskName)
Adds continuum suppression variables.
Definition: ContinuumSuppression.cc:33