JSIRENE Namespace Reference

Detector simulations. More...

Classes
struct	JFoxWolframMoments
	Class for computing Fox-Wolfram moments. More...
struct	JHit_t
	Auxiliary class to set-up Hit. More...
struct	JHitInertiaTensor
	Class for hit inertia tensor calculations. More...
struct	JHits_t
	Auxiliary data structure for list of hits with hit merging capability. More...
struct	JParameters
	Detector simulation parameters. More...
struct	JPoint
	Point along muon trajectory. 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	JSphericityTensor
	Class for sphericity tensor calculations. More...
struct	JTrack
	Muon trajectory. More...
struct	JTrk_t
	Auxiliary class to set-up Trk. More...
struct	JVertex
	Vertex of energy loss of muon. More...
struct	number_of_photo_electrons_type
	Auxiliary data structure for determination of number of photo-electrons. More...

Functions
const JCylinder3D	getMaximumContainmentVolume ()
	Forward function declarations.
double	getVisibleEnergy (const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy of a track.
double	getVisibleEnergy (const Evt &evt, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of an event.
double	getThrust (const std::vector< Trk >::const_iterator __begin, const std::vector< Trk >::const_iterator __end)
	Compute thrust for a given range of tracks.
double	getThrust (const Evt &event)
	Compute thrust for a given event.
bool	operator< (const JPulse &first, const JPulse &second)
	Compare Monte Carlo hit times.
bool	operator< (const JPulse &hit, const double t0)
	Compare Monte Carlo hit times.
JHitType_t	getHitType (const JPDFType_t pdf, const bool shower=false)
	Get hit type corresponding to given PDF type.
Vec	getVisibleEnergyVector (const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of a track.
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.
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.
Vec	getVisibleEnergyVector (const Evt &evt, const JCylinder3D &can=getMaximumContainmentVolume())
	Get the visible energy vector of an event.
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.

Variables
static const JPythia	pythia
	Function object for relative light yield as a function of GEANT particle code.
const struct JSIRENE::number_of_photo_electrons_type	getNumberOfPhotoElectrons

Detailed Description

Detector simulations.

Author

bjung

mdejong

Function Documentation

getMaximumContainmentVolume()

const JCylinder3D JSIRENE::getMaximumContainmentVolume ( )

inline

Forward function declarations.

Auxiliary function to retrieve the maximum cylindrical containment volume.

Returns

maximum cylindrical containment volume

Definition at line 48 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.
Note: The optional parameter can is used only when the given track
corresponds to a muon and if this track does not contain the
mc_usr_keys::energy_lost_in_can information, generated by JSirene.

Parameters

track	track
can	detector can

Returns

visible energy [GeV]

Definition at line 73 of file JVisibleEnergyToolkit.hh.

  {
    using namespace std;
    using namespace JPP;
 
    double Evis = 0.0;
    
    if (is_finalstate(track)) {
 
      const bool isMuon = is_muon(track);
      
      if (isMuon && track.haveusr(mc_usr_keys::energy_lost_in_can)) {
 
        Evis = track.getusr(mc_usr_keys::energy_lost_in_can);
        
      } else if (isMuon) {
 
        // 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;
  }

getVisibleEnergy() [2/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 192 of file JVisibleEnergyToolkit.hh.

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

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;      
    }    
  }

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 128 of file JVisibleEnergyToolkit.hh.

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

getVisibleEnergy() [3/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 143 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 168 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;
  }

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 207 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.

getNumberOfPhotoElectrons

const struct JSIRENE::number_of_photo_electrons_type JSIRENE::getNumberOfPhotoElectrons