Belle II Software  release-05-01-25
initialparticles.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/datastore/DataStore.h>
12 #include <framework/database/DBStore.h>
13 #include <framework/dbobjects/BeamParameters.h>
14 #include <framework/utilities/CalcMeanCov.h>
15 #include <generators/utilities/InitialParticleGeneration.h>
16 
17 #include <gtest/gtest.h>
18 
19 using namespace std;
20 using namespace Belle2;
21 
22 namespace {
23 
26  class InitialParticleGenerationTests : public ::testing::Test {
27  protected:
28 
30  virtual void SetUp()
31  {
32  DataStore::Instance().setInitializeActive(true);
33  generator.initialize();
34  DataStore::Instance().setInitializeActive(false);
35 
36  // have some useful defaults for the beam parameters
37  beamparams.setHER(7.004, 0.0415, {2.63169e-05, 1e-5, 1e-5});
38  beamparams.setLER(4.002, -0.0415, {5.64063e-06, 1e-5, 1e-5});
39  beamparams.setVertex({0, 1, 2}, {4.10916e-07, 0, -2.64802e-06, 0, 1.7405e-11, 0, -2.64802e-06, 0, 0.000237962});
40  }
41 
43  virtual void TearDown()
44  {
45  DataStore::Instance().reset();
46  DBStore::Instance().reset();
47  }
48 
49  InitialParticleGeneration generator{BeamParameters::c_smearALL};
50  BeamParameters beamparams;
51  };
52 
54  TEST_F(InitialParticleGenerationTests, TestCMSGeneration)
55  {
56  beamparams.setGenerationFlags(BeamParameters::c_generateCMS);
57  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
58  auto initialCMS = generator.generate();
59  // transformation should be identity
60  EXPECT_EQ(initialCMS.getLabToCMS(), TLorentzRotation());
61  // her should be along z
62  EXPECT_NEAR(initialCMS.getHER().Theta(), 0, 1e-15);
63  // ler should be opposite z
64  EXPECT_NEAR(initialCMS.getLER().Theta(), M_PI, 1e-15);
65  // both should have the same energy
66  EXPECT_NEAR(initialCMS.getHER().E(), initialCMS.getLER().E(), 1e-15);
67  // and the invariant mass should of course be correct
68  EXPECT_EQ(initialCMS.getMass(), beamparams.getMass());
69 
70  // and now compare it to generation in Lab
71  beamparams.setGenerationFlags(0);
72  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
73  auto initialLAB = generator.generate();
74  // same invariant mass
75  EXPECT_EQ(initialLAB.getMass(), beamparams.getMass());
76  // no smearing so LER and HER need to be identical to beam parameters
77  EXPECT_EQ(initialLAB.getHER(), beamparams.getHER());
78  EXPECT_EQ(initialLAB.getLER(), beamparams.getLER());
79  // so check that the CMS beams are identical to the transformed lab ones
80  EXPECT_EQ(initialCMS.getHER(), initialLAB.getLabToCMS() * initialLAB.getHER());
81  EXPECT_EQ(initialCMS.getLER(), initialLAB.getLabToCMS() * initialLAB.getLER());
82  }
83 
85  TEST_F(InitialParticleGenerationTests, TestEnergySmear)
86  {
87  // two beams completely head on, one with energy uncertainty so invariant
88  // mass uncertainty should be exactly the same as single uncertainty. Here
89  // we give variance so we expect 0.1
90  beamparams.setHER(10, 0, {0.01});
91  beamparams.setLER(10, M_PI, {0});
92  beamparams.setGenerationFlags(BeamParameters::c_smearBeamEnergy);
93  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
94 
95  // so let's generate events with this and store all invariant masses
96  {
97  CalcMeanCov<1> mean;
98  for (int i = 0; i < 100000; ++i) {
99  auto& initial = generator.generate();
100  ASSERT_TRUE(initial.hasGenerationFlags(BeamParameters::c_smearBeamEnergy));
101  mean.add(initial.getMass());
102  }
103  // and compare with expectation
104  EXPECT_NEAR(mean.getStddev(), 0.1, 0.0005);
105  EXPECT_NEAR(mean.getMean(), 20, 0.001);
106  }
107 
108  // uncertainty on both so result should be sqrt(a^2 + b^2) = sqrt(2) * 0.1
109  beamparams.setLER(10, M_PI, {0.01});
110  beamparams.setGenerationFlags(BeamParameters::c_smearBeamEnergy);
111  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
112  {
113  CalcMeanCov<1> mean;
114  for (int i = 0; i < 100000; ++i) {
115  auto& initial = generator.generate();
116  ASSERT_TRUE(initial.hasGenerationFlags(BeamParameters::c_smearBeamEnergy));
117  mean.add(initial.getMass());
118  }
119  EXPECT_NEAR(mean.getStddev(), std::sqrt(2) * 0.1, 0.0005);
120  EXPECT_NEAR(mean.getMean(), 20, 0.001);
121  }
122  }
123 
124 
128  TEST_F(InitialParticleGenerationTests, TestVertexSmear)
129  {
130  beamparams.setGenerationFlags(BeamParameters::c_smearVertex);
131  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
132 
133  CalcMeanCov<3> mean;
134  for (int i = 0; i < 100000; ++i) {
135  auto& initial = generator.generate();
136  mean.add(initial.getVertex().X(), initial.getVertex().Y(), initial.getVertex().Z());
137  }
138  auto cov = beamparams.getCovVertex();
139  auto pos = beamparams.getVertex();
140  for (int i = 0; i < 3; ++i) {
141  EXPECT_NEAR(mean.getMean(i), pos(i), 1e-4);
142  for (int j = 0; j < 3; ++j) {
143  EXPECT_NEAR(mean.getCovariance(i, j), cov(i, j), 1e-6);
144  }
145  }
146  }
147 
156  TEST_F(InitialParticleGenerationTests, TestFlags)
157  {
158  // so loop over all flags
159  for (int flag = 0; flag <= BeamParameters::c_smearALL; ++flag) {
160  // set the flag and overwrite dbstore
161  beamparams.setGenerationFlags(flag);
162  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
163  // rememeber last event and set it to the settings for initialization
164  MCInitialParticles last = beamparams;
165  // no generate a few events and check everything
166  for (int i = 0; i < 5; ++i) {
167  auto& initial = generator.generate();
168  EXPECT_EQ(flag, initial.getGenerationFlags());
169  // create text representation of flags
170  const std::string flags = initial.getGenerationFlagString();
171  if (flag & BeamParameters::c_smearBeam) {
172  // smearing of beam is active, check that the beams are actually
173  // different to the settings and the previous event
174  EXPECT_NE(initial.getHER(), beamparams.getHER()) << flags << " " << i;
175  EXPECT_NE(initial.getLER(), beamparams.getLER()) << flags << " " << i;
176  EXPECT_NE(initial.getHER(), last.getHER()) << flags << " " << i;
177  EXPECT_NE(initial.getLER(), last.getLER()) << flags << " " << i;
178  } else if (!(flag & BeamParameters::c_generateCMS)) {
179  // no smearing no cms generation, so everything should stay fixed.
180  EXPECT_EQ(initial.getHER(), beamparams.getHER()) << flags << " " << i;
181  EXPECT_EQ(initial.getLER(), beamparams.getLER()) << flags << " " << i;
182  EXPECT_EQ(initial.getHER(), last.getHER()) << flags << " " << i;
183  EXPECT_EQ(initial.getLER(), last.getLER()) << flags << " " << i;
184  } else {
185  // we want to compare CMS initial particles to lab beam energies
186  // but transformation is identity if c_generateCMS is set. Reset
187  // it so that the can boost beams to CMS for comparison
188  beamparams.setGenerationFlags(0);
189  EXPECT_EQ(initial.getHER(), beamparams.getLabToCMS() * beamparams.getHER()) << flags << " " << i;
190  EXPECT_EQ(initial.getLER(), beamparams.getLabToCMS() * beamparams.getLER()) << flags << " " << i;
191  // comparison to beamparameters fails in generateCMS so we just
192  // skip the first round here
193  if (i > 0) {
194  EXPECT_EQ(initial.getHER(), last.getHER()) << flags << " " << i;
195  EXPECT_EQ(initial.getLER(), last.getLER()) << flags << " " << i;
196  }
197  }
198  if (flag & BeamParameters::c_smearVertex) {
199  //smearing of the vertex is enabled, make sure the vertex changes
200  EXPECT_NE(initial.getVertex(), beamparams.getVertex());
201  EXPECT_NE(initial.getVertex(), last.getVertex());
202  } else {
203  //otherwise make sure it doesn't change.
204  EXPECT_EQ(initial.getVertex(), beamparams.getVertex());
205  EXPECT_EQ(initial.getVertex(), last.getVertex());
206  }
207  last = initial;
208  }
209  }
210  }
211 
213  TEST_F(InitialParticleGenerationTests, TestValidFlag)
214  {
215  beamparams.setGenerationFlags(0);
216  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
217 
218  MCInitialParticles& initial = generator.generate();
219  TLorentzVector her = initial.getHER();
220  TLorentzVector ler = initial.getLER();
221  TVector3 vertex = initial.getVertex();
222  for (int i = 0; i < 10; ++i) {
223  initial = generator.generate();
224  EXPECT_EQ(her, initial.getHER());
225  EXPECT_EQ(ler, initial.getLER());
226  EXPECT_EQ(vertex, initial.getVertex());
227  }
228  beamparams.setGenerationFlags(BeamParameters::c_smearALL);
229  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
230  for (int i = 0; i < 10; ++i) {
231  initial = generator.generate();
232  EXPECT_NE(her, initial.getHER());
233  EXPECT_NE(ler, initial.getLER());
234  EXPECT_NE(vertex, initial.getVertex());
235  her = initial.getHER();
236  ler = initial.getLER();
237  vertex = initial.getVertex();
238  }
239  }
240 
241  TEST_F(InitialParticleGenerationTests, UpdateVertex)
242  {
243  beamparams.setGenerationFlags(BeamParameters::c_smearBeam);
244  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
245  // first run but no smearing allowed, should return the nominal vertex
246  TVector3 shift = generator.updateVertex();
247  EXPECT_EQ(shift, TVector3(0, 1, 2));
248  // create a new initial particle. Particle exists now, no smearing allowed so no change in shift
249  const MCInitialParticles& initial = generator.generate();
250  auto nominal = initial.getVertex();
251  shift = generator.updateVertex();
252  EXPECT_EQ(shift, TVector3(0, 0, 0));
253  // ok, allow smearing, now we expect shift
254  beamparams.setGenerationFlags(BeamParameters::c_smearALL);
255  DBStore::Instance().addConstantOverride("BeamParameters", new BeamParameters(beamparams));
256  shift = generator.updateVertex();
257  EXPECT_NE(shift, TVector3(0, 0, 0));
258  EXPECT_EQ(nominal + shift, initial.getVertex());
259  // but running again should not shift again
260  shift = generator.updateVertex();
261  EXPECT_EQ(shift, TVector3(0, 0, 0));
262  // unless we force regeneration
263  auto previous = initial.getVertex();
264  shift = generator.updateVertex(true);
265  EXPECT_NE(shift, TVector3(0, 0, 0));
266  EXPECT_EQ(previous + shift, initial.getVertex());
267  }
268 }
Belle2::MCInitialParticles::getLER
const TLorentzVector & getLER() const
Get 4vector of the low energy beam.
Definition: MCInitialParticles.h:140
prepareAsicCrosstalkSimDB.e
e
aux.
Definition: prepareAsicCrosstalkSimDB.py:53
Belle2::MCInitialParticles::getHER
const TLorentzVector & getHER() const
Get 4vector of the high energy beam.
Definition: MCInitialParticles.h:137
Belle2::InitialParticleGeneration
Generate Collision.
Definition: InitialParticleGeneration.h:35
Belle2::MCInitialParticles
This class contains the initial state for the given event.
Definition: MCInitialParticles.h:35
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TEST_F
TEST_F(GlobalLabelTest, LargeNumberOfTimeDependentParameters)
Test large number of time-dep params for registration and retrieval.
Definition: globalLabel.cc:65
Belle2::MCInitialParticles::getVertex
const TVector3 & getVertex() const
Get the position of the collision.
Definition: MCInitialParticles.h:143
Belle2::CalcMeanCov
Class to calculate mean and and covariance between a number of parameters on running data without sto...
Definition: CalcMeanCov.h:45
Belle2::BeamParameters
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
Definition: BeamParameters.h:33