Belle II Software light-2406-ragdoll
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
44{
45}
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 necesary because in basf the
199 * long-lived particles (K_S0, K_L0, Lambda, neutron, pi, K, mu)
200 * are decayed by GEANT3.
201 */
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
245{
246}
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.
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
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:230
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition: MCParticle.h:353
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
Base class for Modules.
Definition: Module.h:72
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
@ 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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:600
void clear()
Reset particles and decay information to make the class reusable.
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24