Belle II Software  release-05-01-25
AAFHInterface.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - 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/logging/Logger.h>
12 #include <framework/gearbox/Unit.h>
13 #include <generators/aafh/AAFHInterface.h>
14 #include <TRandom3.h>
15 #include <TDatabasePDG.h>
16 #include <algorithm>
17 #include <array>
18 #include <iomanip>
19 
20 using namespace Belle2;
21 
22 extern "C" {
24  double aafhdiag36_rndm_(int*)
25  {
26  return gRandom->Rndm();
27  }
29  double aafhdiag36_rnf100_(int*)
30  {
31  return gRandom->Rndm();
32  }
34  void aafhdiag36_esft_maxed_(double* weight)
35  {
36  B2ERROR("AAFH: Maximum weight to small, increase maxFinalWeight to at least "
37  << *weight);
38  }
40  void aafhdiag36_eswe_maxed_(double* weight)
41  {
42  B2ERROR("AAFH: Maximum weight to small, increase maxSubGeneratorWeight to at least "
43  << *weight);
44  }
45 
47  extern void aafhdiag36_init_();
49  extern void aafhdiag36_event_(int*);
51  extern void aafhdiag36_finish_();
52 
55  extern struct {
56  double eb;
57  double eswe;
58  double esft;
59  double wap[4];
60  double wbp[4];
61  double vap[4];
62  } input_;
63 
66  extern struct {
67  double w2min;
68  double w2max;
69  } bound_;
70 
72  extern struct {
73  int iproc;
74  int irejec;
75  int idump[4];
76  int info;
77  int itot;
78  } input2_;
79 
81  extern struct {
82  double face;
83  double facl;
84  double facm;
85  double proc;
86  } factor_;
87 
91  extern struct {
92  double q1[5];
93  double q2[5];
94  double q3[5];
95  double q4[5];
96  double q5[5];
97  double q6[5];
98  } momenz_;
99 
102  extern struct {
103  double crosec;
104  double ercros;
105  } output_;
106 
108  extern struct {
109  double xm;
110  double xmu;
111  double xml;
112  double xm2;
113  double xmu2;
114  double xml2;
115  } masses_;
116 
118  extern struct {
119  double qcharg;
120  double qchrg2;
121  double qchrg3;
122  double qchrg4;
123  } charge_;
124 
126  extern struct {
127  double swe[4];
128  double swek[4];
129  double mwe[4];
130  double sum;
131  double sumk;
132  double maxwe;
133  int iwe[4];
134  int igen;
135  } westat_;
136 
138  extern struct {
139  double sumft;
140  double sumft2;
141  double ftmax;
142  int ieen;
143  } ftstat_;
144 };
145 
146 
148 {
149  input2_.itot = tries;
150 }
151 
152 void AAFHInterface::setGeneratorWeights(const std::vector<double>& weights)
153 {
154  if (weights.size() > 8) {
155  B2WARNING("AAFH: more than 8 generator weights supplied, "
156  "ignoring the extra weights");
157  } else if (weights.size() > 4 && weights.size() < 8) {
158  B2WARNING("AAFH: more than 4 but less than 8 generator weights supplied, "
159  "ignoring everything afer the first four");
160  } else if (weights.size() < 4) {
161  B2ERROR("AAFH: not enough generator weights, need exactly four or eight values");
162  }
163  std::copy_n(weights.begin(), 4, input_.wap);
164  if (weights.size() >= 8) {
165  std::copy_n(weights.begin() + 4, 4, input_.wbp);
166  }
167 }
168 
169 void AAFHInterface::setMaxWeights(double subgeneratorWeight, double finalWeight)
170 {
171  input_.eswe = subgeneratorWeight;
172  input_.esft = finalWeight;
173 }
174 
175 void AAFHInterface::setSupressionLimits(std::vector<double> limits)
176 {
177  if (limits.size() == 1) {
178  B2WARNING("AAFH: Only one suppression limit supplied, using same value for all limits");
179  limits.resize(4, limits[0]);
180  } else if (limits.size() != 4) {
181  B2ERROR("AAFH: Wrong number of suppression limits: supply either one or 4 values");
182  }
183  factor_.face = limits[0];
184  factor_.facm = limits[1];
185  factor_.facl = limits[2];
186  factor_.proc = limits[3];
187 }
188 
190 {
191  return input2_.itot;
192 }
193 
194 std::vector<double> AAFHInterface::getGeneratorWeights() const
195 {
196  std::vector<double> weights(8);
197  std::copy_n(input_.wap, 4, weights.begin());
198  std::copy_n(input_.wbp, 4, weights.begin() + 4);
199  return weights;
200 }
201 
203 {
204  return input_.eswe;
205 }
206 
208 {
209  return input_.esft;
210 }
211 
212 std::vector<double> AAFHInterface::getSuppressionLimits() const
213 {
214  return {factor_.face, factor_.facm, factor_.facl, factor_.proc};
215 }
216 
218 {
219  return output_.crosec * Unit::nb;
220 }
221 
223 {
224  return output_.ercros * Unit::nb;
225 }
226 
227 void AAFHInterface::initialize(double beamEnergy, EMode mode, ERejection rejection)
228 {
229  input_.eb = beamEnergy;
230  //minimum invariant mass is in units of beam energy and needs to be squared
231  bound_.w2min = (m_minMass / beamEnergy) * (m_minMass / beamEnergy);
232  //set mode and rejection scheme
233  input2_.iproc = mode;
234  input2_.irejec = rejection;
235 
236  //initialize the L
237  auto* part = TDatabasePDG::Instance()->GetParticle(m_particle.c_str());
238  if (!part) {
239  B2ERROR("AAFH: Particle '" << m_particle << "' does not exist");
240  return;
241  }
242  //set the mass in units of beam energy
243  masses_.xml = part->Mass() / beamEnergy;
244  //set the charge but with opposite sign
245  charge_.qcharg = -part->Charge() / 3.0;
246  //and remember pdg code
247  m_particlePDG = part->PdgCode();
248 
249  aafhdiag36_init_();
250 }
251 
252 void AAFHInterface::updateParticle(MCParticleGraph::GraphParticle& p, double* q, int pdg, bool isInitial)
253 {
254  // all values in q are in units of beam energy
255  // z-direction inversed since AAFH uses different convention
256  TLorentzVector vec(-q[0]*input_.eb, -q[1]*input_.eb, -q[2]*input_.eb, q[3]*input_.eb);
257  p.set4Vector(vec);
258 
259  // mass is multiplied by charge sign so take the absolute
260  p.setMass(fabs(q[4])*input_.eb);
261  p.setPDG(pdg);
263 
264  if (isInitial) {
265  p.addStatus(MCParticle::c_Initial);
266  }
267 }
268 
270 {
271  int status {0};
272  aafhdiag36_event_(&status);
273  if (status != 0) {
274  B2ERROR("Failed to generate event after " << status << " tries, giving up");
275  return;
276  }
277 
278  //initial beam particles
279  auto& p1 = mpg.addParticle();
280  auto& p2 = mpg.addParticle();
281 
282  //generated particles
283  auto& p3 = mpg.addParticle();
284  auto& p4 = mpg.addParticle();
285  auto& p5 = mpg.addParticle();
286  auto& p6 = mpg.addParticle();
287 
288  //array of generated particle pdg codes depends on generated mode
289  std::array<int, 4> pdg {{0}};
290  switch (input2_.iproc) {
291  case c_MuonParticle:
292  pdg = {{m_particlePDG, -m_particlePDG, 13, -13}};
293  break;
294  case c_MuonMuon:
295  pdg = {{13, -13, 13, -13}};
296  break;
297  case c_ElectonMuon:
298  pdg = {{11, -11, 13, -13}};
299  break;
300  case c_ElectronParticle:
301  pdg = {{11, -11, m_particlePDG, -m_particlePDG}};
302  break;
303  case c_ElectonElectron:
304  pdg = {{11, -11, 11, -11}};
305  break;
306  }
307 
308  updateParticle(p1, momenz_.q1, 11, true);
309  updateParticle(p2, momenz_.q2, -11, true);
310  updateParticle(p3, momenz_.q3, pdg[0]);
311  updateParticle(p4, momenz_.q4, pdg[1]);
312  updateParticle(p5, momenz_.q5, pdg[2]);
313  updateParticle(p6, momenz_.q6, pdg[3]);
314 }
315 
317 {
318  aafhdiag36_finish_();
319  // The following code just prints the results. If B2RESULT has been compiled
320  // out this is unneccesary and cppcheck complains so we enclose this in an
321  // ifdef to only be executed if B2RESULT might actually be shown
322 #ifndef LOG_NO_B2RESULT
323  B2RESULT("AAFH Generation finished.");
324  B2RESULT("Final cross section: "
325  << std::setprecision(3) << output_.crosec << "+-"
326  << std::setprecision(3) << output_.ercros << " nb");
327  B2RESULT("Minimum invariant mass: "
328  << sqrt(bound_.w2min)*input_.eb << " GeV");
329  B2RESULT("Maximum invariant mass: "
330  << sqrt(bound_.w2max)*input_.eb << " GeV");
331  B2RESULT("Maximum subgenerator weight: " << westat_.maxwe << " ("
332  << westat_.mwe[0] << ", "
333  << westat_.mwe[1] << ", "
334  << westat_.mwe[2] << ", "
335  << westat_.mwe[3] << ")");
336  B2RESULT("Maximum final weight: " << ftstat_.ftmax);
337  //The fortran source says that it's advisable if all subgenerators have the
338  //same probability to account for all peaks in the differental cross section
339  //and that it runs with the highest efficiency when the maximum weights in
340  //all sub generators are equal. So let's calculate appropriate subgenerator
341  //weights by looking how many events were generated per sub generator or the
342  //maximum weights
343  std::array<double, 4> idealWAP;
344  std::array<double, 4> idealWBP;
345  //find the first non-zero weight
346  int reference = 0;
347  for (int i = 0; i < 4; ++i) {
348  if (input_.wap[i] != 0) {
349  reference = i;
350  break;
351  }
352  }
353  //and use that es reference
354  const double w0 = westat_.mwe[reference] * input_.wbp[reference];
355  const double n0 = westat_.iwe[reference] / input_.wap[reference];
356  //now calculate the relative weights we need to get
357  //1) same maximum weight which is inversely proportional to to wbp
358  //2) same number of generated events which is proportional to wap
359  for (int i = 0; i < 4; ++i) {
360  if (input_.wap[i] == 0) {
361  idealWBP[i] = input_.wbp[i];
362  idealWAP[i] = 0;
363  } else {
364  idealWBP[i] = (westat_.mwe[i] * input_.wbp[i]) / w0;
365  idealWAP[i] = n0 * input_.wap[i] / westat_.iwe[i];
366  //and normalize to 1.0
367  idealWAP[i] /= idealWAP[reference];
368  }
369  }
370 
371  B2RESULT("To get the same maximum event weight and chance for each sub "
372  "generator use these parameters:");
373  B2RESULT(" 'maxFinalWeight': " << std::setprecision(3)
374  << ftstat_.ftmax * 1.1 << ",");
375  B2RESULT(" 'maxSubgeneratorWeight': " << std::setprecision(3)
376  << westat_.maxwe * 1.1 << ",");
377  B2RESULT(" 'subgeneratorWeights': [" << std::scientific
378  << std::setprecision(3) << idealWAP[0] << ", "
379  << std::setprecision(3) << idealWAP[1] << ", "
380  << std::setprecision(3) << idealWAP[2] << ", "
381  << std::setprecision(3) << idealWAP[3] << ", ");
382  B2RESULT(" " << std::scientific
383  << std::setprecision(3) << idealWBP[0] << ", "
384  << std::setprecision(3) << idealWBP[1] << ", "
385  << std::setprecision(3) << idealWBP[2] << ", "
386  << std::setprecision(3) << idealWBP[3] << "],");
387 #endif
388 }
Belle2::AAFHInterface::getTotalCrossSectionError
double getTotalCrossSectionError() const
Return error on the total cross section.
Definition: AAFHInterface.cc:222
Belle2::AAFHInterface::m_minMass
double m_minMass
minimum invariant mass
Definition: AAFHInterface.h:110
Belle2::AAFHInterface::c_ElectonMuon
@ c_ElectonMuon
E+E- -> E+E-Mu+Mu-.
Definition: AAFHInterface.h:40
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::Unit::nb
static const double nb
[nanobarn]
Definition: Unit.h:94
Belle2::AAFHInterface::getSuppressionLimits
std::vector< double > getSuppressionLimits() const
Get suppression limits.
Definition: AAFHInterface.cc:212
Belle2::AAFHInterface::updateParticle
void updateParticle(MCParticleGraph::GraphParticle &p, double *q, int pdg, bool isInitial=false)
update particle with generated values
Definition: AAFHInterface.cc:252
Belle2::AAFHInterface::setSupressionLimits
void setSupressionLimits(std::vector< double > limits)
Set the suppression limits when calculation the matrix elements.
Definition: AAFHInterface.cc:175
Belle2::AAFHInterface::m_particle
std::string m_particle
name of the particle for modes c_MuonParticle and c_ElectronParticle
Definition: AAFHInterface.h:112
Belle2::AAFHInterface::setMaxTries
void setMaxTries(int tries)
Set the maximum number of tries to generate an event.
Definition: AAFHInterface.cc:147
Belle2::AAFHInterface::finish
void finish()
calculate total cross section
Definition: AAFHInterface.cc:316
Belle2::AAFHInterface::getTotalCrossSection
double getTotalCrossSection() const
Return total cross section.
Definition: AAFHInterface.cc:217
Belle2::AAFHInterface::c_ElectonElectron
@ c_ElectonElectron
E+E- -> E+E-E+E-.
Definition: AAFHInterface.h:44
Belle2::AAFHInterface::getGeneratorWeights
std::vector< double > getGeneratorWeights() const
Get the relative weights for the sub generators.
Definition: AAFHInterface.cc:194
Belle2::AAFHInterface::c_MuonMuon
@ c_MuonMuon
E+E- -> Mu+Mu-Mu+Mu-.
Definition: AAFHInterface.h:38
Belle2::AAFHInterface::c_MuonParticle
@ c_MuonParticle
E+E- -> Mu+Mu-L+L- where L can be anything, defaults to tau.
Definition: AAFHInterface.h:36
Belle2::AAFHInterface::setGeneratorWeights
void setGeneratorWeights(const std::vector< double > &weights)
Set the relative weights for the sub generators.
Definition: AAFHInterface.cc:152
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::AAFHInterface::initialize
void initialize(double beamEnergy, EMode mode, ERejection rejection)
initialize the generator
Definition: AAFHInterface.cc:227
Belle2::AAFHInterface::m_particlePDG
int m_particlePDG
pdg of the particle for modes c_MuonParticle and c_ElectronParticle, gets set in initialize()
Definition: AAFHInterface.h:115
Belle2::AAFHInterface::ERejection
ERejection
Rejection mode.
Definition: AAFHInterface.h:48
Belle2::MCParticle::c_Initial
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition: MCParticle.h:68
Belle2::AAFHInterface::EMode
EMode
Generation mode.
Definition: AAFHInterface.h:34
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::AAFHInterface::getMaxFinalWeight
double getMaxFinalWeight() const
Get the maximum expected subgenerator weight for the rejection method.
Definition: AAFHInterface.cc:207
Belle2::AAFHInterface::getMaxTries
int getMaxTries() const
Get the maximum number of tries to generate an event.
Definition: AAFHInterface.cc:189
Belle2::AAFHInterface::c_ElectronParticle
@ c_ElectronParticle
E+E- -> E+E-L+L- where L can be anything, defaults to tau.
Definition: AAFHInterface.h:42
Belle2::MCParticle::c_StableInGenerator
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:60
Belle2::AAFHInterface::generateEvent
void generateEvent(MCParticleGraph &mpg)
generate one event and add it to the graph in CMS
Definition: AAFHInterface.cc:269
Belle2::AAFHInterface::getMaxSubGeneratorWeight
double getMaxSubGeneratorWeight() const
Get the maximum expected final weight for the rejection method.
Definition: AAFHInterface.cc:202
Belle2::AAFHInterface::setMaxWeights
void setMaxWeights(double subgeneratorWeight, double finalWeight)
Set the maximum expected weights for the rejection method.
Definition: AAFHInterface.cc:169
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58