Belle II Software development
BelleMCOutputModule.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/* Own header. */
10#include <b2bii/modules/BelleMCOutput/BelleMCOutputModule.h>
11
12/* Basf2 headers. */
13#include <framework/gearbox/Unit.h>
14#include <framework/logging/Logger.h>
15
16/* Belle headers. */
17#include <belle_legacy/tables/belletdf.h>
18#include <belle_legacy/tables/hepevt.h>
19#include <belle_legacy/tables/filespec.h>
20
21/* ROOT headers. */
22#include <TMatrixD.h>
23#include <Math/RotationY.h>
24#include <Math/Rotation3D.h>
25
26using namespace Belle2;
27
28REG_MODULE(BelleMCOutput);
29
31 Module()
32{
33 setDescription("Output of MC particle list in Belle format.");
35 addParam("outputFileName", m_OutputFileName,
36 "Output file name.", std::string("belle_mc.mdst"));
37 addParam("decayKsInGenerator", m_DecayKsInGenerator,
38 "Decay K_S0 in generator.", false);
39 addParam("decayLambdaInGenerator", m_DecayLambdaInGenerator,
40 "Decay Lambda in generator.", false);
41}
42
46
48{
49 m_BelleFile = new Belle::Panther_FileIO(m_OutputFileName.c_str(), BBS_WRITE);
50 m_BelleFile->write(BBS_FORMAT, 0);
51 m_BelleFile->save_br("BELLE_FILE_SPECIFICATION");
52 m_BelleFile->save_br("BELLE_RUNHEAD");
53 m_BelleFile->save_br("BELLE_NOMINAL_BEAM");
54 m_BelleFile->save("GEN_HEPEVT");
55 B2WARNING(
56 "************ ATTENTION! ************\n"
57 "\n"
58 " Belle MC generation differs from Belle II. In order to generate MC "
59 "correctly, you must do the following:\n"
60 "\n"
61 " 1. Use global tag with smearing of beam parameters and generation "
62 "flags: b2bii_beamParameters_with_smearing (there are no smearing data "
63 "or generation flags in the tags B2BII or B2BII_MC).\n"
64 " 2. Disable smearing in Belle simulation by removing the module "
65 "\"bpsmear\" from basf gsim scripts.\n"
66 "\n"
67 "It is also necessary to consider the following differences:\n"
68 "\n"
69 " 1. When converting the simulation result back to Belle II format, "
70 "it is recommended to disable conversion of beam parameters "
71 "(convertBeamParameters = False). Since the format of the table "
72 "BELLE_NOMINAL_BEAM is not sufficient to store an arbitrary covariance "
73 "matrix, the covariance matrix of IP position (BeamSpot) is not guaranteed "
74 "to be exactly the same after conversion.\n"
75 " 2. By default, decays of long-lived particles are removed and particles "
76 "are declared to be stable in generator (ISTHEP == 1 in basf) because "
77 "such decays are simulated by GEANT3 in basf. However, you may choose "
78 "to decay such particles via module parameters. This results in direct "
79 "passing of the decay products to simulation, i.e. material effects are "
80 "ignored. In this case, you may need to perform an additional study "
81 "of the difference between data and MC.\n"
82 "\n"
83 "************************************\n"
84 );
85}
86
88{
89 Belle::Belle_file_specification_Manager& fileManager =
90 Belle::Belle_file_specification_Manager::get_manager();
91 fileManager.remove();
92 int position = 0;
93 int length = m_OutputFileName.size();
94 do {
95 Belle::Belle_file_specification& file = fileManager.add();
96 file.Spec(m_OutputFileName.substr(position, std::min(length, 8)).c_str());
97 position += 8;
98 length -= 8;
99 } while (length > 0);
100 Belle::Belle_runhead_Manager& runheadManager =
101 Belle::Belle_runhead_Manager::get_manager();
102 runheadManager.remove();
103 Belle::Belle_runhead& runhead = runheadManager.add();
104 runhead.ExpMC(2);
105 runhead.ExpNo(m_EventMetaData->getExperiment());
106 runhead.RunNo(m_EventMetaData->getRun());
107 runhead.Time(time(nullptr));
108 runhead.Type(0);
109 ROOT::Math::PxPyPzEVector momentumLER = m_BeamParameters->getLER();
110 ROOT::Math::PxPyPzEVector momentumHER = m_BeamParameters->getHER();
111 runhead.ELER(momentumLER.E());
112 runhead.EHER(momentumHER.E());
113 Belle::Belle_nominal_beam_Manager& beamManager =
114 Belle::Belle_nominal_beam_Manager::get_manager();
115 beamManager.remove();
116 Belle::Belle_nominal_beam& beam = beamManager.add();
117 beam.px_high(momentumHER.X());
118 beam.py_high(momentumHER.Y());
119 beam.pz_high(momentumHER.Z());
120 beam.px_low(momentumLER.X());
121 beam.py_low(momentumLER.Y());
122 beam.pz_low(momentumLER.Z());
123 TMatrixDSym herCovariance = m_BeamParameters->getCovHER();
124 beam.sigma_p_high(sqrt(herCovariance[0][0]));
125 TMatrixDSym lerCovariance = m_BeamParameters->getCovLER();
126 beam.sigma_p_low(sqrt(lerCovariance[0][0]));
127 /*
128 * The vertex parameters are in cm, as in basf2.
129 * The unit is not the same as for particles in GEN_HEPEVT.
130 */
131 ROOT::Math::XYZVector vertex = m_BeamParameters->getVertex();
132 beam.ip_x(vertex.X());
133 beam.ip_y(vertex.Y());
134 beam.ip_z(vertex.Z());
135 TMatrixDSym vertexCovariance = m_BeamParameters->getCovVertex();
136 beam.cang_high(momentumHER.Theta());
137 beam.cang_low(M_PI - momentumLER.Theta());
138 double angleIPZX = momentumHER.Theta() / 2;
139 beam.angle_ip_zx(angleIPZX);
140 /*
141 * Transformation of error matrix. It is inverse to the transformation in
142 * belle_legacy/ip/IpProfile.cc.
143 */
144 ROOT::Math::RotationY rotationY(-angleIPZX);
145 ROOT::Math::Rotation3D rotation(rotationY);
146 TMatrixD rotationMatrix(3, 3);
147 rotation.GetRotationMatrix(rotationMatrix);
148 TMatrixDSym vertexCovariance2 = vertexCovariance.Similarity(rotationMatrix);
149 beam.sigma_ip_x(sqrt(vertexCovariance2[0][0]));
150 beam.sigma_ip_y(sqrt(vertexCovariance2[1][1]));
151 beam.sigma_ip_z(sqrt(vertexCovariance2[2][2]));
152 m_BelleFile->write(BBS_BEGIN_RUN, 0);
153}
154
156 const MCParticle* particle, MCParticleGraph::GraphParticle* mother)
157{
159 part = *particle;
160 if (mother != nullptr)
161 part.comesFrom(*mother);
162 int pdg = abs(particle->getPDG());
163 if ((pdg == Const::muon.getPDGCode()) ||
164 (pdg == Const::pion.getPDGCode()) ||
165 (pdg == Const::kaon.getPDGCode()) ||
166 ((pdg == Const::Kshort.getPDGCode()) && !m_DecayKsInGenerator) ||
167 (pdg == Const::Klong.getPDGCode()) ||
168 (pdg == Const::neutron.getPDGCode()) ||
169 ((pdg == Const::Lambda.getPDGCode()) && !m_DecayLambdaInGenerator) ||
170 (pdg == 3222) || // Sigma+
171 (pdg == 3112) || // Sigma-
172 (pdg == 3322) || // Xi0
173 (pdg == 3312) || // Xi-
174 (pdg == 3334)) { // Omega-
176 return;
177 }
178 std::vector<MCParticle*> daughters = particle->getDaughters();
179 for (const MCParticle* daughter : daughters)
180 addParticle(daughter, &part);
181}
182
184{
185 Belle::Gen_hepevt_Manager& hepevtManager =
186 Belle::Gen_hepevt_Manager::get_manager();
187 hepevtManager.remove();
188 /*
189 * The time shift applied by basf module "evtgen" (file beam.cc) is
190 * vertex z coordinate [mm] / (2.0 * 2.99792458).
191 * The vertex coordinate is calculated relative to the IP position.
192 * The position correction happens at the simulation stage (bpsmear) in basf.
193 */
194 double timeShift = (m_MCInitialParticles->getVertex().Z() -
195 m_BeamParameters->getVertex().Z()) /
196 Unit::mm / (2.0 * 0.1 * Const::speedOfLight);
197 /*
198 * Regeneration of MCParticle array. It is necessary because in basf the
199 * long-lived particles (K_S0, K_L0, Lambda, neutron, pi, K, mu)
200 * are decayed by GEANT3.
201 */
202 m_MCParticleGraph.clear();
203 for (const MCParticle& particle : m_MCParticles) {
204 if (particle.getMother() == nullptr)
205 addParticle(&particle, nullptr);
206 }
210 for (const MCParticle& particle : m_MCParticles) {
211 Belle::Gen_hepevt& hepevt = hepevtManager.add();
212 if (particle.hasStatus(MCParticle::c_Initial))
213 hepevt.isthep(3);
214 else if (particle.hasStatus(MCParticle::c_StableInGenerator))
215 hepevt.isthep(1);
216 else
217 hepevt.isthep(2);
218 hepevt.idhep(particle.getPDG());
219 hepevt.reset_mother();
220 const MCParticle* mother = particle.getMother();
221 int motherIndex = 0;
222 if (mother != nullptr)
223 motherIndex = mother->getIndex();
224 hepevt.moFirst(motherIndex);
225 hepevt.moLast(motherIndex);
226 hepevt.daFirst(particle.getFirstDaughter());
227 hepevt.daLast(particle.getLastDaughter());
228 ROOT::Math::PxPyPzEVector momentum = particle.get4Vector();
229 hepevt.PX(momentum.Px());
230 hepevt.PY(momentum.Py());
231 hepevt.PZ(momentum.Pz());
232 hepevt.E(momentum.E());
233 hepevt.M(particle.getMass());
234 ROOT::Math::XYZVector vertex = particle.getVertex();
235 hepevt.VX(vertex.X() / Unit::mm);
236 hepevt.VY(vertex.Y() / Unit::mm);
237 hepevt.VZ(vertex.Z() / Unit::mm);
238 hepevt.T(particle.getProductionTime() / Unit::mm * Const::speedOfLight +
239 timeShift);
240 }
241 m_BelleFile->write(BBS_EVENT, m_EventMetaData->getEvent());
242}
243
247
249{
250 delete m_BelleFile;
251}
void addParticle(const MCParticle *particle, MCParticleGraph::GraphParticle *mother)
Add particle to graph.
StoreObjPtr< MCInitialParticles > m_MCInitialParticles
Initial particles.
virtual void initialize() override
Initializer.
bool m_DecayKsInGenerator
Decay K_S0 in generator.
virtual void event() override
This method is called for each event.
Belle::Panther_FileIO * m_BelleFile
Belle file input-output handler.
virtual void endRun() override
This method is called if the current run ends.
virtual void terminate() override
This method is called at the end of the event processing.
virtual void beginRun() override
Called when entering a new run.
DBObjPtr< BeamParameters > m_BeamParameters
Beam parameters.
virtual ~BelleMCOutputModule()
Destructor.
bool m_DecayLambdaInGenerator
Decay Lambda in generator.
std::string m_OutputFileName
Output file name.
StoreObjPtr< EventMetaData > m_EventMetaData
Event metadata.
StoreArray< MCParticle > m_MCParticles
MC particles.
MCParticleGraph m_MCParticleGraph
MC particle graph.
static const ParticleType neutron
neutron particle
Definition Const.h:675
static const ParticleType Lambda
Lambda particle.
Definition Const.h:679
static const ChargedStable muon
muon particle
Definition Const.h:660
static const ChargedStable pion
charged pion particle
Definition Const.h:661
static const ParticleType Klong
K^0_L particle.
Definition Const.h:678
static const double speedOfLight
[cm/ns]
Definition Const.h:695
static const ParticleType Kshort
K^0_S particle.
Definition Const.h:677
static const ChargedStable kaon
charged kaon particle
Definition Const.h:662
Class to represent Particle data in graph.
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
@ c_checkCyclic
Check for cyclic dependencies.
@ c_clearParticles
Clear the particle list before adding the graph.
@ c_setDecayInfo
Set decay time and vertex.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition MCParticle.h:57
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition MCParticle.h:49
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition MCParticle.h:219
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition MCParticle.h:342
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_Output
This module is an output module (writes data).
Definition Module.h:79
static const double mm
[millimeters]
Definition Unit.h:70
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
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.