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>
17 #include <TDatabasePDG.h>
18 #include <TLorentzVector.h>
21 #include <generators/utilities/InitialParticleGeneration.h>
46 double rescrossphot[40];
47 double rescrossphoterr[40];
48 double rescrossphotfrac[40];
59 double hnmmaxtriallimit;
68 double biasneghitmiss;
83 double babayaganlo_rndm_(
int*)
85 double r = gRandom->Rndm();
90 double babayaganlo_getrandomcmsenergy_()
92 B2FATAL(
"babayaganlo_getrandomcmsenergy() is not implmented");
97 void main_belle2_(
int* mode,
double* ecm,
double* xpar,
int* npar);
100 void babayaga_setvpolnsk_(
const char* vpolnskfilename,
size_t* length);
103 void babayaganlo_warning_overweight_(
double* weight,
double* max)
105 B2WARNING(
"Babayaga.NLO: Maximum weight " << *max <<
" to small, increase fmax to at least " << *weight);
109 void babayaganlo_error_rejection_(
double* ratio)
111 B2ERROR(
"Babayaga.NLO: Event rejected! Ratio of cross section error increase too large: " << *ratio);
115 void babayaganlo_error_negative_(
double* weight)
117 B2ERROR(
"Babayaga.NLO: Event has negative weight: " << *weight);
121 void babayaganlo_result_nominalcmsenergy_(
double* energy)
123 B2RESULT(
"Babayaga.NLO: Nominal CMS energy (GeV): " << *energy);
126 void babayaganlo_result_weightedxsec_(
double* xsec,
double* xsecerr)
128 B2RESULT(
"Babayaga.NLO: Weighted cross section (nb): " << *xsec <<
" +/- " << *xsecerr);
131 void babayaganlo_result_unweightedxsec_(
double* xsec,
double* xsecerr)
133 B2RESULT(
"Babayaga.NLO: Unweighted cross section (nb): " << *xsec <<
" +/- " << *xsecerr);
136 void babayaganlo_result_unweightedxsec_biascorrected_(
double* xsec,
double* xsecerr)
138 B2RESULT(
"Babayaga.NLO: Unweighted cross section, bias corrected (nb): " << *xsec <<
" +/- " << *xsecerr);
141 void babayaganlo_result_unweightedxsec_overweight_(
double* xsec,
double* xsecerr)
143 B2RESULT(
"Babayaga.NLO: Unweighted cross section overweight (nb): " << *xsec <<
" +/- " << *xsecerr);
146 void babayaganlo_result_unweightedxsec_underweight_(
double* xsec,
double* xsecerr)
148 B2RESULT(
"Babayaga.NLO: Unweighted cross section underweight (nb): " << -*xsec <<
" +/- " << *xsecerr);
151 void babayaganlo_result_hitormisseff_(
double* eff)
153 B2RESULT(
"Babayaga.NLO: Hit or miss efficiency: " << *eff * 100. <<
" % ");
156 void babayaganlo_result_nover_(
int* nover)
158 B2RESULT(
"Babayaga.NLO: Points with w > fmax (bias): " << *nover);
161 void babayaganlo_result_biashit_(
double* biashit)
163 B2RESULT(
"Babayaga.NLO: Bias/hit: " << *biashit * 100. <<
" % ");
166 void babayaganlo_result_biashitpmiss_(
double* biashitpmiss)
168 B2RESULT(
"Babayaga.NLO: Bias/(hit+missed): " << *biashitpmiss * 100. <<
" % ");
171 void babayaganlo_result_nneg_(
int* nneg)
173 B2RESULT(
"Babayaga.NLO: Points with w > fmax (bias): " << *nneg);
176 void babayaganlo_result_biasneghit_(
double* biasneghit)
178 B2RESULT(
"Babayaga.NLO: Neg. bias/hit: " << *biasneghit * 100. <<
" % ");
181 void babayaganlo_result_biasneghitmiss_(
double* biasneghitmiss)
183 B2RESULT(
"Babayaga.NLO: Neg. bias/(hit+missed): " << *biasneghitmiss * 100. <<
" % ");
186 void babayaganlo_result_maxweight_(
double* sdifmax,
double* fmax)
188 B2RESULT(
"Babayaga.NLO: Maximum weight (sdifmax): " << *sdifmax <<
", user maximum weight(fmax): " << *fmax);
191 void babayaganlo_result_vpmin_(
double* xsec,
double* xsecerr,
double* frac)
193 B2RESULT(
"Babayaga.NLO: VP Uncertainty, minimum (nb): " << *xsec <<
" +/- " << *xsecerr <<
" (" << *frac <<
" %)");
196 void babayaganlo_result_vpcentral_(
double* xsec,
double* xsecerr)
198 B2RESULT(
"Babayaga.NLO: VP Uncertainty, central (nb): " << *xsec <<
" +/- " << *xsecerr);
201 void babayaganlo_result_vpmax_(
double* xsec,
double* xsecerr,
double* frac)
203 B2RESULT(
"Babayaga.NLO: VP Uncertainty, maximum (nb): " << *xsec <<
" +/- " << *xsecerr <<
" ( " << *frac <<
" %)");
206 void babayaganlo_finishedmaxsearch_(
double* fmax)
208 B2INFO(
"Babayaga.NLO: Finished maximum search: " << *fmax <<
", starting now unweighted generation");
211 void babayaganlo_fatal_usercuts_()
213 B2FATAL(
"Babayaga.NLO: UserMode is selected but wrong user cuts are provided!");
216 void babayaganlo_fatal_usercutsfs_()
218 B2FATAL(
"Babayaga.NLO: UserMode is only possible for ee final states!");
221 void babayaganlo_fatal_usercutsprescale_()
223 B2FATAL(
"Babayaga.NLO: Prescale value must be larger than 0 (zero), hint: ~100 is reasonable!");
226 void babayaganlo_fatal_weightedprescale_()
228 B2FATAL(
"Babayaga.NLO: Generator prescale is not allowed with usermode = weighted");
231 void babayaganlo_result_weightsum_(
double* wsum)
233 B2RESULT(
"Babayaga.NLO: Sum of weights: " << *wsum);
236 void babayaganlo_result_intlum_(
double* lum,
double* lumerr)
238 B2RESULT(
"Babayaga.NLO: Luminosity equivalent (using calc. xsec) (fb-1): " << *lum <<
" +/- " << *lumerr);
241 void babayaganlo_warning_prescaleweight_(
double* weight)
243 B2WARNING(
"Babayaga.NLO: Prescale of less than one, increase prescale safety margin " << *weight);
246 void babayaganlo_error_isnan1_(
double* phsp,
double* w)
248 B2ERROR(
"Babayaga.NLO: phsp (" << *phsp <<
") or w (" << *w <<
") are NAN, skipping event (nan1 type)");
251 void babayaganlo_error_isnan2_(
double* sdif)
253 B2ERROR(
"Babayaga.NLO: sdif (" << *sdif <<
") is NAN, skipping event (nan2 type)");
260 BabayagaNLO::BabayagaNLO()
262 for (
int i = 0; i < 100; ++i) {
269 setDefaultSettings();
272 BabayagaNLO::~BabayagaNLO()
277 void BabayagaNLO::setDefaultSettings()
280 m_cmsEnergyNominal = -1.;
288 m_mode =
"unweighted";
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;
301 m_maxAcollinearity = 180.0;
303 m_ScatteringAngleRange = make_pair(15.0, 165.0);
304 m_ScatteringAngleRangePhoton = make_pair(15.0, 165.0);
307 m_nSearchMax = 50000;
309 m_EnergySpread = 5e-3;
310 m_VPUncertainty =
false;
322 void BabayagaNLO::initExtraInfo()
325 extrainfo.registerInDataStore();
328 void BabayagaNLO::init()
334 void BabayagaNLO::generateEvent(
MCParticleGraph& mcGraph,
double ecm, TVector3 vertex, TLorentzRotation boost)
338 main_belle2_(&mode, &ecm, m_xpar, m_npar);
341 storeParticle(mcGraph, momset_.bq1, 11, vertex, boost,
true);
342 storeParticle(mcGraph, momset_.bp1, -11, vertex, boost,
true);
348 if (m_finalState ==
"gg") {
351 }
else if (m_finalState ==
"mm") {
356 storeParticle(mcGraph, momset_.bq2, pdg, vertex, boost,
false,
false,
false);
357 storeParticle(mcGraph, momset_.bp2, antipdg, vertex, boost,
false,
false,
false);
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);
367 eventMetaDataPtr->setGeneratedWeight(momset_.bweight);
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!");
376 float targetValue = prescale_.bprescale;
377 eventExtraInfo->addExtraInfo(
"GeneratedPrescale", targetValue);
383 void BabayagaNLO::term()
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);
397 main_belle2_(&mode, &ecm, m_xpar, m_npar);
405 void BabayagaNLO::applySettings()
411 m_npar[1] = m_nSearchMax;
412 if (m_VPUncertainty) m_npar[5] = 1;
418 m_xpar[0] = m_cmsEnergyNominal;
420 m_xpar[2] = m_conversionFactor;
421 m_xpar[3] = m_alphaQED0;
422 m_xpar[4] = m_massElectron;
423 m_xpar[5] = m_massMuon;
426 m_xpar[8] = m_widthZ;
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;
433 m_xpar[40] = m_EnergySpread;
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;
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);
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);
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);
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);
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);
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);
480 size_t fileLength = m_NSKDataFile.size();
481 babayaga_setvpolnsk_(m_NSKDataFile.c_str(), &fileLength);
486 main_belle2_(&mode, &ecm, m_xpar, m_npar);
490 void BabayagaNLO::storeParticle(
MCParticleGraph& mcGraph,
const double* mom,
int pdg, TVector3 vertex, TLorentzRotation boost,
491 bool isVirtual,
bool isInitial,
bool isISRFSR)
497 part.setStatus(MCParticle::c_IsVirtual);
498 }
else if (isInitial) {
499 part.setStatus(MCParticle::c_Initial);
503 part.addStatus(MCParticle::c_PrimaryParticle);
506 part.addStatus(MCParticle::c_StableInGenerator);
510 part.addStatus(MCParticle::c_IsISRPhoton);
511 part.addStatus(MCParticle::c_IsFSRPhoton);
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]);
522 TLorentzVector p4 = part.get4Vector();
528 TVector3 v3 = part.getProductionVertex();
530 part.setProductionVertex(v3);
531 part.setValidVertex(
true);