bsim Namespace Reference

bsim namespace for beam simulation classes and functions More...

Classes
class	Ancestor
	More...
class	Decay
	More...
class	Dk2Nu
	More...
class	DkMeta
	More...
class	Location
	More...
class	NuChoice
	More...
class	NuRay
	More...
class	TgtExit
	More...
class	Traj
	More...

Typedefs
typedef enum bsim::flgbitval	flgbitval_t
typedef enum bsim::dkproc	dkproc_t

Enumerations
enum	flgbitval {   kFlgOverflow = 0x00000001, kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000, kFlgOverflow = 0x00000001,   kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000 }
enum	dkproc {   dkp_unknown = 0, dkp_k0l_nuepimep = 1, dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3,   dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5, dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7,   dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9, dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11,   dkp_mum_nusep = 12, dk_pip_numumup = 13, dk_pim_numubmum = 14, dk_mum_capture = 15,   dkp_maximum, dkp_other = 999, dkp_unknown = 0, dkp_k0l_nuepimep = 1,   dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3, dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5,   dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7, dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9,   dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11, dkp_mum_nusep = 12, dk_pip_numumup = 13,   dk_pim_numubmum = 14, dk_mum_capture = 15, dkp_maximum, dkp_other = 999 }
enum	flgbitval {   kFlgOverflow = 0x00000001, kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000, kFlgOverflow = 0x00000001,   kMaskReserved = 0x0000FFFF, kMaskUser = 0xFFFF0000 }
enum	dkproc {   dkp_unknown = 0, dkp_k0l_nuepimep = 1, dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3,   dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5, dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7,   dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9, dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11,   dkp_mum_nusep = 12, dk_pip_numumup = 13, dk_pim_numubmum = 14, dk_mum_capture = 15,   dkp_maximum, dkp_other = 999, dkp_unknown = 0, dkp_k0l_nuepimep = 1,   dkp_k0l_nuebpipem = 2, dkp_k0l_numupimmup = 3, dkp_k0l_numubpipmum = 4, dkp_kp_numumup = 5,   dkp_kp_nuepi0ep = 6, dkp_kp_numupi0mup = 7, dkp_kp_numubmum = 8, dkp_kp_nuebpi0em = 9,   dkp_kp_numubpi0mum = 10, dkp_mup_nusep = 11, dkp_mum_nusep = 12, dk_pip_numumup = 13,   dk_pim_numubmum = 14, dk_mum_capture = 15, dkp_maximum, dkp_other = 999 }

Functions
int	calcEnuWgt (const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
	workhorse routine More...
int	calcEnuWgt (const bsim::Dk2Nu *dk2nu, const TVector3 &xyz, double &enu, double &wgt_xy)
	alternate interface More...
void	calcLocationWeights (const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu)
	user interface More...
Bool_t	IsDefault (Float_t intval)
Bool_t	IsDefault (Int_t intval)
Bool_t	IsDefault (UInt_t intval)
Bool_t	IsDefault (Double_t intval)
Bool_t	IsDefault (Bool_t intval)
void	readWeightLocations (std::string locfilename, std::vector< bsim::Location > &locations)
void	readWeightLocations (std::string locfilename, bsim::DkMeta *dkmeta)
	a variant that will fill the dkmeta object More...
void	printWeightLocations (std::vector< bsim::Location > &locations)
	print the locations More...
void	printWeightLocations (bsim::DkMeta *dkmeta)
	a variant that prints the locations in the dkmeta object More...

Variables
static const Float_t	kDfltFloat = -9999.99
static const Int_t	kDfltInt = -9999
static const UInt_t	kDfltUInt = 9999
static const Double_t	kDfltDouble = -9999.99
static const Bool_t	kDfltBool = false
static const Float_t	kDfltFloat = -9999.99
static const Int_t	kDfltInt = -9999
static const UInt_t	kDfltUInt = 9999
static const Double_t	kDfltDouble = -9999.99
static const Bool_t	kDfltBool = false

Detailed Description

bsim namespace for beam simulation classes and functions

---

Typedef Documentation

typedef enum bsim::dkproc bsim::dkproc_t

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

typedef enum bsim::flgbitval bsim::flgbitval_t

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumeration Type Documentation

enum bsim::dkproc

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

Enumerator
dkp_unknown
dkp_k0l_nuepimep	k0long => nu_e + pi- + e+
dkp_k0l_nuebpipem	k0long => nu_e_bar + p+ + e-
dkp_k0l_numupimmup	k0long => nu_mu + pi- + mu+
dkp_k0l_numubpipmum	k0long => nu_mu_bar + pi+ + mu-
dkp_kp_numumup	k+ => nu_mu + mu+
dkp_kp_nuepi0ep	k+ => nu_e + pi0 + e+
dkp_kp_numupi0mup	k+ => nu_mu + pi0 + mu+
dkp_kp_numubmum	k- => nu_mu_bar + mu-
dkp_kp_nuebpi0em	k- => nu_e_bar + pi0 + e-
dkp_kp_numubpi0mum	k- => nu_mu_bar + pi0 + mu-
dkp_mup_nusep	mu+ => nu_mu_bar + nu_e + e+
dkp_mum_nusep	mu- => nu_mu + nu_e_bar + e-
dk_pip_numumup	pi+ => nu_mu + mu+
dk_pim_numubmum	pi- => nu_mu_bar + mu-
dk_mum_capture	mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)
dkp_maximum	one-beyond end for iterating
dkp_other	flag for unusual cases
dkp_unknown
dkp_k0l_nuepimep	k0long => nu_e + pi- + e+
dkp_k0l_nuebpipem	k0long => nu_e_bar + p+ + e-
dkp_k0l_numupimmup	k0long => nu_mu + pi- + mu+
dkp_k0l_numubpipmum	k0long => nu_mu_bar + pi+ + mu-
dkp_kp_numumup	k+ => nu_mu + mu+
dkp_kp_nuepi0ep	k+ => nu_e + pi0 + e+
dkp_kp_numupi0mup	k+ => nu_mu + pi0 + mu+
dkp_kp_numubmum	k- => nu_mu_bar + mu-
dkp_kp_nuebpi0em	k- => nu_e_bar + pi0 + e-
dkp_kp_numubpi0mum	k- => nu_mu_bar + pi0 + mu-
dkp_mup_nusep	mu+ => nu_mu_bar + nu_e + e+
dkp_mum_nusep	mu- => nu_mu + nu_e_bar + e-
dk_pip_numumup	pi+ => nu_mu + mu+
dk_pim_numubmum	pi- => nu_mu_bar + mu-
dk_mum_capture	mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)
dkp_maximum	one-beyond end for iterating
dkp_other	flag for unusual cases

Definition at line 69 of file dk2nu.h.

                       {
     dkp_unknown         =  0,
     dkp_k0l_nuepimep    =  1,  ///< k0long => nu_e + pi- + e+
     dkp_k0l_nuebpipem   =  2,  ///< k0long => nu_e_bar + p+ + e-
     dkp_k0l_numupimmup  =  3,  ///< k0long => nu_mu + pi- + mu+
     dkp_k0l_numubpipmum =  4,  ///< k0long => nu_mu_bar + pi+ + mu-
     dkp_kp_numumup      =  5,  ///< k+ => nu_mu + mu+
     dkp_kp_nuepi0ep     =  6,  ///< k+ => nu_e + pi0 + e+
     dkp_kp_numupi0mup   =  7,  ///< k+ => nu_mu + pi0 + mu+
     dkp_kp_numubmum     =  8,  ///< k- => nu_mu_bar + mu-
     dkp_kp_nuebpi0em    =  9,  ///< k- => nu_e_bar + pi0 + e-
     dkp_kp_numubpi0mum  = 10,  ///< k- => nu_mu_bar + pi0 + mu-
     dkp_mup_nusep       = 11,  ///< mu+ => nu_mu_bar + nu_e + e+
     dkp_mum_nusep       = 12,  ///< mu- => nu_mu + nu_e_bar + e-
     dk_pip_numumup      = 13,  ///< pi+ => nu_mu + mu+
     dk_pim_numubmum     = 14,  ///< pi- => nu_mu_bar + mu-
     dk_mum_capture      = 15,  ///< mu- + (A,Z) => nu
                                ///   (different kinematics from dkp_mum_nusep)
     dkp_maximum,               ///< one-beyond end for iterating
     dkp_other           = 999  ///< flag for unusual cases
   } dkproc_t;

enum bsim::dkproc

Enumeration of decay processes, stored in "ndecay" stored as integer; these are for reference DK2NU: should there be an associated AsString() method that returns a text (optionally formatted for latex?)?

Enumerator
dkp_unknown
dkp_k0l_nuepimep	k0long => nu_e + pi- + e+
dkp_k0l_nuebpipem	k0long => nu_e_bar + p+ + e-
dkp_k0l_numupimmup	k0long => nu_mu + pi- + mu+
dkp_k0l_numubpipmum	k0long => nu_mu_bar + pi+ + mu-
dkp_kp_numumup	k+ => nu_mu + mu+
dkp_kp_nuepi0ep	k+ => nu_e + pi0 + e+
dkp_kp_numupi0mup	k+ => nu_mu + pi0 + mu+
dkp_kp_numubmum	k- => nu_mu_bar + mu-
dkp_kp_nuebpi0em	k- => nu_e_bar + pi0 + e-
dkp_kp_numubpi0mum	k- => nu_mu_bar + pi0 + mu-
dkp_mup_nusep	mu+ => nu_mu_bar + nu_e + e+
dkp_mum_nusep	mu- => nu_mu + nu_e_bar + e-
dk_pip_numumup	pi+ => nu_mu + mu+
dk_pim_numubmum	pi- => nu_mu_bar + mu-
dk_mum_capture	mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)
dkp_maximum	one-beyond end for iterating
dkp_other	flag for unusual cases
dkp_unknown
dkp_k0l_nuepimep	k0long => nu_e + pi- + e+
dkp_k0l_nuebpipem	k0long => nu_e_bar + p+ + e-
dkp_k0l_numupimmup	k0long => nu_mu + pi- + mu+
dkp_k0l_numubpipmum	k0long => nu_mu_bar + pi+ + mu-
dkp_kp_numumup	k+ => nu_mu + mu+
dkp_kp_nuepi0ep	k+ => nu_e + pi0 + e+
dkp_kp_numupi0mup	k+ => nu_mu + pi0 + mu+
dkp_kp_numubmum	k- => nu_mu_bar + mu-
dkp_kp_nuebpi0em	k- => nu_e_bar + pi0 + e-
dkp_kp_numubpi0mum	k- => nu_mu_bar + pi0 + mu-
dkp_mup_nusep	mu+ => nu_mu_bar + nu_e + e+
dkp_mum_nusep	mu- => nu_mu + nu_e_bar + e-
dk_pip_numumup	pi+ => nu_mu + mu+
dk_pim_numubmum	pi- => nu_mu_bar + mu-
dk_mum_capture	mu- + (A,Z) => nu (different kinematics from dkp_mum_nusep)
dkp_maximum	one-beyond end for iterating
dkp_other	flag for unusual cases

Definition at line 69 of file dk2nu.h.

                       {
     dkp_unknown         =  0,
     dkp_k0l_nuepimep    =  1,  ///< k0long => nu_e + pi- + e+
     dkp_k0l_nuebpipem   =  2,  ///< k0long => nu_e_bar + p+ + e-
     dkp_k0l_numupimmup  =  3,  ///< k0long => nu_mu + pi- + mu+
     dkp_k0l_numubpipmum =  4,  ///< k0long => nu_mu_bar + pi+ + mu-
     dkp_kp_numumup      =  5,  ///< k+ => nu_mu + mu+
     dkp_kp_nuepi0ep     =  6,  ///< k+ => nu_e + pi0 + e+
     dkp_kp_numupi0mup   =  7,  ///< k+ => nu_mu + pi0 + mu+
     dkp_kp_numubmum     =  8,  ///< k- => nu_mu_bar + mu-
     dkp_kp_nuebpi0em    =  9,  ///< k- => nu_e_bar + pi0 + e-
     dkp_kp_numubpi0mum  = 10,  ///< k- => nu_mu_bar + pi0 + mu-
     dkp_mup_nusep       = 11,  ///< mu+ => nu_mu_bar + nu_e + e+
     dkp_mum_nusep       = 12,  ///< mu- => nu_mu + nu_e_bar + e-
     dk_pip_numumup      = 13,  ///< pi+ => nu_mu + mu+
     dk_pim_numubmum     = 14,  ///< pi- => nu_mu_bar + mu-
     dk_mum_capture      = 15,  ///< mu- + (A,Z) => nu
                                ///   (different kinematics from dkp_mum_nusep)
     dkp_maximum,               ///< one-beyond end for iterating
     dkp_other           = 999  ///< flag for unusual cases
   } dkproc_t;

enum bsim::flgbitval

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumerator
kFlgOverflow
kMaskReserved
kMaskUser
kFlgOverflow
kMaskReserved
kMaskUser

Definition at line 57 of file dk2nu.h.

                          {
     kFlgOverflow    = 0x00000001,
     kMaskReserved   = 0x0000FFFF,
     kMaskUser       = 0xFFFF0000
   } flgbitval_t;

enum bsim::flgbitval

All the data members are public as these classes are used as generalized structs. As they will be branches of a TTree no specialized naming indicators signifying that they are member data of a class will be used, nor will any fancy capitalization schemes.

All classes must implement a clear() method that resets their values to an identifiably invalid state or clears any vectors. Additionally classes should provide a AsString() method for formatting themselves

for use output.

Specialized enumerations Proposed flag bits:

Enumerator
kFlgOverflow
kMaskReserved
kMaskUser
kFlgOverflow
kMaskReserved
kMaskUser

Definition at line 57 of file dk2nu.h.

                          {
     kFlgOverflow    = 0x00000001,
     kMaskReserved   = 0x0000FFFF,
     kMaskUser       = 0xFFFF0000
   } flgbitval_t;

Function Documentation

int bsim::calcEnuWgt	(	const bsim::Decay &	decay,
		const TVector3 &	xyz,
		double &	enu,
		double &	wgt_xy
	)

workhorse routine

Definition at line 48 of file calcLocationWeights.cxx.

References beta, om::cerr, dkp_mum_nusep, dkp_mup_nusep, allTimeWatchdog::endl, bsim::Decay::mupare, bsim::Decay::muparpx, bsim::Decay::muparpy, bsim::Decay::muparpz, bsim::Decay::ndecay, bsim::Decay::necm, bsim::Decay::norig, bsim::Decay::ntype, partial, bsim::Decay::pdpx, bsim::Decay::pdpy, bsim::Decay::pdpz, bsim::Decay::ppdxdz, bsim::Decay::ppdydz, bsim::Decay::ppenergy, bsim::Decay::pppz, bsim::Decay::ptype, Munits::rad, ana::Sqrt(), bsim::Decay::vx, bsim::Decay::vy, and bsim::Decay::vz.

Referenced by calcEnuWgt(), calcLocationWeights(), genie::flux::GDk2NuFlux::GenerateNext_weighted(), and genie::flux::GDk2NuFluxXMLHelper::~GDk2NuFluxXMLHelper().

{
  // Neutrino Energy and Weight at arbitrary point
  // Based on:
  //   NuMI-NOTE-BEAM-0109 (MINOS DocDB # 109)
  //   Title:   Neutrino Beam Simulation using PAW with Weighted Monte Carlos
  //   Author:  Rick Milburn
  //   Date:    1995-10-01

  // History:
  // jzh  3/21/96 grab R.H.Milburn's weighing routine
  // jzh  5/ 9/96 substantially modify the weighting function use dot product
  //              instead of rotation vecs to get theta get all info except
  //              det from ADAMO banks neutrino parent is in Particle.inc
  //              Add weighting factor for polarized muon decay
  // jzh  4/17/97 convert more code to double precision because of problems
  //              with Enu>30 GeV
  // rwh 10/ 9/08 transliterate function from f77 to C++

  // Original function description:
  //   Real function for use with PAW Ntuple To transform from destination
  //   detector geometry to the unit sphere moving with decaying hadron with
  //   velocity v, BETA=v/c, etc..  For (pseudo)scalar hadrons the decays will
  //   be isotropic in this  sphere so the fractional area (out of 4-pi) is the
  //   fraction of decays that hit the target.  For a given target point and
  //   area, and given x-y components of decay transverse location and slope,
  //   and given decay distance from target ans given decay GAMMA and
  //   rest-frame neutrino energy, the lab energy at the target and the
  //   fractional solid angle in the rest-frame are determined.
  //   For muon decays, correction for non-isotropic nature of decay is done.

  // Arguments:
  //    dk2nu    :: contains current decay information
  //    xyz      :: 3-vector of position to evaluate
  //                in *beam* frame coordinates  (cm units)
  //    enu      :: resulting energy
  //    wgt_xy   :: resulting weight
  // Return:
  //    (int)    :: error code
  // Assumptions:
  //    Energies given in GeV
  //    Particle codes have been translated from GEANT into PDG codes

  // for now ... these masses _should_ come from TDatabasePDG
  // but use these hard-coded values to "exactly" reproduce old code
  //
  // old mass values are v01_07_* and before
  // new mass values (v01_08_* and after) are Geant4 v4.10.3 values
  //
#ifdef HISTORIC_MASS
  const double kPIMASS      = 0.13957;
  const double kKMASS       = 0.49368;
  const double kK0MASS      = 0.49767;
  const double kMUMASS      = 0.105658389;
  const double kOMEGAMASS   = 1.67245;
#else
  const double kPIMASS      = 0.1395701;     // 0.13957;
  const double kKMASS       = 0.493677;      // 0.49368;
  const double kK0MASS      = 0.497614;      // 0.49767;
  const double kMUMASS      = 0.1056583715;  // 0.105658389;
  const double kOMEGAMASS   = 1.67245;       // 1.67245;
#endif
  // from CLHEP/Units/PhysicalConstants.h
  // used by Geant as CLHEP::neutron_mass_c2
  const double kNEUTRONMASS = 0.93956536;

  const int kpdg_nue       =   12;  // extended Geant 53
  const int kpdg_nuebar    =  -12;  // extended Geant 52
  const int kpdg_numu      =   14;  // extended Geant 56
  const int kpdg_numubar   =  -14;  // extended Geant 55

  const int kpdg_muplus      =   -13;  // Geant  5
  const int kpdg_muminus     =    13;  // Geant  6
  const int kpdg_pionplus    =   211;  // Geant  8
  const int kpdg_pionminus   =  -211;  // Geant  9
  const int kpdg_k0long      =   130;  // Geant 10  ( K0=311, K0S=310 )
  const int kpdg_k0short     =   310;  // Geant 16
  const int kpdg_k0mix       =   311;
  const int kpdg_kaonplus    =   321;  // Geant 11
  const int kpdg_kaonminus   =  -321;  // Geant 12
  const int kpdg_omegaminus  =  3334;  // Geant 24
  const int kpdg_omegaplus   = -3334;  // Geant 32
  const int kpdg_neutron     =  2112;
  const int kpdg_antineutron = -2112;

  const double kRDET = 100.0;   // set to flux per 100 cm radius

  double xpos = xyz.X();
  double ypos = xyz.Y();
  double zpos = xyz.Z();

  enu    = 0.0;  // don't know what the final value is
  wgt_xy = 0.0;  // but set these in case we return early due to error


  // in principle we should get these from the particle DB
  // but for consistency testing use the hardcoded values
  double parent_mass = kPIMASS;
  switch ( decay.ptype ) {
  case kpdg_pionplus:
  case kpdg_pionminus:
    parent_mass = kPIMASS;
    break;
  case kpdg_kaonplus:
  case kpdg_kaonminus:
    parent_mass = kKMASS;
    break;
  case kpdg_k0long:
  case kpdg_k0short:
  case kpdg_k0mix:
    parent_mass = kK0MASS;
    break;
  case kpdg_muplus:
  case kpdg_muminus:
    parent_mass = kMUMASS;
    break;
  case kpdg_omegaminus:
  case kpdg_omegaplus:
    parent_mass = kOMEGAMASS;
    break;
  case kpdg_neutron:
  case kpdg_antineutron:
    parent_mass = kNEUTRONMASS;
    break;
  default:
    std::cerr << "bsim::calcEnuWgt unknown particle type " << decay.ptype
              << std::endl << std::flush;
    enu    = 0.0;
    wgt_xy = 0.0;
    return 1;
  }

  double parentp2 = ( decay.pdpx*decay.pdpx +
                      decay.pdpy*decay.pdpy +
                      decay.pdpz*decay.pdpz );
  double parent_energy = TMath::Sqrt( parentp2 +
                                     parent_mass*parent_mass);
  double parentp = TMath::Sqrt( parentp2 );

  double gamma     = parent_energy / parent_mass;
  double gamma_sqr = gamma * gamma;
  double beta_mag  = TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );

  // Get the neutrino energy in the parent decay CM
  double enuzr = decay.necm;
  // Get angle from parent line of flight to chosen point in beam frame
  double rad = TMath::Sqrt( (xpos-decay.vx)*(xpos-decay.vx) +
                            (ypos-decay.vy)*(ypos-decay.vy) +
                            (zpos-decay.vz)*(zpos-decay.vz) );

  double emrat = 1.0;
  double costh_pardet = -999., theta_pardet = -999.;

  // boost correction, but only if parent hasn't stopped
  if ( parentp > 0. ) {
    costh_pardet = ( decay.pdpx*(xpos-decay.vx) +
                     decay.pdpy*(ypos-decay.vy) +
                     decay.pdpz*(zpos-decay.vz) )
                     / ( parentp * rad);
    if ( costh_pardet >  1.0 ) costh_pardet =  1.0;
    if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
    theta_pardet = TMath::ACos(costh_pardet);

    // Weighted neutrino energy in beam, approx, good for small theta
    emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
  }

  enu = emrat * enuzr;  // the energy ... normally

  // Get solid angle/4pi for detector element
  // small angle approximation, fixed by Alex Radovic
  //SAA//  double sangdet = ( kRDET*kRDET /
  //SAA//                   ( (zpos-decay.vz)*(zpos-decay.vz) ) ) / 4.0;
  double sanddetcomp = TMath::Sqrt( ( (xpos-decay.vx)*(xpos-decay.vx) ) +
                                    ( (ypos-decay.vy)*(ypos-decay.vy) ) +
                                    ( (zpos-decay.vz)*(zpos-decay.vz) )   );
  double sangdet = (1.0-TMath::Cos(TMath::ATan( kRDET / sanddetcomp )))/2.0;

  // Weight for solid angle and lorentz boost
  wgt_xy = sangdet * ( emrat * emrat );  // ! the weight ... normally

  // Done for all except polarized muon decay
  // in which case need to modify weight
  // (must be done in double precision)
  // BUT do this only for case of muon decay, not muon capture
  //   until beamline simulation code gets updated these generally show up as
  //   decay.ndecay == 0, but certainly not dkp_mup_nusep or dkp_mum_nusep
  // so was:
  // if ( decay.ptype  == kpdg_muplus      ||
  //      decay.ptype  == kpdg_muminus        ) {
  // now:
  if ( decay.ndecay == bsim::dkp_mup_nusep ||
       decay.ndecay == bsim::dkp_mum_nusep    ) {
    double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3], partial;

    // Boost neu neutrino to mu decay CM
    beta[0] = decay.pdpx / parent_energy;
    beta[1] = decay.pdpy / parent_energy;
    beta[2] = decay.pdpz / parent_energy;
    p_nu[0] = (xpos-decay.vx)*enu/rad;
    p_nu[1] = (ypos-decay.vy)*enu/rad;
    p_nu[2] = (zpos-decay.vz)*enu/rad;
    partial = gamma *
      (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
    partial = enu - partial/(gamma+1.0);
    // the following calculation is numerically imprecise
    // especially p_dcm_nu[2] leads to taking the difference of numbers
    //  of order ~10's and getting results of order ~0.02's
    // for g3numi we're starting with floats (ie. good to ~1 part in 10^7)
    p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*partial;
    p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*partial;
    p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*partial;
    p_dcm_nu[3] = TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
                               p_dcm_nu[1]*p_dcm_nu[1] +
                               p_dcm_nu[2]*p_dcm_nu[2] );

    // Boost parent of mu to mu production CM
    double particle_energy = decay.ppenergy;
    gamma = particle_energy/parent_mass;
    beta[0] = decay.ppdxdz * decay.pppz / particle_energy;
    beta[1] = decay.ppdydz * decay.pppz / particle_energy;
    beta[2] =                    decay.pppz / particle_energy;
    partial = gamma * ( beta[0]*decay.muparpx +
                        beta[1]*decay.muparpy +
                        beta[2]*decay.muparpz );
    partial = decay.mupare - partial/(gamma+1.0);
    p_pcm_mp[0] = decay.muparpx - beta[0]*gamma*partial;
    p_pcm_mp[1] = decay.muparpy - beta[1]*gamma*partial;
    p_pcm_mp[2] = decay.muparpz - beta[2]*gamma*partial;
    double p_pcm = TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
                                 p_pcm_mp[1]*p_pcm_mp[1] +
                                 p_pcm_mp[2]*p_pcm_mp[2] );

    const double eps = 1.0e-30;  // ? what value to use
    if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
      return 3; // mu missing parent info?
    }
    // Calc new decay angle w.r.t. (anti)spin direction
    double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
                     p_dcm_nu[1]*p_pcm_mp[1] +
                     p_dcm_nu[2]*p_pcm_mp[2] ) /
                   ( p_dcm_nu[3]*p_pcm );
    // protect against small excursions
    if ( costh >  1.0 ) costh =  1.0;
    if ( costh < -1.0 ) costh = -1.0;
    // Calc relative weight due to angle difference
    double wgt_ratio = 0.0;
    switch ( decay.ntype ) {
    case kpdg_nue:
    case kpdg_nuebar:
      wgt_ratio = 1.0 - costh;
      break;
    case kpdg_numu:
    case kpdg_numubar:
    {
      double xnu = 2.0 * enuzr / kMUMASS;
      wgt_ratio = ( (3.0-2.0*xnu )  - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);

      if ( wgt_ratio < 0.0 ) {
        std::cerr << "bsim::calcEnuWgt encountered serious problem: "
                  << " wgt_ratio " << wgt_ratio
                  << " enu " << enu << " costh " << costh << " xnu " << xnu
                  << " enuzr=decay.necm " << enuzr << " kMUMASS " << kMUMASS
                  << " norig " << decay.norig
                  << " ndecay " << decay.ndecay
                  << " ntype " << decay.ntype
                  << " ptype " << decay.ptype
                  << std::endl;
          enu    = 0;
          wgt_xy = 0;
          return 4; // bad, bad, bad calculation
      }
      break;
    }
    default:
      enu    = 0.0;
      wgt_xy = 0.0;
      return 2; // bad neutrino type for muon decay
    }

    wgt_xy = wgt_xy * wgt_ratio;

  } // ptype is muon

  return 0;
}

int bsim::calcEnuWgt	(	const bsim::Dk2Nu *	dk2nu,
		const TVector3 &	xyz,
		double &	enu,
		double &	wgt_xy
	)

alternate interface

Definition at line 337 of file calcLocationWeights.cxx.

References calcEnuWgt(), and bsim::Dk2Nu::decay.

{
  return bsim::calcEnuWgt(dk2nu->decay,xyz,enu,wgt_xy);
}

void bsim::calcLocationWeights	(	const bsim::DkMeta *	dkmeta,
		bsim::Dk2Nu *	dk2nu
	)

user interface

Definition at line 10 of file calcLocationWeights.cxx.

References calcEnuWgt(), om::cerr, bsim::Dk2Nu::decay, allTimeWatchdog::endl, bsim::DkMeta::location, bsim::Dk2Nu::nuray, make_associated_cosmic_defs::p3, fabricate::status, string, Unit(), bsim::Decay::vx, bsim::Decay::vy, and bsim::Decay::vz.

Referenced by convert_flugg(), convert_g4lbne(), convert_g4minerva(), and test_fill_dk2nu().

{
  size_t nloc = dkmeta->location.size();
  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
    // skip calculation for random location ... should already be filled
    const std::string rkey = "random decay";
    if ( dkmeta->location[iloc].name == rkey ) {
      if ( iloc != 0 ) {
        std::cerr << "calcLocationWeights \"" << rkey << "\""
                  << " isn't the 0-th entry" << std::endl;
        return; // assert(0);
      }
      if ( dk2nu->nuray.size() != 1 ) {
        std::cerr << "calcLocationWeights \"" << rkey << "\""
                  << " nuenergy[" << iloc << "] not filled" << std::endl;
        return; //assert(0);
      }
      continue;
    }
    TVector3 xyzDet(dkmeta->location[iloc].x,
                    dkmeta->location[iloc].y,
                    dkmeta->location[iloc].z);  // position to evaluate
    double enu_xy = 0;  // give a default value
    double wgt_xy = 0;  // give a default value
    int status = bsim::calcEnuWgt(dk2nu,xyzDet,enu_xy,wgt_xy);
    if ( status != 0 ) {
      std::cerr << "bsim::calcEnuWgt returned " << status << " for "
                << dkmeta->location[iloc].name << std::endl;
    }
    // with the recalculated energy compute the momentum components
    TVector3 xyzDk(dk2nu->decay.vx,dk2nu->decay.vy,dk2nu->decay.vz);  // origin of decay
    TVector3 p3 = enu_xy * (xyzDet - xyzDk).Unit();
    bsim::NuRay anuray(p3.x(), p3.y(), p3.z(), enu_xy, wgt_xy);
    dk2nu->nuray.push_back(anuray);
  }
}

Bool_t bsim::IsDefault

(

Float_t

intval

)

Definition at line 4 of file dflt.cxx.

References e, and kDfltFloat.

{
   return (TMath::Abs(fval - kDfltFloat) < 1e-4);
}

Bool_t bsim::IsDefault

(

Int_t

intval

)

Definition at line 9 of file dflt.cxx.

References kDfltInt.

{
   return (ival == kDfltInt);
}

Bool_t bsim::IsDefault

(

UInt_t

intval

)

Definition at line 14 of file dflt.cxx.

References kDfltUInt.

{
   return (uival == kDfltUInt);
}

Bool_t bsim::IsDefault

(

Double_t

intval

)

Definition at line 19 of file dflt.cxx.

References e, and kDfltDouble.

{
   return (TMath::Abs(dval - kDfltDouble) < 1e-4);
}

Bool_t bsim::IsDefault

(

Bool_t

intval

)

Definition at line 24 of file dflt.cxx.

References kDfltBool.

{
   return (bval == kDfltBool);
}

void bsim::printWeightLocations

(

std::vector< bsim::Location > &

locations

)

print the locations

Definition at line 83 of file readWeightLocations.cxx.

References om::cout, submit_hadd::l, and fhicl::detail::nl().

Referenced by ConvertInit(), and printWeightLocations().

{
  size_t nl = locations.size();
  std::cout << nl << " locations:\n";
  for ( size_t l = 0; l < nl; ++l ) {
    std::cout << " [" << std::setw(2) << l << "] " << locations[l] << "\n";
  }
}

void bsim::printWeightLocations

(

bsim::DkMeta *

dkmeta

)

a variant that prints the locations in the dkmeta object

Definition at line 92 of file readWeightLocations.cxx.

References bsim::DkMeta::location, and printWeightLocations().

{
  bsim::printWeightLocations(dkmeta->location);
}

void bsim::readWeightLocations	(	std::string	locfilename,
		std::vector< bsim::Location > &	locations
	)

Read a text file that contains a header line followed by quartets of "<xpos> <ypos> <zpos> <text string>" on separate lines. Fill the supplied vectors. Trim off leading/trailing blanks and quotes (single/double) from the string. Convention has it that positions are given in (cm).

Definition at line 17 of file readWeightLocations.cxx.

References plot_validation_datamc::c, MECModelEnuComparisons::i, path, string, tmp, submit_syst::x, submit_syst::y, and test::z.

Referenced by ConvertInit(), readWeightLocations(), test_fill_dk2nu(), and test_read_locations().

{

  const char* path = gSystem->ExpandPathName(locfilename.c_str());
  std::ifstream locfile(path);

  int iline=0;

  char comment_buffer[1000];

  // read lines
  char tmp[1001];
  size_t tmplen = sizeof(tmp);
  while ( ! locfile.eof() ) {
    char c;
    while ( ( c = locfile.get() ) == ' ' ) {}; // eat leading spaces
    if ( c == '#' ) {
      // comment_buffer line
      locfile.getline(comment_buffer,sizeof(comment_buffer));
      continue;
    }
    locfile.putback(c);
    double x, y, z;
    locfile >> x >> y >> z;
    locfile.getline(tmp,tmplen-1);
    size_t i = locfile.gcount();
    // make sure the c-string is null terminated
    size_t inull = i;
    //if ( inull < 0 ) inull = 0;
    if ( inull > tmplen-1 ) inull = tmplen-1;
    tmp[inull] = '\0';
    std::string name(tmp);
    // ignore leading & trailing blanks (and any single/double quotes)
    size_t ilast  = name.find_last_not_of(" \t\n'\"");
    name.erase(ilast+1,std::string::npos);  // trim tail
    size_t ifirst = name.find_first_not_of(" \t\n'\"");
    name.erase(0,ifirst);  // trim head

    ++iline;
    if ( ! locfile.good() ) {
      //if ( verbose)
      //  std::cout << "stopped reading on line " << iline << std::endl;
      break;
    }
    bsim::Location alocation(x,y,z,name);
    locations.push_back(alocation);
  }

}

void bsim::readWeightLocations	(	std::string	locfilename,
		bsim::DkMeta *	dkmeta
	)

a variant that will fill the dkmeta object

read & print the locations where weights are to be calculated

make an entry for the random decay

read and parse the text file for additional positions use the vector version

Definition at line 69 of file readWeightLocations.cxx.

References bsim::DkMeta::location, and readWeightLocations().

{
  ///  read & print the locations where weights are to be calculated

  std::vector<bsim::Location>& locations = dkmeta->location;
  /// make an entry for the random decay
  bsim::Location rndmloc(0,0,0,"random decay");
  locations.push_back(rndmloc);

  /// read and parse the text file for additional positions
  /// use the vector version
  readWeightLocations(locfilename, locations);
}

Variable Documentation

const Bool_t bsim::kDfltBool = false

static

Definition at line 16 of file dflt.h.

Referenced by IsDefault().

const Bool_t bsim::kDfltBool = false

static

Definition at line 16 of file dflt.h.

const Double_t bsim::kDfltDouble = -9999.99

static

Definition at line 15 of file dflt.h.

const Double_t bsim::kDfltDouble = -9999.99

static

Definition at line 15 of file dflt.h.

Referenced by IsDefault().

const Float_t bsim::kDfltFloat = -9999.99

static

Definition at line 12 of file dflt.h.

Referenced by IsDefault().

const Float_t bsim::kDfltFloat = -9999.99

static

Definition at line 12 of file dflt.h.

const Int_t bsim::kDfltInt = -9999

static

Definition at line 13 of file dflt.h.

Referenced by IsDefault().

const Int_t bsim::kDfltInt = -9999

static

Definition at line 13 of file dflt.h.

const UInt_t bsim::kDfltUInt = 9999

static

Definition at line 14 of file dflt.h.

Referenced by IsDefault().

const UInt_t bsim::kDfltUInt = 9999

static

Definition at line 14 of file dflt.h.