genie::flux::GCylindTH1Flux Class Reference

A generic GENIE flux driver. Generates a 'cylindrical' neutrino beam along the input direction, with the input transverse radius and centered at the input position. The energies are generated from the input energy spectrum (TH1D). Multiple neutrino species can be generated (you will need to supply an energy spectrum for each). More...

#include <GCylindTH1Flux.h>

Inheritance diagram for genie::flux::GCylindTH1Flux:

Public Member Functions
	GCylindTH1Flux ()
	~GCylindTH1Flux ()
void	SetNuDirection (const TVector3 &direction)
void	SetBeamSpot (const TVector3 &spot)
void	SetTransverseRadius (double Rt)
void	AddEnergySpectrum (int nu_pdgc, TH1D *spectrum)
void	SetRtDependence (string rdep)
const PDGCodeList &	FluxParticles (void)
	declare list of flux neutrinos that can be generated (for init. purposes) More...
double	MaxEnergy (void)
	declare the max flux neutrino energy that can be generated (for init. purposes) More...
bool	GenerateNext (void)
	generate the next flux neutrino (return false in err) More...
int	PdgCode (void)
	returns the flux neutrino pdg code More...
double	Weight (void)
	returns the flux neutrino weight (if any) More...
const TLorentzVector &	Momentum (void)
	returns the flux neutrino 4-momentum More...
const TLorentzVector &	Position (void)
	returns the flux neutrino 4-position (note: expect SI rather than physical units) More...
bool	End (void)
	true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples) More...
long int	Index (void)
	returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number) More...
void	Clear (Option_t *opt)
	reset state variables based on opt More...
void	GenerateWeighted (bool gen_weighted)
	set whether to generate weighted or unweighted neutrinos More...
Public Member Functions inherited from genie::GFluxI
virtual	~GFluxI ()

Private Member Functions
void	Initialize (void)
void	CleanUp (void)
void	ResetSelection (void)
void	AddAllFluxes (void)
int	SelectNeutrino (double Ev)
double	GeneratePhi (void) const
double	GenerateRt (void) const

Private Attributes
double	fMaxEv
	maximum energy More...
PDGCodeList *	fPdgCList
	list of neutrino pdg-codes More...
int	fgPdgC
	running generated nu pdg-code More...
TLorentzVector	fgP4
	running generated nu 4-momentum More...
TLorentzVector	fgX4
	running generated nu 4-position More...
vector< TH1D * >	fSpectrum
	flux = f(Ev), 1/neutrino species More...
TH1D *	fTotSpectrum
	combined flux = f(Ev) More...
TVector3 *	fDirVec
	neutrino direction More...
TVector3 *	fBeamSpot
	beam spot position More...
double	fRt
	transverse size of neutrino beam More...
TF1 *	fRtDep
	transverse radius dependence More...

Additional Inherited Members
Protected Member Functions inherited from genie::GFluxI
	GFluxI ()

Detailed Description

A generic GENIE flux driver. Generates a 'cylindrical' neutrino beam along the input direction, with the input transverse radius and centered at the input position. The energies are generated from the input energy spectrum (TH1D). Multiple neutrino species can be generated (you will need to supply an energy spectrum for each).

Author

Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

July 4, 2005

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 43 of file GCylindTH1Flux.h.

Constructor & Destructor Documentation

GCylindTH1Flux::GCylindTH1Flux

(

)

Definition at line 35 of file GCylindTH1Flux.cxx.

{
  this->Initialize();
}

GCylindTH1Flux::~GCylindTH1Flux

(

)

Definition at line 40 of file GCylindTH1Flux.cxx.

{
  this->CleanUp();
}

Member Function Documentation

void GCylindTH1Flux::AddAllFluxes ( void  )

private

Definition at line 218 of file GCylindTH1Flux.cxx.

{
  LOG("Flux", pNOTICE) << "Computing combined flux";

  if(fTotSpectrum) delete fTotSpectrum;

  vector<TH1D *>::const_iterator spectrum_iter;

  unsigned int inu=0;
  for(spectrum_iter = fSpectrum.begin();
                       spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
     TH1D * spectrum = *spectrum_iter;

     if(inu==0) { fTotSpectrum = new TH1D(*spectrum); }
     else       { fTotSpectrum->Add(spectrum);        }
     inu++;
  }
}

void GCylindTH1Flux::AddEnergySpectrum	(	int	nu_pdgc,
		TH1D *	spectrum
	)

Definition at line 183 of file GCylindTH1Flux.cxx.

{
  LOG("Flux", pNOTICE) << "Adding flux spectrum for pdg = " << nu_pdgc;

  fPdgCList->push_back(nu_pdgc);

  bool accepted = (count(fPdgCList->begin(),fPdgCList->end(),nu_pdgc) == 1);
  if(!accepted) {
     LOG ("Flux", pWARN)
            << "The pdg-code isn't recognized and the spectrum was ignored";
  } else {
     fSpectrum.push_back(spectrum);

     int    nb  = spectrum->GetNbinsX();
     Axis_t max = spectrum->GetBinLowEdge(nb)+spectrum->GetBinWidth(nb);
     fMaxEv = TMath::Max(fMaxEv, (double)max);

     LOG("Flux", pNOTICE)
          << "Updating maximum energy of flux particles to: " << fMaxEv;

     this->AddAllFluxes(); // update combined flux
  }
}

void GCylindTH1Flux::CleanUp ( void  )

private

Definition at line 145 of file GCylindTH1Flux.cxx.

{
  LOG("Flux", pNOTICE) << "Cleaning up...";

  if (fDirVec     ) delete fDirVec;
  if (fBeamSpot   ) delete fBeamSpot;
  if (fPdgCList   ) delete fPdgCList;
  if (fTotSpectrum) delete fTotSpectrum;
  if (fRtDep      ) delete fRtDep;

  unsigned int nspectra = fSpectrum.size();
  for(unsigned int i = 0; i < nspectra; i++) {
     TH1D * spectrum = fSpectrum[i];
     delete spectrum;
     spectrum = 0;
  }
}

void GCylindTH1Flux::Clear ( Option_t *  opt )

virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 103 of file GCylindTH1Flux.cxx.

{
// Dummy clear method needed to conform to GFluxI interface
//
  LOG("Flux", pERROR) <<
      "No Clear(Option_t * opt) method implemented for opt: "<< opt;
}

bool genie::flux::GCylindTH1Flux::End ( void  )

inlinevirtual

true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)

Implements genie::GFluxI.

Definition at line 64 of file GCylindTH1Flux.h.

{ return  false;     }

const PDGCodeList& genie::flux::GCylindTH1Flux::FluxParticles ( void  )

inlinevirtual

declare list of flux neutrinos that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 57 of file GCylindTH1Flux.h.

{ return *fPdgCList; }

bool GCylindTH1Flux::GenerateNext ( void  )

virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 45 of file GCylindTH1Flux.cxx.

{
  //-- Reset previously generated neutrino code / 4-p / 4-x
  this->ResetSelection();

  //-- Generate an energy from the 'combined' spectrum histogram
  //   and compute the momentum vector
  double Ev = (double) fTotSpectrum->GetRandom();

  TVector3 p3(*fDirVec); // momentum along the neutrino direction
  p3.SetMag(Ev);         // with |p|=Ev
  // EDIT: Check if we're running with DM beam
  if (fPdgCList->ExistsInPDGCodeList(kPdgDarkMatter) || fPdgCList->ExistsInPDGCodeList(kPdgAntiDarkMatter)) {
    double Md = PDGLibrary::Instance()->Find(kPdgDarkMatter)->Mass();
    p3.SetMag(TMath::Sqrt(Ev*Ev - Md*Md));
  }

  fgP4.SetPxPyPzE(p3.Px(), p3.Py(), p3.Pz(), Ev);

  //-- Select a neutrino species from the flux fractions at the
  //   selected energy
  fgPdgC = (*fPdgCList)[this->SelectNeutrino(Ev)];

  //-- Compute neutrino 4-x

  if(fRt <= 0) {
    fgX4.SetXYZT(0.,0.,0.,0.);
  }
  else {
    // Create a vector (vec) that points to a random position at a disk
    // of radius Rt passing through the origin, perpendicular to the
    // input direction.
    TVector3 vec0(*fDirVec);           // vector along the input direction
    TVector3 vec = vec0.Orthogonal();  // orthogonal vector

    double psi = this->GeneratePhi();  // rndm angle [0,2pi]
    double Rt  = this->GenerateRt();   // rndm R [0,Rtransverse]

    vec.Rotate(psi,vec0); // rotate around original vector
    vec.SetMag(Rt);       // set new norm

    // Set the neutrino position as beam_spot + vec
    double x = fBeamSpot->X() + vec.X();
    double y = fBeamSpot->Y() + vec.Y();
    double z = fBeamSpot->Z() + vec.Z();

    fgX4.SetXYZT(x,y,z,0.);
  }

  LOG("Flux", pINFO) << "Generated neutrino pdg-code: " << fgPdgC;
  LOG("Flux", pINFO)
        << "Generated neutrino p4: " << utils::print::P4AsShortString(&fgP4);
  LOG("Flux", pINFO)
             << "Generated neutrino x4: " << utils::print::X4AsString(&fgX4);

  return true;
}

double GCylindTH1Flux::GeneratePhi ( void  ) const

private

Definition at line 271 of file GCylindTH1Flux.cxx.

{
  RandomGen * rnd = RandomGen::Instance();
  double phi = 2.*kPi * rnd->RndFlux().Rndm(); // [0,2pi]
  return phi;
}

double GCylindTH1Flux::GenerateRt ( void  ) const

private

Definition at line 278 of file GCylindTH1Flux.cxx.

{
  double Rt = fRtDep->GetRandom(); // rndm R [0,Rtransverse]
  return Rt;
}

void GCylindTH1Flux::GenerateWeighted ( bool  gen_weighted )

virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 111 of file GCylindTH1Flux.cxx.

{
// Dummy implementation needed to conform to GFluxI interface
//
  LOG("Flux", pERROR) <<
      "No GenerateWeighted(bool gen_weighted) method implemented for " <<
      "gen_weighted: " << gen_weighted;
}

long int genie::flux::GCylindTH1Flux::Index ( void  )

inlinevirtual

returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)

Implements genie::GFluxI.

Definition at line 65 of file GCylindTH1Flux.h.

{ return -1;         }

void GCylindTH1Flux::Initialize ( void  )

private

Definition at line 120 of file GCylindTH1Flux.cxx.

{
  LOG("Flux", pNOTICE) << "Initializing GCylindTH1Flux driver";

  fMaxEv       = 0;
  fPdgCList    = new PDGCodeList;
  fTotSpectrum = 0;
  fDirVec      = 0;
  fBeamSpot    = 0;
  fRt          =-1;
  fRtDep       = 0;

  this->ResetSelection();
  this->SetRtDependence("x");
  //eg, other example: this->SetRtDependence("pow(x,2)");
}

double genie::flux::GCylindTH1Flux::MaxEnergy ( void  )

inlinevirtual

declare the max flux neutrino energy that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 58 of file GCylindTH1Flux.h.

{ return  fMaxEv;    }

const TLorentzVector& genie::flux::GCylindTH1Flux::Momentum ( void  )

inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 62 of file GCylindTH1Flux.h.

{ return  fgP4;      }

int genie::flux::GCylindTH1Flux::PdgCode ( void  )

inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 60 of file GCylindTH1Flux.h.

{ return  fgPdgC;    }

const TLorentzVector& genie::flux::GCylindTH1Flux::Position ( void  )

inlinevirtual

returns the flux neutrino 4-position (note: expect SI rather than physical units)

Implements genie::GFluxI.

Definition at line 63 of file GCylindTH1Flux.h.

{ return  fgX4;      }

void GCylindTH1Flux::ResetSelection ( void  )

private

Definition at line 137 of file GCylindTH1Flux.cxx.

{
// initializing running neutrino pdg-code, 4-position, 4-momentum
  fgPdgC = 0;
  fgP4.SetPxPyPzE (0.,0.,0.,0.);
  fgX4.SetXYZT    (0.,0.,0.,0.);
}

int GCylindTH1Flux::SelectNeutrino ( double  Ev )

private

Definition at line 237 of file GCylindTH1Flux.cxx.

{
  const unsigned int n = fPdgCList->size();
  double fraction[n];

  vector<TH1D *>::const_iterator spectrum_iter;

  unsigned int inu=0;
  for(spectrum_iter = fSpectrum.begin();
                       spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
     TH1D * spectrum = *spectrum_iter;
     fraction[inu++] = spectrum->GetBinContent(spectrum->FindBin(Ev));
  }

  double sum = 0;
  for(inu = 0; inu < n; inu++) {
     sum += fraction[inu];
     fraction[inu] = sum;
     LOG("Flux", pDEBUG) << "SUM-FRACTION(0->" << inu <<") = " << sum;
  }

  RandomGen * rnd = RandomGen::Instance();
  double R = sum * rnd->RndFlux().Rndm();

  LOG("Flux", pDEBUG) << "R e [0,SUM] = " << R;

  for(inu = 0; inu < n; inu++) {if ( R < fraction[inu] ) return inu;}

  LOG("Flux", pERROR) << "Could not select a neutrino species";
  assert(false);

  return -1;
}

void GCylindTH1Flux::SetBeamSpot

(

const TVector3 &

spot

)

Definition at line 169 of file GCylindTH1Flux.cxx.

{
  if(fBeamSpot) delete fBeamSpot;
  fBeamSpot = new TVector3(spot);
}

void GCylindTH1Flux::SetNuDirection

(

const TVector3 &

direction

)

Definition at line 163 of file GCylindTH1Flux.cxx.

{
  if(fDirVec) delete fDirVec;
  fDirVec = new TVector3(direction);
}

void GCylindTH1Flux::SetRtDependence

(

string

rdep

)

Definition at line 207 of file GCylindTH1Flux.cxx.

{
// Set the (functional form of) Rt dependence as string, eg "x*x+sin(x)"
// You do not need to set this method. The default behaviour is to generate
// flux neutrinos uniformly over the area of the cylinder's cross section.

  if(fRtDep) delete fRtDep;

  fRtDep = new TF1("rdep", rdep.c_str(), 0,fRt);
}

void GCylindTH1Flux::SetTransverseRadius

(

double

Rt

)

Definition at line 175 of file GCylindTH1Flux.cxx.

{
  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rt;
  fRt = Rt;

  if(fRtDep) fRtDep->SetRange(0,Rt);
}

double genie::flux::GCylindTH1Flux::Weight ( void  )

inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 61 of file GCylindTH1Flux.h.

{ return  1.0;       }

Member Data Documentation

TVector3* genie::flux::GCylindTH1Flux::fBeamSpot

private

beam spot position

Definition at line 89 of file GCylindTH1Flux.h.

TVector3* genie::flux::GCylindTH1Flux::fDirVec

private

neutrino direction

Definition at line 88 of file GCylindTH1Flux.h.

TLorentzVector genie::flux::GCylindTH1Flux::fgP4

private

running generated nu 4-momentum

Definition at line 84 of file GCylindTH1Flux.h.

int genie::flux::GCylindTH1Flux::fgPdgC

private

running generated nu pdg-code

Definition at line 83 of file GCylindTH1Flux.h.

TLorentzVector genie::flux::GCylindTH1Flux::fgX4

private

running generated nu 4-position

Definition at line 85 of file GCylindTH1Flux.h.

double genie::flux::GCylindTH1Flux::fMaxEv

private

maximum energy

Definition at line 81 of file GCylindTH1Flux.h.

PDGCodeList* genie::flux::GCylindTH1Flux::fPdgCList

private

list of neutrino pdg-codes

Definition at line 82 of file GCylindTH1Flux.h.

double genie::flux::GCylindTH1Flux::fRt

private

transverse size of neutrino beam

Definition at line 90 of file GCylindTH1Flux.h.

TF1* genie::flux::GCylindTH1Flux::fRtDep

private

transverse radius dependence

Definition at line 91 of file GCylindTH1Flux.h.

vector<TH1D *> genie::flux::GCylindTH1Flux::fSpectrum

private

flux = f(Ev), 1/neutrino species

Definition at line 86 of file GCylindTH1Flux.h.

TH1D* genie::flux::GCylindTH1Flux::fTotSpectrum

private

combined flux = f(Ev)

Definition at line 87 of file GCylindTH1Flux.h.

---

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

• Generator/src/Tools/Flux/GCylindTH1Flux.h
• Generator/src/Tools/Flux/GCylindTH1Flux.cxx