GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
genie::hnl::FluxCreator Class Reference

Calculates HNL production kinematics & production vertex. Is a concrete implementation of the FluxRecordVisitorI interface. More...

#include <HNLFluxCreator.h>

Inheritance diagram for genie::hnl::FluxCreator:
Inheritance graph
[legend]
Collaboration diagram for genie::hnl::FluxCreator:
Collaboration graph
[legend]

Public Member Functions

 FluxCreator ()
 
 FluxCreator (string name)
 
 FluxCreator (string name, string config)
 
 ~FluxCreator ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
void SetInputFluxPath (std::string finpath) const
 
void SetGeomFile (string geomfile) const
 
int GetNFluxEntries () const
 
void SetFirstFluxEntry (int iFirst) const
 
FluxContainer RetrieveFluxInfo () const
 
std::vector< double > GetB2UTranslation () const
 
std::vector< double > GetB2URotation () const
 
std::vector< double > GetDetOffset () const
 
std::vector< double > GetDetRotation () const
 
- Public Member Functions inherited from genie::hnl::FluxRecordVisitorI
virtual ~FluxRecordVisitorI ()
 
virtual void SetGeomFile (std::string geomfile) const =0
 
- Public Member Functions inherited from genie::hnl::GeomRecordVisitorI
virtual ~GeomRecordVisitorI ()
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
void SetUsingRootGeom (bool IsUsingRootGeom) const
 
void ImportBoundingBox (TGeoBBox *box) const
 
void SetCurrentEntry (int iCurr) const
 
FluxContainer MakeTupleFluxEntry (int iEntry, std::string finpath) const
 
void FillNonsense (int iEntry, genie::hnl::FluxContainer &gnmf) const
 
void OpenFluxInput (std::string finpath) const
 
void InitialiseTree () const
 
void InitialiseMeta () const
 
TLorentzVector HNLEnergy (genie::hnl::HNLProd_t hnldm, TLorentzVector p4par) const
 
TVector3 PointToRandomPointInBBox () const
 
void ReadBRs () const
 
std::map
< genie::hnl::HNLProd_t,
double > 		GetProductionProbs (int parPDG) const
 
void MakeBBox () const
 
TVector3 ApplyUserRotation (TVector3 vec, bool doBackwards=false) const
 
TVector3 ApplyUserRotation (TVector3 vec, TVector3 oriVec, std::vector< double > rotVec, bool doBackwards=false) const
 
double CalculateDetectorAcceptanceSAA (TVector3 detO) const
 
double CalculateAcceptanceCorrection (TLorentzVector p4par, TLorentzVector p4HNL, double SMECM, double zm, double zp) const
 
double GetAngDeviation (TLorentzVector p4par, TVector3 detO, bool seekingMax) const
 
void GetAngDeviation (TLorentzVector p4par, TVector3 detO, double &zm, double &zp) const
 
double CalculateAreaNormalisation ()
 
std::string CheckGeomPoint (Double_t x, Double_t y, Double_t z) const
 

Static Private Member Functions

static double labangle (double *x, double *par)
 

Private Attributes

std::string fCurrPath = ""
 
bool fPathLoaded = false
 
int iCurrEntry = 0
 
int fFirstEntry = 0
 
int fNEntries = 0
 
std::map
< genie::hnl::HNLProd_t,
double > 		dynamicScores
 
std::map
< genie::hnl::HNLProd_t,
double > 		dynamicScores_pion
 
std::map
< genie::hnl::HNLProd_t,
double > 		dynamicScores_kaon
 
std::map
< genie::hnl::HNLProd_t,
double > 		dynamicScores_muon
 
std::map
< genie::hnl::HNLProd_t,
double > 		dynamicScores_neuk
 
double BR_pi2mu
 
double BR_pi2e
 
double BR_K2mu
 
double BR_K2e
 
double BR_K3mu
 
double BR_K3e
 
double BR_K03mu
 
double BR_K03e
 
bool isParentOnAxis = true
 
TGeoVolume * fTopVol = 0
 
string fGeomFile = ""
 
bool fIsUsingRootGeom = false
 
TChain * ctree = 0
 
TChain * cmeta = 0
 
double fMass
 
std::vector< double > fU4l2s
 
bool fIsMajorana
 
double fLx
 
double fLy
 
double fLz
 
double fLxR
 
double fLyR
 
double fLzR
 
std::vector< double > fB2UTranslation
 
std::vector< double > fB2URotation
 
std::vector< double > fDetRotation
 
std::vector< double > fDetOffset
 
double fCx
 
double fCy
 
double fCz
 
double fAx1
 
double fAz
 
double fAx2
 
double fBx1
 
double fBz
 
double fBx2
 
double fDx
 
double fDy
 
double fDz
 
bool doPol = true
 
bool fixPol = false
 
double fLPx
 
double fLPy
 
double fLPz
 
double fLPE
 
std::vector< double > fFixedPolarisation
 
int fLepPdg
 
int fNuPdg
 
double parentMass
 
double parentMomentum
 
double parentEnergy
 
double fECM
 
double fSMECM
 
double fZm
 
double fZp
 
int fProdChan
 
int fNuProdChan
 
TVector3 fTargetPoint
 
TVector3 fTargetPointUser
 
double potnum
 N POT for this SM-v. More...
 
int decay_ptype
 PDG code of parent. More...
 
double decay_vx
 
double decay_vy
 
double decay_vz
 coordinates of prod vtx [cm] More...
 
double decay_pdpx
 
double decay_pdpy
 
double decay_pdpz
 final parent momentum [GeV] More...
 
double decay_necm
 SM v CM energy [GeV]. More...
 
double decay_nimpwt
 Importance weight from beamsim. More...
 
int arSize
 
int anArSize
 
int trArSize
 
int djob
 
double ppvx
 
double ppvy
 
double ppvz
 
int decay_norig
 
int decay_ndecay
 
int decay_ntype
 
double decay_ppdxdz
 
double decay_ppdydz
 
double decay_pppz
 
double decay_ppenergy
 
int decay_ppmedium
 
double decay_muparpx
 
double decay_muparpy
 
double decay_muparpz
 
double decay_mupare
 
double nuray_px [maxArray]
 
double nuray_py [maxArray]
 
double nuray_pz [maxArray]
 
double nuray_E [maxArray]
 
double nuray_wgt [maxArray]
 
int ancestor_pdg [maxArray]
 
double ancestor_startx [maxArray]
 
double ancestor_starty [maxArray]
 
double ancestor_startz [maxArray]
 
double ancestor_startt [maxArray]
 
double ancestor_startpx [maxArray]
 
double ancestor_startpy [maxArray]
 
double ancestor_startpz [maxArray]
 
double ancestor_stoppx [maxArray]
 
double ancestor_stoppy [maxArray]
 
double ancestor_stoppz [maxArray]
 
double ancestor_polx [maxArray]
 
double ancestor_poly [maxArray]
 
double ancestor_polz [maxArray]
 
double ancestor_pprodpx [maxArray]
 
double ancestor_pprodpy [maxArray]
 
double ancestor_pprodpz [maxArray]
 
int ancestor_nucleus [maxArray]
 
char ancestor_proc [maxArray *maxC]
 
char ancestor_ivol [maxArray *maxC]
 
char ancestor_imat [maxArray *maxC]
 
double tgtexit_tvx
 
double tgtexit_tvy
 
double tgtexit_tvz
 
double tgtexit_tpx
 
double tgtexit_tpy
 
double tgtexit_tpz
 
int tgtexit_tptype
 
int tgtexit_tgen
 
double traj_trkx [maxArray]
 
double traj_trky [maxArray]
 
double traj_trkz [maxArray]
 
double traj_trkpx [maxArray]
 
double traj_trkpy [maxArray]
 
double traj_trkpz [maxArray]
 
int job
 beamsim MC job number More...
 
double pots
 how many pot in this job? More...
 
int mArSize
 
char beamsim [maxC]
 
char physics [maxC]
 
char physcuts [maxC]
 
char tgtcfg [maxC]
 
char horncfg [maxC]
 
char dkvolcfg [maxC]
 
double beam0x
 
double beam0y
 
double beam0z
 
double beamhwidth
 
double beamvwidth
 
double beamdxdz
 
double beamdydz
 
double location_x [maxArray]
 
double location_y [maxArray]
 
double location_z [maxArray]
 
char location_name [maxArray *maxC]
 
genie::hnl::FluxContainer fGnmf
 
double POTScaleWeight
 
std::vector< double > fScales
 
bool fDoingOldFluxCalc = false
 
bool fRerollPoints = false
 
double fRadius
 
bool fSupplyingBEAM = false
 
bool fIsConfigLoaded = false
 
bool fUsingDk2nu = true
 
string fFinPath
 
string fProdHist
 
TH1D * fSpectrum = 0
 
TH1D * fIntegrals = 0
 

Static Private Attributes

static const int maxArray = 30
 
static const int maxC = 100
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::hnl::FluxRecordVisitorI
 FluxRecordVisitorI ()
 
 FluxRecordVisitorI (string name)
 
 FluxRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::hnl::GeomRecordVisitorI
 GeomRecordVisitorI ()
 
 GeomRecordVisitorI (string name)
 
 GeomRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
Registry * ExtractLocalConfig (const Registry &in) const
 
Registry * ExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Calculates HNL production kinematics & production vertex. Is a concrete implementation of the FluxRecordVisitorI interface.

This is a module for GENIE to read in hadron flux ntuples and construct HNL fluxes on the fly.

Core loop: + Open flux entry

Author
John Plows komni.nosp@m.nos-.nosp@m.john..nosp@m.plow.nosp@m.s@phy.nosp@m.sics.nosp@m..ox.a.nosp@m.c.uk
Created:
April 25th, 2022
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 94 of file HNLFluxCreator.h.

Constructor & Destructor Documentation

FluxCreator::FluxCreator ( )

Definition at line 18 of file HNLFluxCreator.cxx.

18  :
19  FluxRecordVisitorI("genie::hnl::FluxCreator")
20 {
21 
22 }
FluxCreator::FluxCreator ( string  name)

Definition at line 24 of file HNLFluxCreator.cxx.

24  :
25  FluxRecordVisitorI(name)
26 {
27 
28 }
tests::name
const char * name
Definition: gtestGAtmoFlux.cxx:195
FluxCreator::FluxCreator ( string  name,
string  config 
)

Definition at line 30 of file HNLFluxCreator.cxx.

30  :
31  FluxRecordVisitorI(name, config)
32 {
33 
34 }
tests::name
const char * name
Definition: gtestGAtmoFlux.cxx:195
FluxCreator::~FluxCreator ( )

Definition at line 36 of file HNLFluxCreator.cxx.

37 {
38 
39 }

Member Function Documentation

TVector3 FluxCreator::ApplyUserRotation ( TVector3  vec,
bool  doBackwards = false 
) const
private

Definition at line 1959 of file HNLFluxCreator.cxx.

References fAx1, fAx2, and fAz.

Referenced by GetAngDeviation(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

1960 {
1961  double vx = vec.X(), vy = vec.Y(), vz = vec.Z();
1962 
1963  double Ax2 = ( doBackwards ) ? -fAx2 : fAx2;
1964  double Az = ( doBackwards ) ? -fAz : fAz;
1965  double Ax1 = ( doBackwards ) ? -fAx1 : fAx1;
1966 
1967  // Ax2 first
1968  double x = vx, y = vy, z = vz;
1969  vy = y * std::cos( Ax2 ) - z * std::sin( Ax2 );
1970  vz = y * std::sin( Ax2 ) + z * std::cos( Ax2 );
1971  y = vy; z = vz;
1972  // then Az
1973  vx = x * std::cos( Az ) - y * std::sin( Az );
1974  vy = x * std::sin( Az ) + y * std::cos( Az );
1975  x = vx; y = vy;
1976  // Ax1 last
1977  vy = y * std::cos( Ax1 ) - z * std::sin( Ax1 );
1978  vz = y * std::sin( Ax1 ) + z * std::cos( Ax1 );
1979 
1980  TVector3 nvec( vx, vy, vz );
1981  return nvec;
1982 }
TVector3 FluxCreator::ApplyUserRotation ( TVector3  vec,
TVector3  oriVec,
std::vector< double >  rotVec,
bool  doBackwards = false 
) const
private

Definition at line 1984 of file HNLFluxCreator.cxx.

1985 {
1986  double vx = vec.X(), vy = vec.Y(), vz = vec.Z();
1987  double ox = oriVec.X(), oy = oriVec.Y(), oz = oriVec.Z();
1988 
1989  vx -= ox; vy -= oy; vz -= oz; // make this rotation about detector origin
1990 
1991  assert( rotVec.size() == 3 ); // want 3 Euler angles, otherwise this is unphysical.
1992  double Ax2 = ( doBackwards ) ? -rotVec.at(2) : rotVec.at(2);
1993  double Az = ( doBackwards ) ? -rotVec.at(1) : rotVec.at(1);
1994  double Ax1 = ( doBackwards ) ? -rotVec.at(0) : rotVec.at(0);
1995 
1996  // Ax2 first
1997  double x = vx, y = vy, z = vz;
1998  vy = y * std::cos( Ax2 ) - z * std::sin( Ax2 );
1999  vz = y * std::sin( Ax2 ) + z * std::cos( Ax2 );
2000  y = vy; z = vz;
2001  // then Az
2002  vx = x * std::cos( Az ) - y * std::sin( Az );
2003  vy = x * std::sin( Az ) + y * std::cos( Az );
2004  x = vx; y = vy;
2005  // Ax1 last
2006  vy = y * std::cos( Ax1 ) - z * std::sin( Ax1 );
2007  vz = y * std::sin( Ax1 ) + z * std::cos( Ax1 );
2008 
2009  // back to beam frame
2010  vx += ox; vy += oy; vz += oz;
2011  TVector3 nvec( vx, vy, vz );
2012  return nvec;
2013 }
double FluxCreator::CalculateAcceptanceCorrection ( TLorentzVector  p4par,
TLorentzVector  p4HNL,
double  SMECM,
double  zm,
double  zp 
) const
private

Definition at line 1803 of file HNLFluxCreator.cxx.

References e, and labangle().

Referenced by MakeTupleFluxEntry().

1805 {
1806  /*
1807  * This method calculates HNL acceptance by taking into account the collimation effect
1808  * HNL are massive so Lorentz boost from parent CM ==> lab is more effective
1809  * This means that, given a desired range of lab-frame emission angles, there are
1810  * more rest-frame emission angles that map into this range.
1811  * Find the measure of the rest-frame that maps onto the allowed lab-frame angles
1812  * and return the ratio over the relevant measure for a SM neutrino
1813  */
1814 
1815  assert( zm >= 0.0 && zp >= zm );
1816  if( zp == zm ) return 1.0;
1817 
1818  double M = p4HNL.M();
1819  if( M < 1.0e-3 ) return 1.0;
1820 
1821  TF1 * fHNL = ( TF1* ) gROOT->GetListOfFunctions()->FindObject( "fHNL" );
1822  if( !fHNL ){ fHNL = new TF1( "fHNL", labangle, 0.0, 180.0, 6 ); }
1823  fHNL->SetParameter( 0, p4par.E() );
1824  fHNL->SetParameter( 1, p4par.Px() );
1825  fHNL->SetParameter( 2, p4par.Py() );
1826  fHNL->SetParameter( 3, p4par.Pz() );
1827  fHNL->SetParameter( 4, p4HNL.P() );
1828  fHNL->SetParameter( 5, p4HNL.E() );
1829 
1830  double ymax = fHNL->GetMaximum(), xmax = fHNL->GetMaximumX();
1831  double range1 = 0.0;
1832 
1833  if( fHNL->GetMinimum() >= fHNL->GetMaximum() ) return 1.0; // bail on constant function
1834 
1835  if( zm < fHNL->GetMinimum() ){ // really good collimation, ignore checks on zm
1836  double z0 = fHNL->GetMinimum();
1837  if( ymax > zp && xmax < 180.0 ){ // >=2 pre-images, add them together
1838  int nPreim = 0;
1839 
1840  // RETHERE: Make this more sophisticated! Assumes 2 preimages, 1 before and 1 after max
1841  double xl1 = fHNL->GetX( z0, 0.0, xmax );
1842  double xh1 = fHNL->GetX( zp, 0.0, xmax );
1843  double xl2 = fHNL->GetX( z0, xmax, 180.0 );
1844  double xh2 = fHNL->GetX( zp, xmax, 180.0 );
1845 
1846  range1 += std::abs( xl1 - xh1 ) + std::abs( xh2 - xl2 ); nPreim = 2;
1847  } else if( ymax > zp && xmax == 180.0 ){ // 1 pre-image, SMv-like case
1848  double xl = fHNL->GetX( z0 ), xh = fHNL->GetX( zp );
1849  range1 = std::abs( xh - xl );
1850 
1851  } else if( ymax <= zp ){ // 1 pre-image but all emissions reach detector
1852  range1 = 180.0;
1853 
1854  }
1855  } else { // not so good collimation, enforce checks on zm
1856  if( ymax <= zm ){ // 0 pre-images
1857  return 0.0;
1858  } else if( ymax > zp && xmax < 180.0 ){ // >=2 pre-images, add them together
1859  int nPreim = 0;
1860 
1861  // RETHERE: Make this more sophisticated! Assumes 2 preimages, 1 before and 1 after max
1862  double xl1 = fHNL->GetX( zm, 0.0, xmax );
1863  double xh1 = fHNL->GetX( zp, 0.0, xmax );
1864  double xl2 = fHNL->GetX( zm, xmax, 180.0 );
1865  double xh2 = fHNL->GetX( zp, xmax, 180.0 );
1866 
1867  range1 += std::abs( xl1 - xh1 ) + std::abs( xh2 - xl2 ); nPreim = 2;
1868  } else if( ymax > zp && xmax == 180.0 ){ // 1 pre-image, SMv-like case
1869  double xl = fHNL->GetX( zm ), xh = fHNL->GetX( zp );
1870  range1 = std::abs( xh - xl );
1871 
1872  } else if( zm < ymax && ymax <= zp ){ // 1 pre-image
1873  double xl = fHNL->GetX( zm, 0., xmax ), xh = fHNL->GetX( zm, xmax, 180.0 );
1874  range1 = std::abs( xh - xl );
1875  }
1876  }
1877 
1878  TF1 * fSMv = ( TF1* ) gROOT->GetListOfFunctions()->FindObject( "fSMv" );
1879  if( !fSMv ){
1880  fSMv = new TF1( "fSMv", labangle, 0.0, 180.0, 6 );
1881  }
1882  fSMv->SetParameter( 0, p4par.E() );
1883  fSMv->SetParameter( 1, p4par.Px() );
1884  fSMv->SetParameter( 2, p4par.Py() );
1885  fSMv->SetParameter( 3, p4par.Pz() );
1886  fSMv->SetParameter( 4, SMECM );
1887  fSMv->SetParameter( 5, SMECM );
1888  double range2 = -1.0;
1889 
1890  if( fSMv->GetMaximum() == fSMv->GetMinimum() ) return 1.0; // bail
1891 
1892  // SMv deviates more from parent than HNL due to masslessness. This means a larger minimum of labangle
1893  // Sometimes the angle is so small, that this calculation fails as there is not SMv preimage to compare
1894  // with. Default to estimating dx/dz * dy/dz ratio in that case.
1895 
1896  if( fSMv->GetMinimum() < zp ){
1897  if( fSMv->GetMinimum() < zm ){
1898  range2 = std::abs( fSMv->GetX( zp ) - fSMv->GetX( zm ) );
1899  } else { // due to monotonicity all of [0.0, fSMv->GetX(zp)] is good
1900  range2 = fSMv->GetX( zp );
1901  }
1902  if( range2 < 1.0e-6 || range1 / range2 > 10.0 ) return 1.0; // sometimes this happens, gotta bail!
1903  } else { // can't decide based on SMv analytically.
1904  TVector3 bv = p4par.BoostVector();
1905 
1906  TLorentzVector vcx( p4HNL.P(), 0.0, 0.0, p4HNL.E() ),
1907  vcy( 0.0, p4HNL.P(), 0.0, p4HNL.E() ),
1908  vcz( 0.0, 0.0, p4HNL.P(), p4HNL.E() );
1909  vcx.Boost( bv ); vcy.Boost( bv ); vcz.Boost( bv );
1910 
1911  TLorentzVector vsx( SMECM, 0.0, 0.0, SMECM ),
1912  vsy( 0.0, SMECM, 0.0, SMECM ),
1913  vsz( 0.0, 0.0, SMECM, SMECM );
1914  vsx.Boost( bv ); vsy.Boost( bv ); vsz.Boost( bv );
1915 
1916  double xpart = std::abs( ( vcx.X() / vcz.Z() ) / ( vsx.X() / vsz.Z() ) );
1917  double ypart = std::abs( ( vcy.Y() / vcz.Z() ) / ( vsy.Y() / vsz.Z() ) );
1918 
1919  return 1.0 / ( xpart * ypart );
1920  }
1921 
1922  assert( range2 > 0.0 );
1923 
1924  return range1 / range2;
1925 
1926 }
genie::hnl::FluxCreator::labangle
static double labangle(double *x, double *par)
Definition: HNLFluxCreator.cxx:1928
e
const double e
Definition: gUpMuFluxGen.cxx:165
double FluxCreator::CalculateAreaNormalisation ( )
private

Definition at line 2023 of file HNLFluxCreator.cxx.

References kRDET.

2024 {
2025  // for now this is just a square of length kRDET
2026  // returns 1 / area
2027  return 1.0 / ( kRDET * kRDET );
2028 }
kRDET
const double kRDET
Definition: HNLFluxCreator.h:85
double FluxCreator::CalculateDetectorAcceptanceSAA ( TVector3  detO) const
private

Definition at line 2015 of file HNLFluxCreator.cxx.

References kRDET, and genie::units::rad.

Referenced by MakeTupleFluxEntry().

2016 {
2017  // sang is solid-angle / 4pi
2018  double rad = std::sqrt( detO.X() * detO.X() + detO.Y() * detO.Y() + detO.Z() * detO.Z() );
2019  double sang = 1.0 - TMath::Cos( TMath::ATan( kRDET / rad ) ); sang *= 0.5;
2020  return sang;
2021 }
genie::units::rad
static constexpr double rad
Definition: Units.h:164
kRDET
const double kRDET
Definition: HNLFluxCreator.h:85
std::string FluxCreator::CheckGeomPoint ( Double_t  x,
Double_t  y,
Double_t  z 
) const
private

Definition at line 2136 of file HNLFluxCreator.cxx.

Referenced by GetAngDeviation(), and PointToRandomPointInBBox().

2137 {
2138  Double_t point[3];
2139  Double_t local[3];
2140  point[0] = x;
2141  point[1] = y;
2142  point[2] = z;
2143  TGeoVolume *vol = gGeoManager->GetTopVolume();
2144  TGeoNode *node = gGeoManager->FindNode(point[0], point[1], point[2]);
2145  gGeoManager->MasterToLocal(point, local);
2146  return gGeoManager->GetPath();
2147 }
void FluxCreator::Configure ( const Registry &  config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 2030 of file HNLFluxCreator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

2031 {
2032  Algorithm::Configure(config);
2033  this->LoadConfig();
2034 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void FluxCreator::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 2036 of file HNLFluxCreator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

2037 {
2038  Algorithm::Configure(config);
2039  this->LoadConfig();
2040 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void FluxCreator::FillNonsense ( int  iEntry,
genie::hnl::FluxContainer &  gnmf 
) const
private

Definition at line 703 of file HNLFluxCreator.cxx.

References genie::hnl::FluxContainer::accCorr, genie::hnl::FluxContainer::acceptance, genie::hnl::FluxContainer::boostCorr, genie::hnl::FluxContainer::Ecm, genie::hnl::FluxContainer::evtno, genie::hnl::FluxContainer::lepPdg, genie::hnl::FluxContainer::nimpwt, genie::hnl::FluxContainer::nuEcm, genie::hnl::FluxContainer::nuPdg, genie::hnl::FluxContainer::nuProdChan, genie::hnl::FluxContainer::p4, genie::hnl::FluxContainer::p4User, genie::hnl::FluxContainer::parp4, genie::hnl::FluxContainer::parp4User, genie::hnl::FluxContainer::parPdg, genie::hnl::FluxContainer::pdg, genie::hnl::FluxContainer::polz, genie::hnl::FluxContainer::prodChan, genie::hnl::FluxContainer::startPoint, genie::hnl::FluxContainer::startPointUser, genie::hnl::FluxContainer::targetPoint, genie::hnl::FluxContainer::targetPointUser, genie::hnl::FluxContainer::XYWgt, genie::hnl::FluxContainer::zetaMinus, and genie::hnl::FluxContainer::zetaPlus.

Referenced by MakeTupleFluxEntry().

704 {
705  gnmf.evtno = iEntry;
706 
707  gnmf.pdg = -9999;
708  gnmf.parPdg = -9999;
709  gnmf.lepPdg = -9999;
710  gnmf.nuPdg = -9999;
711 
712  gnmf.prodChan = -9999;
713  gnmf.nuProdChan = -9999;
714 
715  gnmf.startPoint.SetXYZ(-9999.9, -9999.9, -9999.9);
716  gnmf.targetPoint.SetXYZ(-9999.9, -9999.9, -9999.9);
717  gnmf.startPointUser.SetXYZ(-9999.9, -9999.9, -9999.9);
718  gnmf.targetPointUser.SetXYZ(-9999.9, -9999.9, -9999.9);
719 
720  gnmf.polz.SetXYZ(-9999.9, -9999.9, -9999.9);
721 
722  gnmf.p4.SetPxPyPzE(-9999.9, -9999.9, -9999.9, -9999.9);
723  gnmf.parp4.SetPxPyPzE(-9999.9, -9999.9, -9999.9, -9999.9);
724  gnmf.p4User.SetPxPyPzE(-9999.9, -9999.9, -9999.9, -9999.9);
725  gnmf.parp4User.SetPxPyPzE(-9999.9, -9999.9, -9999.9, -9999.9);
726 
727  gnmf.Ecm = -9999.9;
728  gnmf.nuEcm = -9999.9;
729  gnmf.XYWgt = -9999.9;
730  gnmf.boostCorr = -9999.9;
731  gnmf.accCorr = -9999.9;
732  gnmf.zetaMinus = -9999.9;
733  gnmf.zetaPlus = -9999.9;
734  gnmf.acceptance = -9999.9;
735 
736  gnmf.nimpwt = -9999.9;
737 
738  return;
739 }
genie::hnl::FluxContainer::parp4
TLorentzVector parp4
parent momentum at HNL production in NEAR coords [GeV/c]
Definition: HNLFluxContainer.h:94
genie::hnl::FluxContainer::Ecm
double Ecm
Parent rest-frame energy of HNL [GeV].
Definition: HNLFluxContainer.h:98
genie::hnl::FluxContainer::zetaMinus
double zetaMinus
minimum angular deviation from parent momentum to reach detector [deg]
Definition: HNLFluxContainer.h:105
genie::hnl::FluxContainer::nuEcm
double nuEcm
Parent rest-frame energy of equivalent SM neutrino [GeV].
Definition: HNLFluxContainer.h:99
genie::hnl::FluxContainer::startPoint
TVector3 startPoint
parent decay vertex in NEAR coords [m]
Definition: HNLFluxContainer.h:85
genie::hnl::FluxContainer::lepPdg
int lepPdg
PDG code of lepton produced with HNL on parent decay.
Definition: HNLFluxContainer.h:79
genie::hnl::FluxContainer::polz
TVector3 polz
HNL polarisation vector, in HNL rest frame, in NEAR coords.
Definition: HNLFluxContainer.h:91
genie::hnl::FluxContainer::parPdg
int parPdg
parent PDG code
Definition: HNLFluxContainer.h:78
genie::hnl::FluxContainer::zetaPlus
double zetaPlus
maximum angular deviation from parent momentum to reach detector [deg]
Definition: HNLFluxContainer.h:106
genie::hnl::FluxContainer::startPointUser
TVector3 startPointUser
parent decay vertex in USER coords [m]
Definition: HNLFluxContainer.h:87
genie::hnl::FluxContainer::parp4User
TLorentzVector parp4User
parent momentum at HNL production in USER coords [GeV/c]
Definition: HNLFluxContainer.h:96
genie::hnl::FluxContainer::p4
TLorentzVector p4
HNL momentum in NEAR coords [GeV/c].
Definition: HNLFluxContainer.h:93
genie::hnl::FluxContainer::p4User
TLorentzVector p4User
HNL momentum in USER coords [GeV/c].
Definition: HNLFluxContainer.h:95
genie::hnl::FluxContainer::prodChan
int prodChan
Decay mode that produced HNL.
Definition: HNLFluxContainer.h:82
genie::hnl::FluxContainer::acceptance
double acceptance
full acceptance == XYWgt * boostCorr^2 * accCorr
Definition: HNLFluxContainer.h:108
genie::hnl::FluxContainer::targetPointUser
TVector3 targetPointUser
point in detector HNL is forced towards in USER coords [m]
Definition: HNLFluxContainer.h:88
genie::hnl::FluxContainer::targetPoint
TVector3 targetPoint
point in detector HNL is forced towards in NEAR coords [m]
Definition: HNLFluxContainer.h:86
genie::hnl::FluxContainer::nimpwt
double nimpwt
Weight of parent.
Definition: HNLFluxContainer.h:110
genie::hnl::FluxContainer::accCorr
double accCorr
acceptance correction (collimation effect. SM v == 1)
Definition: HNLFluxContainer.h:104
genie::hnl::FluxContainer::nuPdg
int nuPdg
PDG code of SM neutrino that would have been produced.
Definition: HNLFluxContainer.h:80
genie::hnl::FluxContainer::XYWgt
double XYWgt
geometric acceptance (angular size of detector in parent rest frame)
Definition: HNLFluxContainer.h:101
genie::hnl::FluxContainer::nuProdChan
int nuProdChan
Decay mode that would have produced SM neutrino.
Definition: HNLFluxContainer.h:83
genie::hnl::FluxContainer::boostCorr
double boostCorr
boost correction wrt parent rest-frame (ELAB = ECM * boostCorr)
Definition: HNLFluxContainer.h:102
double FluxCreator::GetAngDeviation ( TLorentzVector  p4par,
TVector3  detO,
bool  seekingMax 
) const
private

Definition at line 1348 of file HNLFluxCreator.cxx.

References fLx, fLy, genie::RandomGen::Instance(), genie::controls::kASmallNum, genie::constants::kPi, LOG, pERROR, and genie::RandomGen::RndGen().

Referenced by MakeTupleFluxEntry().

1349 {
1350  TVector3 ppar = p4par.Vect(); assert( ppar.Mag() > 0.0 );
1351  TVector3 pparUnit = ppar.Unit();
1352  // let face be planar and perpendicular to vector Q
1353  // assuming Q = ( 0, 0, 1 ) == face perpendicular to z
1354  double Qx = 0.0, Qy = 0.0, Qz = 1.0;
1355  // plane: Qx . (x-xC) + Qy . (y-yC) + Qz . (z-zC) = 0
1356  // line: r(t) - r(D) = t * ppar
1357  // V0 \in plane && line
1358  // ==> Qx * ( x(t) - x(C) ) + Qy * ( y(t) - y(C) ) + Qz * ( z(t) - z(C) ) = 0
1359  // ==> t = ( \vec{Q} \cdot \vec{detO} ) / ( \vec{Q} \cdot \vec{ppar} )
1360  double nterm = Qx * detO.X() + Qy * detO.Y() + Qz * detO.Z();
1361  double dterm = Qx * ppar.X() + Qy * ppar.Y() + Qz * ppar.Z();
1362  double t = nterm / dterm;
1363  double x_incp = t * pparUnit.X(), y_incp = t * pparUnit.Y(), z_incp = t * pparUnit.Z();
1364 
1365  // sweep over plane perp to ppar, passing through centre, and calc intersection point
1366  // special case: parent is perfectly on axis so hits detector centre
1367  TVector3 IPdev( detO.X() - x_incp, detO.Y() - y_incp, detO.Z() - z_incp );
1368  bool parentHitsCentre = ( IPdev.Mag() < controls::kASmallNum );
1369 
1370  // see assumption about Q
1371  // to fix probably with a rotation of fLx, fLy by Euler angles onto Q-plane?
1372  // line: r(t) - r(incp) = t * IPdev
1373  // assume square face
1374  double ttx = ( IPdev.X() != 0.0 ) ? fLx / std::abs( IPdev.X() ) : 99999.9;
1375  double tty = ( IPdev.Y() != 0.0 ) ? fLy / std::abs( IPdev.Y() ) : 99999.9;
1376  double tt = std::max( ttx, tty ); // this defines how much the least sweep goes
1377  TVector3 atilde( tt * IPdev.X(), tt * IPdev.Y(), tt * IPdev.Z() );
1378 
1379  double fLT = IPdev.Mag();
1380  double dist = atilde.Mag();
1381 
1382  assert( fLT > 0.0 );
1383  double detRadius = std::max( fLx, fLy ) / 2.0;
1384 
1385  if( parentHitsCentre ){
1386  // calculate angles for four points and return largest (smallest) of them
1387 
1388  // randomly select a phi to go with, make 2 perpendicular vectors from it
1389  double phi = ( RandomGen::Instance()->RndGen() ).Uniform( 0., 2.0 * constants::kPi );
1390  TVector3 r1VecPrim( -pparUnit.Y(), pparUnit.X(), 0.0 ); // perp to ppar == on plane
1391  TVector3 r1Vec = r1VecPrim.Unit();
1392  r1Vec.Rotate( phi, pparUnit );
1393  TVector3 r2Vec( r1Vec ); r2Vec.Rotate( 0.5 * constants::kPi, pparUnit );
1394 
1395  double rprod = r1Vec.X() * r2Vec.X() + r1Vec.Y() * r2Vec.Y() + r1Vec.Z() * r2Vec.Z();
1396 
1397  assert( std::abs( rprod ) < controls::kASmallNum );
1398 
1399  // four IP with det. All have distance detRadius from centre.
1400  TVector3 p1( detO.X() + detRadius*r1Vec.X(),
1401  detO.Y() + detRadius*r1Vec.Y(),
1402  detO.Z() + detRadius*r1Vec.Z() );
1403  TVector3 p2( detO.X() - detRadius*r1Vec.X(),
1404  detO.Y() - detRadius*r1Vec.Y(),
1405  detO.Z() - detRadius*r1Vec.Z() );
1406  TVector3 p3( detO.X() + detRadius*r2Vec.X(),
1407  detO.Y() + detRadius*r2Vec.Y(),
1408  detO.Z() + detRadius*r2Vec.Z() );
1409  TVector3 p4( detO.X() - detRadius*r2Vec.X(),
1410  detO.Y() - detRadius*r2Vec.Y(),
1411  detO.Z() - detRadius*r2Vec.Z() );
1412 
1413  // Return largest(smallest) angle using inner product magic
1414  double thLarge = -999.9; double thSmall = 999.9;
1415  double th1 = TMath::ACos( ( p1.X()*pparUnit.X() + p1.Y()*pparUnit.Y() + p1.Z()*pparUnit.Z() ) / p1.Mag() ); if( th1 > thLarge ){ thLarge = th1; } else if( th1 < thSmall ){ thSmall = th1; }
1416  double th2 = TMath::ACos( ( p2.X()*pparUnit.X() + p2.Y()*pparUnit.Y() + p2.Z()*pparUnit.Z() ) / p2.Mag() ); if( th2 > thLarge ){ thLarge = th2; } else if( th2 < thSmall ){ thSmall = th2; }
1417  double th3 = TMath::ACos( ( p3.X()*pparUnit.X() + p3.Y()*pparUnit.Y() + p3.Z()*pparUnit.Z() ) / p3.Mag() ); if( th3 > thLarge ){ thLarge = th3; } else if( th3 < thSmall ){ thSmall = th3; }
1418  double th4 = TMath::ACos( ( p4.X()*pparUnit.X() + p4.Y()*pparUnit.Y() + p4.Z()*pparUnit.Z() ) / p4.Mag() ); if( th4 > thLarge ){ thLarge = th4; } else if( th4 < thSmall ){ thSmall = th4; }
1419 
1420  return ( seekingMax ) ? thLarge * 180.0 / constants::kPi : thSmall * 180.0 / constants::kPi;
1421  } else {
1422  // find direction from IP to det centre.
1423  TVector3 rVec = IPdev.Unit();
1424  // two IP with det. Closer(farther) has distance detRadius -(+) d( IP, detO )
1425  // actually, if IPdev endpoint lies inside detector have to go other way
1426  double dh = fLT + dist, dl = fLT - dist;
1427  // get those vectors and do inner product magic
1428  TVector3 ph( x_incp + dh * rVec.X(), y_incp + dh * rVec.Y(), z_incp + dh * rVec.Z() );
1429  TVector3 pl( x_incp - dl * rVec.X(), y_incp - dl * rVec.Y(), z_incp - dl * rVec.Z() );
1430  double thh = TMath::ACos( ( ph.X()*pparUnit.X() + ph.Y()*pparUnit.Y() + ph.Z()*pparUnit.Z() ) / ph.Mag() );
1431  double thl = TMath::ACos( ( pl.X()*pparUnit.X() + pl.Y()*pparUnit.Y() + pl.Z()*pparUnit.Z() ) / pl.Mag() );
1432 
1433  return ( seekingMax ) ? thh * 180.0 / constants::kPi : thl * 180.0 / constants::kPi;
1434  }
1435 
1436  LOG( "HNL", pERROR )
1437  << "Could not calculate the angle range for detector at ( " << detO.X() << ", "
1438  << detO.Y() << ", " << detO.Z() << " ) [m] from HNL production point with parent momentum = ( "
1439  << ppar.X() << ", " << ppar.Y() << ", " << ppar.Z() << " ) [GeV]. Returning zero.";
1440  return 0.0;
1441 }
pERROR
#define pERROR
Definition: Messenger.h:59
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::controls::kASmallNum
static const double kASmallNum
Definition: Controls.h:40
genie::RandomGen::RndGen
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
void FluxCreator::GetAngDeviation ( TLorentzVector  p4par,
TVector3  detO,
double &  zm,
double &  zp 
) const
private

Definition at line 1443 of file HNLFluxCreator.cxx.

References ApplyUserRotation(), CheckGeomPoint(), genie::units::cm, fBx1, fBx2, fBz, fCx, fCy, fCz, fDetOffset, fDetRotation, fDx, fDy, fDz, fLx, fLxR, fLy, fLyR, fLz, fLzR, LOG, genie::units::m, pDEBUG, and genie::utils::print::Vec3AsString().

1444 {
1445  // implementation of GetAngDeviation that uses ROOT geometry. More robust than analytical geom
1446  // (fewer assumptions about detector position)
1447 
1448  TVector3 ppar = p4par.Vect(); assert( ppar.Mag() > 0.0 );
1449  TVector3 pparUnit = ppar.Unit();
1450 
1451  const double sx1 = TMath::Sin(fBx1), cx1 = TMath::Cos(fBx1);
1452  const double sz = TMath::Sin(fBz), cz = TMath::Cos(fBz);
1453  const double sx2 = TMath::Sin(fBx2), cx2 = TMath::Cos(fBx2);
1454 
1455  const double xun[3] = { cz, -cx1*sz, sx1*sz };
1456  const double yun[3] = { sz*cx2, cx1*cz*cx2 - sx1*sx2, -sx1*cz*cx2 - cx1*sx2 };
1457  const double zun[3] = { sz*sx2, cx1*cz*sx2 + sx1*cx2, -sx1*cz*sx2 + cx1*cx2 };
1458 
1459  LOG("HNL", pDEBUG)
1460  << "\nxun = ( " << xun[0] << ", " << xun[1] << ", " << xun[2] << " )"
1461  << "\nyun = ( " << yun[0] << ", " << yun[1] << ", " << yun[2] << " )"
1462  << "\nzun = ( " << zun[0] << ", " << zun[1] << ", " << zun[2] << " )";
1463 
1464  /*
1465  TVector3 detO_cm( (detO.X() + fDetOffset.at(0)) * units::m / units::cm,
1466  (detO.Y() + fDetOffset.at(1)) * units::m / units::cm,
1467  (detO.Z() + fDetOffset.at(2)) * units::m / units::cm );
1468  */
1469  TVector3 detO_cm( (detO.X()) * units::m / units::cm,
1470  (detO.Y()) * units::m / units::cm,
1471  (detO.Z()) * units::m / units::cm );
1472  double inProd = zun[0] * detO_cm.X() + zun[1] * detO_cm.Y() + zun[2] * detO_cm.Z(); // cm
1473  assert( pparUnit.X() * zun[0] + pparUnit.Y() * zun[1] + pparUnit.Z() * zun[2] != 0.0 );
1474  inProd /= ( pparUnit.X() * zun[0] + pparUnit.Y() * zun[1] + pparUnit.Z() * zun[2] );
1475 
1476  // using vector formulation, find point of closest approach between parent momentum from
1477  // parent decay point, and detector centre.
1478 
1479  TVector3 dumori(0.0, 0.0, 0.0); // tgt-hall frame origin is 0
1480  TVector3 detori( (fDetOffset.at(0)) * units::m / units::cm,
1481  (fDetOffset.at(1)) * units::m / units::cm,
1482  (fDetOffset.at(2)) * units::m / units::cm ); // for rotations of the detector
1483 
1484  // Do all this in NEAR coords and m to avoid ambiguity.
1485 
1486  TVector3 fCvec_near( fCx, fCy, fCz ); // NEAR m
1487  TVector3 fDvec( fDx, fDy, fDz ); // in BEAM coords, m
1488  TVector3 fDvec_beam = this->ApplyUserRotation( fDvec, true ); // in NEAR coords, m
1489  TVector3 detO_near( fCvec_near.X() - fDvec_beam.X(),
1490  fCvec_near.Y() - fDvec_beam.Y(),
1491  fCvec_near.Z() - fDvec_beam.Z() );
1492  TVector3 detO_user = this->ApplyUserRotation( detO_near, dumori, fDetRotation, false ); // tgt-hall --> det
1493 
1494  const double aConst[3] = { fDvec_beam.X(), fDvec_beam.Y(), fDvec_beam.Z() };
1495  const double aConstUser[3] = { -detO_user.X(), -detO_user.Y(), -detO_user.Z() };
1496  /*
1497  const double dConst[3] = { fCx + fDetOffset.at(0), fCy + fDetOffset.at(1), fCz + fDetOffset.at(2) };
1498  */
1499  const double dConst[3] = { fCx, fCy, fCz };
1500  const double nConst[3] = { pparUnit.X(), pparUnit.Y(), pparUnit.Z() };
1501 
1502  // formula for POCA is \vec{a} + \vec{n} * ( (\vec{d} - \vec{a}) \cdot \vec{n} )
1503  const double nMult = nConst[0] * (dConst[0] - aConst[0]) +
1504  nConst[1] * (dConst[1] - aConst[1]) +
1505  nConst[2] * (dConst[2] - aConst[2]);
1506 
1507  // don't use the actual POCA, force calculation from that point V0 such that z(V0) = z(C)
1508  // and V0 lies on line joining decay point and POCA
1509 
1510  const double POCA_m[3] = { aConst[0] + nMult * nConst[0],
1511  aConst[1] + nMult * nConst[1],
1512  aConst[2] + nMult * nConst[2] }; // NEAR
1513  /*
1514  const double zConstMult = ((fCz + fDetOffset.at(2)) - aConst[2]) / (POCA_m[2] - aConst[2]);
1515  */
1516  const double zConstMult = ((fCz) - aConst[2]) / (POCA_m[2] - aConst[2]);
1517  const double startPoint_m[3] = { aConst[0] + nMult * zConstMult * nConst[0],
1518  aConst[1] + nMult * zConstMult * nConst[1],
1519  aConst[2] + nMult * zConstMult * nConst[2] }; // NEAR
1520 
1521  LOG( "HNL", pDEBUG )
1522  << "\ndetO_cm = " << utils::print::Vec3AsString( &detO_cm )
1523  << "\npparUnit = " << utils::print::Vec3AsString( &pparUnit );
1524 
1525  const double startPoint[3] = { startPoint_m[0] * units::m / units::cm,
1526  startPoint_m[1] * units::m / units::cm,
1527  startPoint_m[2] * units::m / units::cm }; // NEAR, cm
1528 
1529  /*
1530  const double sweepVect[3] = { (fCx + fDetOffset.at(0)) * units::m / units::cm - startPoint[0],
1531  (fCy + fDetOffset.at(1)) * units::m / units::cm - startPoint[1],
1532  (fCz + fDetOffset.at(2)) * units::m / units::cm - startPoint[2] }; // NEAR, cm
1533  */
1534  const double sweepVect[3] = { (fCx) * units::m / units::cm - startPoint[0],
1535  (fCy) * units::m / units::cm - startPoint[1],
1536  (fCz) * units::m / units::cm - startPoint[2] }; // NEAR, cm
1537  const double swvMag = std::sqrt( sweepVect[0]*sweepVect[0] + sweepVect[1]*sweepVect[1] + sweepVect[2]*sweepVect[2] ); assert( swvMag > 0.0 );
1538 
1539  // Note the geometry manager works in the *detector frame*. Transform to that.
1540  /*
1541  TVector3 detStartPoint( startPoint[0] - (fCx + fDetOffset.at(0)) * units::m / units::cm,
1542  startPoint[1] - (fCy + fDetOffset.at(1)) * units::m / units::cm,
1543  startPoint[2] - (fCz + fDetOffset.at(2)) * units::m / units::cm );
1544  */
1545  TVector3 detStartPoint( startPoint[0] - (fCx) * units::m / units::cm,
1546  startPoint[1] - (fCy) * units::m / units::cm,
1547  startPoint[2] - (fCz) * units::m / units::cm );
1548  TVector3 detSweepVect( sweepVect[0], sweepVect[1], sweepVect[2] );
1549 
1550  TVector3 detPpar = ppar;
1551 
1552  LOG( "HNL", pDEBUG )
1553  << "\nStartPoint = " << utils::print::Vec3AsString( &detStartPoint )
1554  << "\nSweepVect = " << utils::print::Vec3AsString( &detSweepVect )
1555  << "\nparent p3 = " << utils::print::Vec3AsString( &detPpar );
1556 
1557  detStartPoint = this->ApplyUserRotation( detStartPoint, detori, fDetRotation, false ); // passive transformation
1558  detSweepVect = this->ApplyUserRotation( detSweepVect, dumori, fDetRotation, true );
1559  detPpar = this->ApplyUserRotation( detPpar, dumori, fDetRotation, true );
1560 
1561  // first check that detStartPoint is not already in the detector! If it is, we should flag this now.
1562  std::string detPathString = this->CheckGeomPoint( detStartPoint.X(), detStartPoint.Y(), detStartPoint.Z() ); int iDNode = 1; // 1 past beginning
1563  bool startsInsideDet = ( detPathString.find("/", iDNode) != string::npos );
1564 
1565  TLorentzVector detPpar_4v( detPpar.X(), detPpar.Y(), detPpar.Z(), p4par.E() );
1566 
1567  // now sweep along sweepVect until we hit either side of the detector.
1568  // This will give us two points in space
1569 
1570  double minusPoint[3] = { 0.0, 0.0, 0.0 }; double plusPoint[3] = { 0.0, 0.0, 0.0 }; // USER, cm
1571  if( startsInsideDet ){
1572  minusPoint[0] = (gGeoManager->GetCurrentPoint())[0];
1573  minusPoint[1] = (gGeoManager->GetCurrentPoint())[1];
1574  minusPoint[2] = (gGeoManager->GetCurrentPoint())[2];
1575  }
1576 
1577  gGeoManager->SetCurrentPoint( detStartPoint.X(), detStartPoint.Y(), detStartPoint.Z() );
1578  gGeoManager->SetCurrentDirection( detSweepVect.X() / swvMag, detSweepVect.Y() / swvMag, detSweepVect.Z() / swvMag );
1579  const double sStepSize = 0.05 * std::min( std::min( fLx, fLy ), fLz );
1580 
1581  // start stepping. Let's do this manually cause FindNextBoundaryAndStep() can be finicky.
1582  // if start inside detector, then only find exit point, otherwise find entry point.
1583  if( startsInsideDet ){ // only find exit
1584 
1585  // to avoid large time inefficiencies just go to the edge of the box and step back by 10% of the size
1586  // this is *roughly correct* if not hugely correct
1587 
1588  double currx = (gGeoManager->GetCurrentPoint())[0];
1589  double curry = (gGeoManager->GetCurrentPoint())[1];
1590  double currz = (gGeoManager->GetCurrentPoint())[2];
1591  double currDist = std::sqrt( currx*currx + curry*curry + currz*currz ); // cm
1592 
1593  double curdx = (gGeoManager->GetCurrentDirection())[0];
1594  double curdy = (gGeoManager->GetCurrentDirection())[1];
1595  double curdz = (gGeoManager->GetCurrentDirection())[2];
1596 
1597  double stepMod = 0.99;
1598  double boxSize = std::sqrt( fLxR*fLxR + fLyR*fLyR + fLzR*fLzR )/2.0; // m
1599  double halfBoxSize = boxSize/2.0;
1600  double desiredDist = stepMod * halfBoxSize * units::m / units::cm; // cm
1601 
1602  // now we're at distance d on the one side of our detector centre
1603  // we want to get to distance D on the other side of our detector centre
1604  // meaning we should step by d + D
1605  double largeStep = currDist + desiredDist;
1606 
1607  gGeoManager->SetCurrentPoint( currx + largeStep * curdx,
1608  curry + largeStep * curdy,
1609  currz + largeStep * curdz );
1610 
1611  /* // this is the very correct, very slow way of doing it
1612  while( detPathString.find("/", iDNode) != string::npos ){
1613  gGeoManager->SetCurrentPoint( (gGeoManager->GetCurrentPoint())[0] + (gGeoManager->GetCurrentDirection())[0] * sStepSize,
1614  (gGeoManager->GetCurrentPoint())[1] + (gGeoManager->GetCurrentDirection())[1] * sStepSize,
1615  (gGeoManager->GetCurrentPoint())[2] + (gGeoManager->GetCurrentDirection())[2] * sStepSize );
1616  detPathString = this->CheckGeomPoint( (gGeoManager->GetCurrentPoint())[0], (gGeoManager->GetCurrentPoint())[1], (gGeoManager->GetCurrentPoint())[2] );
1617  }
1618  */
1619 
1620  plusPoint[0] = (gGeoManager->GetCurrentPoint())[0];
1621  plusPoint[1] = (gGeoManager->GetCurrentPoint())[1];
1622  plusPoint[2] = (gGeoManager->GetCurrentPoint())[2];
1623  } else { // find entry and exit
1624  // first, entry
1625 
1626  // to avoid large time inefficiencies just go to the edge of the box and step back by 10% of the size
1627  // this is *roughly correct* if not hugely correct
1628 
1629  // find that point that lies on the line along direction from current point that is 90% of box size from box centre
1630  double currx = (gGeoManager->GetCurrentPoint())[0];
1631  double curry = (gGeoManager->GetCurrentPoint())[1];
1632  double currz = (gGeoManager->GetCurrentPoint())[2];
1633  double currDist = std::sqrt( currx*currx + curry*curry + currz*currz );
1634 
1635  double stepMod = 0.99;
1636  double boxSize = std::sqrt( fLxR*fLxR + fLyR*fLyR + fLzR*fLzR )/2.0; // m
1637  double halfBoxSize = boxSize / 2.0;
1638  double desiredDist = stepMod * halfBoxSize * units::m / units::cm;
1639 
1640  // we're at some distance D and want to get to distance d, which means we step forward by
1641  // R = D-d = D(1-d/D)
1642  double largeStep = currDist - desiredDist;
1643 
1644  double curdx = (gGeoManager->GetCurrentDirection())[0];
1645  double curdy = (gGeoManager->GetCurrentDirection())[1];
1646  double curdz = (gGeoManager->GetCurrentDirection())[2];
1647 
1648  gGeoManager->SetCurrentPoint( currx + largeStep * curdx,
1649  curry + largeStep * curdy,
1650  currz + largeStep * curdz );
1651 
1652  /* // slow but correct
1653  while( detPathString.find("/", iDNode) == string::npos ){
1654  gGeoManager->SetCurrentPoint( (gGeoManager->GetCurrentPoint())[0] + (gGeoManager->GetCurrentDirection())[0] * sStepSize,
1655  (gGeoManager->GetCurrentPoint())[1] + (gGeoManager->GetCurrentDirection())[1] * sStepSize,
1656  (gGeoManager->GetCurrentPoint())[2] + (gGeoManager->GetCurrentDirection())[2] * sStepSize );
1657  detPathString = this->CheckGeomPoint( (gGeoManager->GetCurrentPoint())[0], (gGeoManager->GetCurrentPoint())[1], (gGeoManager->GetCurrentPoint())[2] );
1658  }
1659  */
1660 
1661  minusPoint[0] = (gGeoManager->GetCurrentPoint())[0];
1662  minusPoint[1] = (gGeoManager->GetCurrentPoint())[1];
1663  minusPoint[2] = (gGeoManager->GetCurrentPoint())[2];
1664 
1665  // then, exit
1666 
1667  // quick and dirty way: reflect about the origin
1668 
1669  currx = (gGeoManager->GetCurrentPoint())[0];
1670  curry = (gGeoManager->GetCurrentPoint())[1];
1671  currz = (gGeoManager->GetCurrentPoint())[2];
1672  /*
1673  double cdx = fDetOffset.at(0) * units::m / units::cm - currx; // cm
1674  double cdy = fDetOffset.at(1) * units::m / units::cm - curry; // cm
1675  double cdz = fDetOffset.at(2) * units::m / units::cm - currz; // cm
1676 
1677  gGeoManager->SetCurrentPoint( fDetOffset.at(0) * units::m / units::cm + cdx,
1678  fDetOffset.at(1) * units::m / units::cm + cdy,
1679  fDetOffset.at(2) * units::m / units::cm + cdz );
1680  */
1681 
1682  gGeoManager->SetCurrentPoint( -currx, -curry, -currz );
1683 
1684  /* // slow but correct
1685  while( detPathString.find("/", iDNode) != string::npos ){
1686  gGeoManager->SetCurrentPoint( (gGeoManager->GetCurrentPoint())[0] + (gGeoManager->GetCurrentDirection())[0] * sStepSize,
1687  (gGeoManager->GetCurrentPoint())[1] + (gGeoManager->GetCurrentDirection())[1] * sStepSize,
1688  (gGeoManager->GetCurrentPoint())[2] + (gGeoManager->GetCurrentDirection())[2] * sStepSize );
1689  detPathString = this->CheckGeomPoint( (gGeoManager->GetCurrentPoint())[0], (gGeoManager->GetCurrentPoint())[1], (gGeoManager->GetCurrentPoint())[2] );
1690  }
1691  */
1692 
1693  plusPoint[0] = (gGeoManager->GetCurrentPoint())[0];
1694  plusPoint[1] = (gGeoManager->GetCurrentPoint())[1];
1695  plusPoint[2] = (gGeoManager->GetCurrentPoint())[2];
1696  }
1697 
1698  /*
1699  int bdIdx = 0; const int bdIdxMax = 1e+4;
1700  if(!startsInsideDet){
1701  while( gGeoManager->FindNextBoundaryAndStep() && bdIdx < bdIdxMax ){
1702  bdIdx++;
1703  if( bdIdx % 100 == 0 )
1704  LOG( "HNL", pDEBUG ) << "bdIdx = " << bdIdx;
1705  if( std::abs( (gGeoManager->GetCurrentPoint())[0] ) != fLxR/2.0 * units::m / units::cm &&
1706  std::abs( (gGeoManager->GetCurrentPoint())[1] ) != fLyR/2.0 * units::m / units::cm &&
1707  std::abs( (gGeoManager->GetCurrentPoint())[2] ) != fLzR/2.0 * units::m / units::cm ){
1708  plusPoint[0] = (gGeoManager->GetCurrentPoint())[0];
1709  plusPoint[1] = (gGeoManager->GetCurrentPoint())[1];
1710  plusPoint[2] = (gGeoManager->GetCurrentPoint())[2];
1711  }
1712  }
1713  }
1714  */
1715 
1716  TVector3 originPoint( -(fCx + fDetOffset.at(0)), -(fCy + fDetOffset.at(1)), -(fCz + fDetOffset.at(2)) );
1717 
1718  TVector3 vMUser( minusPoint[0] * units::cm / units::m,
1719  minusPoint[1] * units::cm / units::m,
1720  minusPoint[2] * units::cm / units::m ); // USER, m
1721  TVector3 vPUser( plusPoint[0] * units::cm / units::m,
1722  plusPoint[1] * units::cm / units::m,
1723  plusPoint[2] * units::cm / units::m ); // USER, m
1724 
1725  // rotate to NEAR frame
1726  TVector3 vMNear = this->ApplyUserRotation( vMUser, originPoint, fDetRotation, true ); // det --> tgt-hall
1727  /*
1728  vMNear.SetXYZ( vMNear.X() + (fCx + fDetOffset.at(0)),
1729  vMNear.Y() + (fCy + fDetOffset.at(1)),
1730  vMNear.Z() + (fCz + fDetOffset.at(2)) );
1731  */
1732  vMNear.SetXYZ( vMNear.X() + (fCx),
1733  vMNear.Y() + (fCy),
1734  vMNear.Z() + (fCz) );
1735  TVector3 vPNear = this->ApplyUserRotation( vPUser, originPoint, fDetRotation, true ); // det --> tgt-hall
1736  /*
1737  vPNear.SetXYZ( vPNear.X() + (fCx + fDetOffset.at(0)),
1738  vPNear.Y() + (fCy + fDetOffset.at(1)),
1739  vPNear.Z() + (fCz + fDetOffset.at(2)) );
1740  */
1741  vPNear.SetXYZ( vPNear.X() + (fCx),
1742  vPNear.Y() + (fCy),
1743  vPNear.Z() + (fCz) );
1744 
1745  // with 3 points and 1 vector we calculate the angles.
1746  // Points: D(decay), E(entry), X(exit) [all in local, cm]. Vector: detPpar [local, GeV/GeV]
1747  // angles are <DE, detPpar> and <DX, detPpar>
1748 
1749  // now obtain the angles themselves and return in deg.
1750  /*
1751  TVector3 decayPoint_user( (fDx - (fCx + fDetOffset.at(0))) * units::m / units::cm,
1752  (fDy - (fCy + fDetOffset.at(1))) * units::m / units::cm,
1753  (fDz - (fCz + fDetOffset.at(2))) * units::m / units::cm ); // USER, cm
1754  */
1755  TVector3 decayPoint_user( aConstUser[0] * units::m / units::cm,
1756  aConstUser[1] * units::m / units::cm,
1757  aConstUser[2] * units::m / units::cm );
1758  TVector3 decayPoint_near = this->ApplyUserRotation( decayPoint_user, detori, fDetRotation, false ); // NEAR, cm
1759  /*
1760  decayPoint_near.SetXYZ( decayPoint_near.X() + (fCx + fDetOffset.at(0)) * units::m / units::cm,
1761  decayPoint_near.Y() + (fCy + fDetOffset.at(1)) * units::m / units::cm,
1762  decayPoint_near.Z() + (fCz + fDetOffset.at(2)) * units::m / units::cm );
1763  */
1764  decayPoint_near.SetXYZ( decayPoint_near.X() + (fCx) * units::m / units::cm,
1765  decayPoint_near.Y() + (fCy) * units::m / units::cm,
1766  decayPoint_near.Z() + (fCz) * units::m / units::cm );
1767 
1768  TVector3 minusVec( minusPoint[0] - decayPoint_user.X(), minusPoint[1] - decayPoint_user.Y(), minusPoint[2] - decayPoint_user[2] ); // USER, cm
1769  TVector3 plusVec( plusPoint[0] - decayPoint_user.X(), plusPoint[1] - decayPoint_user.Y(), plusPoint[2] - decayPoint_user[2] ); // USER, cm
1770  TVector3 startVec( detStartPoint[0] - decayPoint_user.X(), detStartPoint[1] - decayPoint_user.Y(), detStartPoint[2] - decayPoint_user.Z() ); // USER, cm
1771 
1772  TVector3 minusVec_near = this->ApplyUserRotation( minusVec, detori, fDetRotation, false ); // NEAR, cm
1773  TVector3 plusVec_near = this->ApplyUserRotation( plusVec, detori, fDetRotation, false ); // NEAR, cm
1774  TVector3 startVec_near = this->ApplyUserRotation( startVec, detori, fDetRotation, false ); // NEAR, cm
1775 
1776  double minusNum = startVec.X() * minusVec.X() + startVec.Y() * minusVec.Y() + startVec.Z() * minusVec.Z(); // USER AND USER
1777  double minusDen = startVec.Mag() * minusVec.Mag(); assert( minusDen > 0.0 ); // USER AND USER
1778 
1779  zm = TMath::ACos( minusNum / minusDen ) * TMath::RadToDeg();
1780 
1781  double plusNum = startVec.X() * plusVec.X() + startVec.Y() * plusVec.Y() + startVec.Z() * plusVec.Z(); // USER AND USER
1782  double plusDen = startVec.Mag() * plusVec.Mag(); assert( plusDen > 0.0 ); // USER AND USER
1783 
1784  zp = TMath::ACos( plusNum / plusDen ) * TMath::RadToDeg();
1785 
1786  if( zm > zp ){
1787  double tmpzp = zp;
1788  zp = zm;
1789  zm = tmpzp;
1790  }
1791 
1792  /*
1793  LOG( "HNL", pDEBUG )
1794  << "\nIn DETECTOR coordinates:"
1795  << "\nentered at ( " << minusPoint[0] << ", " << minusPoint[1] << ", " << minusPoint[2] << " ) [cm]"
1796  << "\nexited at ( " << plusPoint[0] << ", " << plusPoint[1] << ", " << plusPoint[2] << " ) [cm]"
1797  << "\nstarted at ( " << decVec.X() << ", " << decVec.Y() << ", " << decVec.Z() << " ) [cm]"
1798  << "\nmomentum ( " << detPpar.X() << ", " << detPpar.Y() << ", " << detPpar.Z() << " )"
1799  << "\nmeaning zm = " << zm << ", zp = " << zp << " [deg]";
1800  */
1801 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
genie::hnl::FluxCreator::fDetOffset
std::vector< double > fDetOffset
Definition: HNLFluxCreator.h:214
genie::hnl::FluxCreator::fDetRotation
std::vector< double > fDetRotation
Definition: HNLFluxCreator.h:213
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::FluxCreator::CheckGeomPoint
std::string CheckGeomPoint(Double_t x, Double_t y, Double_t z) const
Definition: HNLFluxCreator.cxx:2136
genie::hnl::FluxCreator::ApplyUserRotation
TVector3 ApplyUserRotation(TVector3 vec, bool doBackwards=false) const
Definition: HNLFluxCreator.cxx:1959
genie::units::m
static constexpr double m
Definition: Units.h:71
genie::utils::print::Vec3AsString
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
std::vector< double > genie::hnl::FluxCreator::GetB2URotation ( ) const
inlinevirtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 125 of file HNLFluxCreator.h.

References fB2URotation.

Referenced by ReadInConfig().

125 { return fB2URotation; }
genie::hnl::FluxCreator::fB2URotation
std::vector< double > fB2URotation
Definition: HNLFluxCreator.h:212
std::vector< double > genie::hnl::FluxCreator::GetB2UTranslation ( ) const
inlinevirtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 124 of file HNLFluxCreator.h.

References fB2UTranslation.

Referenced by ReadInConfig().

124 { return fB2UTranslation; }
genie::hnl::FluxCreator::fB2UTranslation
std::vector< double > fB2UTranslation
Definition: HNLFluxCreator.h:212
std::vector< double > genie::hnl::FluxCreator::GetDetOffset ( ) const
inlinevirtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 126 of file HNLFluxCreator.h.

References fDetOffset.

126 { return fDetOffset; }
genie::hnl::FluxCreator::fDetOffset
std::vector< double > fDetOffset
Definition: HNLFluxCreator.h:214
std::vector< double > genie::hnl::FluxCreator::GetDetRotation ( ) const
inlinevirtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 127 of file HNLFluxCreator.h.

References fDetRotation.

127 { return fDetRotation; }
genie::hnl::FluxCreator::fDetRotation
std::vector< double > fDetRotation
Definition: HNLFluxCreator.h:213
int FluxCreator::GetNFluxEntries ( ) const
virtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 134 of file HNLFluxCreator.cxx.

References fCurrPath, fNEntries, and OpenFluxInput().

Referenced by main().

135 {
136  if( fNEntries <= 0 ){
137  this->OpenFluxInput( fCurrPath );
138  }
139  return fNEntries;
140 }
genie::hnl::FluxCreator::OpenFluxInput
void OpenFluxInput(std::string finpath) const
Definition: HNLFluxCreator.cxx:741
std::map< HNLProd_t, double > FluxCreator::GetProductionProbs ( int  parPDG) const
private

Definition at line 1027 of file HNLFluxCreator.cxx.

References BR_K03e, BR_K03mu, BR_K2e, BR_K2mu, BR_K3e, BR_K3mu, BR_pi2e, BR_pi2mu, dynamicScores_kaon, dynamicScores_muon, dynamicScores_neuk, dynamicScores_pion, fU4l2s, genie::AlgFactory::GetAlgorithm(), genie::AlgFactory::Instance(), genie::hnl::kHNLProdKaon2Electron, genie::hnl::kHNLProdKaon2Muon, genie::hnl::kHNLProdKaon3Electron, genie::hnl::kHNLProdKaon3Muon, genie::hnl::kHNLProdMuon3Nue, genie::hnl::kHNLProdMuon3Numu, genie::hnl::kHNLProdMuon3Nutau, genie::hnl::kHNLProdNeuk3Electron, genie::hnl::kHNLProdNeuk3Muon, genie::hnl::kHNLProdNull, genie::hnl::kHNLProdPion2Electron, genie::hnl::kHNLProdPion2Muon, genie::hnl::BRCalculator::KinematicScaling(), genie::kPdgK0L, genie::kPdgKP, genie::kPdgMuon, genie::kPdgPiP, LOG, pDEBUG, pERROR, pWARN, and ReadBRs().

Referenced by MakeTupleFluxEntry().

1028 {
1029  // check if we've calculated scores before
1030  switch( std::abs( parPDG ) ){
1031  case kPdgPiP : if( dynamicScores_pion.size() > 0 ) return dynamicScores_pion; break;
1032  case kPdgKP : if( dynamicScores_kaon.size() > 0 ) return dynamicScores_kaon; break;
1033  case kPdgMuon: if( dynamicScores_muon.size() > 0 ) return dynamicScores_muon; break;
1034  case kPdgK0L : if( dynamicScores_neuk.size() > 0 ) return dynamicScores_neuk; break;
1035  default: LOG( "HNL", pWARN ) << "Unknown parent. Proceeding, but results may be unphysical"; break;
1036  }
1037 
1038  std::map< HNLProd_t, double > dynScores;
1039 
1040  // first get branching ratios to SM
1041  ReadBRs();
1042  // then get HNL parameter space
1043 
1044  double Ue42 = fU4l2s.at(0);
1045  double Um42 = fU4l2s.at(1);
1046  double Ut42 = fU4l2s.at(2);
1047 
1048  // also, construct an BRCalculator * object to handle the scalings.
1049  const Algorithm * algBRCalc = AlgFactory::Instance()->GetAlgorithm("genie::hnl::BRCalculator", "Default");
1050  const BRCalculator * BRCalc = dynamic_cast< const BRCalculator * >( algBRCalc );
1051 
1052  // first get pure kinematic part of the BRs
1053  double KScale[4] = { -1.0, -1.0, -1.0, -1.0 }, mixScale[4] = { -1.0, -1.0, -1.0, -1.0 };
1054  double totalMix = 0.0;
1055  switch( std::abs( parPDG ) ){
1056  case genie::kPdgMuon:
1057  KScale[0] = BRCalc->KinematicScaling( kHNLProdMuon3Numu );
1058  KScale[1] = BRCalc->KinematicScaling( kHNLProdMuon3Nue ); // same, convenience for later
1059  KScale[2] = BRCalc->KinematicScaling( kHNLProdMuon3Nutau ); // same, convenience for later
1060  mixScale[0] = 1.0 * Um42 * KScale[0]; totalMix += mixScale[0];
1061  mixScale[1] = 1.0 * Ue42 * KScale[1]; totalMix += mixScale[1];
1062  mixScale[2] = 1.0 * Ut42 * KScale[2]; totalMix += mixScale[2];
1063 
1064  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdMuon3Numu, mixScale[0] / totalMix } ) );
1065  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdMuon3Nue, mixScale[1] / totalMix } ) );
1066  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdMuon3Nutau, mixScale[2] / totalMix } ) );
1067 
1068  // it can happen that HNL is not coupled to the only kinematically available channel.
1069  // Return bogus map if that happens
1070  if( totalMix <= 0.0 ){
1071  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNull, -999.9 } ) );
1072  dynamicScores_muon = dynScores;
1073  return dynScores;
1074  }
1075 
1076  dynamicScores_muon = dynScores;
1077  break;
1078  case genie::kPdgKP:
1079  KScale[0] = BRCalc->KinematicScaling( kHNLProdKaon2Muon );
1080  KScale[1] = BRCalc->KinematicScaling( kHNLProdKaon2Electron );
1081  KScale[2] = BRCalc->KinematicScaling( kHNLProdKaon3Muon );
1082  KScale[3] = BRCalc->KinematicScaling( kHNLProdKaon3Electron );
1083  mixScale[0] = BR_K2mu * Um42 * KScale[0]; totalMix += mixScale[0];
1084  mixScale[1] = BR_K2e * Ue42 * KScale[1]; totalMix += mixScale[1];
1085  mixScale[2] = BR_K3mu * Um42 * KScale[2]; totalMix += mixScale[2];
1086  mixScale[3] = BR_K3e * Ue42 * KScale[3]; totalMix += mixScale[3];
1087 
1088  if( totalMix <= 0.0 ){
1089  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNull, -999.9 } ) );
1090  dynamicScores_pion = dynScores;
1091  return dynScores;
1092  }
1093 
1094  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdKaon2Muon, mixScale[0] / totalMix } ) );
1095  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdKaon2Electron, mixScale[1] / totalMix } ) );
1096  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdKaon3Muon, mixScale[2] / totalMix } ) );
1097  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdKaon3Electron, mixScale[3] / totalMix } ) );
1098 
1099  dynamicScores_kaon = dynScores;
1100  break;
1101  case genie::kPdgPiP:
1102 
1103  KScale[0] = BRCalc->KinematicScaling( kHNLProdPion2Muon );
1104  KScale[1] = BRCalc->KinematicScaling( kHNLProdPion2Electron );
1105  mixScale[0] = BR_pi2mu * Um42 * KScale[0]; totalMix += mixScale[0];
1106  mixScale[1] = BR_pi2e * Ue42 * KScale[1]; totalMix += mixScale[1];
1107 
1108  if( totalMix <= 0.0 ){
1109  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNull, -999.9 } ) );
1110  dynamicScores_pion = dynScores;
1111  return dynScores;
1112  }
1113 
1114  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdPion2Muon, mixScale[0] / totalMix } ) );
1115  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdPion2Electron, mixScale[1] / totalMix } ) );
1116 
1117  dynamicScores_pion = dynScores;
1118  break;
1119  case genie::kPdgK0L:
1120 
1121  KScale[0] = BRCalc->KinematicScaling( kHNLProdNeuk3Muon );
1122  KScale[1] = BRCalc->KinematicScaling( kHNLProdNeuk3Electron );
1123  mixScale[0] = BR_K03mu * Um42 * KScale[0]; totalMix += mixScale[0];
1124  mixScale[1] = BR_K03e * Ue42 * KScale[1]; totalMix += mixScale[1];
1125 
1126  if( totalMix <= 0.0 ){
1127  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNull, -999.9 } ) );
1128  dynamicScores_neuk = dynScores;
1129  return dynScores;
1130  }
1131 
1132  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNeuk3Muon, mixScale[0] / totalMix } ) );
1133  dynScores.insert( std::pair< HNLProd_t, double >( { kHNLProdNeuk3Electron, mixScale[1] / totalMix } ) );
1134 
1135  dynamicScores_neuk = dynScores;
1136  break;
1137  default:
1138  LOG( "HNL", pERROR )
1139  << "Unknown parent particle. Cannot make scales, exiting."; exit(1);
1140  }
1141 
1142  LOG( "HNL", pDEBUG )
1143  << "Score map now has " << dynScores.size() << " elements. Returning.";
1144  return dynScores;
1145 
1146 }
genie::hnl::FluxCreator::dynamicScores_neuk
std::map< genie::hnl::HNLProd_t, double > dynamicScores_neuk
Definition: HNLFluxCreator.h:192
pERROR
#define pERROR
Definition: Messenger.h:59
genie::hnl::BRCalculator
Manages HNL BR (prod and decay)
Definition: HNLBRCalculator.h:53
genie::Algorithm
Algorithm abstract base class.
Definition: Algorithm.h:54
genie::hnl::BRCalculator::KinematicScaling
double KinematicScaling(genie::hnl::HNLProd_t hnlprod) const
Definition: HNLBRCalculator.cxx:126
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::FluxCreator::dynamicScores_muon
std::map< genie::hnl::HNLProd_t, double > dynamicScores_muon
Definition: HNLFluxCreator.h:191
genie::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::hnl::FluxCreator::dynamicScores_pion
std::map< genie::hnl::HNLProd_t, double > dynamicScores_pion
Definition: HNLFluxCreator.h:189
genie::kPdgK0L
const int kPdgK0L
Definition: PDGCodes.h:176
genie::hnl::FluxCreator::dynamicScores_kaon
std::map< genie::hnl::HNLProd_t, double > dynamicScores_kaon
Definition: HNLFluxCreator.h:190
pWARN
#define pWARN
Definition: Messenger.h:60
genie::AlgFactory::Instance
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
genie::hnl::FluxCreator::fU4l2s
std::vector< double > fU4l2s
Definition: HNLFluxCreator.h:204
genie::kPdgMuon
const int kPdgMuon
Definition: PDGCodes.h:37
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
TLorentzVector FluxCreator::HNLEnergy ( genie::hnl::HNLProd_t  hnldm,
TLorentzVector  p4par 
) const
private

accept_decay

Definition at line 1148 of file HNLFluxCreator.cxx.

References doPol, fFixedPolarisation, genie::PDGLibrary::Find(), fixPol, fLepPdg, fLPE, fLPx, fLPy, fLPz, genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::controls::kMaxUnweightDecayIterations, genie::kPdgElectron, genie::kPdgHNL, genie::kPdgMuon, genie::kPdgTau, LOG, genie::units::m, genie::utils::print::P4AsString(), pDEBUG, pERROR, pINFO, pNOTICE, genie::utils::hnl::ProdAsString(), genie::utils::hnl::ProductionProductList(), genie::PDGCodeList::push_back(), pWARN, genie::RandomGen::RndHadro(), genie::exceptions::EVGThreadException::SetReason(), and genie::exceptions::EVGThreadException::SwitchOnFastForward().

Referenced by MakeTupleFluxEntry().

1149 {
1150  // first boost to parent rest frame
1151  TLorentzVector p4par_rest = p4par;
1152  TVector3 boost_beta = p4par.BoostVector();
1153  p4par_rest.Boost( -boost_beta );
1154 
1155  LOG( "HNL", pDEBUG )
1156  << "Attempting to decay rest-system p4 = " << utils::print::P4AsString(&p4par_rest)
1157  << " as " << utils::hnl::ProdAsString( hnldm );
1158 
1159  // get PDGCodeList and truncate 1st member
1160  PDGCodeList fullList = utils::hnl::ProductionProductList( hnldm );
1161  bool allow_duplicate = true;
1162  PDGCodeList decayList( allow_duplicate );
1163  double * mass = new double[decayList.size()];
1164  double sum = 0.0;
1165 
1166  for( std::vector<int>::const_iterator pdg_iter = fullList.begin(); pdg_iter != fullList.end(); ++pdg_iter )
1167  {
1168  if( pdg_iter != fullList.begin() ){
1169  int pdgc = *pdg_iter;
1170  decayList.push_back( pdgc );
1171  }
1172  }
1173 
1174  int iv = 0;
1175  for( std::vector<int>::const_iterator pdg_iter = decayList.begin(); pdg_iter != decayList.end(); ++pdg_iter )
1176  {
1177  int pdgc = *pdg_iter;
1178  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
1179  mass[iv++] = m; sum += m;
1180  }
1181 
1182  // Set the decay
1183  TGenPhaseSpace fPhaseSpaceGenerator;
1184  bool permitted = fPhaseSpaceGenerator.SetDecay( p4par_rest, decayList.size(), mass );
1185  if(!permitted) {
1186  LOG("HNL", pERROR)
1187  << " *** Phase space decay is not permitted \n"
1188  << " Total particle mass = " << sum << "\n"
1189  << " Decaying system p4 = " << utils::print::P4AsString(&p4par_rest);
1190  // clean-up
1191  delete [] mass;
1192  // throw exception
1193  genie::exceptions::EVGThreadException exception;
1194  exception.SetReason("Decay not permitted kinematically");
1195  exception.SwitchOnFastForward();
1196  throw exception;
1197  }
1198 
1199  // Get the maximum weight
1200  double wmax = -1;
1201  for(int idec=0; idec<200; idec++) {
1202  double w = fPhaseSpaceGenerator.Generate();
1203  wmax = TMath::Max(wmax,w);
1204  }
1205  assert(wmax>0);
1206  wmax *= 2;
1207 
1208  LOG("HNL", pNOTICE)
1209  << "Max phase space gen. weight @ current HNL system: " << wmax;
1210 
1211  // Generate an unweighted decay
1212  RandomGen * rnd = RandomGen::Instance();
1213 
1214  bool accept_decay=false;
1215  unsigned int itry=0;
1216  while(!accept_decay)
1217  {
1218  itry++;
1219 
1220  if(itry > controls::kMaxUnweightDecayIterations) {
1221  // report, clean-up and return
1222  LOG("HNL", pWARN)
1223  << "Couldn't generate an unweighted phase space decay after "
1224  << itry << " attempts";
1225  // clean up
1226  delete [] mass;
1227  // throw exception
1228  genie::exceptions::EVGThreadException exception;
1229  exception.SetReason("Couldn't select decay after N attempts");
1230  exception.SwitchOnFastForward();
1231  throw exception;
1232  }
1233  double w = fPhaseSpaceGenerator.Generate();
1234  if(w > wmax) {
1235  LOG("HNL", pWARN)
1236  << "Decay weight = " << w << " > max decay weight = " << wmax;
1237  }
1238  double gw = wmax * rnd->RndHadro().Rndm();
1239  accept_decay = (gw<=w);
1240 
1241  LOG("HNL", pINFO)
1242  << "Decay weight = " << w << " / R = " << gw
1243  << " - accepted: " << accept_decay;
1244 
1245  } //!accept_decay
1246 
1247  // Grab 0th entry energy and return that
1248  int idp = 0; TLorentzVector p4HNL, p4HNL_rest;
1249  // search for the charged lepton in this stack, get the 4-vector in parent rest frame
1250  TLorentzVector p4Lep, p4Lep_HNLRest;
1251 
1252  p4HNL.SetPxPyPzE( 0.0, 0.0, 0.0, 0.0 );
1253  p4HNL_rest.SetPxPyPzE( 0.0, 0.0, 0.0, 0.0 );
1254  p4Lep.SetPxPyPzE( 0.0, 0.0, 0.0, 0.0 );
1255  p4Lep_HNLRest.SetPxPyPzE( 0.0, 0.0, 0.0, 0.0 );
1256 
1257  for(std::vector<int>::const_iterator pdg_iter = decayList.begin(); pdg_iter != decayList.end(); ++pdg_iter) {
1258  int pdgc = *pdg_iter;
1259  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
1260 
1261  if( std::abs( pdgc ) == kPdgHNL ) p4HNL = *p4fin;
1262  if( std::abs( pdgc ) == kPdgElectron ||
1263  std::abs( pdgc ) == kPdgMuon ||
1264  std::abs( pdgc ) == kPdgTau ){
1265  p4Lep = *p4fin;
1266  fLepPdg = pdgc;
1267  }
1268  idp++;
1269  }
1270 
1271  if( doPol ){
1272  // boost this to HNL rest frame.
1273  TVector3 boostHNL = p4HNL.BoostVector();
1274  p4Lep_HNLRest = p4Lep; fLPE = p4Lep_HNLRest.E(); // still in parent rest frame here
1275  p4Lep_HNLRest.Boost( boostHNL );
1276  // save this
1277  fLPx = ( fixPol ) ? fFixedPolarisation.at(0) : p4Lep.Px() / p4Lep.P(); // note that this is for a true random decay.
1278  fLPy = ( fixPol ) ? fFixedPolarisation.at(1) : p4Lep.Py() / p4Lep.P(); // We still need to take the geometrical
1279  fLPz = ( fixPol ) ? fFixedPolarisation.at(2) : p4Lep.Pz() / p4Lep.P(); // constraint into account.
1280  } else {
1281  fLPx = 0.0;
1282  fLPy = 0.0;
1283  fLPz = 0.0;
1284  }
1285 
1286  delete [] mass;
1287  return p4HNL; // rest frame momentum!
1288 }
genie::hnl::FluxCreator::fFixedPolarisation
std::vector< double > fFixedPolarisation
Definition: HNLFluxCreator.h:226
pERROR
#define pERROR
Definition: Messenger.h:59
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::utils::hnl::ProductionProductList
PDGCodeList ProductionProductList(genie::hnl::HNLProd_t hnldm)
genie::utils::print::P4AsString
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
genie::controls::kMaxUnweightDecayIterations
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
genie::kPdgElectron
const int kPdgElectron
Definition: PDGCodes.h:35
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::exceptions::EVGThreadException
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
Definition: EVGThreadException.h:35
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::kPdgTau
const int kPdgTau
Definition: PDGCodes.h:39
genie::kPdgHNL
const int kPdgHNL
Definition: PDGCodes.h:224
genie::utils::hnl::ProdAsString
string ProdAsString(genie::hnl::HNLProd_t hnlprod)
pINFO
#define pINFO
Definition: Messenger.h:62
pWARN
#define pWARN
Definition: Messenger.h:60
genie::RandomGen::RndHadro
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::kPdgMuon
const int kPdgMuon
Definition: PDGCodes.h:37
genie::units::m
static constexpr double m
Definition: Units.h:71
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::ImportBoundingBox ( TGeoBBox *  box) const
private

Definition at line 2121 of file HNLFluxCreator.cxx.

References genie::units::cm, fDoingOldFluxCalc, fLx, fLxR, fLy, fLyR, fLz, fLzR, fRadius, LOG, genie::units::m, and pDEBUG.

Referenced by ProcessEventRecord().

2122 {
2123  LOG( "HNL", pDEBUG ) << "Importing bounding box...";
2124  fLxR = 2.0 * box->GetDX() * units::cm / units::m;
2125  fLyR = 2.0 * box->GetDY() * units::cm / units::m;
2126  fLzR = 2.0 * box->GetDZ() * units::cm / units::m;
2127 
2128  double testRadius = fRadius; // m
2129  if( !fDoingOldFluxCalc ){
2130  fLx = std::min( testRadius, fLxR );
2131  fLy = std::min( testRadius, fLyR );
2132  fLz = std::min( testRadius, fLzR );
2133  }
2134 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::units::m
static constexpr double m
Definition: Units.h:71
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::InitialiseMeta ( ) const
private

Definition at line 938 of file HNLFluxCreator.cxx.

References beam0x, beam0y, beam0z, beamdxdz, beamdydz, beamhwidth, beamsim, beamvwidth, cmeta, dkvolcfg, horncfg, job, location_name, location_x, location_y, location_z, mArSize, maxArray, maxC, physcuts, physics, pots, and tgtcfg.

Referenced by MakeTupleFluxEntry().

939 {
940  job = 0;
941  pots = 0.0;
942 
943  beam0x = -9999.9;
944  beam0y = -9999.9;
945  beam0z = -9999.9;
946  beamhwidth = -9999.9;
947  beamvwidth = -9999.9;
948  beamdxdz = -9999.9;
949  beamdydz = -9999.9;
950  mArSize = 0;
951 
952  for( int i = 0; i < maxC; i++ ){
953  beamsim[i] = 0;
954  physics[i] = 0;
955  physcuts[i] = 0;
956  tgtcfg[i] = 0;
957  horncfg[i] = 0;
958  dkvolcfg[i] = 0;
959  }
960 
961  for( int i = 0; i < maxArray; i++ ){
962  location_x[i] = -9999.9;
963  location_y[i] = -9999.9;
964  location_z[i] = -9999.9;
965 
966  for( int j = 0; j < maxC; j++ ){ location_name[i*maxC+j] = 0; }
967  }
968 
969  // necessary branches
970  cmeta->SetBranchAddress( "job", &job );
971  cmeta->SetBranchAddress( "pots", &pots );
972 
973  // extra branches
974  if( cmeta->GetBranch( "beamsim" ) ) cmeta->SetBranchAddress( "beamsim", beamsim );
975  if( cmeta->GetBranch( "physics" ) ) cmeta->SetBranchAddress( "physics", physics );
976  if( cmeta->GetBranch( "physcuts" ) ) cmeta->SetBranchAddress( "physcuts", physcuts );
977  if( cmeta->GetBranch( "tgtcfg" ) ) cmeta->SetBranchAddress( "tgtcfg", tgtcfg );
978  if( cmeta->GetBranch( "horncfg" ) ) cmeta->SetBranchAddress( "horncfg", horncfg );
979  if( cmeta->GetBranch( "dkvolcfg" ) ) cmeta->SetBranchAddress( "dkvolcfg", dkvolcfg );
980  if( cmeta->GetBranch( "beam0x" ) ) cmeta->SetBranchAddress( "beam0x", &beam0x );
981  if( cmeta->GetBranch( "beam0y" ) ) cmeta->SetBranchAddress( "beam0y", &beam0y );
982  if( cmeta->GetBranch( "beam0z" ) ) cmeta->SetBranchAddress( "beam0z", &beam0z );
983  if( cmeta->GetBranch( "beamhwidth" ) ) cmeta->SetBranchAddress( "beamhwidth", &beamhwidth );
984  if( cmeta->GetBranch( "beamvwidth" ) ) cmeta->SetBranchAddress( "beamvwidth", &beamvwidth );
985  if( cmeta->GetBranch( "beamdxdz" ) ) cmeta->SetBranchAddress( "beamdxdz", &beamdxdz );
986  if( cmeta->GetBranch( "beamdydz" ) ) cmeta->SetBranchAddress( "beamdydz", &beamdydz );
987  if( cmeta->GetBranch( "arSize" ) ) cmeta->SetBranchAddress( "arSize", &mArSize );
988  if( cmeta->GetBranch( "location_x" ) ) cmeta->SetBranchAddress( "location_x", location_x );
989  if( cmeta->GetBranch( "location_y" ) ) cmeta->SetBranchAddress( "location_y", location_y );
990  if( cmeta->GetBranch( "location_z" ) ) cmeta->SetBranchAddress( "location_z", location_z );
991  if( cmeta->GetBranch( "location_name" ) ) cmeta->SetBranchAddress( "location_name", location_name );
992 }
genie::hnl::FluxCreator::maxArray
static const int maxArray
Definition: HNLFluxCreator.h:238
genie::hnl::FluxCreator::pots
double pots
how many pot in this job?
Definition: HNLFluxCreator.h:276
genie::hnl::FluxCreator::location_y
double location_y[maxArray]
Definition: HNLFluxCreator.h:284
genie::hnl::FluxCreator::job
int job
beamsim MC job number
Definition: HNLFluxCreator.h:275
genie::hnl::FluxCreator::location_name
char location_name[maxArray *maxC]
Definition: HNLFluxCreator.h:285
genie::hnl::FluxCreator::location_z
double location_z[maxArray]
Definition: HNLFluxCreator.h:284
genie::hnl::FluxCreator::location_x
double location_x[maxArray]
Definition: HNLFluxCreator.h:284
genie::hnl::FluxCreator::maxC
static const int maxC
Definition: HNLFluxCreator.h:238
void FluxCreator::InitialiseTree ( ) const
private

Definition at line 798 of file HNLFluxCreator.cxx.

References anArSize, ancestor_imat, ancestor_ivol, ancestor_nucleus, ancestor_pdg, ancestor_polx, ancestor_poly, ancestor_polz, ancestor_pprodpx, ancestor_pprodpy, ancestor_pprodpz, ancestor_proc, ancestor_startpx, ancestor_startpy, ancestor_startpz, ancestor_startx, ancestor_starty, ancestor_startz, ancestor_stoppx, ancestor_stoppy, ancestor_stoppz, arSize, ctree, decay_mupare, decay_muparpx, decay_muparpy, decay_muparpz, decay_ndecay, decay_necm, decay_nimpwt, decay_norig, decay_ntype, decay_pdpx, decay_pdpy, decay_pdpz, decay_ppdxdz, decay_ppdydz, decay_ppenergy, decay_ppmedium, decay_pppz, decay_ptype, decay_vx, decay_vy, decay_vz, djob, maxArray, maxC, nuray_E, nuray_px, nuray_py, nuray_pz, nuray_wgt, potnum, ppvx, ppvy, ppvz, tgtexit_tgen, tgtexit_tptype, tgtexit_tpx, tgtexit_tpy, tgtexit_tpz, tgtexit_tvx, tgtexit_tvy, tgtexit_tvz, traj_trkpx, traj_trkpy, traj_trkpz, traj_trkx, traj_trky, traj_trkz, and trArSize.

Referenced by MakeTupleFluxEntry().

799 {
800  potnum = 0.0;
801  decay_ptype = 0;
802  decay_vx = 0.0; decay_vy = 0.0; decay_vz = 0.0;
803  decay_pdpx = 0.0; decay_pdpy = 0.0; decay_pdpz = 0.0;
804  decay_nimpwt = 0.0;
805 
806  arSize = 0, anArSize = 0, trArSize = 0;
807  djob = -9999;
808  ppvx = -9999.9, ppvy = -9999.9, ppvz = -9999.9;
809  decay_norig = -9999, decay_ndecay = -9999, decay_ntype = -9999;
810  decay_ppdxdz = -9999.9, decay_ppdydz = -9999.9, decay_pppz = -9999.9, decay_ppenergy = -9999.9;
811  decay_ppmedium = -9999;
812  decay_muparpx = -9999.9, decay_muparpy = -9999.9, decay_muparpz = -9999.9, decay_mupare = -9999.9;
813 
814  tgtexit_tvx = -9999.9, tgtexit_tvy = -9999.9, tgtexit_tvz = -9999.9;
815  tgtexit_tpx = -9999.9, tgtexit_tpy = -9999.9, tgtexit_tpz = -9999.9;
816  tgtexit_tptype = -9999, tgtexit_tgen = -9999;
817 
818  for( int i = 0; i < maxArray; i++ ){
819  nuray_px[i] = -9999.9;
820  nuray_py[i] = -9999.9;
821  nuray_pz[i] = -9999.9;
822  nuray_E[i] = -9999.9;
823  nuray_wgt[i] = -9999.9;
824 
825  ancestor_pdg[i] = -9999;
826  ancestor_startx[i] = -9999.9;
827  ancestor_starty[i] = -9999.9;
828  ancestor_startz[i] = -9999.9;
829  ancestor_startpx[i] = -9999.9;
830  ancestor_startpy[i] = -9999.9;
831  ancestor_startpz[i] = -9999.9;
832  ancestor_stoppx[i] = -9999.9;
833  ancestor_stoppy[i] = -9999.9;
834  ancestor_stoppz[i] = -9999.9;
835  ancestor_polx[i] = -9999.9;
836  ancestor_poly[i] = -9999.9;
837  ancestor_polz[i] = -9999.9;
838  ancestor_pprodpx[i] = -9999.9;
839  ancestor_pprodpy[i] = -9999.9;
840  ancestor_pprodpz[i] = -9999.9;
841  ancestor_nucleus[i] = -9999;
842 
843  traj_trkx[i] = -9999.9;
844  traj_trky[i] = -9999.9;
845  traj_trkz[i] = -9999.9;
846  traj_trkpx[i] = -9999.9;
847  traj_trkpy[i] = -9999.9;
848  traj_trkpz[i] = -9999.9;
849 
850  for( int j = 0; j < maxC; j++ ){
851  ancestor_proc[i*maxC+j] = 0;
852  ancestor_ivol[i*maxC+j] = 0;
853  ancestor_imat[i*maxC+j] = 0;
854  }
855  }
856 
857  // necessary branches
858  ctree->SetBranchAddress( "potnum", &potnum );
859  ctree->SetBranchAddress( "decay_ptype", &decay_ptype );
860  ctree->SetBranchAddress( "decay_vx", &decay_vx );
861  ctree->SetBranchAddress( "decay_vy", &decay_vy );
862  ctree->SetBranchAddress( "decay_vz", &decay_vz );
863  ctree->SetBranchAddress( "decay_pdpx", &decay_pdpx );
864  ctree->SetBranchAddress( "decay_pdpy", &decay_pdpy );
865  ctree->SetBranchAddress( "decay_pdpz", &decay_pdpz );
866  ctree->SetBranchAddress( "decay_necm", &decay_necm );
867  ctree->SetBranchAddress( "decay_nimpwt", &decay_nimpwt );
868 
869  // extra branches
870  if( ctree->GetBranch( "job" ) ) ctree->SetBranchAddress( "job", &djob );
871  if( ctree->GetBranch( "decay_norig" ) ) ctree->SetBranchAddress( "decay_norig", &decay_norig );
872  if( ctree->GetBranch( "decay_ndecay" ) ) ctree->SetBranchAddress( "decay_ndecay", &decay_ndecay );
873  if( ctree->GetBranch( "decay_ntype" ) ) ctree->SetBranchAddress( "decay_ntype", &decay_ntype );
874  if( ctree->GetBranch( "decay_ppdxdz" ) ) ctree->SetBranchAddress( "decay_ppdxdz", &decay_ppdxdz );
875  if( ctree->GetBranch( "decay_ppdydz" ) ) ctree->SetBranchAddress( "decay_ppdydz", &decay_ppdydz );
876  if( ctree->GetBranch( "decay_pppz" ) ) ctree->SetBranchAddress( "decay_pppz", &decay_pppz );
877  if( ctree->GetBranch( "decay_ppenergy" ) ) ctree->SetBranchAddress( "decay_ppenergy", &decay_ppenergy );
878  if( ctree->GetBranch( "decay_ppmedium" ) ) ctree->SetBranchAddress( "decay_ppmedium", &decay_ppmedium );
879  if( ctree->GetBranch( "decay_ptype" ) ) ctree->SetBranchAddress( "decay_ptype", &decay_ptype );
880  if( ctree->GetBranch( "decay_muparpx" ) ) ctree->SetBranchAddress( "decay_muparpx", &decay_muparpx );
881  if( ctree->GetBranch( "decay_muparpy" ) ) ctree->SetBranchAddress( "decay_muparpy", &decay_muparpy );
882  if( ctree->GetBranch( "decay_muparpz" ) ) ctree->SetBranchAddress( "decay_muparpz", &decay_muparpz );
883  if( ctree->GetBranch( "decay_mupare" ) ) ctree->SetBranchAddress( "decay_mupare", &decay_mupare );
884 
885  if( ctree->GetBranch( "nuray_size" ) ) ctree->SetBranchAddress( "nuray_size", &arSize );
886  if( ctree->GetBranch( "nuray_px" ) ) ctree->SetBranchAddress( "nuray_px", nuray_px );
887  if( ctree->GetBranch( "nuray_py" ) ) ctree->SetBranchAddress( "nuray_py", nuray_py );
888  if( ctree->GetBranch( "nuray_pz" ) ) ctree->SetBranchAddress( "nuray_pz", nuray_pz );
889  if( ctree->GetBranch( "nuray_E" ) ) ctree->SetBranchAddress( "nuray_E", nuray_E );
890  if( ctree->GetBranch( "nuray_wgt" ) ) ctree->SetBranchAddress( "nuray_wgt", nuray_wgt );
891 
892  if( ctree->GetBranch( "ancestor_size" ) ) ctree->SetBranchAddress( "ancestor_size", &anArSize );
893  if( ctree->GetBranch( "ancestor_pdg" ) ) ctree->SetBranchAddress( "ancestor_pdg", ancestor_pdg );
894  if( ctree->GetBranch( "ancestor_startx" ) ) ctree->SetBranchAddress( "ancestor_startx", ancestor_startx );
895  if( ctree->GetBranch( "ancestor_starty" ) ) ctree->SetBranchAddress( "ancestor_starty", ancestor_starty );
896  if( ctree->GetBranch( "ancestor_startz" ) ) ctree->SetBranchAddress( "ancestor_startz", ancestor_startz );
897  if( ctree->GetBranch( "ancestor_startpx" ) ) ctree->SetBranchAddress( "ancestor_startpx", ancestor_startpx );
898  if( ctree->GetBranch( "ancestor_startpy" ) ) ctree->SetBranchAddress( "ancestor_startpy", ancestor_starty );
899  if( ctree->GetBranch( "ancestor_startpz" ) ) ctree->SetBranchAddress( "ancestor_startpz", ancestor_startpz );
900  if( ctree->GetBranch( "ancestor_stoppx" ) ) ctree->SetBranchAddress( "ancestor_stoppx", ancestor_stoppx );
901  if( ctree->GetBranch( "ancestor_stoppy" ) ) ctree->SetBranchAddress( "ancestor_stoppy", ancestor_stoppy );
902  if( ctree->GetBranch( "ancestor_stoppz" ) ) ctree->SetBranchAddress( "ancestor_stoppz", ancestor_stoppz );
903  if( ctree->GetBranch( "ancestor_polx" ) ) ctree->SetBranchAddress( "ancestor_polx", ancestor_polx );
904  if( ctree->GetBranch( "ancestor_poly" ) ) ctree->SetBranchAddress( "ancestor_poly", ancestor_poly );
905  if( ctree->GetBranch( "ancestor_polz" ) ) ctree->SetBranchAddress( "ancestor_polz", ancestor_polz );
906  if( ctree->GetBranch( "ancestor_pprodpx" ) ) ctree->SetBranchAddress( "ancestor_pprodpx", ancestor_pprodpx );
907  if( ctree->GetBranch( "ancestor_pprodpy" ) ) ctree->SetBranchAddress( "ancestor_pprodpy", ancestor_pprodpy );
908  if( ctree->GetBranch( "ancestor_pprodpz" ) ) ctree->SetBranchAddress( "ancestor_pprodpz", ancestor_pprodpz );
909  if( ctree->GetBranch( "ancestor_proc" ) ) ctree->SetBranchAddress( "ancestor_proc", ancestor_proc );
910  if( ctree->GetBranch( "ancestor_ivol" ) ) ctree->SetBranchAddress( "ancestor_ivol", ancestor_ivol );
911  if( ctree->GetBranch( "ancestor_imat" ) ) ctree->SetBranchAddress( "ancestor_imat", ancestor_imat );
912 
913  if( ctree->GetBranch( "ppvx" ) ) ctree->SetBranchAddress( "ppvx", &ppvx );
914  if( ctree->GetBranch( "ppvy" ) ) ctree->SetBranchAddress( "ppvy", &ppvy );
915  if( ctree->GetBranch( "ppvz" ) ) ctree->SetBranchAddress( "ppvz", &ppvz );
916 
917  if( ctree->GetBranch( "tgtexit_tvx" ) ) ctree->SetBranchAddress( "tgtexit_tvx", &tgtexit_tvx );
918  if( ctree->GetBranch( "tgtexit_tvy" ) ) ctree->SetBranchAddress( "tgtexit_tvy", &tgtexit_tvy );
919  if( ctree->GetBranch( "tgtexit_tvz" ) ) ctree->SetBranchAddress( "tgtexit_tvz", &tgtexit_tvz );
920 
921  if( ctree->GetBranch( "tgtexit_tpx" ) ) ctree->SetBranchAddress( "tgtexit_tpx", &tgtexit_tpx );
922  if( ctree->GetBranch( "tgtexit_tpy" ) ) ctree->SetBranchAddress( "tgtexit_tpy", &tgtexit_tpy );
923  if( ctree->GetBranch( "tgtexit_tpz" ) ) ctree->SetBranchAddress( "tgtexit_tpz", &tgtexit_tpz );
924 
925  if( ctree->GetBranch( "tgtexit_tptype" ) ) ctree->SetBranchAddress( "tgtexit_tptype", &tgtexit_tptype );
926  if( ctree->GetBranch( "tgtexit_tgen" ) ) ctree->SetBranchAddress( "tgtexit_tgen", &tgtexit_tgen );
927 
928  if( ctree->GetBranch( "traj_size" ) ) ctree->SetBranchAddress( "traj_size", &trArSize );
929  if( ctree->GetBranch( "traj_trkx" ) ) ctree->SetBranchAddress( "traj_trkx", traj_trkx );
930  if( ctree->GetBranch( "traj_trky" ) ) ctree->SetBranchAddress( "traj_trky", traj_trky );
931  if( ctree->GetBranch( "traj_trkz" ) ) ctree->SetBranchAddress( "traj_trkz", traj_trkz );
932  if( ctree->GetBranch( "traj_trkpx" ) ) ctree->SetBranchAddress( "traj_trkpx", traj_trkpx );
933  if( ctree->GetBranch( "traj_trkpy" ) ) ctree->SetBranchAddress( "traj_trkpy", traj_trkpy );
934  if( ctree->GetBranch( "traj_trkpz" ) ) ctree->SetBranchAddress( "traj_trkpz", traj_trkpz );
935 
936 }
genie::hnl::FluxCreator::decay_ptype
int decay_ptype
PDG code of parent.
Definition: HNLFluxCreator.h:242
genie::hnl::FluxCreator::maxArray
static const int maxArray
Definition: HNLFluxCreator.h:238
genie::hnl::FluxCreator::nuray_px
double nuray_px[maxArray]
Definition: HNLFluxCreator.h:256
genie::hnl::FluxCreator::potnum
double potnum
N POT for this SM-v.
Definition: HNLFluxCreator.h:241
genie::hnl::FluxCreator::ancestor_polz
double ancestor_polz[maxArray]
Definition: HNLFluxCreator.h:262
genie::hnl::FluxCreator::decay_vz
double decay_vz
coordinates of prod vtx [cm]
Definition: HNLFluxCreator.h:243
genie::hnl::FluxCreator::ancestor_startz
double ancestor_startz[maxArray]
Definition: HNLFluxCreator.h:259
genie::hnl::FluxCreator::traj_trkpx
double traj_trkpx[maxArray]
Definition: HNLFluxCreator.h:272
genie::hnl::FluxCreator::nuray_E
double nuray_E[maxArray]
Definition: HNLFluxCreator.h:256
genie::hnl::FluxCreator::decay_necm
double decay_necm
SM v CM energy [GeV].
Definition: HNLFluxCreator.h:245
genie::hnl::FluxCreator::ancestor_pprodpx
double ancestor_pprodpx[maxArray]
Definition: HNLFluxCreator.h:263
genie::hnl::FluxCreator::traj_trkz
double traj_trkz[maxArray]
Definition: HNLFluxCreator.h:271
genie::hnl::FluxCreator::decay_nimpwt
double decay_nimpwt
Importance weight from beamsim.
Definition: HNLFluxCreator.h:246
genie::hnl::FluxCreator::ancestor_startx
double ancestor_startx[maxArray]
Definition: HNLFluxCreator.h:259
genie::hnl::FluxCreator::traj_trkpy
double traj_trkpy[maxArray]
Definition: HNLFluxCreator.h:272
genie::hnl::FluxCreator::ancestor_stoppx
double ancestor_stoppx[maxArray]
Definition: HNLFluxCreator.h:261
genie::hnl::FluxCreator::ancestor_starty
double ancestor_starty[maxArray]
Definition: HNLFluxCreator.h:259
genie::hnl::FluxCreator::traj_trkx
double traj_trkx[maxArray]
Definition: HNLFluxCreator.h:271
genie::hnl::FluxCreator::ancestor_startpx
double ancestor_startpx[maxArray]
Definition: HNLFluxCreator.h:260
genie::hnl::FluxCreator::nuray_py
double nuray_py[maxArray]
Definition: HNLFluxCreator.h:256
genie::hnl::FluxCreator::decay_pdpz
double decay_pdpz
final parent momentum [GeV]
Definition: HNLFluxCreator.h:244
genie::hnl::FluxCreator::ancestor_imat
char ancestor_imat[maxArray *maxC]
Definition: HNLFluxCreator.h:265
genie::hnl::FluxCreator::ancestor_proc
char ancestor_proc[maxArray *maxC]
Definition: HNLFluxCreator.h:265
genie::hnl::FluxCreator::ancestor_startpz
double ancestor_startpz[maxArray]
Definition: HNLFluxCreator.h:260
genie::hnl::FluxCreator::traj_trkpz
double traj_trkpz[maxArray]
Definition: HNLFluxCreator.h:272
genie::hnl::FluxCreator::ancestor_pprodpz
double ancestor_pprodpz[maxArray]
Definition: HNLFluxCreator.h:263
genie::hnl::FluxCreator::ancestor_stoppz
double ancestor_stoppz[maxArray]
Definition: HNLFluxCreator.h:261
genie::hnl::FluxCreator::maxC
static const int maxC
Definition: HNLFluxCreator.h:238
genie::hnl::FluxCreator::nuray_wgt
double nuray_wgt[maxArray]
Definition: HNLFluxCreator.h:256
genie::hnl::FluxCreator::ancestor_stoppy
double ancestor_stoppy[maxArray]
Definition: HNLFluxCreator.h:261
genie::hnl::FluxCreator::ancestor_startpy
double ancestor_startpy[maxArray]
Definition: HNLFluxCreator.h:260
genie::hnl::FluxCreator::ancestor_pprodpy
double ancestor_pprodpy[maxArray]
Definition: HNLFluxCreator.h:263
genie::hnl::FluxCreator::ancestor_nucleus
int ancestor_nucleus[maxArray]
Definition: HNLFluxCreator.h:264
genie::hnl::FluxCreator::nuray_pz
double nuray_pz[maxArray]
Definition: HNLFluxCreator.h:256
genie::hnl::FluxCreator::ancestor_ivol
char ancestor_ivol[maxArray *maxC]
Definition: HNLFluxCreator.h:265
genie::hnl::FluxCreator::ancestor_poly
double ancestor_poly[maxArray]
Definition: HNLFluxCreator.h:262
genie::hnl::FluxCreator::ancestor_pdg
int ancestor_pdg[maxArray]
Definition: HNLFluxCreator.h:258
genie::hnl::FluxCreator::traj_trky
double traj_trky[maxArray]
Definition: HNLFluxCreator.h:271
genie::hnl::FluxCreator::ancestor_polx
double ancestor_polx[maxArray]
Definition: HNLFluxCreator.h:262
double FluxCreator::labangle ( double *  x,
double *  par 
)
staticprivate

Definition at line 1928 of file HNLFluxCreator.cxx.

References genie::constants::kPi.

Referenced by CalculateAcceptanceCorrection().

1929 {
1930  double xrad = x[0] * TMath::DegToRad();
1931  double Ehad = par[0], pxhad = par[1], pyhad = par[2], pzhad = par[3];
1932  double phnl = par[4], Ehnl = par[5];
1933 
1934  TLorentzVector p4had( pxhad, pyhad, pzhad, Ehad );
1935  TVector3 boost_vec = p4had.BoostVector(); // beta of parent in lab frame
1936 
1937  // assume phi invariance so create HNL rest-frame momentum along y'z' plane
1938  TLorentzVector pncm( 0.0, phnl * TMath::Sin( xrad ), phnl * TMath::Cos( xrad ), Ehnl );
1939 
1940  // boost into lab frame
1941  pncm.Boost( boost_vec );
1942 
1943  // return lab frame theta wrt parent momentum in deg
1944  double num = pxhad * pncm.X() + pyhad * pncm.Y() + pzhad * pncm.Z();
1945  double den = p4had.P() * pncm.P();
1946  double theta = TMath::ACos( num / den ) * 180.0 / constants::kPi;
1947  return theta;
1948 }
void FluxCreator::LoadConfig ( void  )
private

Definition at line 2042 of file HNLFluxCreator.cxx.

References doPol, fAx1, fAx2, fAz, fB2URotation, fB2UTranslation, fBx1, fBx2, fBz, fCx, fCy, fCz, fDetOffset, fDetRotation, fDoingOldFluxCalc, fFixedPolarisation, fIsConfigLoaded, fIsMajorana, fixPol, fMass, fRadius, fRerollPoints, fScales, fSupplyingBEAM, fU4l2s, genie::Algorithm::GetParam(), genie::Algorithm::GetParamVect(), isParentOnAxis, genie::utils::hnl::IsProdKinematicallyAllowed(), genie::hnl::kHNLProdKaon2Electron, genie::hnl::kHNLProdKaon2Muon, genie::hnl::kHNLProdKaon3Electron, genie::hnl::kHNLProdKaon3Muon, genie::hnl::kHNLProdMuon3Nue, genie::hnl::kHNLProdNeuk3Electron, genie::hnl::kHNLProdNeuk3Muon, genie::hnl::kHNLProdPion2Electron, genie::hnl::kHNLProdPion2Muon, LOG, pDEBUG, and POTScaleWeight.

Referenced by Configure().

2043 {
2044  if( fIsConfigLoaded ) return;
2045 
2046  LOG("HNL", pDEBUG)
2047  << "Loading flux-creation parameters from file...";
2048 
2049  this->GetParam( "HNL-Mass", fMass );
2050  this->GetParamVect( "HNL-LeptonMixing", fU4l2s );
2051  this->GetParam( "HNL-Majorana", fIsMajorana );
2052 
2053  this->GetParamVect( "Near2User_T", fB2UTranslation );
2054  this->GetParamVect( "Near2User_R", fDetRotation );
2055  this->GetParamVect( "Near2Beam_R", fB2URotation );
2056  this->GetParamVect( "DetCentre_User", fDetOffset );
2057 
2058  this->GetParamVect( "ParentPOTScalings", fScales );
2059  this->GetParam( "DoOldFluxCalculation", fDoingOldFluxCalc );
2060  this->GetParam( "RerollPoints", fRerollPoints );
2061  this->GetParam( "CollectionRadius", fRadius );
2062  this->GetParam( "InputFluxesInBEAM", fSupplyingBEAM );
2063  this->GetParam( "IncludePolarisation", doPol );
2064  this->GetParam( "FixPolarisationDirection", fixPol );
2065  this->GetParamVect( "HNL-PolDir", fFixedPolarisation );
2066 
2067  this->GetParam( "IsParentOnAxis", isParentOnAxis );
2068 
2069  fCx = fB2UTranslation.at(0);
2070  fCy = fB2UTranslation.at(1);
2071  fCz = fB2UTranslation.at(2);
2072 
2073  fAx1 = fB2URotation.at(0);
2074  fAz = fB2URotation.at(1);
2075  fAx2 = fB2URotation.at(2);
2076 
2077  fBx1 = fDetRotation.at(0);
2078  fBz = fDetRotation.at(1);
2079  fBx2 = fDetRotation.at(2);
2080 
2081  POTScaleWeight = 1.0;
2082  if( utils::hnl::IsProdKinematicallyAllowed( kHNLProdMuon3Nue ) )
2083  POTScaleWeight = fScales[0]; // all POT contribute
2084  else if( utils::hnl::IsProdKinematicallyAllowed( kHNLProdPion2Muon ) ||
2085  utils::hnl::IsProdKinematicallyAllowed( kHNLProdPion2Electron ) )
2086  POTScaleWeight = fScales[1]; // muons do not contribute
2087  else if( utils::hnl::IsProdKinematicallyAllowed( kHNLProdNeuk3Muon ) ||
2088  utils::hnl::IsProdKinematicallyAllowed( kHNLProdNeuk3Electron ) )
2089  POTScaleWeight = fScales[2]; // only kaons contribute
2090  else if( utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon2Muon ) ||
2091  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon2Electron ) ||
2092  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon3Muon ) ||
2093  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon3Electron ) )
2094  POTScaleWeight = fScales[3]; // only charged kaons contribute
2095 
2096  /*
2097  LOG( "HNL", pDEBUG )
2098  << "Read the following parameters :"
2099  << "\n Mass = " << fMass
2100  << "\n couplings = " << fU4l2s.at(0) << " : " << fU4l2s.at(1) << " : " << fU4l2s.at(2)
2101  << "\n translation = " << fB2UTranslation.at(0) << ", " << fB2UTranslation.at(1) << ", " << fB2UTranslation.at(2)
2102  << "\n rotation = " << fB2URotation.at(0) << ", " << fB2URotation.at(1) << ", " << fB2URotation.at(2)
2103  << "\n isParentOnAxis = " << isParentOnAxis
2104  << "\n POTScaleWeight = " << POTScaleWeight;
2105  */
2106 
2107  fIsConfigLoaded = true;
2108 }
genie::hnl::FluxCreator::fFixedPolarisation
std::vector< double > fFixedPolarisation
Definition: HNLFluxCreator.h:226
genie::hnl::FluxCreator::fB2UTranslation
std::vector< double > fB2UTranslation
Definition: HNLFluxCreator.h:212
genie::hnl::FluxCreator::fDetOffset
std::vector< double > fDetOffset
Definition: HNLFluxCreator.h:214
genie::Algorithm::GetParamVect
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
genie::hnl::FluxCreator::fB2URotation
std::vector< double > fB2URotation
Definition: HNLFluxCreator.h:212
genie::hnl::FluxCreator::fDetRotation
std::vector< double > fDetRotation
Definition: HNLFluxCreator.h:213
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::hnl::FluxCreator::fU4l2s
std::vector< double > fU4l2s
Definition: HNLFluxCreator.h:204
genie::hnl::FluxCreator::fScales
std::vector< double > fScales
Definition: HNLFluxCreator.h:290
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::MakeBBox ( ) const
private

Definition at line 1950 of file HNLFluxCreator.cxx.

References fLx, fLy, fLz, fRadius, LOG, and pFATAL.

Referenced by MakeTupleFluxEntry().

1951 {
1952  if( fRadius < 0.0 ) fRadius = 1.0;
1953  LOG( "HNL", pFATAL )
1954  << "WARNING: This is a bounding box centred at config-given point and radius " << fRadius << " m";
1955 
1956  fLx = fRadius; fLy = fRadius; fLz = fRadius;
1957 }
pFATAL
#define pFATAL
Definition: Messenger.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
FluxContainer FluxCreator::MakeTupleFluxEntry ( int  iEntry,
std::string  finpath 
) const
private

Definition at line 157 of file HNLFluxCreator.cxx.

References genie::hnl::FluxContainer::accCorr, genie::hnl::FluxContainer::acceptance, ApplyUserRotation(), genie::hnl::FluxContainer::boostCorr, CalculateAcceptanceCorrection(), CalculateDetectorAcceptanceSAA(), genie::units::cm, ctree, decay_necm, decay_nimpwt, decay_pdpx, decay_pdpy, decay_pdpz, decay_ptype, decay_vx, decay_vy, decay_vz, genie::hnl::FluxContainer::delay, dynamicScores, e, genie::hnl::FluxContainer::Ecm, genie::hnl::FluxContainer::evtno, fCx, fCy, fCz, fDetOffset, fDetRotation, fDoingOldFluxCalc, fDx, fDy, fDz, fECM, fFirstEntry, fFixedPolarisation, FillNonsense(), genie::PDGLibrary::Find(), fIsMajorana, fIsUsingRootGeom, fixPol, fLepPdg, fLPE, fLPx, fLPy, fLPz, fNuPdg, fNuProdChan, fProdChan, fRerollPoints, fSMECM, fSupplyingBEAM, fTargetPoint, fZm, fZp, GetAngDeviation(), GetProductionProbs(), HNLEnergy(), InitialiseMeta(), InitialiseTree(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), isParentOnAxis, genie::utils::hnl::IsProdKinematicallyAllowed(), genie::hnl::kHNLProdKaon2Electron, genie::hnl::kHNLProdKaon2Muon, genie::hnl::kHNLProdKaon3Electron, genie::hnl::kHNLProdKaon3Muon, genie::hnl::kHNLProdMuon3Nue, genie::hnl::kHNLProdMuon3Numu, genie::hnl::kHNLProdMuon3Nutau, genie::hnl::kHNLProdNeuk3Electron, genie::hnl::kHNLProdNeuk3Muon, genie::hnl::kHNLProdNull, genie::hnl::kHNLProdPion2Electron, genie::hnl::kHNLProdPion2Muon, genie::kPdgHNL, genie::kPdgK0L, genie::kPdgKP, genie::kPdgMuon, genie::kPdgNuE, genie::kPdgNuMu, genie::kPdgPiP, genie::units::kSpeedOfLight, genie::hnl::FluxContainer::lepPdg, LOG, genie::units::m, MakeBBox(), genie::hnl::FluxContainer::nimpwt, genie::units::ns, genie::hnl::FluxContainer::nuEcm, genie::hnl::FluxContainer::nuPdg, genie::hnl::FluxContainer::nuProdChan, OpenFluxInput(), genie::hnl::FluxContainer::p4, genie::utils::print::P4AsString(), genie::hnl::FluxContainer::p4User, parentEnergy, parentMass, parentMomentum, genie::hnl::FluxContainer::parp4, genie::hnl::FluxContainer::parp4User, genie::hnl::FluxContainer::parPdg, pDEBUG, genie::hnl::FluxContainer::pdg, pERROR, pNOTICE, PointToRandomPointInBBox(), genie::hnl::FluxContainer::polz, genie::utils::hnl::ProdAsString(), genie::hnl::FluxContainer::prodChan, genie::RandomGen::RndGen(), genie::units::s, genie::hnl::FluxContainer::startPoint, genie::hnl::FluxContainer::startPointUser, genie::hnl::FluxContainer::targetPoint, genie::hnl::FluxContainer::targetPointUser, genie::utils::print::Vec3AsString(), genie::hnl::FluxContainer::XYWgt, genie::hnl::FluxContainer::zetaMinus, and genie::hnl::FluxContainer::zetaPlus.

Referenced by ProcessEventRecord().

158 {
159  // This method creates 1 HNL from the flux info and saves the information
160  // Essentially, it replaces a SMv with an HNL
161  FluxContainer * tmpGnmf = new FluxContainer();
162  FluxContainer gnmf = *tmpGnmf;
163  delete tmpGnmf;
164 
165  // Open flux input and initialise trees
166  if( iEntry == fFirstEntry ){
167  this->OpenFluxInput( finpath );
168  this->InitialiseTree();
169  this->InitialiseMeta();
170  if(fDoingOldFluxCalc) this->MakeBBox();
171  } else if( iEntry < fFirstEntry ){
172  this->FillNonsense( iEntry, gnmf );
173  return gnmf;
174  }
175 
176  TVector3 dumori( 0.0, 0.0, 0.0 ); // use to rotate VECTORS
177  TVector3 originPoint( -(fCx + fDetOffset.at(0)),
178  -(fCy + fDetOffset.at(1)),
179  -(fCz + fDetOffset.at(2)) ); // use to rotate POINTS. m
180 
181  // All these in m
182  TVector3 fCvec_beam( fCx, fCy, fCz );
183  TVector3 fCvec = this->ApplyUserRotation( fCvec_beam ); // in NEAR coords
184  fLepPdg = 0;
185 
186  ctree->GetEntry(iEntry);
187 
188  // explicitly check if there are any allowed decays for this parent
189  bool canGoForward = true;
190  switch( std::abs( decay_ptype ) ){
191  case kPdgPiP:
192  canGoForward =
193  utils::hnl::IsProdKinematicallyAllowed( kHNLProdPion2Muon ) ||
194  utils::hnl::IsProdKinematicallyAllowed( kHNLProdPion2Electron ); break;
195  case kPdgKP:
196  canGoForward =
197  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon2Muon ) ||
198  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon2Electron ) ||
199  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon3Muon ) ||
200  utils::hnl::IsProdKinematicallyAllowed( kHNLProdKaon3Electron ); break;
201  case kPdgMuon:
202  canGoForward =
203  utils::hnl::IsProdKinematicallyAllowed( kHNLProdMuon3Nue ); break; // SM nus are massless so doesn't matter
204  case kPdgK0L:
205  canGoForward =
206  utils::hnl::IsProdKinematicallyAllowed( kHNLProdNeuk3Muon ) ||
207  utils::hnl::IsProdKinematicallyAllowed( kHNLProdNeuk3Electron ); break;
208  }
209 
210  if( !canGoForward ){
211  this->FillNonsense( iEntry, gnmf ); return gnmf;
212  }
213 
214  // turn cm to m and make origin wrt detector
215  fDx = decay_vx * units::cm / units::m; // BEAM, m
216  fDy = decay_vy * units::cm / units::m;
217  fDz = decay_vz * units::cm / units::m;
218 
219  if( !fSupplyingBEAM ){
220  TVector3 tmpVec( fDx, fDy, fDz ); // NEAR
221  tmpVec = this->ApplyUserRotation( tmpVec, false ); // BEAM
222  fDx = tmpVec.X(); fDy = tmpVec.Y(); fDz = tmpVec.Z();
223  }
224 
225  TVector3 fDvec( fDx, fDy, fDz ); // in BEAM coords
226  TVector3 fDvec_beam = this->ApplyUserRotation( fDvec, true ); // in NEAR coords
227 
228  TVector3 fDvec_user( fDvec_beam.X() - fCx, fDvec_beam.Y() - fCy, fDvec_beam.Z() - fCz ); // in USER coords
229  fDvec_user = this->ApplyUserRotation( fDvec_user, originPoint, fDetRotation, false );
230 
231  LOG( "HNL", pDEBUG )
232  << "\nIn BEAM coords, fDvec = " << utils::print::Vec3AsString( &fDvec )
233  << "\nIn NEAR coords, fDvec = " << utils::print::Vec3AsString( &fDvec_beam );
234 
235  TVector3 detO_beam( fCvec_beam.X() - fDvec_beam.X(),
236  fCvec_beam.Y() - fDvec_beam.Y(),
237  fCvec_beam.Z() - fDvec_beam.Z() ); // separation in NEAR coords
238  TVector3 detO( fCvec.X() - fDvec.X(),
239  fCvec.Y() - fDvec.Y(),
240  fCvec.Z() - fDvec.Z() ); // separation in BEAM coords
241  TVector3 detO_user( detO_beam.X(), detO_beam.Y(), detO_beam.Z() );
242 
243  detO_user = this->ApplyUserRotation( detO_user, dumori, fDetRotation, false ); // tgt-hall --> det
244 
245  double acc_saa = this->CalculateDetectorAcceptanceSAA( detO_user );
246 
247  // set parent mass
248 
249  double dpdpx = decay_pdpx, dpdpy = decay_pdpy, dpdpz = decay_pdpz; // BEAM GeV
250  if( !fSupplyingBEAM ){
251  TVector3 tmpVec( dpdpx, dpdpy, dpdpz ); // NEAR
252  tmpVec = this->ApplyUserRotation( tmpVec, false ); // BEAM
253  dpdpx = tmpVec.X(); dpdpy = tmpVec.Y(); dpdpz = tmpVec.Z();
254  }
255 
256  switch( std::abs( decay_ptype ) ){
257  case kPdgPiP: case kPdgKP: case kPdgMuon: case kPdgK0L:
258  parentMass = PDGLibrary::Instance()->Find(decay_ptype)->Mass(); break;
259  default:
260  LOG( "HNL", pERROR ) << "Parent with PDG code " << decay_ptype << " not handled!"
261  << "\n\tProceeding, but results are possibly unphysical.";
262  parentMass = PDGLibrary::Instance()->Find(decay_ptype)->Mass(); break;
263  }
264  parentMomentum = std::sqrt( dpdpx*dpdpx + dpdpy*dpdpy + dpdpz*dpdpz );
265  parentEnergy = std::sqrt( parentMass*parentMass + parentMomentum*parentMomentum );
266 
267  TLorentzVector p4par = ( isParentOnAxis ) ?
268  TLorentzVector( parentMomentum * (detO.Unit()).X(),
269  parentMomentum * (detO.Unit()).Y(),
270  parentMomentum * (detO.Unit()).Z(),
271  parentEnergy ) :
272  TLorentzVector( dpdpx, dpdpy, dpdpz, parentEnergy ); // in BEAM coords
273 
274  TLorentzVector p4par_beam( dpdpx, dpdpy, dpdpz, parentEnergy ); // in BEAM coords
275  TVector3 p3par_beam = p4par_beam.Vect();
276  TVector3 p3par_near = this->ApplyUserRotation( p3par_beam, true );
277  TLorentzVector p4par_near( p3par_near.X(), p3par_near.Y(), p3par_near.Z(), parentEnergy ); // in NEAR coords
278  LOG( "HNL", pDEBUG )
279  << "\nIn BEAM coords: p3par_beam = " << utils::print::Vec3AsString( &p3par_beam )
280  << "\nIn NEAR coords: p3par_near = " << utils::print::Vec3AsString( &p3par_near );
281  if( !isParentOnAxis ){
282  // rotate p4par to NEAR coordinates
283  TVector3 tmpv3 = ApplyUserRotation( p4par.Vect(), true );
284  p4par.SetPxPyPzE( tmpv3.Px(), tmpv3.Py(), tmpv3.Pz(), p4par.E() );
285  }
286 
287  TLorentzVector p4par_user = p4par_near;
288  // rotate it to user coords
289  TVector3 ppar_user = p4par_user.Vect();
290  ppar_user = this->ApplyUserRotation( ppar_user, dumori, fDetRotation, false ); // tgt-hall --> det
291  p4par_user.SetPxPyPzE( ppar_user.Px(), ppar_user.Py(), ppar_user.Pz(), p4par_user.E() );
292 
293  TVector3 boost_beta = p4par_beam.BoostVector(); // in BEAM coords
294 
295  // now calculate which decay channel produces the HNL.
296  dynamicScores = this->GetProductionProbs( decay_ptype );
297  assert( dynamicScores.size() > 0 );
298 
299  if( dynamicScores.find( kHNLProdNull ) != dynamicScores.end() ){ // exists kin allowed channel but 0 coupling
300  this->FillNonsense( iEntry, gnmf ); return gnmf;
301  }
302 
303  RandomGen * rnd = RandomGen::Instance();
304  double score = rnd->RndGen().Uniform( 0.0, 1.0 );
305  HNLProd_t prodChan;
306  // compare with cumulative prob. If < 1st in map, pick 1st chan. If >= 1st and < (1st+2nd), pick 2nd, etc
307 
308  unsigned int imap = 0; double s1 = 0.0;
309  std::map< HNLProd_t, double >::iterator pdit = dynamicScores.begin();
310  while( score >= s1 && pdit != dynamicScores.end() ){
311  s1 += (*pdit).second;
312  if( parentMass > 0.495 ){
313  LOG( "HNL", pDEBUG )
314  << "(*pdit).first = " << utils::hnl::ProdAsString( (*pdit).first )
315  << " : (*pdit).second = " << (*pdit).second;
316  }
317  if( score >= s1 ){
318  imap++; pdit++;
319  }
320  }
321  assert( imap < dynamicScores.size() ); // should have decayed to *some* HNL
322  prodChan = (*pdit).first;
323 
324  // bookkeep this
325  fProdChan = static_cast<int>(prodChan);
326  switch( prodChan ){
327  case kHNLProdNeuk3Electron: fNuProdChan = 1; fNuPdg = kPdgNuE; break;
328  case kHNLProdNeuk3Muon: fNuProdChan = 2; fNuPdg = kPdgNuMu; break;
329  case kHNLProdKaon2Electron: fNuProdChan = 4; fNuPdg = kPdgNuE; break;
330  case kHNLProdKaon2Muon: fNuProdChan = 3; fNuPdg = kPdgNuMu; break;
331  case kHNLProdKaon3Electron: fNuProdChan = 6; fNuPdg = kPdgNuE; break;
332  case kHNLProdKaon3Muon: fNuProdChan = 5; fNuPdg = kPdgNuMu; break;
333  case kHNLProdMuon3Nue:
334  case kHNLProdMuon3Numu:
335  case kHNLProdMuon3Nutau: fNuProdChan = 9; fNuPdg = kPdgNuE; break;
336  case kHNLProdPion2Electron: fNuProdChan = 8; fNuPdg = kPdgNuE; break;
337  case kHNLProdPion2Muon: fNuProdChan = 7; fNuPdg = kPdgNuMu; break;
338  default: fNuProdChan = -999; fNuPdg = -999; break;
339  }
340 
341  LOG( "HNL", pDEBUG )
342  << "Selected channel: " << utils::hnl::ProdAsString( prodChan );
343 
344  // decay channel specified, now time to make kinematics
345  TLorentzVector p4HNL_rest = HNLEnergy( prodChan, p4par );
346  // this is a random direction rest-frame HNL.
347  fECM = p4HNL_rest.E();
348 
349  // we will now boost detO into rest frame, force rest to point to the new direction, boost the result, and compare the boost corrections
350  double boost_correction_two = 0.0;
351 
352  // 17-Jun-22: Notice the time component needs to be nonzero to get this to work!
353 
354  // first guess: betaHNL ~= 1 . Do the Lorentz boosts knocking betaHNL downwards until we hit det centre
355  double betaMag = boost_beta.Mag();
356  double gamma = std::sqrt( 1.0 / ( 1.0 - betaMag * betaMag ) );
357  double betaLab = 1.0; // first guess
358 
359  // now make a TLorentzVector in lab frame to boost back to rest.
360  double timeBit = detO.Mag(); // / units::kSpeedOfLight ; // s
361  TLorentzVector detO_4v( detO.X(), detO.Y(), detO.Z(), timeBit ); detO_4v.Boost( -boost_beta ); // BEAM with BEAM
362  TVector3 detO_rest_unit = (detO_4v.Vect()).Unit();
363  TLorentzVector p4HNL_rest_good( p4HNL_rest.P() * detO_rest_unit.X(),
364  p4HNL_rest.P() * detO_rest_unit.Y(),
365  p4HNL_rest.P() * detO_rest_unit.Z(),
366  p4HNL_rest.E() );
367 
368  double pLep_rest = std::sqrt( fLPx*fLPx + fLPy*fLPy + fLPz*fLPz );
369  TLorentzVector p4Lep_rest_good( -1.0 * pLep_rest * detO_rest_unit.X(),
370  -1.0 * pLep_rest * detO_rest_unit.Y(),
371  -1.0 * pLep_rest * detO_rest_unit.Z(),
372  fLPE );
373 
374  // boost HNL into lab frame!
375 
376  TLorentzVector p4HNL_good = p4HNL_rest_good;
377  p4HNL_good.Boost( boost_beta );
378  boost_correction_two = p4HNL_good.E() / p4HNL_rest.E();
379 
380  TVector3 detO_unit = detO.Unit(); // BEAM
381 
382  TVector3 p4HNL_good_vect = p4HNL_good.Vect();
383  TVector3 p4HNL_good_unit = p4HNL_good_vect.Unit();
384 
385  // now calculate how far away from target point we are.
386  // dist = || detO x p4HNL || / || p4HNL_good || where x == cross product
387  TVector3 distNum = detO.Cross( p4HNL_good_unit );
388  double dist = distNum.Mag(); // m
389 
390  double prevDist = 2.0 * dist;
391 
392  while( betaLab > 0.0 && ( dist > 1.0e-3 && dist < prevDist &&
393  std::abs(dist - prevDist) > 1.0e-3 * prevDist ) ){ // 1mm tolerance
394 
395  // that didn't work. Knock betaLab down a little bit and try again.
396  prevDist = dist;
397  betaLab -= 1.0e-4;
398  timeBit = detO.Mag() / ( betaLab );
399  detO_4v.SetXYZT( detO.X(), detO.Y(), detO.Z(), timeBit );
400  detO_4v.Boost( -boost_beta );
401  detO_rest_unit = (detO_4v.Vect()).Unit();
402  p4HNL_rest_good.SetPxPyPzE( p4HNL_rest.P() * detO_rest_unit.X(),
403  p4HNL_rest.P() * detO_rest_unit.Y(),
404  p4HNL_rest.P() * detO_rest_unit.Z(),
405  p4HNL_rest.E() );
406 
407  // boost into lab frame
408  p4HNL_good = p4HNL_rest_good;
409  p4HNL_good.Boost( boost_beta );
410 
411  detO_unit = detO.Unit();
412  p4HNL_good_vect = p4HNL_good.Vect();
413  p4HNL_good_unit = p4HNL_good_vect.Unit();
414 
415  distNum = detO.Cross( p4HNL_good_unit );
416  dist = distNum.Mag(); // m
417  }
418 
419  // but we don't care about that. We just want to obtain a proxy for betaHNL in lab frame.
420  // Then we can use the dk2nu-style formula modified for betaHNL!
421 
422  /*
423  * it is NOT sufficient to boost this into lab frame!
424  * Only a small portion of the CM decays can possibly reach the detector,
425  * imposing a constraint on the allowed directions of p4HNL_rest.
426  * You will miscalculate the HNL energy if you just Boost here.
427  */
428  // explicitly calculate the boost correction to lab-frame energy
429  // in a dk2nu-like fashion. See bsim::CalcEnuWgt()
430  //double betaHNL = p4HNL_rest.P() / p4HNL_rest.E();
431  double betaHNL = p4HNL_good.P() / p4HNL_good.E();
432  double costh_pardet = 0.0;
433  double boost_correction = 0.0;
434  if( parentMomentum > 0.0 ){
435  costh_pardet = ( p4par_beam.X() * detO.X() +
436  p4par_beam.Y() * detO.Y() +
437  p4par_beam.Z() * detO.Z() ) / ( parentMomentum * detO.Mag() );
438  if( costh_pardet < -1.0 ) costh_pardet = -1.0;
439  if( costh_pardet > 1.0 ) costh_pardet = 1.0;
440  boost_correction = 1.0 / ( gamma * ( 1.0 - betaMag * betaHNL * costh_pardet ) );
441  // assume boost is on z' direction where z' = parent momentum direction, subbing betaMag ==> betaMag * costh_pardet
442  if( true && boost_correction * p4HNL_rest.E() > p4HNL_rest.M() ) {
443  boost_correction = 1.0 / ( gamma * ( 1.0 - betaMag * betaHNL * costh_pardet ) );
444  } else {
445  boost_correction = p4HNL_good.E() / p4HNL_rest_good.E();
446  }
447  }
448 
449  assert( boost_correction > 0.0 && boost_correction_two > 0.0 );
450 
451  // so now we have the random decay. Direction = parent direction, energy = what we calculated
452  double EHNL = p4HNL_rest.E() * boost_correction;
453  double MHNL = p4HNL_rest.M();
454  double PHNL = std::sqrt( EHNL * EHNL - MHNL * MHNL );
455  TVector3 pdu = ( p4par.Vect() ).Unit(); // in BEAM coords
456  TLorentzVector p4HNL_rand( PHNL * pdu.X(), PHNL * pdu.Y(), PHNL * pdu.Z(), EHNL );
457 
458  // find random point in BBox and force momentum to point to that point
459 
460  double FDx = fDvec_beam.X();
461  double FDy = fDvec_beam.Y();
462  double FDz = fDvec_beam.Z();
463 
464  TVector3 absolutePoint = this->PointToRandomPointInBBox( ); // in NEAR coords, m
465  TVector3 fRVec_beam( absolutePoint.X() - FDx, absolutePoint.Y() - FDy, absolutePoint.Z() - FDz ); // NEAR, m
466  // rotate it and get unit
467  TVector3 fRVec_unit = (this->ApplyUserRotation( fRVec_beam )).Unit(); // BEAM, m/m
468  TVector3 fRVec_actualBeam = this->ApplyUserRotation( fRVec_beam ); // BEAM, m
469  TVector3 fRVec_user = this->ApplyUserRotation( fRVec_beam, dumori, fDetRotation, false ); // USER m
470 
471 
472  TVector3 rRVec_near( fRVec_beam.X() + FDx, fRVec_beam.Y() + FDy, fRVec_beam.Z() + FDz );
473  TVector3 rRVec_beam = this->ApplyUserRotation( rRVec_near );
474  TVector3 rRVec_user = this->ApplyUserRotation( rRVec_near, dumori, fDetRotation, false );
475  /*
476  rRVec_user.SetXYZ( rRVec_user.X() + (fCx + fDetOffset.at(0)),
477  rRVec_user.Y() + (fCy + fDetOffset.at(1)),
478  rRVec_user.Z() + (fCz + fDetOffset.at(2)) );
479  */
480  rRVec_user.SetXYZ( rRVec_user.X() + (fCx),
481  rRVec_user.Y() + (fCy),
482  rRVec_user.Z() + (fCz) );
483 
484  // force HNL to point along this direction
485  TLorentzVector p4HNL( p4HNL_rand.P() * fRVec_unit.X(),
486  p4HNL_rand.P() * fRVec_unit.Y(),
487  p4HNL_rand.P() * fRVec_unit.Z(),
488  p4HNL_rand.E() ); // in BEAM coords
489 
490  TVector3 pHNL_near = this->ApplyUserRotation( p4HNL.Vect(), true ); // BEAM --> NEAR coords
491  TLorentzVector p4HNL_near( pHNL_near.X(), pHNL_near.Y(), pHNL_near.Z(), p4HNL.E() );
492 
493  LOG( "HNL", pDEBUG )
494  << "\nRandom: " << utils::print::P4AsString( &p4HNL_rand )
495  << "\nPointed [NEAR]: " << utils::print::P4AsString( &p4HNL_near )
496  << "\nRest: " << utils::print::P4AsString( &p4HNL_rest );
497 
498  // update polarisation
499  TLorentzVector p4Lep_good = p4Lep_rest_good; // in parent rest frame
500  p4Lep_good.Boost( boost_beta ); // in lab frame
501  TVector3 boost_beta_HNL = p4HNL_near.BoostVector();
502  p4Lep_good.Boost( -boost_beta_HNL ); // in HNL rest frame
503 
504  fLPx = ( fixPol ) ? fFixedPolarisation.at(0) : p4Lep_good.Px() / p4Lep_good.P();
505  fLPy = ( fixPol ) ? fFixedPolarisation.at(1) : p4Lep_good.Py() / p4Lep_good.P();
506  fLPz = ( fixPol ) ? fFixedPolarisation.at(2) : p4Lep_good.Pz() / p4Lep_good.P();
507 
508  // calculate acceptance correction
509  // first, get minimum and maximum deviation from parent momentum to hit detector in degrees
510  double zm = 0.0, zp = 0.0;
511  if( fIsUsingRootGeom ){
512  this->GetAngDeviation( p4par_near, detO_beam, zm, zp ); // using NEAR and NEAR
513  } else { // !fIsUsingRootGeom
514  zm = ( isParentOnAxis ) ? 0.0 : this->GetAngDeviation( p4par_near, detO_beam, false );
515  zp = this->GetAngDeviation( p4par_near, detO_beam, true );
516  }
517 
518  if( zm == -999.9 && zp == 999.9 ){
519  this->FillNonsense( iEntry, gnmf ); return gnmf;
520  }
521 
522  if( isParentOnAxis ){
523  double tzm = zm, tzp = zp;
524  zm = 0.0;
525  zp = (tzp - tzm)/2.0; // 1/2 * angular opening
526  }
527 
528  fSMECM = decay_necm;
529  fZm = zm; fZp = zp;
530  double accCorr = this->CalculateAcceptanceCorrection( p4par, p4HNL_rest, decay_necm, zm, zp );
531  //if( !fDoingOldFluxCalc ){
532  if( fRerollPoints ){
533  // if accCorr == 0 then we must ~bail and find the next event.~
534  // We must reroll the point a bunch of times. Then we can skip.
535  // Ideally we'd be able to tell how much detector lives within the reachable region [zm, zp]
536  int iAccFail = 0; const int iAccFailBail = 10;
537  while( iAccFail < iAccFailBail && accCorr == 0.0 ){
538  LOG( "HNL", pNOTICE )
539  << "Point with separation " << utils::print::Vec3AsString( &fRVec_beam ) << " is unreachable "
540  << "by HNL from parent with momentum " << utils::print::P4AsString( &p4par ) << " !"
541  << "\nRerolling point. This is the " << iAccFail << "th try out of " << iAccFailBail;
542 
543  // find random point in BBox and force momentum to point to that point
544  // first, separation in beam frame
545  absolutePoint = this->PointToRandomPointInBBox( ); // always NEAR, m
546  /*
547  fRVec_beam.SetXYZ( absolutePoint.X() - (fCx + fDetOffset.at(0)),
548  absolutePoint.Y() - (fCy + fDetOffset.at(1)),
549  absolutePoint.Z() - (fCz + fDetOffset.at(2)) );
550  */
551  fRVec_beam.SetXYZ( absolutePoint.X() - (fCx),
552  absolutePoint.Y() - (fCy),
553  absolutePoint.Z() - (fCz) );
554  // rotate it and get unit
555  fRVec_unit = (this->ApplyUserRotation( fRVec_beam )).Unit(); // BEAM
556  // force HNL to point along this direction
557  p4HNL.SetPxPyPzE( p4HNL_rand.P() * fRVec_unit.X(),
558  p4HNL_rand.P() * fRVec_unit.Y(),
559  p4HNL_rand.P() * fRVec_unit.Z(),
560  p4HNL_rand.E() ); // BEAM
561 
562  pHNL_near = this->ApplyUserRotation( p4HNL.Vect(), true ); // NEAR
563  p4HNL_near.SetPxPyPzE( pHNL_near.X(), pHNL_near.Y(), pHNL_near.Z(), p4HNL.E() );
564 
565  LOG( "HNL", pDEBUG )
566  << "\nRandom: " << utils::print::P4AsString( &p4HNL_rand )
567  << "\nPointed: " << utils::print::P4AsString( &p4HNL )
568  << "\nRest: " << utils::print::P4AsString( &p4HNL_rest );
569 
570  // update polarisation
571  p4Lep_good = p4Lep_rest_good; // in parent rest frame
572  p4Lep_good.Boost( boost_beta ); // in lab frame
573  boost_beta_HNL = p4HNL_near.BoostVector(); // NEAR coords
574  p4Lep_good.Boost( -boost_beta_HNL ); // in HNL rest frame
575 
576  fLPx = ( fixPol ) ? fFixedPolarisation.at(0) : p4Lep_good.Px() / p4Lep_good.P();
577  fLPy = ( fixPol ) ? fFixedPolarisation.at(1) : p4Lep_good.Py() / p4Lep_good.P();
578  fLPz = ( fixPol ) ? fFixedPolarisation.at(2) : p4Lep_good.Pz() / p4Lep_good.P();
579 
580  // calculate acceptance correction
581  // first, get minimum and maximum deviation from parent momentum to hit detector in degrees
582  zm = 0.0; zp = 0.0;
583  if( fIsUsingRootGeom ){
584  this->GetAngDeviation( p4par_near, detO_beam, zm, zp ); // NEAR and NEAR
585  } else { // !fIsUsingRootGeom
586  zm = ( isParentOnAxis ) ? 0.0 : this->GetAngDeviation( p4par_near, detO_beam, false );
587  zp = this->GetAngDeviation( p4par_near, detO_beam, true );
588  }
589 
590  if( zm == -999.9 && zp == 999.9 ){
591  this->FillNonsense( iEntry, gnmf ); return gnmf;
592  }
593 
594  if( isParentOnAxis ){
595  double tzm = zm, tzp = zp;
596  zm = 0.0;
597  zp = (tzp - tzm)/2.0; // 1/2 * angular opening
598  }
599 
600  accCorr = this->CalculateAcceptanceCorrection( p4par_near, p4HNL_rest, decay_necm, zm, zp );
601  iAccFail++;
602  }
603  }
604  if( accCorr == 0.0 ){ // NOW we can give up and return.
605  this->FillNonsense( iEntry, gnmf ); return gnmf;
606  }
607 
608  // also have to factor in boost correction itself... that's same as energy boost correction squared
609  // which means a true acceptance of...
610  double acceptance = acc_saa * boost_correction * boost_correction * accCorr;
611 
612  // finally, a delay calculation
613  // if SMv arrives at t=0, then HNL arrives at t = c * ( 1 - beta_HNL ) / L
614 
615  double detDist = std::sqrt( detO.X() * detO.X() +
616  detO.Y() * detO.Y() +
617  detO.Z() * detO.Z() ); // m
618  const double kSpeedOfLightNs = units::kSpeedOfLight * units::ns / units::s; // m / ns
619  double delay = detDist / kSpeedOfLightNs * ( 1.0 / betaHNL - 1.0 );
620  delay *= units::ns / units::s;
621 
622  /*
623  LOG( "HNL", pDEBUG )
624  << "\ndetDist = " << detDist << " [m]"
625  << "\nbetaHNL = " << betaHNL
626  << "\ndelay = " << delay << " [ns]";
627  */
628 
629  // write 4-position all this happens at
630 
631  double dxvx = decay_vx, dxvy = decay_vy, dxvz = decay_vz;
632 
633  if( !fSupplyingBEAM ){
634  TVector3 tmpVec( dxvx, dxvy, dxvz ); // NEAR
635  tmpVec = this->ApplyUserRotation( tmpVec, false ); // BEAM
636  dxvx = tmpVec.X(); dxvy = tmpVec.Y(); dxvz = tmpVec.Z();
637  }
638 
639  TLorentzVector x4HNL_beam( dxvx, dxvy, dxvz, delay ); // in cm, ns, BEAM coords
640  TVector3 x3HNL_beam = x4HNL_beam.Vect();
641  TVector3 x3HNL_near = this->ApplyUserRotation( x3HNL_beam, true );
642  TLorentzVector x4HNL_near( x3HNL_near.X(), x3HNL_near.Y(), x3HNL_near.Z(), delay );
643 
644  TLorentzVector x4HNL( fDvec_user.X(), fDvec_user.Y(), fDvec_user.Z(), delay ); // in m, ns, USER coords
645  TLorentzVector x4HNL_cm( units::m / units::cm * x4HNL.X(),
646  units::m / units::cm * x4HNL.Y(),
647  units::m / units::cm * x4HNL.Z(), delay ); // in cm, ns, USER
648 
649  LOG( "HNL", pDEBUG ) << "Filling gnmf...";
650 
651  // fill all the flux stuff now!
652 
653  int typeMod = 1;
654  if( !fIsMajorana ) typeMod = ( decay_ptype > 0 ) ? 1 : -1;
655  // fix for muons which are backwards...
656  int parPdg = decay_ptype;
657  if( std::abs(parPdg) == kPdgMuon ) typeMod *= -1;
658 
659  gnmf.evtno = iEntry;
660 
661  gnmf.pdg = typeMod * kPdgHNL;
662  gnmf.parPdg = parPdg;
663  gnmf.lepPdg = fLepPdg;
664  gnmf.nuPdg = typeMod * fNuPdg;
665 
666  gnmf.prodChan = fProdChan;
667  gnmf.nuProdChan = fNuProdChan;
668 
669  gnmf.startPoint.SetXYZ( fDvec_beam.X(), fDvec_beam.Y(), fDvec_beam.Z() ); // NEAR m
670  gnmf.targetPoint.SetXYZ( fTargetPoint.X(), fTargetPoint.Y(), fTargetPoint.Z() ); // NEAR m
671  gnmf.startPointUser.SetXYZ( fDvec_user.X() - fCx, fDvec_user.Y() - fCy, fDvec_user.Z() - fCz ); // USER m
672  gnmf.targetPointUser.SetXYZ( fTargetPoint.X() - fCx, fTargetPoint.Y() - fCy, fTargetPoint.Z() - fCz ); // USER m
673  gnmf.delay = delay; // ns
674 
675  gnmf.polz.SetXYZ( fLPx, fLPy, fLPz );
676 
677  TLorentzVector p4HNL_user = p4HNL_near;
678  // rotate it to user coords
679  TVector3 pHNL_user = p4HNL_user.Vect();
680  pHNL_user = this->ApplyUserRotation( pHNL_user, dumori, fDetRotation, false ); // tgt-hall --> det
681  p4HNL_user.SetPxPyPzE( pHNL_user.Px(), pHNL_user.Py(), pHNL_user.Pz(), p4HNL_user.E() );
682 
683  gnmf.p4 = p4HNL_near;
684  gnmf.parp4 = p4par_near;
685  gnmf.p4User = p4HNL_user;
686  gnmf.parp4User = p4par_user;
687 
688  gnmf.Ecm = fECM;
689  gnmf.nuEcm = fSMECM;
690 
691  gnmf.XYWgt = acc_saa;
692  gnmf.boostCorr = p4HNL.E() / fECM;
693  gnmf.accCorr = accCorr;
694  gnmf.zetaMinus = fZm;
695  gnmf.zetaPlus = fZp;
696  gnmf.acceptance = acceptance;
697  gnmf.nimpwt = decay_nimpwt;
698 
699  return gnmf;
700 
701 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
genie::hnl::HNLProd_t
enum genie::hnl::t_HNLProd HNLProd_t
genie::hnl::FluxContainer::parp4
TLorentzVector parp4
parent momentum at HNL production in NEAR coords [GeV/c]
Definition: HNLFluxContainer.h:94
genie::hnl::FluxCreator::fFixedPolarisation
std::vector< double > fFixedPolarisation
Definition: HNLFluxCreator.h:226
genie::hnl::FluxCreator::decay_ptype
int decay_ptype
PDG code of parent.
Definition: HNLFluxCreator.h:242
genie::hnl::FluxContainer::Ecm
double Ecm
Parent rest-frame energy of HNL [GeV].
Definition: HNLFluxContainer.h:98
genie::kPdgNuE
const int kPdgNuE
Definition: PDGCodes.h:28
genie::hnl::FluxContainer::zetaMinus
double zetaMinus
minimum angular deviation from parent momentum to reach detector [deg]
Definition: HNLFluxContainer.h:105
genie::hnl::FluxContainer::nuEcm
double nuEcm
Parent rest-frame energy of equivalent SM neutrino [GeV].
Definition: HNLFluxContainer.h:99
genie::hnl::FluxCreator::PointToRandomPointInBBox
TVector3 PointToRandomPointInBBox() const
Definition: HNLFluxCreator.cxx:1290
pERROR
#define pERROR
Definition: Messenger.h:59
genie::hnl::FluxContainer::delay
double delay
delay HNL would have wrt SMv [ns]
Definition: HNLFluxContainer.h:89
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::hnl::FluxCreator::fDetOffset
std::vector< double > fDetOffset
Definition: HNLFluxCreator.h:214
genie::hnl::FluxCreator::decay_vz
double decay_vz
coordinates of prod vtx [cm]
Definition: HNLFluxCreator.h:243
genie::hnl::FluxCreator::GetAngDeviation
double GetAngDeviation(TLorentzVector p4par, TVector3 detO, bool seekingMax) const
Definition: HNLFluxCreator.cxx:1348
genie::hnl::FluxContainer::startPoint
TVector3 startPoint
parent decay vertex in NEAR coords [m]
Definition: HNLFluxContainer.h:85
genie::hnl::FluxContainer::lepPdg
int lepPdg
PDG code of lepton produced with HNL on parent decay.
Definition: HNLFluxContainer.h:79
genie::hnl::FluxContainer::polz
TVector3 polz
HNL polarisation vector, in HNL rest frame, in NEAR coords.
Definition: HNLFluxContainer.h:91
genie::hnl::FluxContainer::parPdg
int parPdg
parent PDG code
Definition: HNLFluxContainer.h:78
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::hnl::FluxCreator::CalculateDetectorAcceptanceSAA
double CalculateDetectorAcceptanceSAA(TVector3 detO) const
Definition: HNLFluxCreator.cxx:2015
genie::hnl::FluxContainer
A GENIE flux container specific for HNL containers. Based on the dk2nu flux paradigm and genie::flux:...
Definition: HNLFluxContainer.h:62
genie::utils::print::P4AsString
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
genie::hnl::FluxContainer::zetaPlus
double zetaPlus
maximum angular deviation from parent momentum to reach detector [deg]
Definition: HNLFluxContainer.h:106
genie::units::ns
static constexpr double ns
Definition: Units.h:98
genie::hnl::FluxContainer::startPointUser
TVector3 startPointUser
parent decay vertex in USER coords [m]
Definition: HNLFluxContainer.h:87
genie::hnl::FluxCreator::decay_necm
double decay_necm
SM v CM energy [GeV].
Definition: HNLFluxCreator.h:245
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::hnl::FluxContainer::parp4User
TLorentzVector parp4User
parent momentum at HNL production in USER coords [GeV/c]
Definition: HNLFluxContainer.h:96
genie::hnl::FluxCreator::decay_nimpwt
double decay_nimpwt
Importance weight from beamsim.
Definition: HNLFluxCreator.h:246
genie::hnl::FluxCreator::CalculateAcceptanceCorrection
double CalculateAcceptanceCorrection(TLorentzVector p4par, TLorentzVector p4HNL, double SMECM, double zm, double zp) const
Definition: HNLFluxCreator.cxx:1803
genie::hnl::FluxCreator::fDetRotation
std::vector< double > fDetRotation
Definition: HNLFluxCreator.h:213
genie::hnl::FluxCreator::dynamicScores
std::map< genie::hnl::HNLProd_t, double > dynamicScores
Definition: HNLFluxCreator.h:188
e
const double e
Definition: gUpMuFluxGen.cxx:165
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::FluxContainer::p4
TLorentzVector p4
HNL momentum in NEAR coords [GeV/c].
Definition: HNLFluxContainer.h:93
genie::hnl::FluxContainer::p4User
TLorentzVector p4User
HNL momentum in USER coords [GeV/c].
Definition: HNLFluxContainer.h:95
genie::hnl::FluxContainer::prodChan
int prodChan
Decay mode that produced HNL.
Definition: HNLFluxContainer.h:82
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgHNL
const int kPdgHNL
Definition: PDGCodes.h:224
genie::hnl::FluxCreator::ApplyUserRotation
TVector3 ApplyUserRotation(TVector3 vec, bool doBackwards=false) const
Definition: HNLFluxCreator.cxx:1959
genie::hnl::FluxCreator::OpenFluxInput
void OpenFluxInput(std::string finpath) const
Definition: HNLFluxCreator.cxx:741
genie::hnl::FluxCreator::GetProductionProbs
std::map< genie::hnl::HNLProd_t, double > GetProductionProbs(int parPDG) const
Definition: HNLFluxCreator.cxx:1027
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::utils::hnl::ProdAsString
string ProdAsString(genie::hnl::HNLProd_t hnlprod)
genie::kPdgK0L
const int kPdgK0L
Definition: PDGCodes.h:176
genie::hnl::FluxContainer::acceptance
double acceptance
full acceptance == XYWgt * boostCorr^2 * accCorr
Definition: HNLFluxContainer.h:108
genie::hnl::FluxContainer::targetPointUser
TVector3 targetPointUser
point in detector HNL is forced towards in USER coords [m]
Definition: HNLFluxContainer.h:88
genie::hnl::FluxCreator::decay_pdpz
double decay_pdpz
final parent momentum [GeV]
Definition: HNLFluxCreator.h:244
genie::hnl::FluxContainer::targetPoint
TVector3 targetPoint
point in detector HNL is forced towards in NEAR coords [m]
Definition: HNLFluxContainer.h:86
genie::hnl::FluxContainer::nimpwt
double nimpwt
Weight of parent.
Definition: HNLFluxContainer.h:110
genie::hnl::FluxCreator::FillNonsense
void FillNonsense(int iEntry, genie::hnl::FluxContainer &gnmf) const
Definition: HNLFluxCreator.cxx:703
genie::hnl::FluxContainer::accCorr
double accCorr
acceptance correction (collimation effect. SM v == 1)
Definition: HNLFluxContainer.h:104
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::hnl::FluxContainer::nuPdg
int nuPdg
PDG code of SM neutrino that would have been produced.
Definition: HNLFluxContainer.h:80
genie::hnl::FluxContainer::XYWgt
double XYWgt
geometric acceptance (angular size of detector in parent rest frame)
Definition: HNLFluxContainer.h:101
genie::RandomGen::RndGen
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
genie::hnl::FluxContainer::nuProdChan
int nuProdChan
Decay mode that would have produced SM neutrino.
Definition: HNLFluxContainer.h:83
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::hnl::FluxContainer::boostCorr
double boostCorr
boost correction wrt parent rest-frame (ELAB = ECM * boostCorr)
Definition: HNLFluxContainer.h:102
genie::kPdgMuon
const int kPdgMuon
Definition: PDGCodes.h:37
genie::units::kSpeedOfLight
static constexpr double kSpeedOfLight
Definition: Units.h:32
genie::hnl::FluxCreator::HNLEnergy
TLorentzVector HNLEnergy(genie::hnl::HNLProd_t hnldm, TLorentzVector p4par) const
Definition: HNLFluxCreator.cxx:1148
genie::units::m
static constexpr double m
Definition: Units.h:71
genie::utils::print::Vec3AsString
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::OpenFluxInput ( std::string  finpath) const
private

Definition at line 741 of file HNLFluxCreator.cxx.

References cmeta, ctree, dir, fCurrPath, fFirstEntry, fNEntries, fPathLoaded, iCurrEntry, LOG, pDEBUG, and pFATAL.

Referenced by GetNFluxEntries(), and MakeTupleFluxEntry().

742 {
743  //if( std::strcmp( finpath.c_str(), fCurrPath.c_str() ) == 0 ) return;
744 
745  iCurrEntry = fFirstEntry;
746  fCurrPath = finpath;
747  finpath.append("/");
748 
749  LOG( "HNL", pDEBUG )
750  << "Getting flux input from finpath = " << finpath.c_str();
751 
752  // recurse over files in this directory and add to chain
753  if(!ctree){
754  ctree = new TChain( "dkTree" );
755  cmeta = new TChain( "dkMeta" );
756  }
757 
758  if( fPathLoaded ) return;
759 
760  TSystemDirectory dir( finpath.c_str(), finpath.c_str() );
761  TList * files = dir.GetListOfFiles(); int nFiles = 0;
762  assert( files );
763  files->Sort();
764 
765  TSystemFile * file;
766  TString fname;
767  TIter next(files);
768 
769  while( (file=( TSystemFile * ) next()) && !fPathLoaded ){
770  fname = file->GetName();
771  if( !file->IsDirectory() ){
772  TString fullpath = TString( finpath.c_str() ) + fname;
773  nFiles++;
774  ctree->Add( fullpath );
775  cmeta->Add( fullpath );
776  }
777  }
778 
779  if( !ctree ){ LOG( "HNL", pFATAL ) << "Could not open flux tree!"; }
780  if( !cmeta ){ LOG( "HNL", pFATAL ) << "Could not open meta tree!"; }
781  assert( ctree && cmeta );
782 
783  const int nEntriesInMeta = cmeta->GetEntries();
784  int nEntries = ctree->GetEntries();
785 
786  fNEntries = nEntries;
787 
788  LOG( "HNL", pDEBUG )
789  << "\nThere were " << nEntriesInMeta << " entries in meta with " << nEntries << " total nus"
790  << "\n got from " << nFiles << " files";
791 
792  fPathLoaded = true;
793 
794  delete file;
795  delete files;
796 }
pFATAL
#define pFATAL
Definition: Messenger.h:56
dir
string dir
Definition: INukeNucleonCorr.cxx:61
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
TVector3 FluxCreator::PointToRandomPointInBBox ( ) const
private

Definition at line 1290 of file HNLFluxCreator.cxx.

References ApplyUserRotation(), CheckGeomPoint(), genie::units::cm, fCx, fCy, fCz, fDetOffset, fDetRotation, fDoingOldFluxCalc, fLx, fLy, fLz, fTargetPoint, genie::RandomGen::Instance(), LOG, genie::units::m, pDEBUG, genie::RandomGen::RndGen(), and genie::utils::print::Vec3AsString().

Referenced by MakeTupleFluxEntry().

1291 {
1292  RandomGen * rnd = RandomGen::Instance();
1293  /*
1294  double ox = fCx + fDetOffset.at(0), oy = fCy + fDetOffset.at(1), oz = fCz + fDetOffset.at(2); // NEAR, m
1295  */
1296  double ox = fCx, oy = fCy, oz = fCz; // NEAR, m
1297 
1298  double rx = (rnd->RndGen()).Uniform( -fLx/2.0, fLx/2.0 ),
1299  ry = (rnd->RndGen()).Uniform( -fLy/2.0, fLy/2.0 ),
1300  rz = (rnd->RndGen()).Uniform( -fLz/2.0, fLz/2.0 ); // USER, m
1301 
1302  double ux = (rx + fDetOffset.at(0)) * units::m / units::cm;
1303  double uy = (ry + fDetOffset.at(1)) * units::m / units::cm;
1304  double uz = (rz + fDetOffset.at(2)) * units::m / units::cm;
1305  TVector3 checkPoint( ux, uy, uz ); // USER, cm
1306 
1307  TVector3 originPoint( -ox, -oy, -oz );
1308  TVector3 dumori( 0.0, 0.0, 0.0 );
1309  if( !fDoingOldFluxCalc ){
1310  // user-coordinates of this point. [cm]
1311  //double ux = rx - ox, uy = ry - oy, uz = rz - oz;
1312 
1313  LOG( "HNL", pDEBUG )
1314  << "\nChecking point " << utils::print::Vec3AsString(&checkPoint) << " [m, user]";
1315 
1316  // check if the point is inside the geometry, otherwise do it again
1317  std::string pathString = this->CheckGeomPoint( ux, uy, uz ); int iNode = 1; // 1 past beginning
1318  int iBad = 0;
1319  while( pathString.find( "/", iNode ) == string::npos && iBad < 10 ){
1320  rx = (rnd->RndGen()).Uniform( -fLx/2.0, fLx/2.0 ); ux = (rx + fDetOffset.at(0)) * units::m / units::cm;
1321  ry = (rnd->RndGen()).Uniform( -fLy/2.0, fLy/2.0 ); uy = (ry + fDetOffset.at(1)) * units::m / units::cm;
1322  rz = (rnd->RndGen()).Uniform( -fLz/2.0, fLz/2.0 ); uz = (rz + fDetOffset.at(2)) * units::m / units::cm;
1323  checkPoint.SetXYZ( ux, uy, uz );
1324  pathString = this->CheckGeomPoint( ux, uy, uz ); iNode = 1;
1325  iBad++;
1326  }
1327  assert( pathString.find( "/", iNode ) != string::npos );
1328  }
1329 
1330  // turn u back into [m] from [cm]
1331  ux *= units::cm / units::m; uy *= units::cm / units::m; uz *= units::cm / units::m;
1332  // return the absolute point in space [NEAR, m] that we're pointing to!
1333  checkPoint.SetXYZ( ux, uy, uz );
1334  checkPoint = this->ApplyUserRotation( checkPoint, dumori, fDetRotation, true ); // det --> tgt-hall
1335  checkPoint.SetXYZ( checkPoint.X() + ox, checkPoint.Y() + oy, checkPoint.Z() + oz );
1336 
1337  TVector3 vec( ux, uy, uz ); // USER m
1338  LOG( "HNL", pDEBUG )
1339  << "\nPointing to this point in BBox (USER coords): " << utils::print::Vec3AsString( &vec ) << "[m]"
1340  << "\nIn NEAR coords this is " << utils::print::Vec3AsString( &checkPoint ) << "[m]";
1341 
1342  // update bookkeeping
1343  fTargetPoint = checkPoint;
1344 
1345  return checkPoint;
1346 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::hnl::FluxCreator::fDetOffset
std::vector< double > fDetOffset
Definition: HNLFluxCreator.h:214
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::hnl::FluxCreator::fDetRotation
std::vector< double > fDetRotation
Definition: HNLFluxCreator.h:213
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::FluxCreator::CheckGeomPoint
std::string CheckGeomPoint(Double_t x, Double_t y, Double_t z) const
Definition: HNLFluxCreator.cxx:2136
genie::hnl::FluxCreator::ApplyUserRotation
TVector3 ApplyUserRotation(TVector3 vec, bool doBackwards=false) const
Definition: HNLFluxCreator.cxx:1959
genie::RandomGen::RndGen
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
genie::units::m
static constexpr double m
Definition: Units.h:71
genie::utils::print::Vec3AsString
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::ProcessEventRecord ( GHepRecord *  event_rec) const
virtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 41 of file HNLFluxCreator.cxx.

References genie::hnl::FluxContainer::acceptance, genie::GHepRecord::AddParticle(), genie::units::cm, genie::hnl::FluxContainer::delay, fCurrPath, fFirstEntry, fGeomFile, fGnmf, fIsUsingRootGeom, fLPx, fLPy, fLPz, fUsingDk2nu, genie::GHepParticle::GetX4(), iCurrEntry, ImportBoundingBox(), genie::kIStInitialState, genie::kPdgHNL, genie::hnl::FluxContainer::lepPdg, LOG, genie::units::m, MakeTupleFluxEntry(), genie::hnl::FluxContainer::nimpwt, genie::hnl::FluxContainer::p4User, genie::hnl::FluxContainer::parPdg, genie::GHepRecord::Particle(), pDEBUG, genie::hnl::FluxContainer::pdg, pFATAL, POTScaleWeight, SetCurrentEntry(), genie::GHepParticle::SetPosition(), SetUsingRootGeom(), genie::GHepRecord::SetVertex(), genie::GHepRecord::SetWeight(), genie::GHepRecord::SetXSec(), genie::hnl::FluxContainer::startPoint, genie::GHepRecord::Weight(), and genie::GHepRecord::XSec().

Referenced by main().

42 {
43  // Adds the inital state HNL at the event record.
44  // Also assigns the production vertex to evrec (this will be overwritten by subsequent modules)
45  // Also adds (acceptance*nimpwt)^(-1) component of weight
46 
47  this->SetCurrentEntry( evrec->XSec() );
48 
49  if( !fIsUsingRootGeom ){
50  this->SetUsingRootGeom(true); // must always be true
51 
52  gGeoManager = TGeoManager::Import( fGeomFile.c_str() );
53 
54  TGeoVolume * top_volume = gGeoManager->GetTopVolume();
55  assert( top_volume );
56  TGeoShape * ts = top_volume->GetShape();
57  TGeoBBox * box = (TGeoBBox *) ts;
58 
59  this->ImportBoundingBox( box );
60  }
61 
62  if( fUsingDk2nu ){
63 
64  if( iCurrEntry < fFirstEntry ) iCurrEntry = fFirstEntry;
65 
66  if( iCurrEntry >= fFirstEntry ) {
67 
68  if( iCurrEntry == fFirstEntry ){
69  FluxContainer * pfGnmf = new FluxContainer();
70  fGnmf = *pfGnmf;
71  delete pfGnmf;
72  }
73 
74  fGnmf = this->MakeTupleFluxEntry( iCurrEntry, fCurrPath );
75 
76  if( std::abs(fGnmf.pdg) == genie::kPdgHNL ){ // only add particle if parent is valid
77 
78  LOG( "HNL", pDEBUG ) << fGnmf;
79 
80  double invAccWeight = fGnmf.nimpwt * fGnmf.acceptance;
81  evrec->SetWeight( evrec->Weight() / invAccWeight );
82 
83  // scale by how many POT it takes to make the appropriate parent
84  /*
85  * To incorporate populations of parents, we take the cumulative multiplicity
86  * i.e. HNL light enough to be made by every parent get scaled by
87  * n1 = \sigma(p + target) / \sigma(p + target ; parent-producing)
88  * For HNL that are heavier than a muon, we don't take muons into account. So
89  * we up the scaling to incorporate their dropping out as
90  * n2 = \sigma(p + target) / \sigma(p + target ; parent-producing ; no muon) - etc.
91  */
92  evrec->SetWeight( evrec->Weight() * POTScaleWeight );
93 
94  // set prod-vertex in cm, ns, NEAR coords
95  TVector3 xVtxNear = fGnmf.startPoint; // in m, NEAR
96  double tVtx = fGnmf.delay; // ns
97  TLorentzVector pVtx( xVtxNear.X() * units::m / units::cm,
98  xVtxNear.Y() * units::m / units::cm,
99  xVtxNear.Z() * units::m / units::cm, tVtx ); // cm ns NEAR
100  evrec->SetVertex( pVtx ); // HNL production vertex. NOT where HNL decays to visible FS.
101 
102  // construct Particle(0). Don't worry about daughter links at this stage.
103  // this must be in USER coords
104  TLorentzVector probeP4 = fGnmf.p4User; // USER
105  GHepParticle ptHNL( fGnmf.pdg, kIStInitialState, -1, -1, -1, -1, probeP4, pVtx );
106  evrec->AddParticle( ptHNL );
107  }
108 
109  if( fGnmf.acceptance >= 0.0 ){
110  // set some event information where subsequent events can see it.
111  // Use the x4 position of the HNL. First, ensure the Vertex() is correctly set.
112  TLorentzVector * vx4 = evrec->Particle(0)->GetX4();
113  evrec->SetVertex( *vx4 );
114  TLorentzVector tmpx4( fLPx, fLPy, fGnmf.parPdg, fGnmf.lepPdg );
115  if( fLPz >= 0.0 ) evrec->SetXSec( 1.0 );
116  else evrec->SetXSec( -1.0 );
117  evrec->Particle(0)->SetPosition( tmpx4 );
118 
119  } // if( fGnmf.fgXYWgt >= 0 )
120  } // if( iCurrEntry > fFirstEntry )
121  } else {
122  LOG( "HNL", pFATAL )
123  << "No input dk2nu flux detected. Cannot proceed.";
124  exit(1);
125  }
126 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
genie::hnl::FluxContainer::delay
double delay
delay HNL would have wrt SMv [ns]
Definition: HNLFluxContainer.h:89
genie::hnl::FluxCreator::ImportBoundingBox
void ImportBoundingBox(TGeoBBox *box) const
Definition: HNLFluxCreator.cxx:2121
genie::hnl::FluxContainer::startPoint
TVector3 startPoint
parent decay vertex in NEAR coords [m]
Definition: HNLFluxContainer.h:85
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::hnl::FluxContainer::lepPdg
int lepPdg
PDG code of lepton produced with HNL on parent decay.
Definition: HNLFluxContainer.h:79
genie::hnl::FluxContainer::parPdg
int parPdg
parent PDG code
Definition: HNLFluxContainer.h:78
genie::hnl::FluxContainer
A GENIE flux container specific for HNL containers. Based on the dk2nu flux paradigm and genie::flux:...
Definition: HNLFluxContainer.h:62
genie::hnl::FluxCreator::fGnmf
genie::hnl::FluxContainer fGnmf
Definition: HNLFluxCreator.h:287
genie::hnl::FluxCreator::SetUsingRootGeom
void SetUsingRootGeom(bool IsUsingRootGeom) const
Definition: HNLFluxCreator.cxx:152
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::FluxContainer::p4User
TLorentzVector p4User
HNL momentum in USER coords [GeV/c].
Definition: HNLFluxContainer.h:95
genie::kPdgHNL
const int kPdgHNL
Definition: PDGCodes.h:224
genie::hnl::FluxContainer::acceptance
double acceptance
full acceptance == XYWgt * boostCorr^2 * accCorr
Definition: HNLFluxContainer.h:108
genie::hnl::FluxContainer::nimpwt
double nimpwt
Weight of parent.
Definition: HNLFluxContainer.h:110
genie::hnl::FluxCreator::MakeTupleFluxEntry
FluxContainer MakeTupleFluxEntry(int iEntry, std::string finpath) const
Definition: HNLFluxCreator.cxx:157
genie::units::m
static constexpr double m
Definition: Units.h:71
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::hnl::FluxCreator::SetCurrentEntry
void SetCurrentEntry(int iCurr) const
Definition: HNLFluxCreator.cxx:142
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::ReadBRs ( ) const
private

Definition at line 994 of file HNLFluxCreator.cxx.

References BR_K03e, BR_K03mu, BR_K2e, BR_K2mu, BR_K3e, BR_K3mu, BR_pi2e, BR_pi2mu, genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), genie::kPdgK0L, genie::kPdgKP, and genie::kPdgPiP.

Referenced by GetProductionProbs().

995 {
996  TParticlePDG * pionParticle = PDGLibrary::Instance()->Find( kPdgPiP );
997  TParticlePDG * kaonParticle = PDGLibrary::Instance()->Find( kPdgKP );
998  TParticlePDG * neukParticle = PDGLibrary::Instance()->Find( kPdgK0L );
999 
1000  TObjArray * pionDecayList = pionParticle->DecayList();
1001  TObjArray * kaonDecayList = kaonParticle->DecayList();
1002  TObjArray * neukDecayList = neukParticle->DecayList();
1003 
1004  TDecayChannel * pion2muChannel = ( TDecayChannel * ) pionDecayList->At(0);
1005  TDecayChannel * pion2elChannel = ( TDecayChannel * ) pionDecayList->At(1);
1006 
1007  TDecayChannel * kaon2muChannel = ( TDecayChannel * ) kaonDecayList->At(0);
1008  //TDecayChannel * kaon2elChannel = 0; // tiny BR, not in genie_pdg_table.txt
1009  TDecayChannel * kaon3muChannel = ( TDecayChannel * ) kaonDecayList->At(5);
1010  TDecayChannel * kaon3elChannel = ( TDecayChannel * ) kaonDecayList->At(4);
1011 
1012  TDecayChannel * neuk3muChannel = ( TDecayChannel * ) neukDecayList->At(4);
1013  TDecayChannel * neuk3elChannel = ( TDecayChannel * ) neukDecayList->At(2);
1014 
1015  BR_pi2mu = pion2muChannel->BranchingRatio();
1016  BR_pi2e = pion2elChannel->BranchingRatio();
1017 
1018  BR_K2mu = kaon2muChannel->BranchingRatio();
1019  BR_K2e = 1.6e-5; // From PDG 2021
1020  BR_K3mu = kaon3muChannel->BranchingRatio();
1021  BR_K3e = kaon3elChannel->BranchingRatio();
1022 
1023  BR_K03mu = 2.0 * neuk3muChannel->BranchingRatio(); // one from K0L-->mu+ and one from -->mu-
1024  BR_K03e = 2.0 * neuk3elChannel->BranchingRatio();
1025 }
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgK0L
const int kPdgK0L
Definition: PDGCodes.h:176
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
FluxContainer FluxCreator::RetrieveFluxInfo ( ) const
virtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 147 of file HNLFluxCreator.cxx.

References fGnmf.

Referenced by main().

148 {
149  return fGnmf;
150 }
genie::hnl::FluxCreator::fGnmf
genie::hnl::FluxContainer fGnmf
Definition: HNLFluxCreator.h:287
void FluxCreator::SetCurrentEntry ( int  iCurr) const
private

Definition at line 142 of file HNLFluxCreator.cxx.

References iCurrEntry.

Referenced by ProcessEventRecord().

143 {
144  iCurrEntry = iCurr;
145 }
void FluxCreator::SetFirstFluxEntry ( int  iFirst) const
virtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 2116 of file HNLFluxCreator.cxx.

References fFirstEntry.

Referenced by main().

2117 {
2118  fFirstEntry = iFirst;
2119 }
void FluxCreator::SetGeomFile ( string  geomfile) const

Definition at line 2110 of file HNLFluxCreator.cxx.

References fGeomFile, LOG, and pDEBUG.

Referenced by main().

2111 {
2112  LOG( "HNL", pDEBUG ) << "Setting geometry file to " << geomfile;
2113  fGeomFile = geomfile;
2114 }
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::SetInputFluxPath ( std::string  finpath) const
virtual

Implements genie::hnl::FluxRecordVisitorI.

Definition at line 128 of file HNLFluxCreator.cxx.

References fCurrPath, LOG, and pDEBUG.

Referenced by main().

129 {
130  LOG( "HNL", pDEBUG ) << "Setting input path to " << finpath;
131  fCurrPath = finpath;
132 }
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void FluxCreator::SetUsingRootGeom ( bool  IsUsingRootGeom) const
private

Definition at line 152 of file HNLFluxCreator.cxx.

References fIsUsingRootGeom.

Referenced by ProcessEventRecord().

153 {
154  fIsUsingRootGeom = IsUsingRootGeom;
155 }

Member Data Documentation

int genie::hnl::FluxCreator::anArSize
mutableprivate

Definition at line 248 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

char genie::hnl::FluxCreator::ancestor_imat[maxArray *maxC]
mutableprivate

Definition at line 265 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

char genie::hnl::FluxCreator::ancestor_ivol[maxArray *maxC]
mutableprivate

Definition at line 265 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::ancestor_nucleus[maxArray]
mutableprivate

Definition at line 264 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::ancestor_pdg[maxArray]
mutableprivate

Definition at line 258 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_polx[maxArray]
mutableprivate

Definition at line 262 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_poly[maxArray]
mutableprivate

Definition at line 262 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_polz[maxArray]
mutableprivate

Definition at line 262 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_pprodpx[maxArray]
mutableprivate

Definition at line 263 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_pprodpy[maxArray]
mutableprivate

Definition at line 263 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_pprodpz[maxArray]
mutableprivate

Definition at line 263 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

char genie::hnl::FluxCreator::ancestor_proc[maxArray *maxC]
mutableprivate

Definition at line 265 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_startpx[maxArray]
mutableprivate

Definition at line 260 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_startpy[maxArray]
mutableprivate

Definition at line 260 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_startpz[maxArray]
mutableprivate

Definition at line 260 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_startt[maxArray]
mutableprivate

Definition at line 259 of file HNLFluxCreator.h.

double genie::hnl::FluxCreator::ancestor_startx[maxArray]
mutableprivate

Definition at line 259 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_starty[maxArray]
mutableprivate

Definition at line 259 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_startz[maxArray]
mutableprivate

Definition at line 259 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_stoppx[maxArray]
mutableprivate

Definition at line 261 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_stoppy[maxArray]
mutableprivate

Definition at line 261 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ancestor_stoppz[maxArray]
mutableprivate

Definition at line 261 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::arSize
mutableprivate

Definition at line 248 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::beam0x
mutableprivate

Definition at line 281 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beam0y
mutableprivate

Definition at line 281 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beam0z
mutableprivate

Definition at line 281 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beamdxdz
mutableprivate

Definition at line 283 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beamdydz
mutableprivate

Definition at line 283 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beamhwidth
mutableprivate

Definition at line 282 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

char genie::hnl::FluxCreator::beamsim[maxC]
mutableprivate

Definition at line 279 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::beamvwidth
mutableprivate

Definition at line 282 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::BR_K03e
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_K03mu
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_K2e
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_K2mu
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_K3e
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_K3mu
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_pi2e
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

double genie::hnl::FluxCreator::BR_pi2mu
mutableprivate

Definition at line 194 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and ReadBRs().

TChain * genie::hnl::FluxCreator::cmeta = 0
mutableprivate

Definition at line 201 of file HNLFluxCreator.h.

Referenced by InitialiseMeta(), and OpenFluxInput().

TChain* genie::hnl::FluxCreator::ctree = 0
mutableprivate

Definition at line 201 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), MakeTupleFluxEntry(), and OpenFluxInput().

double genie::hnl::FluxCreator::decay_mupare
mutableprivate

Definition at line 254 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_muparpx
mutableprivate

Definition at line 254 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_muparpy
mutableprivate

Definition at line 254 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_muparpz
mutableprivate

Definition at line 254 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::decay_ndecay
mutableprivate

Definition at line 251 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_necm
mutableprivate

SM v CM energy [GeV].

Definition at line 245 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_nimpwt
mutableprivate

Importance weight from beamsim.

Definition at line 246 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

int genie::hnl::FluxCreator::decay_norig
mutableprivate

Definition at line 251 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::decay_ntype
mutableprivate

Definition at line 251 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_pdpx
mutableprivate

Definition at line 244 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_pdpy
mutableprivate

Definition at line 244 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_pdpz
mutableprivate

final parent momentum [GeV]

Definition at line 244 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_ppdxdz
mutableprivate

Definition at line 252 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_ppdydz
mutableprivate

Definition at line 252 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_ppenergy
mutableprivate

Definition at line 252 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::decay_ppmedium
mutableprivate

Definition at line 253 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::decay_pppz
mutableprivate

Definition at line 252 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::decay_ptype
mutableprivate

PDG code of parent.

Definition at line 242 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_vx
mutableprivate

Definition at line 243 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_vy
mutableprivate

Definition at line 243 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::decay_vz
mutableprivate

coordinates of prod vtx [cm]

Definition at line 243 of file HNLFluxCreator.h.

Referenced by InitialiseTree(), and MakeTupleFluxEntry().

int genie::hnl::FluxCreator::djob
mutableprivate

Definition at line 249 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

char genie::hnl::FluxCreator::dkvolcfg[maxC]
mutableprivate

Definition at line 280 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

bool genie::hnl::FluxCreator::doPol = true
mutableprivate

Definition at line 223 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), and LoadConfig().

std::map< genie::hnl::HNLProd_t, double > genie::hnl::FluxCreator::dynamicScores
mutableprivate

Definition at line 188 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

std::map< genie::hnl::HNLProd_t, double > genie::hnl::FluxCreator::dynamicScores_kaon
mutableprivate

Definition at line 190 of file HNLFluxCreator.h.

Referenced by GetProductionProbs().

std::map< genie::hnl::HNLProd_t, double > genie::hnl::FluxCreator::dynamicScores_muon
mutableprivate

Definition at line 191 of file HNLFluxCreator.h.

Referenced by GetProductionProbs().

std::map< genie::hnl::HNLProd_t, double > genie::hnl::FluxCreator::dynamicScores_neuk
mutableprivate

Definition at line 192 of file HNLFluxCreator.h.

Referenced by GetProductionProbs().

std::map< genie::hnl::HNLProd_t, double > genie::hnl::FluxCreator::dynamicScores_pion
mutableprivate

Definition at line 189 of file HNLFluxCreator.h.

Referenced by GetProductionProbs().

double genie::hnl::FluxCreator::fAx1
mutableprivate

Definition at line 216 of file HNLFluxCreator.h.

Referenced by ApplyUserRotation(), and LoadConfig().

double genie::hnl::FluxCreator::fAx2
mutableprivate

Definition at line 216 of file HNLFluxCreator.h.

Referenced by ApplyUserRotation(), and LoadConfig().

double genie::hnl::FluxCreator::fAz
mutableprivate

Definition at line 216 of file HNLFluxCreator.h.

Referenced by ApplyUserRotation(), and LoadConfig().

std::vector< double > genie::hnl::FluxCreator::fB2URotation
mutableprivate

Definition at line 212 of file HNLFluxCreator.h.

Referenced by GetB2URotation(), and LoadConfig().

std::vector< double > genie::hnl::FluxCreator::fB2UTranslation
mutableprivate

Definition at line 212 of file HNLFluxCreator.h.

Referenced by GetB2UTranslation(), and LoadConfig().

double genie::hnl::FluxCreator::fBx1
mutableprivate

Definition at line 217 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and LoadConfig().

double genie::hnl::FluxCreator::fBx2
mutableprivate

Definition at line 217 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and LoadConfig().

double genie::hnl::FluxCreator::fBz
mutableprivate

Definition at line 217 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and LoadConfig().

std::string genie::hnl::FluxCreator::fCurrPath = ""
mutableprivate

Definition at line 179 of file HNLFluxCreator.h.

Referenced by GetNFluxEntries(), OpenFluxInput(), ProcessEventRecord(), and SetInputFluxPath().

double genie::hnl::FluxCreator::fCx
mutableprivate

Definition at line 215 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fCy
mutableprivate

Definition at line 215 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fCz
mutableprivate

Definition at line 215 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

std::vector< double > genie::hnl::FluxCreator::fDetOffset
mutableprivate

Definition at line 214 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), GetDetOffset(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

std::vector< double > genie::hnl::FluxCreator::fDetRotation
mutableprivate

Definition at line 213 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), GetDetRotation(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

bool genie::hnl::FluxCreator::fDoingOldFluxCalc = false
mutableprivate

Definition at line 292 of file HNLFluxCreator.h.

Referenced by ImportBoundingBox(), LoadConfig(), MakeTupleFluxEntry(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fDx
mutableprivate

Definition at line 219 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fDy
mutableprivate

Definition at line 219 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fDz
mutableprivate

Definition at line 219 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fECM
mutableprivate

Definition at line 231 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

string genie::hnl::FluxCreator::fFinPath
mutableprivate

Definition at line 299 of file HNLFluxCreator.h.

int genie::hnl::FluxCreator::fFirstEntry = 0
mutableprivate

Definition at line 183 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry(), OpenFluxInput(), ProcessEventRecord(), and SetFirstFluxEntry().

std::vector< double > genie::hnl::FluxCreator::fFixedPolarisation
mutableprivate

Definition at line 226 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), LoadConfig(), and MakeTupleFluxEntry().

string genie::hnl::FluxCreator::fGeomFile = ""
mutableprivate

Definition at line 198 of file HNLFluxCreator.h.

Referenced by ProcessEventRecord(), and SetGeomFile().

genie::hnl::FluxContainer genie::hnl::FluxCreator::fGnmf
mutableprivate

Definition at line 287 of file HNLFluxCreator.h.

Referenced by ProcessEventRecord(), and RetrieveFluxInfo().

TH1D * genie::hnl::FluxCreator::fIntegrals = 0
mutableprivate

Definition at line 300 of file HNLFluxCreator.h.

bool genie::hnl::FluxCreator::fIsConfigLoaded = false
mutableprivate

Definition at line 296 of file HNLFluxCreator.h.

Referenced by LoadConfig().

bool genie::hnl::FluxCreator::fIsMajorana
mutableprivate

Definition at line 205 of file HNLFluxCreator.h.

Referenced by LoadConfig(), and MakeTupleFluxEntry().

bool genie::hnl::FluxCreator::fIsUsingRootGeom = false
mutableprivate

Definition at line 199 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry(), ProcessEventRecord(), and SetUsingRootGeom().

bool genie::hnl::FluxCreator::fixPol = false
mutableprivate

Definition at line 223 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), LoadConfig(), and MakeTupleFluxEntry().

int genie::hnl::FluxCreator::fLepPdg
mutableprivate

Definition at line 228 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fLPE
mutableprivate

Definition at line 225 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), and MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fLPx
mutableprivate

Definition at line 224 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), MakeTupleFluxEntry(), and ProcessEventRecord().

double genie::hnl::FluxCreator::fLPy
mutableprivate

Definition at line 224 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), MakeTupleFluxEntry(), and ProcessEventRecord().

double genie::hnl::FluxCreator::fLPz
mutableprivate

Definition at line 224 of file HNLFluxCreator.h.

Referenced by HNLEnergy(), MakeTupleFluxEntry(), and ProcessEventRecord().

double genie::hnl::FluxCreator::fLx
mutableprivate

Definition at line 209 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), ImportBoundingBox(), MakeBBox(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fLxR
mutableprivate

Definition at line 210 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and ImportBoundingBox().

double genie::hnl::FluxCreator::fLy
mutableprivate

Definition at line 209 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), ImportBoundingBox(), MakeBBox(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fLyR
mutableprivate

Definition at line 210 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and ImportBoundingBox().

double genie::hnl::FluxCreator::fLz
mutableprivate

Definition at line 209 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), ImportBoundingBox(), MakeBBox(), and PointToRandomPointInBBox().

double genie::hnl::FluxCreator::fLzR
mutableprivate

Definition at line 210 of file HNLFluxCreator.h.

Referenced by GetAngDeviation(), and ImportBoundingBox().

double genie::hnl::FluxCreator::fMass
mutableprivate

Definition at line 203 of file HNLFluxCreator.h.

Referenced by LoadConfig().

int genie::hnl::FluxCreator::fNEntries = 0
mutableprivate

Definition at line 185 of file HNLFluxCreator.h.

Referenced by GetNFluxEntries(), and OpenFluxInput().

int genie::hnl::FluxCreator::fNuPdg
mutableprivate

Definition at line 229 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

int genie::hnl::FluxCreator::fNuProdChan
mutableprivate

Definition at line 234 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

bool genie::hnl::FluxCreator::fPathLoaded = false
mutableprivate

Definition at line 179 of file HNLFluxCreator.h.

Referenced by OpenFluxInput().

int genie::hnl::FluxCreator::fProdChan
mutableprivate

Definition at line 234 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

string genie::hnl::FluxCreator::fProdHist
mutableprivate

Definition at line 299 of file HNLFluxCreator.h.

double genie::hnl::FluxCreator::fRadius
mutableprivate

Definition at line 294 of file HNLFluxCreator.h.

Referenced by ImportBoundingBox(), LoadConfig(), and MakeBBox().

bool genie::hnl::FluxCreator::fRerollPoints = false
mutableprivate

Definition at line 293 of file HNLFluxCreator.h.

Referenced by LoadConfig(), and MakeTupleFluxEntry().

std::vector<double> genie::hnl::FluxCreator::fScales
mutableprivate

Definition at line 290 of file HNLFluxCreator.h.

Referenced by LoadConfig().

double genie::hnl::FluxCreator::fSMECM
mutableprivate

Definition at line 231 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

TH1D* genie::hnl::FluxCreator::fSpectrum = 0
mutableprivate

Definition at line 300 of file HNLFluxCreator.h.

bool genie::hnl::FluxCreator::fSupplyingBEAM = false
mutableprivate

Definition at line 295 of file HNLFluxCreator.h.

Referenced by LoadConfig(), and MakeTupleFluxEntry().

TVector3 genie::hnl::FluxCreator::fTargetPoint
mutableprivate

Definition at line 236 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry(), and PointToRandomPointInBBox().

TVector3 genie::hnl::FluxCreator::fTargetPointUser
mutableprivate

Definition at line 236 of file HNLFluxCreator.h.

TGeoVolume* genie::hnl::FluxCreator::fTopVol = 0
mutableprivate

Definition at line 197 of file HNLFluxCreator.h.

std::vector< double > genie::hnl::FluxCreator::fU4l2s
mutableprivate

Definition at line 204 of file HNLFluxCreator.h.

Referenced by GetProductionProbs(), and LoadConfig().

bool genie::hnl::FluxCreator::fUsingDk2nu = true
mutableprivate

Definition at line 298 of file HNLFluxCreator.h.

Referenced by ProcessEventRecord().

double genie::hnl::FluxCreator::fZm
mutableprivate

Definition at line 232 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

double genie::hnl::FluxCreator::fZp
mutableprivate

Definition at line 232 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

char genie::hnl::FluxCreator::horncfg[maxC]
mutableprivate

Definition at line 280 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

int genie::hnl::FluxCreator::iCurrEntry = 0
mutableprivate

Definition at line 181 of file HNLFluxCreator.h.

Referenced by OpenFluxInput(), ProcessEventRecord(), and SetCurrentEntry().

bool genie::hnl::FluxCreator::isParentOnAxis = true
mutableprivate

Definition at line 196 of file HNLFluxCreator.h.

Referenced by LoadConfig(), and MakeTupleFluxEntry().

int genie::hnl::FluxCreator::job
mutableprivate

beamsim MC job number

Definition at line 275 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

char genie::hnl::FluxCreator::location_name[maxArray *maxC]
mutableprivate

Definition at line 285 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::location_x[maxArray]
mutableprivate

Definition at line 284 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::location_y[maxArray]
mutableprivate

Definition at line 284 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::location_z[maxArray]
mutableprivate

Definition at line 284 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

int genie::hnl::FluxCreator::mArSize
mutableprivate

Definition at line 278 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

const int genie::hnl::FluxCreator::maxArray = 30
staticprivate

Definition at line 238 of file HNLFluxCreator.h.

Referenced by InitialiseMeta(), and InitialiseTree().

const int genie::hnl::FluxCreator::maxC = 100
staticprivate

Definition at line 238 of file HNLFluxCreator.h.

Referenced by InitialiseMeta(), and InitialiseTree().

double genie::hnl::FluxCreator::nuray_E[maxArray]
mutableprivate

Definition at line 256 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::nuray_px[maxArray]
mutableprivate

Definition at line 256 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::nuray_py[maxArray]
mutableprivate

Definition at line 256 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::nuray_pz[maxArray]
mutableprivate

Definition at line 256 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::nuray_wgt[maxArray]
mutableprivate

Definition at line 256 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::parentEnergy
mutableprivate

Definition at line 230 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

double genie::hnl::FluxCreator::parentMass
mutableprivate

Definition at line 230 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

double genie::hnl::FluxCreator::parentMomentum
mutableprivate

Definition at line 230 of file HNLFluxCreator.h.

Referenced by MakeTupleFluxEntry().

char genie::hnl::FluxCreator::physcuts[maxC]
mutableprivate

Definition at line 279 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

char genie::hnl::FluxCreator::physics[maxC]
mutableprivate

Definition at line 279 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::potnum
mutableprivate

N POT for this SM-v.

Definition at line 241 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::pots
mutableprivate

how many pot in this job?

Definition at line 276 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

double genie::hnl::FluxCreator::POTScaleWeight
mutableprivate

Definition at line 289 of file HNLFluxCreator.h.

Referenced by LoadConfig(), and ProcessEventRecord().

double genie::hnl::FluxCreator::ppvx
mutableprivate

Definition at line 250 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ppvy
mutableprivate

Definition at line 250 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::ppvz
mutableprivate

Definition at line 250 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

char genie::hnl::FluxCreator::tgtcfg[maxC]
mutableprivate

Definition at line 280 of file HNLFluxCreator.h.

Referenced by InitialiseMeta().

int genie::hnl::FluxCreator::tgtexit_tgen
mutableprivate

Definition at line 269 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::tgtexit_tptype
mutableprivate

Definition at line 269 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tpx
mutableprivate

Definition at line 268 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tpy
mutableprivate

Definition at line 268 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tpz
mutableprivate

Definition at line 268 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tvx
mutableprivate

Definition at line 267 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tvy
mutableprivate

Definition at line 267 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::tgtexit_tvz
mutableprivate

Definition at line 267 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trkpx[maxArray]
mutableprivate

Definition at line 272 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trkpy[maxArray]
mutableprivate

Definition at line 272 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trkpz[maxArray]
mutableprivate

Definition at line 272 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trkx[maxArray]
mutableprivate

Definition at line 271 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trky[maxArray]
mutableprivate

Definition at line 271 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

double genie::hnl::FluxCreator::traj_trkz[maxArray]
mutableprivate

Definition at line 271 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

int genie::hnl::FluxCreator::trArSize
mutableprivate

Definition at line 248 of file HNLFluxCreator.h.

Referenced by InitialiseTree().

The documentation for this class was generated from the following files:
Generated on Wed Jul 17 2024 14:47:06 for GENIEGenerator by   doxygen 1.8.5