Belle II Software development
InitialParticleGeneration.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 <generators/utilities/InitialParticleGeneration.h>
10#include <framework/logging/Logger.h>
11#include <generators/utilities/beamHelpers.h>
12#include <cmath>
13#include <Eigen/Dense>
14
15
16namespace Belle2 {
21
23 {
24 m_event.registerInDataStore();
25 }
26
27 ROOT::Math::XYZVector InitialParticleGeneration::generateVertex(const ROOT::Math::XYZVector& initial, const TMatrixDSym& cov,
29 {
30 if (m_event->hasGenerationFlags(BeamParameters::c_smearVertex)) {
31 if (!gen.size()) gen.setMeanCov(initial, cov);
32 return gen.generateVec3();
33 }
34 return initial;
35 }
36
37
38 TMatrixDSym InitialParticleGeneration::adjustCovMatrix(TMatrixDSym cov) const
39 {
40 // no smearing
41 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeam))
42 return TMatrixDSym(3);
43
44 // don't smear energy
45 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeamEnergy)) {
46 for (int i = 0; i < 3; ++i)
47 cov(i, 0) = cov(0, i) = 0;
48 }
49
50 // don't smear angles
51 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeamDirection)) {
52 for (int i = 1; i < 3; ++i)
53 for (int j = 1; j < 3; ++j)
54 cov(i, j) = cov(j, i) = 0;
55 }
56
57 return cov;
58
59 }
60
61
63 static Eigen::MatrixXd toEigen(TMatrixDSym m)
64 {
65 Eigen::MatrixXd mEig(m.GetNrows(), m.GetNcols());
66 for (int i = 0; i < m.GetNrows(); ++i)
67 for (int j = 0; j < m.GetNcols(); ++j)
68 mEig(i, j) = m(i, j);
69
70 return mEig;
71 }
72
73
74 // init random number generator which generates random boost based on the beam parameters
76 const ROOT::Math::PxPyPzEVector& pLER, const TMatrixDSym& covHER, const TMatrixDSym& covLER)
77 {
78 //Calculate initial parameters of the HER beam before smearing
79 double E0her = pHER.E();
80 double thX0her = std::atan(pHER.Px() / pHER.Pz());
81 double thY0her = std::atan(pHER.Py() / pHER.Pz());
82
83 //Calculate initial parameters of the LER beam before smearing
84 double E0ler = pLER.E();
85 double thX0ler = std::atan(pLER.Px() / pLER.Pz());
86 double thY0ler = std::atan(pLER.Py() / pLER.Pz());
87
88 // calculate gradient of transformation (EH, thXH, thYH, EL, thXL, thYL) -> (Ecms, bX, bY, bZ, angleCmsX, angleCmsY)
89 Eigen::MatrixXd grad = getGradientMatrix(E0her, thX0her, thY0her, E0ler, thX0ler, thY0ler);
90 Eigen::VectorXd mu = getCentralValues(E0her, thX0her, thY0her, E0ler, thX0ler, thY0ler);
91
92 // calculate covariance smearing matrix in the (Ecms, bX, bY, bZ, angleCmsX, angleCmsY) coordinate system
93 Eigen::MatrixXd covT = transformCov(toEigen(adjustCovMatrix(covHER)), toEigen(adjustCovMatrix(covLER)), grad);
94
95 // return random number generator which allows to generate Lorentz transformation parameters based on the cmsEnergy
96 return ConditionalGaussGenerator(mu, covT);
97
98 }
99
100
101 ROOT::Math::PxPyPzEVector InitialParticleGeneration::generateBeam(const ROOT::Math::PxPyPzEVector& initial, const TMatrixDSym& cov,
103 {
104 //no smearing, nothing to do
105 if (!(m_event->getGenerationFlags() & BeamParameters::c_smearBeam))
106 return initial;
107
108
109 //Calculate initial parameters of the beam before smearing
110 double E0 = initial.E();
111 double thX0 = std::atan(initial.Px() / initial.Pz());
112 double thY0 = std::atan(initial.Py() / initial.Pz());
113
114 //init random generator if there is a change in beam parameters or at the beginning
115 if (gen.size() != 3) gen.setMeanCov(ROOT::Math::XYZVector(E0, thX0, thY0), cov);
116
117 //generate the actual smeared vector (E, thX, thY)
118 Eigen::VectorXd p = gen.generate();
119
120 //if we don't smear beam energy set the E to initial value
121 if (!m_event->hasGenerationFlags(BeamParameters::c_smearBeamEnergy))
122 p[0] = E0;
123
124 //if we don't smear beam direction set thX and thY to initial values
125 if (!m_event->hasGenerationFlags(BeamParameters::c_smearBeamDirection)) {
126 p[1] = thX0;
127 p[2] = thY0;
128 }
129
130 //store smeared values
131 double E = p[0];
132 double thX = p[1];
133 double thY = p[2];
134
135 //Convert values back to 4-vector
136 bool isHER = initial.Z() > 0;
137 return BeamParameters::getFourVector(E, thX, thY, isHER);
138 }
139
144
146 {
147 if (!m_event) {
148 m_event.create();
149 }
150 if (!m_beamParams.isValid()) {
151 B2FATAL("Cannot generate beam without valid BeamParameters");
152 }
153 if (m_beamParams.hasChanged()) {
154 m_generateHER.reset();
155 m_generateLER.reset();
156 m_generateVertex.reset();
157 }
158 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & allowedFlags);
159 ROOT::Math::PxPyPzEVector her = generateBeam(m_beamParams->getHER(), m_beamParams->getCovHER(), m_generateHER);
160 ROOT::Math::PxPyPzEVector ler = generateBeam(m_beamParams->getLER(), m_beamParams->getCovLER(), m_generateLER);
161 ROOT::Math::XYZVector vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
162 m_event->set(her, ler, vtx);
163 //Check if we want to go to CMS, if so boost both
164 if (m_beamParams->hasGenerationFlags(BeamParameters::c_generateCMS)) {
165 m_event->setGenerationFlags(0);
166 her = m_event->getLabToCMS() * her;
167 ler = m_event->getLabToCMS() * ler;
168 m_event->set(her, ler, vtx);
169 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & allowedFlags);
170 }
171 return *m_event;
172 }
173
174
175 // it generates random Lorentz transformation based on the HER/LER smearing covariance matrices
177 {
178 if (!m_event) {
179 m_event.create();
180 }
181 if (!m_beamParams.isValid()) {
182 B2FATAL("Cannot generate beam without valid BeamParameters");
183 }
184 if (m_beamParams.hasChanged()) {
186 m_beamParams->getCovHER(), m_beamParams->getCovLER());
187 m_generateVertex.reset();
188 }
189 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & m_allowedFlags);
190 ROOT::Math::XYZVector vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
191
192 //TODO is m_event of type MCInitialParticles important?
193 //Eigen::VectorXd collSys = m_generateLorentzTransformation.generate(Ecms);
194 //m_event->set(her, ler, vtx);
195 //m_event->setByLorentzTransformation(collSys[0], collSys[1], collSys[2], collSys[3], collSys[4], collSys[5], vtx);
196
197 return vtx;
198 }
199
200
201
202 ROOT::Math::XYZVector InitialParticleGeneration::updateVertex(bool force)
203 {
204 if (!m_beamParams.isValid()) {
205 B2FATAL("Cannot generate beam without valid BeamParameters");
206 }
207 if (!m_event) {
208 // generate a new mc initial particle without smearing except for the vertex
210 return m_event->getVertex();
211 }
212 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearVertex) or
213 (m_event->hasGenerationFlags(BeamParameters::c_smearVertex) and not force)) {
214 // already has a smeared vertex or smearing disallowed. nothing to do
215 return {0, 0, 0};
216 }
217 // Ok, now we need to just update the vertex, make sure the parameters are up to date ...
218 if (m_beamParams.hasChanged()) {
219 m_generateHER.reset();
220 m_generateLER.reset();
221 m_generateVertex.reset();
222 }
223 // Add the smear vertex flag
224 m_event->setGenerationFlags(m_event->getGenerationFlags() | BeamParameters::c_smearVertex);
225 // And generate a vertex ...
226 auto previous = m_event->getVertex();
227 auto vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
228 m_event->setVertex(vtx);
229 // Everything done
230 return vtx - previous;
231 }
232
233} //Belle2 namespace
R E
internal precision of FFTW codelets
static ROOT::Math::PxPyPzEVector getFourVector(double energy, double angleX, double angleY, bool isHER)
Calculate FourVector of a beam from energy and angles in xz and yz planes.
Class implementing n-dimensional random number generator from Gaussian distribution where the first c...
DBObjPtr< BeamParameters > m_beamParams
Datastore object containing the nominal beam parameters.
int m_allowedFlags
Allowed generation flags.
MultivariateNormalGenerator m_generateHER
Generator for HER.
StoreObjPtr< MCInitialParticles > m_event
Datastore object containing the generated event.
MultivariateNormalGenerator m_generateLER
Generator for LER.
MultivariateNormalGenerator m_generateVertex
Generator for Vertex.
ConditionalGaussGenerator m_generateLorentzTransformation
Generator of the Lorentz transformation.
This class contains the initial state for the given event.
@ c_smearBeam
smear the full beam momentum (energy and direction)
@ c_generateCMS
generate initial event in CMS instead of lab
@ c_smearBeamEnergy
smear energy of HER and LER (but not direction)
@ c_smearBeamDirection
smear direction of HER and LER (but not energy)
Class to generate normal distributed, correlated random numbers given the mean values and the covaria...
size_t size() const
Return the number of elements to be generated on generate()
ROOT::Math::XYZVector generateVec3() const
Generate a set of correlated random numbers with the previously set mean and covariance and return a ...
bool setMeanCov(int n, const double *mean, const double *cov)
set the mean and covariance for the distribution with array interface: mean and covariance are passed...
Eigen::VectorXd generate() const
Generate a set of correlated random numbers with the previously set mean and covariance.
ROOT::Math::XYZVector getVertexConditional()
Generate vertex position and possibly update the generator of Lorentz transformation.
void initialize()
function to be executed on initialize()
ROOT::Math::XYZVector updateVertex(bool force=false)
Update the vertex position:
static Eigen::MatrixXd toEigen(TMatrixDSym m)
transform matrix from ROOT to Eigen format
MCInitialParticles & generate()
Generate a new event.
Eigen::MatrixXd transformCov(const Eigen::MatrixXd &covHER, const Eigen::MatrixXd &covLER, const Eigen::MatrixXd &grad)
transform covariance matrix to the variables (Ecms, boostX, boostY, boostZ, angleXZ,...
ROOT::Math::XYZVector generateVertex(const ROOT::Math::XYZVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
generate the vertex
ConditionalGaussGenerator initConditionalGenerator(const ROOT::Math::PxPyPzEVector &pHER, const ROOT::Math::PxPyPzEVector &pLER, const TMatrixDSym &covHER, const TMatrixDSym &covLER)
Initialize the conditional generator using HER & LER 4-vectors and HER & LER covariance matrices desc...
Eigen::VectorXd getCentralValues(double eH, double thXH, double thYH, double eL, double thXL, double thYL)
get vector including (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ) from beam energies and angles
TMatrixDSym adjustCovMatrix(TMatrixDSym cov) const
adjust smearing covariance matrix based on the generation flags
Eigen::MatrixXd getGradientMatrix(double eH, double thXH, double thYH, double eL, double thXL, double thYL)
get gradient matrix of (eH, thXH, thYH, eL, thXL, thYL) -> (eCMS, boostX, boostY, boostY,...
ROOT::Math::PxPyPzEVector generateBeam(const ROOT::Math::PxPyPzEVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
generate 4 vector for one beam
Abstract base class for different kinds of events.