11 #include <framework/logging/Logger.h>
12 #include <framework/gearbox/Unit.h>
13 #include <generators/aafh/AAFHInterface.h>
15 #include <TDatabasePDG.h>
24 double aafhdiag36_rndm_(
int*)
26 return gRandom->Rndm();
29 double aafhdiag36_rnf100_(
int*)
31 return gRandom->Rndm();
34 void aafhdiag36_esft_maxed_(
double* weight)
36 B2ERROR(
"AAFH: Maximum weight to small, increase maxFinalWeight to at least "
40 void aafhdiag36_eswe_maxed_(
double* weight)
42 B2ERROR(
"AAFH: Maximum weight to small, increase maxSubGeneratorWeight to at least "
47 extern void aafhdiag36_init_();
49 extern void aafhdiag36_event_(
int*);
51 extern void aafhdiag36_finish_();
149 input2_.itot = tries;
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");
163 std::copy_n(weights.begin(), 4, input_.wap);
164 if (weights.size() >= 8) {
165 std::copy_n(weights.begin() + 4, 4, input_.wbp);
171 input_.eswe = subgeneratorWeight;
172 input_.esft = finalWeight;
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");
183 factor_.face = limits[0];
184 factor_.facm = limits[1];
185 factor_.facl = limits[2];
186 factor_.proc = limits[3];
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);
214 return {factor_.face, factor_.facm, factor_.facl, factor_.proc};
229 input_.eb = beamEnergy;
233 input2_.iproc = mode;
234 input2_.irejec = rejection;
237 auto* part = TDatabasePDG::Instance()->GetParticle(
m_particle.c_str());
239 B2ERROR(
"AAFH: Particle '" <<
m_particle <<
"' does not exist");
243 masses_.xml = part->Mass() / beamEnergy;
245 charge_.qcharg = -part->Charge() / 3.0;
256 TLorentzVector vec(-q[0]*input_.eb, -q[1]*input_.eb, -q[2]*input_.eb, q[3]*input_.eb);
260 p.setMass(fabs(q[4])*input_.eb);
272 aafhdiag36_event_(&status);
274 B2ERROR(
"Failed to generate event after " << status <<
" tries, giving up");
289 std::array<int, 4> pdg {{0}};
290 switch (input2_.iproc) {
295 pdg = {{13, -13, 13, -13}};
298 pdg = {{11, -11, 13, -13}};
304 pdg = {{11, -11, 11, -11}};
318 aafhdiag36_finish_();
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);
343 std::array<double, 4> idealWAP;
344 std::array<double, 4> idealWBP;
347 for (
int i = 0; i < 4; ++i) {
348 if (input_.wap[i] != 0) {
354 const double w0 = westat_.mwe[reference] * input_.wbp[reference];
355 const double n0 = westat_.iwe[reference] / input_.wap[reference];
359 for (
int i = 0; i < 4; ++i) {
360 if (input_.wap[i] == 0) {
361 idealWBP[i] = input_.wbp[i];
364 idealWBP[i] = (westat_.mwe[i] * input_.wbp[i]) / w0;
365 idealWAP[i] = n0 * input_.wap[i] / westat_.iwe[i];
367 idealWAP[i] /= idealWAP[reference];
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] <<
"],");