simb::MCFlux Class Reference

#include <MCFlux.h>

Public Member Functions
	MCFlux ()
void	Reset ()
double	Flux (int pdgcode, int which=0) const
void	SetFluxPos (double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
void	SetFluxNeg (double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
void	SetFluxGen (double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
void	ReDecay (double &newE, double &newW, double x, double y, double z)

Public Attributes
int	frun
int	fevtno
double	fndxdz
double	fndydz
double	fnpz
double	fnenergy
double	fndxdznea
double	fndydznea
double	fnenergyn
double	fnwtnear
double	fndxdzfar
double	fndydzfar
double	fnenergyf
double	fnwtfar
int	fnorig
int	fndecay
int	fntype
double	fvx
double	fvy
double	fvz
double	fpdpx
double	fpdpy
double	fpdpz
double	fppdxdz
double	fppdydz
double	fpppz
double	fppenergy
int	fppmedium
int	fptype
double	fppvx
double	fppvy
double	fppvz
double	fmuparpx
double	fmuparpy
double	fmuparpz
double	fmupare
double	fnecm
double	fnimpwt
double	fxpoint
double	fypoint
double	fzpoint
double	ftvx
double	ftvy
double	ftvz
double	ftpx
double	ftpy
double	ftpz
int	ftptype
int	ftgen
int	ftgptype
double	ftgppx
double	ftgppy
double	ftgppz
double	ftprivx
double	ftprivy
double	ftprivz
double	fbeamx
double	fbeamy
double	fbeamz
double	fbeampx
double	fbeampy
double	fbeampz
simb::flux_code_	fFluxType
double	fgenx
	origin of ray from flux generator More...
double	fgeny
double	fgenz
double	fdk2gen
	distance from decay to ray origin More...
double	fgen2vtx
	distance from ray origin to event vtx More...

Private Attributes
float	fFluxPos [6]
	e,ebar,mu,mubar,tau,taubar flux, +horn focus More...
float	fFluxNeg [6]
	Fluxes as aboce, for negative horn focus. More...
float	fFluxGen [6]
	Fluxes as above, assumed by generator. More...

Friends
std::ostream &	operator<< (std::ostream &output, const simb::MCFlux &mcflux)

Detailed Description

Definition at line 26 of file MCFlux.h.

Constructor & Destructor Documentation

simb::MCFlux::MCFlux

(

)

Definition at line 18 of file MCFlux.cxx.

    : frun(-999)
    , fevtno(-999)
    , fndxdz(-999.)
    , fndydz(-999.)
    , fnpz(-999.)
    , fnenergy(-999.)
    , fndxdznea(-999.)
    , fndydznea(-999.)
    , fnenergyn(-999.)
    , fnwtnear(-999.)
    , fndxdzfar(-999.)
    , fndydzfar(-999.)
    , fnenergyf(-999.)
    , fnwtfar(-999.)
    , fnorig(-999)
    , fndecay(-999)
    , fntype(-999)
    , fvx(-999.)
    , fvy(-999.)
    , fvz(-999.)
    , fpdpx(-999.)
    , fpdpy(-999.)
    , fpdpz(-999.)
    , fppdxdz(-999.)
    , fppdydz(-999.)
    , fpppz(-999.)
    , fppenergy(-999.)
    , fppmedium(-999)
    , fptype(-999)     // converted to PDG
    , fppvx(-999.)
    , fppvy(-999.)
    , fppvz(-999.)
    , fmuparpx(-999.)
    , fmuparpy(-999.)
    , fmuparpz(-999.)
    , fmupare(-999.)
    , fnecm(-999.)
    , fnimpwt(-999.)
    , fxpoint(-999.)
    , fypoint(-999.)
    , fzpoint(-999.)
    , ftvx(-999.)
    , ftvy(-999.)
    , ftvz(-999.)
    , ftpx(-999.)
    , ftpy(-999.)
    , ftpz(-999.)
    , ftptype(-999)   // converted to PDG
    , ftgen(-999)
    , ftgptype(-999)  // converted to PDG
    , ftgppx(-999.)
    , ftgppy(-999.)
    , ftgppz(-999.)
    , ftprivx(-999.)
    , ftprivy(-999.)
    , ftprivz(-999.)
    , fbeamx(-999.)
    , fbeamy(-999.)
    , fbeamz(-999.)
    , fbeampx(-999.)
    , fbeampy(-999.)
    , fbeampz(-999.)
    , fFluxType(simb::kGenerator)
    , fgenx(-999.)
    , fgeny(-999.)
    , fgenz(-999.)
    , fdk2gen(-999.)
    , fgen2vtx(-999.)
  {

    for (int i=0; i<6; ++i) fFluxGen[i]=fFluxPos[i]=fFluxNeg[i]= 0;

  }

Member Function Documentation

double simb::MCFlux::Flux	(	int	pdgcode,
		int	which = 0
	)		const

Flux of neutrinos at the same energy

Parameters

pdgcode	: PDG code of neutrino flavor
which	: Which flux (generated, plus focus, minus focus)

Returns

: #neutrinos/cm^2/10^20 POT

Definition at line 176 of file MCFlux.cxx.

  {
    const float* flux = 0;
    if (which==kGenerator)      flux = fFluxGen;
    if (which==kHistPlusFocus)  flux = fFluxPos;
    if (which==kHistMinusFocus) flux = fFluxNeg;
    if (pdgcode== 12) return flux[0];
    if (pdgcode==-12) return flux[1];
    if (pdgcode== 14) return flux[2];
    if (pdgcode==-14) return flux[3];
    if (pdgcode== 16) return flux[4];
    if (pdgcode==-16) return flux[5];
    return 0.0;
  }

void simb::MCFlux::ReDecay	(	double &	newE,
		double &	newW,
		double	x,
		double	y,
		double	z
	)

Definition at line 224 of file MCFlux.cxx.

  {
    //note x,y,z are assumed to be in cm
    //x,y,z are also assumed to be in the beam reference frame
    //i.e. 0,0,0 are at the target,
    //z points along the beamline
    //right handed coords, x points to the left if you are looking down the beam

    //these constants should probably be defined elsewhere.
    //I'll put them here until I figure out where they should properly go
    const double pimass=.13957;
    const double mumass=0.105658389;
    const double kmass=0.49368;
    const double k0mass=0.49767;
    double mass=pimass;

    //determine parent particle mass based on fptype PDG code
    if(fptype==13||fptype==-13){
      mass = mumass;
    }
    else if(fptype==211||fptype==-211){
      mass = pimass;
    }
    else if(fptype==321||fptype==-321){
      mass = kmass;
    }
    else if(fptype==130){
      mass = k0mass;
    }
    else{
      std::cerr<<"dont know ptype "<<fptype<<"can not compute new decay weight"<<std::endl;
      return;
    }

    //compute kinematics of parent particle at decay point
    double p=sqrt(1.*fpdpx*fpdpx+1.*fpdpy*fpdpy+1.*fpdpz*fpdpz);
    double Eplab=sqrt(1.*fpdpx*fpdpx+1.*fpdpy*fpdpy+1.*fpdpz*fpdpz+1.*mass*mass);
    double gamma = Eplab/mass;
    double beta = sqrt((gamma*gamma-1)/(gamma*gamma));

    //compute components of vector between decay point
    //and the point you're aiming at
    double rnx=1.*(x-fvx);
    double rny=1.*(y-fvy);
    double rnz=1.*(z-fvz);
    double rn=sqrt(rnx*rnx+rny*rny+rnz*rnz);

    //compute angle between parent momentum
    //and where we want the neutrino to go
    double rndotp = (rnx*fpdpx+rny*fpdpy+rnz*fpdpz);
    double costhetan = rndotp/(rn*p);

    //do some checking of the calculation
    if(std::abs(costhetan)>1){
      if(costhetan>0){
        costhetan = 1;
      }
      else{
        costhetan=-1;
      }
    }

    //now compute the weights
    double MN=1.;
    if(p>0){//if it didn't decay at rest
      //boost
      MN=1./(gamma*(1-beta*costhetan));
    }


    newE = MN*fnecm;
    //  std::cout<<"New E "<<newE<<" old ecm "<<fnecm<<" mn "<<MN<<" gamma "<<gamma<<" beta "<<beta<<" costhetan "<<costhetan<<std::endl;

    //solid angle
    // small angle approximation: // double san = 10000./(4*rn*rn);
    // Alex Radovic's removal of small angle approximation
    const double kRDET = 100.;
    double san = (1.0-cos(atan( kRDET / rn )))/2.0;

    //  std::cout<<"san "<<san<<" fvz-z "<<fvz-z<<std::endl;
    newW = san*MN*MN;
    //  std::cout<<"mn "<<MN<<std::endl;


    //if its a (polarized) muon decay, we have modify the weight
    if(fptype==13||fptype==-13){
      double betav[3]={0.};
      double partialn=0.;
      double partial=0.;
      double P_nun[3]={0.};
      double P_dcm_nun[4]={0.};
      double P_pcm_mp[4]={0.};

      betav[0] = fpdpx/Eplab;
      betav[1] = fpdpy/Eplab;
      betav[2] = fpdpz/Eplab;

      P_nun[0] = rnx*newE/rn;
      P_nun[1] = rny*newE/rn;
      P_nun[2] = rnz*newE/rn;

      partialn =gamma*(betav[0]*P_nun[0]+betav[1]*P_nun[1]+betav[2]*P_nun[2]);
      partialn = newE - partialn /(gamma+1.);

      P_dcm_nun[0] = P_nun[0] - betav[0]*gamma*partialn;
      P_dcm_nun[1] = P_nun[1] - betav[1]*gamma*partialn;
      P_dcm_nun[2] = P_nun[2] - betav[2]*gamma*partialn;
      P_dcm_nun[3] = sqrt(pow(P_dcm_nun[0],2)
                          +pow(P_dcm_nun[1],2)
                          +pow(P_dcm_nun[2],2));


      gamma = fppenergy/mass;
      betav[0] = fppdxdz*fpppz/fppenergy;
      betav[1] = fppdydz*fpppz/fppenergy;
      betav[2] = fpppz/fppenergy;

      partial = gamma*(betav[0]*fmuparpx +
                       +betav[1]*fmuparpy + betav[2]*fmuparpz);
      partial = fmupare-partial/(gamma+1.);
      P_pcm_mp[0] = fmuparpx - betav[0]*gamma*partial;
      P_pcm_mp[1] = fmuparpy - betav[1]*gamma*partial;
      P_pcm_mp[2] = fmuparpz - betav[2]*gamma*partial;
      P_pcm_mp[3] = sqrt(pow(P_pcm_mp[0],2)+
                         pow(P_pcm_mp[1],2)+
                         pow(P_pcm_mp[2],2));

      //calc new  decay angle w.r.t. (anti)spin direction
      double costhn  = 0.;
      if(P_dcm_nun[3]!=0&&P_pcm_mp[3]!=0){
        costhn = ( P_dcm_nun[0]*P_pcm_mp[0]+
                  P_dcm_nun[1]*P_pcm_mp[1]+
                  P_dcm_nun[2]*P_pcm_mp[2])/(P_dcm_nun[3]*P_pcm_mp[3]);
      }

      if(std::abs(costhn)>1){
        if(costhn>0){
          costhn = 1;
        }
        else{
          costhn=-1;
        }
      }

      double wt_ration;
      if(fntype==14||fntype==-14){
        double xnu = 2.*fnecm/mass;
        wt_ration = ( (3.-2.*xnu) - (1.-2.*xnu)*costhn ) / (3.-2.*xnu);
      }
      else{
        wt_ration=1.-costhn;
      }

      newW*=wt_ration;
    }

    return;
  }

void simb::MCFlux::Reset

(

void

)

Definition at line 94 of file MCFlux.cxx.

  {
    frun      = -999;
    fevtno    = -999;
    fndxdz    = -999.;
    fndydz    = -999.;
    fnpz      = -999.;
    fnenergy  = -999.;
    fndxdznea = -999.;
    fndydznea = -999.;
    fnenergyn = -999.;
    fnwtnear  = -999.;
    fndxdzfar = -999.;
    fndydzfar = -999.;
    fnenergyf = -999.;
    fnwtfar   = -999.;
    fnorig    = -999;
    fndecay   = -999;
    fntype    = -999;
    fvx       = -999.;
    fvy       = -999.;
    fvz       = -999.;
    fpdpx     = -999.;
    fpdpy     = -999.;
    fpdpz     = -999.;
    fppdxdz   = -999.;
    fppdydz   = -999.;
    fpppz     = -999.;
    fppenergy = -999.;
    fppmedium = -999;
    fptype    = -999;     // converted to PDG
    fppvx     = -999.;
    fppvy     = -999.;
    fppvz     = -999.;
    fmuparpx  = -999.;
    fmuparpy  = -999.;
    fmuparpz  = -999.;
    fmupare   = -999.;
    fnecm     = -999.;
    fnimpwt   = -999.;
    fxpoint   = -999.;
    fypoint   = -999.;
    fzpoint   = -999.;
    ftvx      = -999.;
    ftvy      = -999.;
    ftvz      = -999.;
    ftpx      = -999.;
    ftpy      = -999.;
    ftpz      = -999.;
    ftptype   = -999;   // converted to PDG
    ftgen     = -999;
    ftgptype  = -999;  // converted to PDG
    ftgppx    = -999.;
    ftgppy    = -999.;
    ftgppz    = -999.;
    ftprivx   = -999.;
    ftprivy   = -999.;
    ftprivz   = -999.;
    fbeamx    = -999.;
    fbeamy    = -999.;
    fbeamz    = -999.;
    fbeampx   = -999.;
    fbeampy   = -999.;
    fbeampz   = -999.;

    fgenx     = -999.;
    fgeny     = -999.;
    fgenz     = -999.;
    fdk2gen   = -999.;
    fgen2vtx  = -999.;

    return;
  }

void simb::MCFlux::SetFluxGen	(	double	nue,
		double	nuebar,
		double	numu,
		double	numubar,
		double	nutau,
		double	nutaubar
	)

Definition at line 214 of file MCFlux.cxx.

  {
    fFluxGen[0] = nue;   fFluxGen[1] = nuebar;
    fFluxGen[2] = numu;  fFluxGen[3] = numubar;
    fFluxGen[4] = nutau; fFluxGen[5] = nutaubar;
  }

void simb::MCFlux::SetFluxNeg	(	double	nue,
		double	nuebar,
		double	numu,
		double	numubar,
		double	nutau,
		double	nutaubar
	)

Definition at line 204 of file MCFlux.cxx.

  {
    fFluxNeg[0] = nue;   fFluxNeg[1] = nuebar;
    fFluxNeg[2] = numu;  fFluxNeg[3] = numubar;
    fFluxNeg[4] = nutau; fFluxNeg[5] = nutaubar;
  }

void simb::MCFlux::SetFluxPos	(	double	nue,
		double	nuebar,
		double	numu,
		double	numubar,
		double	nutau,
		double	nutaubar
	)

Definition at line 193 of file MCFlux.cxx.

  {
    fFluxPos[0] = nue;   fFluxPos[1] = nuebar;
    fFluxPos[2] = numu;  fFluxPos[3] = numubar;
    fFluxPos[4] = nutau; fFluxPos[5] = nutaubar;
  }

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  output, const simb::MCFlux &  mcflux  )

friend

Definition at line 384 of file MCFlux.cxx.

  {
    output
      << "MCFlux:" << std::endl
      //  123456789012
      << "  flux job  " << std::setw(11) << mcflux.frun      << " "
      << "  pot #     " << std::setw(11) << mcflux.fevtno    << std::endl
      << "  ntype     " << std::setw(11) << mcflux.fntype    << " "
      << "  ptype     " << std::setw(11) << mcflux.fptype    << " "
      << "  impwt     " << std::setw(11) << mcflux.fnimpwt   << std::endl
      << "  ndecay    " << std::setw(11) << mcflux.fndecay   << " "
      << "  ppmedium  " << std::setw(11) << mcflux.fppmedium << " "
      << "  tptype    " << std::setw(11) << mcflux.ftptype   << std::endl
      << "  vxyz      " << std::setw(11) << mcflux.fvx       << " "
      <<                   std::setw(11) << mcflux.fvy       << " "
      <<                   std::setw(11) << mcflux.fvz       << std::endl
      << "  pdpxyz    " << std::setw(11) << mcflux.fpdpx     << " "
      <<                   std::setw(11) << mcflux.fpdpy     << " "
      <<                   std::setw(11) << mcflux.fpdpz     << std::endl
      << "  tpxyz     " << std::setw(11) << mcflux.ftpx      << " "
      <<                   std::setw(11) << mcflux.ftpy      << " "
      <<                   std::setw(11) << mcflux.ftpz      << std::endl
      << "  ppd[xy]dz " << std::setw(11) << mcflux.fppdxdz   << " "
      <<                   std::setw(11) << mcflux.fppdydz   << std::endl
      << "  pppxyz    " << std::setw(11) << mcflux.fppdxdz*mcflux.fpppz << " "
      <<                   std::setw(11) << mcflux.fppdydz*mcflux.fpppz << " "
      <<                   std::setw(11) << mcflux.fpppz     << std::endl
      << "  muparpxyz " << std::setw(11) << mcflux.fmuparpx  << " "
      <<                   std::setw(11) << mcflux.fmuparpy  << " "
      <<                   std::setw(11) << mcflux.fmuparpz  << " "
      << "  mupare    " << std::setw(11) << mcflux.fmupare   << std::endl
      << "  necm      " << std::setw(11) << mcflux.fnecm     << " "
      << "  dk2gen    " << std::setw(11) << mcflux.fdk2gen   << std::endl
      << "  near E    " << std::setw(11) << mcflux.fnenergyn << " "
      << "  wgt       " << std::setw(11) << mcflux.fnwtnear  << std::endl
      << "  far E     " << std::setw(11) << mcflux.fnenergyf << " "
      << "  wgt       " << std::setw(11) << mcflux.fnwtfar   << std::endl;
    return output;
  }

Member Data Documentation

double simb::MCFlux::fbeampx

Definition at line 94 of file MCFlux.h.

double simb::MCFlux::fbeampy

Definition at line 95 of file MCFlux.h.

double simb::MCFlux::fbeampz

Definition at line 96 of file MCFlux.h.

double simb::MCFlux::fbeamx

Definition at line 91 of file MCFlux.h.

double simb::MCFlux::fbeamy

Definition at line 92 of file MCFlux.h.

double simb::MCFlux::fbeamz

Definition at line 93 of file MCFlux.h.

double simb::MCFlux::fdk2gen

distance from decay to ray origin

Definition at line 103 of file MCFlux.h.

int simb::MCFlux::fevtno

Definition at line 36 of file MCFlux.h.

float simb::MCFlux::fFluxGen[6]

private

Fluxes as above, assumed by generator.

Definition at line 110 of file MCFlux.h.

float simb::MCFlux::fFluxNeg[6]

private

Fluxes as aboce, for negative horn focus.

Definition at line 109 of file MCFlux.h.

float simb::MCFlux::fFluxPos[6]

private

e,ebar,mu,mubar,tau,taubar flux, +horn focus

Definition at line 108 of file MCFlux.h.

simb::flux_code_ simb::MCFlux::fFluxType

Definition at line 98 of file MCFlux.h.

double simb::MCFlux::fgen2vtx

distance from ray origin to event vtx

Definition at line 104 of file MCFlux.h.

double simb::MCFlux::fgenx

origin of ray from flux generator

Definition at line 100 of file MCFlux.h.

double simb::MCFlux::fgeny

Definition at line 101 of file MCFlux.h.

double simb::MCFlux::fgenz

Definition at line 102 of file MCFlux.h.

double simb::MCFlux::fmupare

Definition at line 70 of file MCFlux.h.

double simb::MCFlux::fmuparpx

Definition at line 67 of file MCFlux.h.

double simb::MCFlux::fmuparpy

Definition at line 68 of file MCFlux.h.

double simb::MCFlux::fmuparpz

Definition at line 69 of file MCFlux.h.

int simb::MCFlux::fndecay

Definition at line 50 of file MCFlux.h.

double simb::MCFlux::fndxdz

Definition at line 37 of file MCFlux.h.

double simb::MCFlux::fndxdzfar

Definition at line 45 of file MCFlux.h.

double simb::MCFlux::fndxdznea

Definition at line 41 of file MCFlux.h.

double simb::MCFlux::fndydz

Definition at line 38 of file MCFlux.h.

double simb::MCFlux::fndydzfar

Definition at line 46 of file MCFlux.h.

double simb::MCFlux::fndydznea

Definition at line 42 of file MCFlux.h.

double simb::MCFlux::fnecm

Definition at line 71 of file MCFlux.h.

double simb::MCFlux::fnenergy

Definition at line 40 of file MCFlux.h.

double simb::MCFlux::fnenergyf

Definition at line 47 of file MCFlux.h.

double simb::MCFlux::fnenergyn

Definition at line 43 of file MCFlux.h.

double simb::MCFlux::fnimpwt

Definition at line 72 of file MCFlux.h.

int simb::MCFlux::fnorig

Definition at line 49 of file MCFlux.h.

double simb::MCFlux::fnpz

Definition at line 39 of file MCFlux.h.

int simb::MCFlux::fntype

Definition at line 51 of file MCFlux.h.

double simb::MCFlux::fnwtfar

Definition at line 48 of file MCFlux.h.

double simb::MCFlux::fnwtnear

Definition at line 44 of file MCFlux.h.

double simb::MCFlux::fpdpx

Definition at line 55 of file MCFlux.h.

double simb::MCFlux::fpdpy

Definition at line 56 of file MCFlux.h.

double simb::MCFlux::fpdpz

Definition at line 57 of file MCFlux.h.

double simb::MCFlux::fppdxdz

Definition at line 58 of file MCFlux.h.

double simb::MCFlux::fppdydz

Definition at line 59 of file MCFlux.h.

double simb::MCFlux::fppenergy

Definition at line 61 of file MCFlux.h.

int simb::MCFlux::fppmedium

Definition at line 62 of file MCFlux.h.

double simb::MCFlux::fpppz

Definition at line 60 of file MCFlux.h.

double simb::MCFlux::fppvx

Definition at line 64 of file MCFlux.h.

double simb::MCFlux::fppvy

Definition at line 65 of file MCFlux.h.

double simb::MCFlux::fppvz

Definition at line 66 of file MCFlux.h.

int simb::MCFlux::fptype

Definition at line 63 of file MCFlux.h.

int simb::MCFlux::frun

Definition at line 35 of file MCFlux.h.

int simb::MCFlux::ftgen

Definition at line 83 of file MCFlux.h.

double simb::MCFlux::ftgppx

Definition at line 85 of file MCFlux.h.

double simb::MCFlux::ftgppy

Definition at line 86 of file MCFlux.h.

double simb::MCFlux::ftgppz

Definition at line 87 of file MCFlux.h.

int simb::MCFlux::ftgptype

Definition at line 84 of file MCFlux.h.

double simb::MCFlux::ftprivx

Definition at line 88 of file MCFlux.h.

double simb::MCFlux::ftprivy

Definition at line 89 of file MCFlux.h.

double simb::MCFlux::ftprivz

Definition at line 90 of file MCFlux.h.

int simb::MCFlux::ftptype

Definition at line 82 of file MCFlux.h.

double simb::MCFlux::ftpx

Definition at line 79 of file MCFlux.h.

double simb::MCFlux::ftpy

Definition at line 80 of file MCFlux.h.

double simb::MCFlux::ftpz

Definition at line 81 of file MCFlux.h.

double simb::MCFlux::ftvx

Definition at line 76 of file MCFlux.h.

double simb::MCFlux::ftvy

Definition at line 77 of file MCFlux.h.

double simb::MCFlux::ftvz

Definition at line 78 of file MCFlux.h.

double simb::MCFlux::fvx

Definition at line 52 of file MCFlux.h.

double simb::MCFlux::fvy

Definition at line 53 of file MCFlux.h.

double simb::MCFlux::fvz

Definition at line 54 of file MCFlux.h.

double simb::MCFlux::fxpoint

Definition at line 73 of file MCFlux.h.

double simb::MCFlux::fypoint

Definition at line 74 of file MCFlux.h.

double simb::MCFlux::fzpoint

Definition at line 75 of file MCFlux.h.

---

The documentation for this class was generated from the following files:

• nusimdata/nusimdata/SimulationBase/MCFlux.h
• nusimdata/nusimdata/SimulationBase/MCFlux.cxx