JSIRENE Namespace Reference

Detector simulations. More...

Classes
struct	JFoxWolframMoments
	Class for computing Fox-Wolfram moments. More...
struct	JSphericityTensor
	Class for sphericity tensor calculations. More...
struct	JHitInertiaTensor
	Class for hit inertia tensor calculations. More...
class	JPulse
	Auxiliary class for a time-over-threshold pulse from a PMT. More...
struct	JPythia
	Auxiliary class to determine EM-equivalent energy as a function of PDG particle code and energy. More...
class	JSeaWater
	Sea water composition. More...
struct	JParameters
	Detector simulation parameters. More...
struct	JHit_t
	Auxiliary class to set-up Hit. More...
struct	JHits_t
	Auxiliary data structure for list of hits with hit merging capability. More...
struct	JTrk_t
	Auxiliary class to set-up Trk. More...
struct	JPoint
	Point along muon trajectory. More...
struct	JTrack
	Muon trajectory. More...
struct	JVertex
	Vertex of energy loss of muon. More...

Functions
double	getThrust (const std::vector< Trk >::const_iterator __begin, const std::vector< Trk >::const_iterator __end)
	Compute thrust for a given range of tracks. More...
double	getThrust (const Evt &event)
	Compute thrust for a given event. More...
bool	operator< (const JPulse &first, const JPulse &second)
	Compare Monte Carlo hit times. More...
bool	operator< (const JPulse &hit, const double t0)
	Compare Monte Carlo hit times. More...
JHitType_t	getHitType (const JPDFType_t pdf, const bool shower=false)
	Get hit type corresponding to given PDF type. More...
const JCylinder3D	getMaximumContainmentVolume ()
	Auxiliary function to retrieve the maximum cylindrical containment volume. More...
double	getVisibleEnergy (const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy of a track. More...
Vec	getVisibleEnergyVector (const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of a track. More...
double	getVisibleEnergy (std::vector< Trk >::const_iterator __begin, std::vector< Trk >::const_iterator __end, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy of a given range of tracks. More...
Vec	getVisibleEnergyVector (std::vector< Trk >::const_iterator __begin, std::vector< Trk >::const_iterator __end, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of a given range of tracks. More...
double	getVisibleEnergy (const Evt &evt, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of an event. More...
Vec	getVisibleEnergyVector (const Evt &evt, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of an event. More...
double	opa_weight_high_e (const double ekin)
double	ngamma_elec (const double ekin)
double	weight_pion (const double ekin)
double	weight_kaon (const double ekin)
double	weight_kshort (const double ekin)
double	weight_klong (const double ekin)
double	weight_proton (const double ekin)
double	weight_neutron (const double ekin)
double	opa_efrac (const int ipart, const double ekin)
double	pythia (const int type, const double E)
	Get equivalent EM-energy for given pion energy. More...

Variables
static const JPythia	pythia
	Function object for relative light yield as a function of GEANT particle code. More...
struct {
std::vector< size_t >   JSIRENE::ns
std::vector< double >   JSIRENE::vs
}	getNumberOfPhotoElectrons
	Auxiliary data structure for determination of number of photo-electrons. More...

Detailed Description

Detector simulations.

Author

mdejong

Function Documentation

getThrust() [1/2]

double JSIRENE::getThrust ( const std::vector< Trk >::const_iterator  __begin, const std::vector< Trk >::const_iterator  __end  )

inline

Compute thrust for a given range of tracks.

The definition was taken from the description in section 15.2.2 of arXiv:hep-ph/0603175v2.

Parameters

__begin	beginning of track data
__end	end of track data

Returns

thrust

< Thrust axis

< residual angle

Definition at line 43 of file JEventShapeVariables.hh.

  {
    using namespace std;
    using namespace JPP;
    
    static const int    MAX_STEPS = 20;
    static const double epsilon = 1e-9;
 
 
    Vec axis0(0.0, 0.0, 0.0); //!< Thrust axis
 
    // Set thrust axis start value (= summed final state particle momenta)
    
    for (vector<Trk>::const_iterator i = __begin; i != __end; ++i) {
      if (is_finalstate(*i)) {
        axis0 += getKineticEnergy(*i) * i->dir;
      }
    }
 
    if (axis0.len() < epsilon) { // If total final state momentum is zero, set start value equal to first final state track momentum
      
      vector<Trk>::const_iterator i = find_if(__begin, __end, is_finalstate);
      
      if (i != __end) {
        axis0 = getKineticEnergy(*i) * i->dir;
      } else {
        THROW(JValueOutOfRange, "getThrust(): No final state tracks given.");
      }
    }
    
    axis0.normalize();
 
    // Find axis which maximizes the thrust
 
    double res = numeric_limits<double>::max(); //!< residual angle
    
    for (int k = 0; k < MAX_STEPS && res > epsilon; ++k) {
 
      Vec axis1(0.0, 0.0, 0.0);
      
      for (vector<Trk>::const_iterator i = __begin; i != __end; ++i) {
        
        if (is_finalstate(*i)) {
 
          const Vec p = getKineticEnergy(*i) * i->dir;
          
          axis1 += (p.dot(axis0) > 0.0 ? p : -p);
        }
      }
 
      axis1.normalize();
 
      res = acos( axis0.dot(axis1) );
 
      axis0 = axis1;
    }
 
    // Compute thrust
 
    double pSum = 0.0;
    double norm = 0.0;
 
    for (vector<Trk>::const_iterator i = __begin; i != __end; ++i) {
 
      if (is_finalstate(*i)) {
 
        const Vec p = getKineticEnergy(*i) * i->dir;
        
        pSum += fabs(p.dot(axis0));
        norm += p.len();
      }
    }
 
    return pSum / norm;
  }

getThrust() [2/2]

double JSIRENE::getThrust ( const Evt &  event )

inline

Compute thrust for a given event.

The definition was taken from the description in section 15.2.2 of arXiv:hep-ph/0603175v2.

Parameters

event

event

Returns

thrust

Definition at line 128 of file JEventShapeVariables.hh.

  {
    return getThrust(event.mc_trks.begin(), event.mc_trks.end());
  }

operator<() [1/2]

bool JSIRENE::operator< ( const JPulse &  first, const JPulse &  second  )

inline

Compare Monte Carlo hit times.

Parameters

first	first hit
second	second hit

Returns

true if first hit earlier than second; else false

Definition at line 81 of file JPulse.hh.

  {
    return first.getLowerLimit() < second.getLowerLimit();
  }

operator<() [2/2]

bool JSIRENE::operator< ( const JPulse &  hit, const double  t0  )

inline

Compare Monte Carlo hit times.

Parameters

hit	hit
t0	time [ns]

Returns

true if hit earlier than given time; else false

Definition at line 94 of file JPulse.hh.

  {
    return hit.getLowerLimit() < t0;
  }

getHitType()

JHitType_t JSIRENE::getHitType ( const JPDFType_t  pdf, const bool  shower = false  )

inline

Get hit type corresponding to given PDF type.

Parameters

pdf	PDF type
shower	force origin from neutrino interaction shower

Returns

hit type

Definition at line 151 of file JSireneToolkit.hh.

  {
    using namespace JAANET;
    using namespace JPHYSICS;
    
    switch (pdf) {
 
    case DIRECT_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_DIRECT;
 
    case SCATTERED_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_SCATTERED;
 
    case DIRECT_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_DIRECT;
 
    case SCATTERED_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_SCATTERED;
 
    case DIRECT_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_DIRECT : HIT_TYPE_BREMSSTRAHLUNG_DIRECT );
 
    case SCATTERED_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_SCATTERED : HIT_TYPE_BREMSSTRAHLUNG_SCATTERED);
 
    default:
      return HIT_TYPE_UNKNOWN;      
    }    
  }

getMaximumContainmentVolume()

const JCylinder3D JSIRENE::getMaximumContainmentVolume ( )

inline

Auxiliary function to retrieve the maximum cylindrical containment volume.

Returns

maximum cylindrical containment volume

Definition at line 46 of file JVisibleEnergyToolkit.hh.

  {
    using namespace JPP;
 
    const double R_Earth = R_EARTH_KM * 1e3; // m
    
    const JVector2D center(0.0, 0.0);
    const JCircle2D circle(center, R_Earth);
    
    return JCylinder3D(circle, -R_Earth, R_Earth);
  }

getVisibleEnergy() [1/3]

double JSIRENE::getVisibleEnergy	(	const Trk &	track,
		const JCylinder3D &	can = getMaximumContainmentVolume()
	)

Get the visible energy of a track.

This method accounts for muon radiative energy losses.

Parameters

track	track
can	detector can

Returns

visible energy [GeV]

Definition at line 67 of file JVisibleEnergyToolkit.hh.

  {
    using namespace std;
    using namespace JPP;
 
    double Evis = 0.0;
    
    if (is_finalstate(track)) {
      
      if (is_muon(track)) {
 
        // Determine muon pathlength inside detector [m]
 
        const JCylinder3D::intersection_type& intersection = can.getIntersection(getAxis(track));
        
        const double Lmuon = gWater.getX(track.E, MASS_MUON / getSinThetaC());
        const double Leff  = (min(Lmuon, max(intersection.second, 0.0)) -
                              min(Lmuon, max(intersection.first,  0.0)));
                                  
 
        // Determine visible energy deposition [GeV]
 
        const double Emidpoint = gWater.getE(track.E, Lmuon/2.0);
      
        const double dEb = gWater.getEb(track.E, Leff);
        const double dEc = Leff / geanc();
        const double dEd = Leff * getDeltaRaysFromMuon(Emidpoint);
      
        Evis = dEb + dEc + dEd;
          
      } else if (!is_neutrino(track) && JPDB::getInstance().hasPDG(track.type)) {
 
        Evis = pythia(track.type, getKineticEnergy(track));
      }
    }
 
    return Evis;
  }

getVisibleEnergyVector() [1/3]

Vec JSIRENE::getVisibleEnergyVector ( const Trk &  track, const JCylinder3D &  can = getMaximumContainmentVolume()  )

inline

Get the visible energy vector of a track.

This method accounts for muon radiative energy losses.

Parameters

track	track
can	detector can

Returns

visible energy vector [GeV]

Definition at line 116 of file JVisibleEnergyToolkit.hh.

                                                                                            {
    return getVisibleEnergy(track, can) * track.dir;
  }

getVisibleEnergy() [2/3]

double JSIRENE::getVisibleEnergy ( std::vector< Trk >::const_iterator  __begin, std::vector< Trk >::const_iterator  __end, const JCylinder3D &  can = getMaximumContainmentVolume()  )

inline

Get the visible energy of a given range of tracks.

This method accounts for muon radiative energy losses.

Parameters

__begin	start of track data
__end	end of track data
can	detector can

Returns

visible energy [GeV]

Definition at line 131 of file JVisibleEnergyToolkit.hh.

  {
    using namespace std;
    
    double Evis = 0.0;
    
    for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
      Evis += getVisibleEnergy(*track, can);
    }
    
    return Evis;
  }

getVisibleEnergyVector() [2/3]

Vec JSIRENE::getVisibleEnergyVector ( std::vector< Trk >::const_iterator  __begin, std::vector< Trk >::const_iterator  __end, const JCylinder3D &  can = getMaximumContainmentVolume()  )

inline

Get the visible energy vector of a given range of tracks.

This method accounts for muon radiative energy losses.

Parameters

__begin	start of track data
__end	end of track data
can	detector can

Returns

visible energy vector [GeV]

Definition at line 156 of file JVisibleEnergyToolkit.hh.

  {
    using namespace std;
    
    Vec Evis(0.0, 0.0, 0.0);    
 
    for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
      Evis += getVisibleEnergyVector(*track, can);
    }
 
    return Evis;
  }

getVisibleEnergy() [3/3]

double JSIRENE::getVisibleEnergy ( const Evt &  evt, const JCylinder3D &  can = getMaximumContainmentVolume()  )

inline

Get the visible energy vector of an event.

This method accounts for muon radiative energy losses.

Parameters

evt	event
can	detector can

Returns

visible energy [GeV]

Definition at line 180 of file JVisibleEnergyToolkit.hh.

  {
    return getVisibleEnergy(evt.mc_trks.begin(), evt.mc_trks.end(), can);
  }

getVisibleEnergyVector() [3/3]

Vec JSIRENE::getVisibleEnergyVector ( const Evt &  evt, const JCylinder3D &  can = getMaximumContainmentVolume()  )

inline

Get the visible energy vector of an event.

This method accounts for muon radiative energy losses.

Parameters

evt	event
can	detector can

Returns

visible energy vector [GeV]

Definition at line 195 of file JVisibleEnergyToolkit.hh.

  {
    return getVisibleEnergyVector(evt.mc_trks.begin(), evt.mc_trks.end(), can);
  }

opa_weight_high_e()

double JSIRENE::opa_weight_high_e ( const double  ekin )

inline

Definition at line 20 of file pythia.hh.

  {
    static const double Ma    =  72.425;
    static const double Mb    = -49.417;
    static const double Mc    =   5.858;
    static const double Md    = 207.252;
    static const double Me    = 132.784;
    static const double Mf    = -10.277;
    static const double Mg    = -19.441;
    static const double Mh    =  58.598;
    static const double Mi    =  53.161;
    static const double Mkref =   2.698;
    //   leading_coef Mlc=(Ma-Mf)/Mkref
    static const double Mlc   = (Ma-Mf)/Mkref;
 
    double num, denom, lognp, Ee, logE;
      
    //           implement the energy-fractions according to M. Dentler: ANTARES-SOFT-2012-009. E must be the kinetic energy in GeV. Here, 'everything else' is treated as a pion. Other than charged pions, mostly kaons of all varieties and protons fall into this category. These fits were made mostly at high energies - at low energies, other fits are used.
         
    if (ekin > 0.2) 
      logE=log10(ekin);
    else 
      // nphotons constant below 0.2 GeV (lowest fit point)
      return 0.292;
 
    num   = Me + Md*logE + Mc*pow(logE,2) + Mb*pow(logE,3) + Ma *pow(logE,4) + Mlc*pow(logE,5);
    denom = Mi + Mh*logE + Mg*pow(logE,2) + Mf*pow(logE,3) + Mlc*pow(logE,4);
    
    if (denom <= 0) {
      return 0.;
    }
 
    lognp = num/denom;
    Ee    = pow(10.0, lognp-Mkref);
    
    return Ee/ekin;
  }

ngamma_elec()

double JSIRENE::ngamma_elec ( const double  ekin )

inline

Definition at line 58 of file pythia.hh.

  {
    static const double eleca =  1.33356e5;
    static const double elecb =  1.66113e2;
    static const double elecc = 16.4949;
    static const double elecd =  1.5385e5;
    static const double elece =  6.04871e5;
    
    return (eleca*ekin+elecb) * exp(-ekin/elecc) + (elecd*ekin+elece)* (1.-exp(-ekin/elecc));
  }

weight_pion()

double JSIRENE::weight_pion ( const double  ekin )

inline

Definition at line 69 of file pythia.hh.

  {
    static const double pia =  0.538346;
    static const double pib =  1.32041;
    static const double pic =  0.737415;
    static const double pid = -0.813861;
    static const double pie = -2.22444;
    static const double pif = -2.15795;
    static const double pig = -6.47242;
    static const double pih = -2.7567;
    static const double pix =  8.83498;
      
    if      (ekin < 0.06) 
      return 1e4*pia/ngamma_elec(ekin);
    else if (ekin < 0.15) 
      return pix*1e4*ekin/ngamma_elec(ekin);
    else {
      const double lek=log(ekin);
      return (pib*1e5*ekin + (pow(ekin,(1.-1./pic))) * (pid*1e4 + 1e4*pie*lek + 1e4*pif*pow(lek,2) + 1e3*pig*pow(lek,3) + 1e2*pih*pow(lek,4)))/ngamma_elec(ekin);
    }
  }

weight_kaon()

double JSIRENE::weight_kaon ( const double  ekin )

inline

Definition at line 92 of file pythia.hh.

  {
    static const double ka  = 12.7537;
    static const double kb  = -1.052;
    static const double kc  =  3.49559;
    static const double kd  = 16.7914;
    static const double ke  = -0.355066;
    static const double kf  =  2.77116;
 
    if (ekin > 0.26) 
      return (1e4*ka*(ekin+kb)*(1.-exp(-ekin/kc)) + 1e4*(kd*ekin+ke)*exp(-ekin/kc))/ngamma_elec(ekin);
    else
      return kf*1e4/ngamma_elec(ekin);
  }      

weight_kshort()

double JSIRENE::weight_kshort ( const double  ekin )

inline

Definition at line 108 of file pythia.hh.

  {
    static const double k0sa  = 0.351242;
    static const double k0sb  = 0.613076;
    static const double k0sc  = 6.24682;
    static const double k0sd  = 2.8858;
 
    return (k0sa*1e5 + k0sb*1e5*ekin + ekin*k0sc*1e4*log(k0sd + 1./ekin))/ngamma_elec(ekin);
  }

weight_klong()

double JSIRENE::weight_klong ( const double  ekin )

inline

Definition at line 119 of file pythia.hh.

  {
    static const double k0la  = 2.18152;
    static const double k0lb  = -0.632798;
    static const double k0lc  = 0.999514;
    static const double k0ld  = 1.36052;
    static const double k0le  = 4.22484;
    static const double k0lf  = -0.307859;
      
    if (ekin < 1.5)
      return (1e4*k0la + ekin*1e5*k0lb*log(ekin)  + 1e5*k0lc*pow(ekin,1.5))/ngamma_elec(ekin);
    else
      return (ekin*k0ld*1e5 + pow(ekin,1.-1./k0le)*k0lf*1e5)/ngamma_elec(ekin);
  }

weight_proton()

double JSIRENE::weight_proton ( const double  ekin )

inline

Definition at line 135 of file pythia.hh.

  {     
    static const double pa   =  12.1281;
    static const double pb   = -17.1528;
    static const double pc   =   0.573158;
    static const double pd   =  34.1436;
    static const double pe   =  -0.28944;
    static const double pf   = -13.2619;
    static const double pg   =  24.1357;
      
    return (1e4*(pa*ekin+pb)*(1.-exp(-ekin/pc)) + 1e4*(pd*ekin +pe+ pf*pow(ekin,2) + pg*pow(ekin,3))*exp(-ekin/pc)) / ngamma_elec(ekin);
  }

weight_neutron()

double JSIRENE::weight_neutron ( const double  ekin )

inline

Definition at line 149 of file pythia.hh.

  {
    static const double na   = 1.24605;
    static const double nb   = 0.63819;
    static const double nc   = -0.802822;
    static const double nd   = -0.935327;
    static const double ne   = -6.1126;
    static const double nf   = -1.96894;
    static const double ng   = 0.716954;
    static const double nh   = 2.68246;
    static const double ni   = -9.39464;
    static const double nj   = 15.0737;
      
    if (ekin > 0.5) {
      const double lek=log(ekin);
      return (na*1e5*ekin + 1e3*pow(ekin,1.-1./nb) * (100.*nc+100.*nd*lek + 10.*ne*pow(lek,2) + 11.*nf*pow(lek,3)))/ngamma_elec(ekin);
    } else
      return (1e3*ng + 1e4*nh*ekin + 1e4*ni*pow(ekin,2) + 1e4*nj*pow(ekin,3))/ngamma_elec(ekin);
  }

opa_efrac()

double JSIRENE::opa_efrac ( const int  ipart, const double  ekin  )

inline

Definition at line 173 of file pythia.hh.

  {
    using namespace JPP;
 
    double etemp, weight_part;
 
    // neutrino or muon
    if (ipart == GEANT4_TYPE_NEUTRINO       || 
        ipart == GEANT4_TYPE_ANTIMUON       || 
        ipart == GEANT4_TYPE_MUON           ||
        ipart == GEANT4_TYPE_ANTITAU        || 
        ipart == GEANT4_TYPE_TAU) {
      return 0.0;
    }
      
    // 100% energy to EM showers if electron, positron, pi0, or gamma
    if (ipart == GEANT4_TYPE_PHOTON         || 
        ipart == GEANT4_TYPE_ANTIELECTRON   || 
        ipart == GEANT4_TYPE_ELECTRON       ||
        ipart == GEANT4_TYPE_NEUTRAL_PION   ||
        ipart == GEANT4_TYPE_ETA) {
      return 1.0;
    }
 
    // low-energy fits went up to 40 GeV      
    if (ekin <= 40) 
      etemp=ekin;
    else
      etemp=40.;
      
    // otherwise, calculate fractions the long way. First we get the
    // expected number of photons from the particle:
    if      (ipart == GEANT4_TYPE_PION_PLUS ||
             ipart == GEANT4_TYPE_PION_MINUS)
      weight_part=weight_pion(etemp);
    else if (ipart == GEANT4_TYPE_KAON_LONG)
      weight_part=weight_klong(etemp);
    else if (ipart == GEANT4_TYPE_KAON_SHORT) 
      weight_part=weight_kshort(etemp);
    else if (ipart == GEANT4_TYPE_KAON_PLUS ||
             ipart == GEANT4_TYPE_KAON_MINUS)
      weight_part=weight_kaon(etemp);
    else if (ipart == GEANT4_TYPE_NEUTRON) 
      weight_part=weight_neutron(etemp);
    else if (ipart == GEANT4_TYPE_PROTON    ||
             ipart == GEANT4_TYPE_ANTIPROTON)
      weight_part=weight_proton(etemp);
    else
      //         when in doubt, treat it as a proton
      weight_part=weight_proton(etemp);
        
    if (ekin <= 40)
      return weight_part;
    else if (ekin >= 1e7) 
      return opa_weight_high_e(ekin);
    else {
      // smooth transition between low-energy weights and high-energy pion weight
      const double wp40  = weight_part;
      const double whe40 = opa_weight_high_e(40.);
      const double whex  = opa_weight_high_e(ekin);      
      return whex - (whe40-wp40)*(7.-log10(ekin))/5.398;
    }
  }

pythia()

double JSIRENE::pythia ( const int  type, const double  E  )

inline

Get equivalent EM-energy for given pion energy.

Parameters

type	particle type [PDG]
E	particle energy [GeV]

Returns

EM-equivalent energy [GeV]

Definition at line 245 of file pythia.hh.

  {
    int   ipart = JPDB::getInstance().getPDG(type).geant;
    float ekin  = (float) E;
 
    return opa_efrac(ipart, ekin) * E;
  }

Variable Documentation

pythia

const JPythia JSIRENE::pythia

static

Function object for relative light yield as a function of GEANT particle code.

Definition at line 96 of file JPythia.hh.

ns

std::vector<size_t> JSIRENE::ns

mutableprivate

Definition at line 138 of file JSireneToolkit.hh.

vs

std::vector<double> JSIRENE::vs

mutableprivate

Definition at line 139 of file JSireneToolkit.hh.

const { ... } JSIRENE::getNumberOfPhotoElectrons

Auxiliary data structure for determination of number of photo-electrons.