BBBrem.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <mdst/dataobjects/MCParticleGraph.h>
 
#include <Math/LorentzRotation.h>
#include <Math/Vector3D.h>
 
namespace Belle2 {
  class BBBrem {
 
  public:
 
    BBBrem() :
      m_eventCount(0),
      m_unweighted(true),
      m_weightCount(0),
      m_weightCountOver(0),
      m_maxWeight(0.0),
      m_maxWeightDelivered(0.0),
      m_sumWeightDelivered(0.0),
      m_sumWeightDeliveredSqr(0.0),
      m_sumWeightDeliveredOver(0.0),
      m_sumWeightDeliveredSqrOver(0.0),
      m_cmsEnergy(10.58),
      m_photonEFrac(0.000001),
      m_crossSection(0.0),
      m_crossSectionError(0.0),
      m_crossSectionOver(0.0),
      m_crossSectionErrorOver(0.0),
      m_densityCorrectionMode(0),
      m_densityCorrectionParameter(0.0),
      alpha(0.0),
      rme(0.0),
      s(0.0),
      rme2(0.0),
      rme2s(0.0),
      rls(0.0),
      z0(0.0),
      a1(0.0),
      a2(0.0),
      ac(0.0),
      sigapp(0.0),
      eb(0.0),
      pb(0.0),
      rin2pb(0.0),
      weight(0.0)
    {for (int i = 0; i < 4; i++) p1[i] = p2[i] = q1[i] = q2[i] = qk[i] = 0.0;}
 
 
    //BBBrem() : m_applyBoost(false), m_unweighted(true), m_cmsEnergy(10.58), m_photonEFrac(0.000001) {};
 
    ~BBBrem() {};
 
    void init(double cmsEnergy, double minPhotonEFrac, bool unweighted = true, double maxWeight = 2000.0, int densitymode = 1,
              double densityparameter = 1.68e-17);
 
    double generateEvent(MCParticleGraph& mcGraph, ROOT::Math::XYZVector vertex, ROOT::Math::LorentzRotation boost);
 
    double getCrossSection() { return m_crossSection; }
 
    double getCrossSectionError() { return m_crossSectionError; }
 
    double getCrossSectionOver() { return m_crossSectionOver; }
 
    double getCrossSectionErrorOver() { return m_crossSectionErrorOver; }
 
    double getMaxWeightDelivered() { return m_maxWeightDelivered; }
 
    double getSumWeightDelivered() { return m_sumWeightDelivered; }
 
    void term();
 
 
  protected:
 
    int m_eventCount;               
    bool m_unweighted;              
    long m_weightCount;             
    long m_weightCountOver;         
    double m_maxWeight;             
    double m_maxWeightDelivered;    
    double m_sumWeightDelivered;    
    double m_sumWeightDeliveredSqr; 
    double m_sumWeightDeliveredOver;    
    double m_sumWeightDeliveredSqrOver; 
    double m_cmsEnergy;             
    double m_photonEFrac;           
    double m_crossSection;          
    double m_crossSectionError;     
    double m_crossSectionOver;      
    double m_crossSectionErrorOver; 
    int m_densityCorrectionMode;   
    double m_densityCorrectionParameter;   
    void calcOutgoingLeptonsAndWeight();
 
    void storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, ROOT::Math::XYZVector vertex,
                       ROOT::Math::LorentzRotation boost,
                       bool isVirtual = false, bool isInitial = false);
 
  private:
 
    static constexpr double tomb = 3.8937966e5 / 1e6;   
    static constexpr double twopi = 2.0 * M_PI;         
    //Variable names directly taken from the FORTRAN code. Sorry.
    double alpha;  
    double rme;    
    double s;      
    double rme2;   
    double rme2s;  
    double rls;    
    double z0;     
    double a1;  
    double a2;  
    double ac;  
    double sigapp;      
    double eb;     
    double pb;     
    double rin2pb; 
    double p1[4];  
    double p2[4];  
    double q1[4];  
    double q2[4];  
    double qk[4];  
    double weight; 
  };
 
}