Belle II Software  release-05-01-25
BabayagaNLO.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Torben Ferber *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/babayaganlo/BabayagaNLO.h>
12 #include <framework/gearbox/Unit.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/dataobjects/EventMetaData.h>
15 #include <analysis/dataobjects/EventExtraInfo.h>
16 
17 #include <TDatabasePDG.h>
18 #include <TLorentzVector.h>
19 #include <TRandom3.h>
20 
21 #include <generators/utilities/InitialParticleGeneration.h>
22 
23 using namespace std;
24 using namespace Belle2;
25 
26 extern "C" {
27 
28  extern struct {
29  double bweight;
30  double bp1[4];
31  double bq1[4];
32  double bp2[4];
33  double bq2[4];
34  double bphot[4][100];
35  int bnphot;
36  } momset_;
37 
38  extern struct {
39  double bprescale;
40  } prescale_;
41 
42  //results for cross section and photon multiplicity
43  extern struct {
44  double rescross;
45  double rescrosserr;
46  double rescrossphot[40];
47  double rescrossphoterr[40];
48  double rescrossphotfrac[40];
49  double resnphmax;
50  } bresults_;
51 
52  //results and statistics for hit/miss
53  extern struct {
54  double hnmcross;
55  double hnmcrosserr;
56  double hnmeff;
57  double hnmsdifmax;
58  double hnmfmax;
59  double hnmmaxtriallimit;
60  double hnmmaxtrial;
61  } bhitnmiss_;
62 
63  //bias for hit/miss
64  extern struct {
65  double biashit;
66  double biashitpmiss;
67  double biasneghit;
68  double biasneghitmiss;
69  double sezover;
70  double errsezover;
71  double sezneg;
72  double errsezneg;
73  double nover;
74  double nneg;
75  } bbias_;
76 
77  //monitoring
78  extern struct {
79  double monsdif;
80  } bmonitoring_;
81 
83  double babayaganlo_rndm_(int*)
84  {
85  double r = gRandom->Rndm();
86  return r;
87  }
88 
90  double babayaganlo_getrandomcmsenergy_()
91  {
92  B2FATAL("babayaganlo_getrandomcmsenergy() is not implmented");
93  }
94 
96 // void main_belle2_(int* mode, double* xpar, int* npar, double *ecm);
97  void main_belle2_(int* mode, double* ecm, double* xpar, int* npar);
98 
100  void babayaga_setvpolnsk_(const char* vpolnskfilename, size_t* length);
101 
103  void babayaganlo_warning_overweight_(double* weight, double* max)
104  {
105  B2WARNING("Babayaga.NLO: Maximum weight " << *max << " to small, increase fmax to at least " << *weight);
106  }
107 
109  void babayaganlo_error_rejection_(double* ratio)
110  {
111  B2ERROR("Babayaga.NLO: Event rejected! Ratio of cross section error increase too large: " << *ratio);
112  }
113 
115  void babayaganlo_error_negative_(double* weight)
116  {
117  B2ERROR("Babayaga.NLO: Event has negative weight: " << *weight);
118  }
119 
121  void babayaganlo_result_nominalcmsenergy_(double* energy)
122  {
123  B2RESULT("Babayaga.NLO: Nominal CMS energy (GeV): " << *energy);
124  }
125 
126  void babayaganlo_result_weightedxsec_(double* xsec, double* xsecerr)
127  {
128  B2RESULT("Babayaga.NLO: Weighted cross section (nb): " << *xsec << " +/- " << *xsecerr);
129  }
130 
131  void babayaganlo_result_unweightedxsec_(double* xsec, double* xsecerr)
132  {
133  B2RESULT("Babayaga.NLO: Unweighted cross section (nb): " << *xsec << " +/- " << *xsecerr);
134  }
135 
136  void babayaganlo_result_unweightedxsec_biascorrected_(double* xsec, double* xsecerr)
137  {
138  B2RESULT("Babayaga.NLO: Unweighted cross section, bias corrected (nb): " << *xsec << " +/- " << *xsecerr);
139  }
140 
141  void babayaganlo_result_unweightedxsec_overweight_(double* xsec, double* xsecerr)
142  {
143  B2RESULT("Babayaga.NLO: Unweighted cross section overweight (nb): " << *xsec << " +/- " << *xsecerr);
144  }
145 
146  void babayaganlo_result_unweightedxsec_underweight_(double* xsec, double* xsecerr)
147  {
148  B2RESULT("Babayaga.NLO: Unweighted cross section underweight (nb): " << -*xsec << " +/- " << *xsecerr);
149  }
150 
151  void babayaganlo_result_hitormisseff_(double* eff)
152  {
153  B2RESULT("Babayaga.NLO: Hit or miss efficiency: " << *eff * 100. << " % ");
154  }
155 
156  void babayaganlo_result_nover_(int* nover)
157  {
158  B2RESULT("Babayaga.NLO: Points with w > fmax (bias): " << *nover);
159  }
160 
161  void babayaganlo_result_biashit_(double* biashit)
162  {
163  B2RESULT("Babayaga.NLO: Bias/hit: " << *biashit * 100. << " % ");
164  }
165 
166  void babayaganlo_result_biashitpmiss_(double* biashitpmiss)
167  {
168  B2RESULT("Babayaga.NLO: Bias/(hit+missed): " << *biashitpmiss * 100. << " % ");
169  }
170 
171  void babayaganlo_result_nneg_(int* nneg)
172  {
173  B2RESULT("Babayaga.NLO: Points with w > fmax (bias): " << *nneg);
174  }
175 
176  void babayaganlo_result_biasneghit_(double* biasneghit)
177  {
178  B2RESULT("Babayaga.NLO: Neg. bias/hit: " << *biasneghit * 100. << " % ");
179  }
180 
181  void babayaganlo_result_biasneghitmiss_(double* biasneghitmiss)
182  {
183  B2RESULT("Babayaga.NLO: Neg. bias/(hit+missed): " << *biasneghitmiss * 100. << " % ");
184  }
185 
186  void babayaganlo_result_maxweight_(double* sdifmax, double* fmax)
187  {
188  B2RESULT("Babayaga.NLO: Maximum weight (sdifmax): " << *sdifmax << ", user maximum weight(fmax): " << *fmax);
189  }
190 
191  void babayaganlo_result_vpmin_(double* xsec, double* xsecerr, double* frac)
192  {
193  B2RESULT("Babayaga.NLO: VP Uncertainty, minimum (nb): " << *xsec << " +/- " << *xsecerr << " (" << *frac << " %)");
194  }
195 
196  void babayaganlo_result_vpcentral_(double* xsec, double* xsecerr)
197  {
198  B2RESULT("Babayaga.NLO: VP Uncertainty, central (nb): " << *xsec << " +/- " << *xsecerr);
199  }
200 
201  void babayaganlo_result_vpmax_(double* xsec, double* xsecerr, double* frac)
202  {
203  B2RESULT("Babayaga.NLO: VP Uncertainty, maximum (nb): " << *xsec << " +/- " << *xsecerr << " ( " << *frac << " %)");
204  }
205 
206  void babayaganlo_finishedmaxsearch_(double* fmax)
207  {
208  B2INFO("Babayaga.NLO: Finished maximum search: " << *fmax << ", starting now unweighted generation");
209  }
210 
211  void babayaganlo_fatal_usercuts_()
212  {
213  B2FATAL("Babayaga.NLO: UserMode is selected but wrong user cuts are provided!");
214  }
215 
216  void babayaganlo_fatal_usercutsfs_()
217  {
218  B2FATAL("Babayaga.NLO: UserMode is only possible for ee final states!");
219  }
220 
221  void babayaganlo_fatal_usercutsprescale_()
222  {
223  B2FATAL("Babayaga.NLO: Prescale value must be larger than 0 (zero), hint: ~100 is reasonable!");
224  }
225 
226  void babayaganlo_fatal_weightedprescale_()
227  {
228  B2FATAL("Babayaga.NLO: Generator prescale is not allowed with usermode = weighted");
229  }
230 
231  void babayaganlo_result_weightsum_(double* wsum)
232  {
233  B2RESULT("Babayaga.NLO: Sum of weights: " << *wsum);
234  }
235 
236  void babayaganlo_result_intlum_(double* lum, double* lumerr)
237  {
238  B2RESULT("Babayaga.NLO: Luminosity equivalent (using calc. xsec) (fb-1): " << *lum << " +/- " << *lumerr);
239  }
240 
241  void babayaganlo_warning_prescaleweight_(double* weight)
242  {
243  B2WARNING("Babayaga.NLO: Prescale of less than one, increase prescale safety margin " << *weight);
244  }
245 
246  void babayaganlo_error_isnan1_(double* phsp, double* w)
247  {
248  B2ERROR("Babayaga.NLO: phsp (" << *phsp << ") or w (" << *w << ") are NAN, skipping event (nan1 type)");
249  }
250 
251  void babayaganlo_error_isnan2_(double* sdif)
252  {
253  B2ERROR("Babayaga.NLO: sdif (" << *sdif << ") is NAN, skipping event (nan2 type)");
254  }
255 
256 
257 
258 }
259 
260 BabayagaNLO::BabayagaNLO()
261 {
262  for (int i = 0; i < 100; ++i) {
263  m_npar[i] = 0;
264  m_xpar[i] = 0.0;
265  }
266 
267  m_sDif = 0.;
268 
269  setDefaultSettings();
270 }
271 
272 BabayagaNLO::~BabayagaNLO()
273 {
274 
275 }
276 
277 void BabayagaNLO::setDefaultSettings()
278 {
279 // m_cmsEnergy = -1.;
280  m_cmsEnergyNominal = -1.;
281 
282  m_applyBoost = true;
283 
284  m_finalState = "ee";
285  m_vacPol = "hlmnt";
286  m_order = "exp";
287  m_model = "matched";
288  m_mode = "unweighted";
289  m_userMode = "NONE";
290 
291  m_pi = 3.1415926535897932384626433832795029;
292  m_conversionFactor = 0.389379660e6;
293  m_alphaQED0 = 1.0 / 137.0359895;
294  m_massElectron = 0.51099906 * Unit::MeV;
295  m_massMuon = 105.65836900 * Unit::MeV;
296  m_massW = 80.385 * Unit::GeV;
297  m_massZ = 91.1882 * Unit::GeV;
298  m_widthZ = 2.4952 * Unit::GeV;
299  m_eMin = 0.1 * Unit::GeV;
300  m_epsilon = 5.e-7;
301  m_maxAcollinearity = 180.0;// * Unit::Deg;
302 
303  m_ScatteringAngleRange = make_pair(15.0, 165.0); //in [deg]
304  m_ScatteringAngleRangePhoton = make_pair(15.0, 165.0); //in [deg] (unused in original babayaga)
305 
306  m_nPhot = -1;
307  m_nSearchMax = 50000;
308  m_fMax = -1.;
309  m_EnergySpread = 5e-3;
310  m_VPUncertainty = false;
311  m_eemin = 0.0;
312  m_temin = 0.0;
313  m_egmin = 0.0;
314  m_tgmin = 0.0;
315  m_eeveto = 0.0;
316  m_teveto = 0.0;
317  m_egveto = 0.0;
318  m_tgveto = 0.0;
319  m_maxprescale = 1.0;
320 }
321 
322 void BabayagaNLO::initExtraInfo()
323 {
324  StoreObjPtr<EventExtraInfo> extrainfo;
325  extrainfo.registerInDataStore();
326 }
327 
328 void BabayagaNLO::init()
329 {
330  applySettings();
331 }
332 
333 
334 void BabayagaNLO::generateEvent(MCParticleGraph& mcGraph, double ecm, TVector3 vertex, TLorentzRotation boost)
335 {
336  //Generate event
337  int mode = 1;
338  main_belle2_(&mode, &ecm, m_xpar, m_npar);
339 
340  //Store the initial particles as virtual particles into the MCParticleGraph
341  storeParticle(mcGraph, momset_.bq1, 11, vertex, boost, true);
342  storeParticle(mcGraph, momset_.bp1, -11, vertex, boost, true);
343 
344  //Store the final state fermions or photons (for 'gg') as real particle into the MCParticleGraph
345  int pdg = 11;
346  int antipdg = -11;
347 
348  if (m_finalState == "gg") {
349  pdg = 22;
350  antipdg = 22;
351  } else if (m_finalState == "mm") {
352  pdg = 13;
353  antipdg = -13;
354  }
355 
356  storeParticle(mcGraph, momset_.bq2, pdg, vertex, boost, false, false, false);
357  storeParticle(mcGraph, momset_.bp2, antipdg, vertex, boost, false, false, false);
358 
359  //Store the real ISR and FSR photons into the MCParticleGraph
360  for (int iPhot = 0; iPhot < momset_.bnphot; ++iPhot) {
361  double photMom[4] = {momset_.bphot[0][iPhot], momset_.bphot[1][iPhot], momset_.bphot[2][iPhot], momset_.bphot[3][iPhot]};
362  storeParticle(mcGraph, photMom, 22, vertex, boost, false, false, true);
363  }
364 
365  //set event weight
366  StoreObjPtr<EventMetaData> eventMetaDataPtr("EventMetaData", DataStore::c_Event);
367  eventMetaDataPtr->setGeneratedWeight(momset_.bweight);
368 
369  //set event prescale
370  StoreObjPtr<EventExtraInfo> eventExtraInfo;
371  if (not eventExtraInfo.isValid())
372  eventExtraInfo.create();
373  if (eventExtraInfo->hasExtraInfo("GeneratedPrescale")) {
374  B2WARNING("EventExtraInfo with given name is already set! I won't set it again!");
375  } else {
376  float targetValue = prescale_.bprescale;
377  eventExtraInfo->addExtraInfo("GeneratedPrescale", targetValue);
378  }
379 
380 }
381 
382 
383 void BabayagaNLO::term()
384 {
385 
386  B2RESULT("Babayaga.NLO: Final state: " << m_finalState);
387  B2RESULT("Babayaga.NLO: Mode: " << m_mode);
388  B2RESULT("Babayaga.NLO: Order: " << m_order);
389  B2RESULT("Babayaga.NLO: Model: " << m_model);
390  B2RESULT("Babayaga.NLO: Vac. pol. (VP): " << m_vacPol);
391  B2RESULT("Babayaga.NLO: Usercuts: " << m_userMode);
392 
393  // all other results are displayed from fortran using a callback to C extern logging
394 
395  int mode = 2;
396  double ecm = -1.;
397  main_belle2_(&mode, &ecm, m_xpar, m_npar);
398 
399 }
400 
401 //=========================================================================
402 // Protected methods
403 //=========================================================================
404 
405 void BabayagaNLO::applySettings()
406 {
407  //--------------------
408  // Integer parameters
409  //--------------------
410  m_npar[0] = m_nPhot;
411  m_npar[1] = m_nSearchMax;
412  if (m_VPUncertainty) m_npar[5] = 1;
413  else m_npar[5] = 0;
414 
415  //--------------------
416  // Double parameters
417  //--------------------
418  m_xpar[0] = m_cmsEnergyNominal;
419  m_xpar[1] = m_pi;
420  m_xpar[2] = m_conversionFactor;
421  m_xpar[3] = m_alphaQED0;
422  m_xpar[4] = m_massElectron;
423  m_xpar[5] = m_massMuon;
424  m_xpar[6] = m_massW;
425  m_xpar[7] = m_massZ;
426  m_xpar[8] = m_widthZ;
427  m_xpar[9] = m_eMin;
428  m_xpar[10] = m_maxAcollinearity;
429  m_xpar[11] = m_epsilon;
430  m_xpar[20] = m_ScatteringAngleRange.first;
431  m_xpar[21] = m_ScatteringAngleRange.second;
432  m_xpar[30] = m_fMax;
433  m_xpar[40] = m_EnergySpread; // only for some internal checks, not the actual smaering!
434 
435  //user cuts
436  m_xpar[50] = m_eemin;
437  m_xpar[51] = m_temin;
438  m_xpar[52] = m_egmin;
439  m_xpar[53] = m_tgmin;
440  m_xpar[54] = m_eeveto;
441  m_xpar[55] = m_teveto;
442  m_xpar[56] = m_egveto;
443  m_xpar[57] = m_tgveto;
444  m_xpar[58] = m_maxprescale;
445 
446  //--------------------
447  // string parameters, wrapped to integers
448  //--------------------
449  if (m_finalState == "ee") m_npar[20] = 1;
450  else if (m_finalState == "gg") m_npar[20] = 2;
451  else if (m_finalState == "mm") m_npar[20] = 3;
452  else B2FATAL("Invalid final state: " << m_finalState);
453 
454  if (m_model == "matched") m_npar[21] = 1;
455  else if (m_model == "ps") m_npar[21] = 2;
456  else B2FATAL("Invalid matching model: " << m_model);
457 
458  if (m_order == "born") m_npar[22] = 1;
459  else if (m_order == "alpha") m_npar[22] = 2;
460  else if (m_order == "exp") m_npar[22] = 3;
461  else B2FATAL("Invalid QED order: " << m_order);
462 
463  if (m_vacPol == "off") m_npar[23] = 1;
464  else if (m_vacPol == "hadr5") m_npar[23] = 2;
465  else if (m_vacPol == "hlmnt") m_npar[23] = 3;
466  else B2FATAL("Invalid vacuum polarization code: " << m_vacPol);
467 
468  if (m_mode == "unweighted" || m_mode == "uw") m_npar[24] = 1;
469  else if (m_mode == "weighted" || m_mode == "w") m_npar[24] = 2;
470  else B2FATAL("Invalid mode: " << m_mode);
471 
472  if (m_userMode == "NONE") m_npar[25] = 1;
473  else if (m_userMode == "GAMMA") m_npar[25] = 2;
474  else if (m_userMode == "EGAMMA") m_npar[25] = 3;
475  else if (m_userMode == "ETRON") m_npar[25] = 4;
476  else if (m_userMode == "PRESCALE") m_npar[25] = 5;
477  else B2FATAL("Invalid user mode: " << m_userMode);
478 
479  // set the data file for the novosibirsk routine (not supported yet, but needs to be set)
480  size_t fileLength = m_NSKDataFile.size();
481  babayaga_setvpolnsk_(m_NSKDataFile.c_str(), &fileLength);
482 
483  //use mode to control init/generation/finalize in FORTRAN code
484  int mode = -1;
485  double ecm = -1.;
486  main_belle2_(&mode, &ecm, m_xpar, m_npar);
487 }
488 
489 
490 void BabayagaNLO::storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, TVector3 vertex, TLorentzRotation boost,
491  bool isVirtual, bool isInitial, bool isISRFSR)
492 {
493 
494 // Create particle
496  if (isVirtual) {
497  part.setStatus(MCParticle::c_IsVirtual);
498  } else if (isInitial) {
499  part.setStatus(MCParticle::c_Initial);
500  }
501 
502  // all particle of a generator are primary
503  part.addStatus(MCParticle::c_PrimaryParticle);
504 
505  // all particles produced by BABAYAGA are stable
506  part.addStatus(MCParticle::c_StableInGenerator);
507 
508  // add ISR/FSR flags to all photons that are not the primary gg pair
509  if (isISRFSR) {
510  part.addStatus(MCParticle::c_IsISRPhoton);
511  part.addStatus(MCParticle::c_IsFSRPhoton);
512  }
513 
514  part.setPDG(pdg);
515  part.setFirstDaughter(0);
516  part.setLastDaughter(0);
517  part.setMomentum(TVector3(mom[0], mom[1], mom[2]));
518  part.setMass(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
519  part.setEnergy(mom[3]);
520 
521  //boost
522  TLorentzVector p4 = part.get4Vector();
523  p4 = boost * p4;
524  part.set4Vector(p4);
525 
526  //set vertex
527  if (!isInitial) {
528  TVector3 v3 = part.getProductionVertex();
529  v3 = v3 + vertex;
530  part.setProductionVertex(v3);
531  part.setValidVertex(true);
532  }
533 }
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86