AAFHInterface Class Reference

Class to inferface AAFH/DIAG36 Generator written in Fortran. More...

#include <AAFHInterface.h>

Public Types
enum	EMode {   c_MuonParticle = 1 ,   c_MuonMuon = 2 ,   c_ElectonMuon = 3 ,   c_ElectronParticle = 4 ,   c_ElectonElectron = 5 }
	Generation mode. More...
enum	ERejection {   c_Once = 1 ,   c_Twice = 2 }
	Rejection mode. More...

Public Member Functions
	AAFHInterface ()=default
	Default Constructor.
void	setMaxTries (int tries)
	Set the maximum number of tries to generate an event.
void	setGeneratorWeights (const std::vector< double > &weights)
	Set the relative weights for the sub generators.
void	setMaxWeights (double subgeneratorWeight, double finalWeight)
	Set the maximum expected weights for the rejection method.
void	setSupressionLimits (std::vector< double > limits)
	Set the suppression limits when calculation the matrix elements.
void	setMinimumMass (double minMass)
	Set the minimum invariant mass for the generated event.
void	setParticle (const std::string &particle)
	Set the particle type for modes c_MuonParticle and c_ElectronParticle.
int	getMaxTries () const
	Get the maximum number of tries to generate an event.
std::vector< double >	getGeneratorWeights () const
	Get the relative weights for the sub generators.
double	getMaxSubGeneratorWeight () const
	Get the maximum expected final weight for the rejection method.
double	getMaxFinalWeight () const
	Get the maximum expected subgenerator weight for the rejection method.
std::vector< double >	getSuppressionLimits () const
	Get suppression limits.
double	getMinimumMass () const
	Get the minimum invariant mass for the generated event.
std::string	getParticle () const
	Get the particle type for modes c_MuonParticle and c_ElectronParticle.
double	getTotalCrossSection () const
	Return total cross section.
double	getTotalCrossSectionError () const
	Return error on the total cross section.
void	initialize (double beamEnergy, EMode mode, ERejection rejection)
	initialize the generator
void	generateEvent (MCParticleGraph &mpg)
	generate one event and add it to the graph in CMS
void	finish ()
	calculate total cross section

Private Member Functions
void	updateParticle (MCParticleGraph::GraphParticle &p, double *q, int pdg, bool isInitial=false)
	update particle with generated values

Private Attributes
double	m_minMass {0}
	minimum invariant mass
std::string	m_particle {"tau-"}
	name of the particle for modes c_MuonParticle and c_ElectronParticle
int	m_particlePDG {0}
	pdg of the particle for modes c_MuonParticle and c_ElectronParticle, gets set in initialize()

Detailed Description

Class to inferface AAFH/DIAG36 Generator written in Fortran.

Almost all parameters can be set using this class. First all parameters should be set, then call initialize(), an arbitary number of generateEvent() and then finish(). After finish(), the total cross section can be obained using getTotalCrossSection()

Definition at line 29 of file AAFHInterface.h.

Member Enumeration Documentation

EMode

enum EMode

Generation mode.

Enumerator
c_MuonParticle	E+E- -> Mu+Mu-L+L- where L can be anything, defaults to tau.
c_MuonMuon	E+E- -> Mu+Mu-Mu+Mu-.
c_ElectonMuon	E+E- -> E+E-Mu+Mu-.
c_ElectronParticle	E+E- -> E+E-L+L- where L can be anything, defaults to tau.
c_ElectonElectron	E+E- -> E+E-E+E-.

Definition at line 32 of file AAFHInterface.h.

               {
      c_MuonParticle = 1,
      c_MuonMuon = 2,
      c_ElectonMuon = 3,
      c_ElectronParticle = 4,
      c_ElectonElectron = 5,
    };

ERejection

enum ERejection

Rejection mode.

Enumerator
c_Once	Apply rejection only once for the final event.
c_Twice	Apply rejection first for each sub generator and then for the final event.

Definition at line 46 of file AAFHInterface.h.

                    {
      c_Once = 1,
      c_Twice = 2,
    };

Member Function Documentation

finish()

void finish

(

)

calculate total cross section

Definition at line 319 of file AAFHInterface.cc.

{
  aafhdiag36_finish_();
  // The following code just prints the results. If B2RESULT has been compiled
  // out this is unneccesary and cppcheck complains so we enclose this in an
  // ifdef to only be executed if B2RESULT might actually be shown
#ifndef LOG_NO_B2RESULT
  B2RESULT("AAFH Generation finished.");
  B2RESULT("Final cross section:         "
           << std::setprecision(3) << output_.crosec << "+-"
           << std::setprecision(3) << output_.ercros << " nb");
  B2RESULT("Minimum invariant mass:      "
           << sqrt(bound_.w2min)*input_.eb << " GeV");
  B2RESULT("Maximum invariant mass:      "
           << sqrt(bound_.w2max)*input_.eb << " GeV");
  B2RESULT("Maximum subgenerator weight: " << westat_.maxwe << " ("
           << westat_.mwe[0] << ", "
           << westat_.mwe[1] << ", "
           << westat_.mwe[2] << ", "
           << westat_.mwe[3] << ")");
  B2RESULT("Maximum final weight:        " << ftstat_.ftmax);
  //The fortran source says that it's advisable if all subgenerators have the
  //same probability to account for all peaks in the differental cross section
  //and that it runs with the highest efficiency when the maximum weights in
  //all sub generators are equal. So let's calculate appropriate subgenerator
  //weights by looking how many events were generated per sub generator or the
  //maximum weights
  std::array<double, 4> idealWAP;
  std::array<double, 4> idealWBP;
  //find the first non-zero weight
  int reference = 0;
  for (int i = 0; i < 4; ++i) {
    if (input_.wap[i] != 0) {
      reference = i;
      break;
    }
  }
  //and use that es reference
  const double w0 = westat_.mwe[reference] * input_.wbp[reference];
  const double n0 = westat_.iwe[reference] / input_.wap[reference];
  //now calculate the relative weights we need to get
  //1) same maximum weight which is inversely proportional to to wbp
  //2) same number of generated events which is proportional to wap
  for (int i = 0; i < 4; ++i) {
    if (input_.wap[i] == 0) {
      idealWBP[i] = input_.wbp[i];
      idealWAP[i] = 0;
    } else {
      idealWBP[i] = (westat_.mwe[i] * input_.wbp[i]) / w0;
      idealWAP[i] = n0 * input_.wap[i] / westat_.iwe[i];
      //and normalize to 1.0
      idealWAP[i] /= idealWAP[reference];
    }
  }
 
  B2RESULT("To get the same maximum event weight and chance for each sub "
           "generator use these parameters:");
  B2RESULT("   'maxFinalWeight': " << std::setprecision(3)
           << ftstat_.ftmax * 1.1 << ",");
  B2RESULT("   'maxSubgeneratorWeight': " << std::setprecision(3)
           << westat_.maxwe * 1.1 << ",");
  B2RESULT("   'subgeneratorWeights': [" << std::scientific
           << std::setprecision(3) << idealWAP[0] << ", "
           << std::setprecision(3) << idealWAP[1] << ", "
           << std::setprecision(3) << idealWAP[2] << ", "
           << std::setprecision(3) << idealWAP[3] << ", ");
  B2RESULT("                           " << std::scientific
           << std::setprecision(3) << idealWBP[0] << ", "
           << std::setprecision(3) << idealWBP[1] << ", "
           << std::setprecision(3) << idealWBP[2] << ", "
           << std::setprecision(3) << idealWBP[3] << "],");
#endif
}

generateEvent()

void generateEvent

(

MCParticleGraph &

mpg

)

generate one event and add it to the graph in CMS

Definition at line 272 of file AAFHInterface.cc.

{
  int status {0};
  aafhdiag36_event_(&status);
  if (status != 0) {
    B2ERROR("Failed to generate event after " << status << " tries, giving up");
    return;
  }
 
  //initial beam particles
  auto& p1 = mpg.addParticle();
  auto& p2 = mpg.addParticle();
 
  //generated particles
  auto& p3 = mpg.addParticle();
  auto& p4 = mpg.addParticle();
  auto& p5 = mpg.addParticle();
  auto& p6 = mpg.addParticle();
 
  //array of generated particle pdg codes depends on generated mode
  std::array<int, 4> pdg {{0}};
  switch (input2_.iproc) {
    case c_MuonParticle:
      pdg = {{m_particlePDG, -m_particlePDG, 13, -13}};
      break;
    case c_MuonMuon:
      pdg = {{13, -13, 13, -13}};
      break;
    case c_ElectonMuon:
      pdg = {{11, -11, 13, -13}};
      break;
    case c_ElectronParticle:
      pdg = {{11, -11, m_particlePDG, -m_particlePDG}};
      break;
    case c_ElectonElectron:
      pdg = {{11, -11, 11, -11}};
      break;
  }
 
  updateParticle(p1, momenz_.q1, 11, true);
  updateParticle(p2, momenz_.q2, -11, true);
  updateParticle(p3, momenz_.q3, pdg[0]);
  updateParticle(p4, momenz_.q4, pdg[1]);
  updateParticle(p5, momenz_.q5, pdg[2]);
  updateParticle(p6, momenz_.q6, pdg[3]);
}

getGeneratorWeights()

std::vector< double > getGeneratorWeights

(

)

const

Get the relative weights for the sub generators.

Definition at line 197 of file AAFHInterface.cc.

{
  std::vector<double> weights(8);
  std::copy_n(input_.wap, 4, weights.begin());
  std::copy_n(input_.wbp, 4, weights.begin() + 4);
  return weights;
}

getMaxFinalWeight()

double getMaxFinalWeight

(

)

const

Get the maximum expected subgenerator weight for the rejection method.

Definition at line 210 of file AAFHInterface.cc.

{
  return input_.esft;
}

getMaxSubGeneratorWeight()

double getMaxSubGeneratorWeight

(

)

const

Get the maximum expected final weight for the rejection method.

Definition at line 205 of file AAFHInterface.cc.

{
  return input_.eswe;
}

getMaxTries()

int getMaxTries

(

)

const

Get the maximum number of tries to generate an event.

Definition at line 192 of file AAFHInterface.cc.

{
  return input2_.itot;
}

getMinimumMass()

double getMinimumMass ( ) const

inline

Get the minimum invariant mass for the generated event.

Definition at line 85 of file AAFHInterface.h.

{ return m_minMass; }

getParticle()

std::string getParticle ( ) const

inline

Get the particle type for modes c_MuonParticle and c_ElectronParticle.

Definition at line 87 of file AAFHInterface.h.

{ return m_particle; }

getSuppressionLimits()

std::vector< double > getSuppressionLimits

(

)

const

Get suppression limits.

Definition at line 215 of file AAFHInterface.cc.

{
  return {factor_.face, factor_.facm, factor_.facl, factor_.proc};
}

getTotalCrossSection()

double getTotalCrossSection

(

)

const

Return total cross section.

Definition at line 220 of file AAFHInterface.cc.

{
  return output_.crosec * Unit::nb;
}

getTotalCrossSectionError()

double getTotalCrossSectionError

(

)

const

Return error on the total cross section.

Definition at line 225 of file AAFHInterface.cc.

{
  return output_.ercros * Unit::nb;
}

initialize()

void initialize	(	double	beamEnergy,
		EMode	mode,
		ERejection	rejection
	)

initialize the generator

Definition at line 230 of file AAFHInterface.cc.

{
  input_.eb = beamEnergy;
  //minimum invariant mass is in units of beam energy and needs to be squared
  bound_.w2min = (m_minMass / beamEnergy) * (m_minMass / beamEnergy);
  //set mode and rejection scheme
  input2_.iproc = mode;
  input2_.irejec = rejection;
 
  //initialize the L
  auto* part = TDatabasePDG::Instance()->GetParticle(m_particle.c_str());
  if (!part) {
    B2ERROR("AAFH: Particle '" << m_particle << "' does not exist");
    return;
  }
  //set the mass in units of beam energy
  masses_.xml = part->Mass() / beamEnergy;
  //set the charge but with opposite sign
  charge_.qcharg = -part->Charge() / 3.0;
  //and remember pdg code
  m_particlePDG = part->PdgCode();
 
  aafhdiag36_init_();
}

setGeneratorWeights()

void setGeneratorWeights

(

const std::vector< double > &

weights

)

Set the relative weights for the sub generators.

Definition at line 151 of file AAFHInterface.cc.

{
  if (weights.size() > 8) {
    B2WARNING("AAFH: more than 8 generator weights supplied, "
              "ignoring the extra weights");
  } else if (weights.size() > 4 && weights.size() < 8) {
    B2WARNING("AAFH: more than 4 but less than 8 generator weights supplied, "
              "ignoring everything afer the first four");
  } else if (weights.size() < 4) {
    B2ERROR("AAFH: not enough generator weights, need exactly four or eight values");
  }
  std::copy_n(weights.begin(), 4, input_.wap);
  if (weights.size() >= 8) {
    std::copy_n(weights.begin() + 4, 4, input_.wbp);
  }
}

setMaxTries()

void setMaxTries

(

int

tries

)

Set the maximum number of tries to generate an event.

Definition at line 146 of file AAFHInterface.cc.

{
  input2_.itot = tries;
}

setMaxWeights()

void setMaxWeights	(	double	subgeneratorWeight,
		double	finalWeight
	)

Set the maximum expected weights for the rejection method.

If the weights are to large generation will be ineffective, if they are to small it will generate Error messages

Definition at line 168 of file AAFHInterface.cc.

{
  input_.eswe = subgeneratorWeight;
  input_.esft = finalWeight;
}

setMinimumMass()

void setMinimumMass ( double  minMass )

inline

Set the minimum invariant mass for the generated event.

Definition at line 67 of file AAFHInterface.h.

{ m_minMass = minMass; }

setParticle()

void setParticle ( const std::string &  particle )

inline

Set the particle type for modes c_MuonParticle and c_ElectronParticle.

Definition at line 69 of file AAFHInterface.h.

    {
      m_particle = particle;
    }

setSupressionLimits()

void setSupressionLimits

(

std::vector< double >

limits

)

Set the suppression limits when calculation the matrix elements.

Definition at line 174 of file AAFHInterface.cc.

{
  if (limits.size() == 1) {
    B2WARNING("AAFH: Only one suppression limit supplied, using same value for all limits");
    factor_.face = limits[0];
    factor_.facm = limits[0];
    factor_.facl = limits[0];
    factor_.proc = limits[0];
  } else if (limits.size() == 4) {
    factor_.face = limits[0];
    factor_.facm = limits[1];
    factor_.facl = limits[2];
    factor_.proc = limits[3];
  } else {
    B2ERROR("AAFH: Wrong number of suppression limits: supply either one or 4 values");
  }
}

updateParticle()

void updateParticle ( MCParticleGraph::GraphParticle &  p, double *  q, int  pdg, bool  isInitial = false  )

private

update particle with generated values

Parameters

p	MCParticle to be updated
q	4-Vector + mass from the generator
pdg	pdgcode of the particle
isInitial	incoming beam particles

Definition at line 255 of file AAFHInterface.cc.

{
  // all values in q are in units of beam energy
  // z-direction inversed since AAFH uses different convention
  ROOT::Math::PxPyPzEVector vec(-q[0]*input_.eb, -q[1]*input_.eb, -q[2]*input_.eb, q[3]*input_.eb);
  p.set4Vector(vec);
 
  // mass is multiplied by charge sign so take the absolute
  p.setMass(fabs(q[4])*input_.eb);
  p.setPDG(pdg);
  p.setStatus(MCParticle::c_PrimaryParticle | MCParticle::c_StableInGenerator);
 
  if (isInitial) {
    p.addStatus(MCParticle::c_Initial);
  }
}

Member Data Documentation

m_minMass

double m_minMass {0}

private

minimum invariant mass

Definition at line 108 of file AAFHInterface.h.

m_particle

std::string m_particle {"tau-"}

private

name of the particle for modes c_MuonParticle and c_ElectronParticle

Definition at line 110 of file AAFHInterface.h.

m_particlePDG

int m_particlePDG {0}

private

pdg of the particle for modes c_MuonParticle and c_ElectronParticle, gets set in initialize()

Definition at line 113 of file AAFHInterface.h.

---

The documentation for this class was generated from the following files:

• generators/aafh/include/AAFHInterface.h
• generators/aafh/src/AAFHInterface.cc