GAstroFlux.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \class genie::flux::GAstroFlux
5 
6 \brief A base class for the concrete astrophysical neutrino flux drivers.
7 
8 
9  <<<< NOTE: CODE UNDER DEVELOPMENT >>>>
10 
11 
12  COORDINATE SYSTEMS / NEUTRINO GENERATION & PROPAGATION :
13 
14  Neutrinos are generated on a sphere with radius R_{earth}.
15  Especially,
16  - For diffuse fluxes:
17  Neutrinos can be generated anywhere on that surface.
18  - For point sources:
19  Neutrinos are generated at fixed right ascension and declination.
20  Then time is randomized, to account for Earth's rotation, and the
21  Equatorial Coordinates are converted to GEF. So, neutrinos from each
22  point source are generated on circles parallel to the Earth's Equator.
23 
24  Initially, neutrino coordinates are generated in the Geocentric Earth-
25  Fixed (GEF) Coordinate System (later to be converted to the appropriate
26  detector coordinate system - See further below).
27 
28  * Definition:
29  Geocentric Earth-Fixed (GEF) Coordinate System
30  +z: Points to North Pole.
31  xy: Equatorial plane.
32  +x: Points to the Prime Meridian.
33  +y: As needed to make a right-handed coordinate system.
34 
35  Neutrinos are then propagated towards the detector.
36  The Earth opaqueness to ultra high energy neutrinos is taken into
37  account. The Earth density profile is modelled using the PREM
38  (Preliminary Earth Model, The Encyclopedia of Solid Earth Geophysics,
39  David E. James, ed., Van Nostrand Reinhold, New York, 1989, p.331).
40 
41  The detector position is determined in the Spherical/Geographic System
42  by its geographic latitude (angle relative to Equator), its geographic
43  longitude (angle relative to Prime Meridian) and its depth from the
44  surface.
45 
46  The generated flux neutrinos, propagated through the Earth towards the
47  detector) are then positioned on the surface of a sphere with radius Rd
48  which should fully enclose the neutrino detector. The centre of that
49  sphere is taken to be the origin of the detector coordinate system.
50  The transverse coords are appropriately randomized so that neutrinos
51  from any given direction bath the entire sphere enclosing the detector.
52 
53  The final flux neutrino coordinates are given in the detector coordinate
54  system. The default detector coordinate system is the Topocentric Horizontal
55  (THZ) Coordinate System. Alternative user-defined topocentric systems can
56  be defined by specifying the appropriate rotation from THZ.
57 
58  * Definition:
59  Topocentric Horizontal (THZ) Coordinate System
60  (default detector coordinate system)
61  +z: Points towards the local zenith.
62  +x: On same plane as local meridian, pointing south.
63  +y: As needed to make a right-handed coordinate system.
64  origin: detector centre
65 
66  WEIGHTING SCHEMES:
67 
68  For a detector with geometrical cross section ~ 1km^2, the solid
69  angle acceptance changes by ~10 orders of magnitude across the
70  surface of the Earth.
71 
72  The driver supports both weighted and un-weighted flux generation
73  schemes. However, because of the enormous changes in solid angle
74  acceptance and energy, only the weighted scheme is practical.
75 
76  PHYSICS:
77 
78  The relative neutrino population needs to be set by the user using
79  the SetRelNuPopulations() method.
80  If run without arguments, the following relative populations are set:
81  nue:numu:nutau:nuebar:numubar:nutaubar = 1:2:0:1:2:0
82 
83  The energy spectrum is follows a power law. The user needs to
84  specify the power-law index by calling SetEnergyPowLawIdx().
85 
86 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
87  University of Liverpool & STFC Rutherford Appleton Lab
88 
89 \created March 27, 2010
90 
91 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
92  For the full text of the license visit http://copyright.genie-mc.org
93  or see $GENIE/LICENSE
94 */
95 //____________________________________________________________________________
96 
97 #ifndef _GASTRO_FLUX_H_
98 #define _GASTRO_FLUX_H_
99 
100 #include <string>
101 #include <map>
102 
103 #include <TLorentzVector.h>
104 #include <TVector3.h>
105 #include <TRotation.h>
106 
109 
110 class TH1D;
111 class TH2D;
112 
113 using std::string;
114 using std::map;
115 
116 namespace genie {
117 namespace flux {
118 
119 const double kAstroDefMaxEv = 1E+20 * units::GeV; ///<
120 const double kAstroDefMinEv = 1E-3 * units::GeV; ///<
121 const int kAstroNlog10EvBins = 1000; ///<
122 const int kAstroNCosThetaBins = 500; ///<
123 const int kAstroNPhiBins = 500; ///<
124 
125 //
126 //
127 //
128 
129 class GAstroFlux: public GFluxI {
130 
131 public :
132  virtual ~GAstroFlux();
133 
134  //
135  // methods implementing the GENIE GFluxI interface
136  //
137  virtual const PDGCodeList & FluxParticles (void) { return *fPdgCList; }
138  virtual double MaxEnergy (void);
139  virtual bool GenerateNext (void);
140  virtual int PdgCode (void) { return fgPdgC; }
141  virtual double Weight (void) { return fgWeight; }
142  virtual const TLorentzVector & Momentum (void) { return fgP4; }
143  virtual const TLorentzVector & Position (void) { return fgX4; }
144  virtual bool End (void) { return false; }
145  virtual long int Index (void) { return -1; }
146  virtual void Clear (Option_t * opt);
147  virtual void GenerateWeighted (bool gen_weighted);
148 
149  //
150  // configuration methods specific to all astrophysical neutrino flux drivers
151  //
152  void ForceMinEnergy (double emin);
153  void ForceMaxEnergy (double emax);
154  void SetDetectorPosition (double latitude, double longitude, double depth, double size);
155  void SetRelNuPopulations (double nnue=1, double nnumu=2, double nnutau=0, double nnuebar=1, double nnumubar=2, double nnutaubar=0);
156  void SetEnergyPowLawIdx (double n);
157  void SetUserCoordSystem (TRotation & rotation); ///< rotation Topocentric Horizontal -> User-defined Topocentric Coord System
158 
159 protected:
160 
161  class NuGenerator;
162  class NuPropagator;
163 
164  // abstract class, ctor hidden
165  GAstroFlux();
166 
167  // protected methods
168  void Initialize (void);
169  void CleanUp (void);
170  void ResetSelection (void);
171 
172  //
173  // protected data members
174  //
175 
176  PDGCodeList * fPdgCList; ///< declared list of neutrino pdg-codes that can be thrown by current instance
177  int fgPdgC; ///< (current) generated nu pdg-code
178  TLorentzVector fgP4; ///< (current) generated nu 4-momentum
179  TLorentzVector fgX4; ///< (current) generated nu 4-position
180  double fgWeight; ///< (current) generated nu weight
181  // configuration properties set by the user
182  double fMaxEvCut; ///< (config) user-defined maximum energy cut
183  double fMinEvCut; ///< (config) user-defined minimum energy cut
184  bool fGenWeighted; ///< (config) generate a weighted or unweighted flux?
185  double fDetGeoLatitude; ///< (config) detector: geographic latitude
186  double fDetGeoLongitude; ///< (config) detector: geographic longitude
187  double fDetGeoDepth; ///< (config) detector: depth from surface
188  double fDetSize; ///< (config) detector: size (detector should be enclosed in sphere of this radius)
189  map<int,double> fRelNuPopulations; ///< (config) relative neutrino populations
190  TRotation fRotGEF2THz; ///< (config) coord. system rotation: GEF translated to detector centre -> THZ
191  TRotation fRotTHz2User; ///< (config) coord. system rotation: THZ -> Topocentric user-defined
192  // internal flags and utility objects
193  TVector3 fDetCenter; ///<
194  TH1D * fEnergySpectrum; ///<
198 
199  //
200  // utility classes
201  //
202  class NuGenerator {
203  public:
206  bool SelectNuPdg (bool weighted, const map<int,double> & nupdgpdf, int & nupdg, double & wght);
207  bool SelectEnergy(bool weighted, TH1D & log10epdf, double log10emin, double log10emax, double & log10e, double & wght);
208  bool SelectOrigin(bool weighted, TH2D & opdf, double & phi, double & costheta, double & wght);
209  };
210  class NuPropagator {
211  public:
212  NuPropagator(double stepsz) : fStepSize(stepsz/units::km) { }
214  bool Go(double phi_start, double costheta_start, const TVector3 & detector_centre, double detector_sz, int nu_pdg, double Ev);
215  int NuPdgAtDetVolBoundary (void) { return fNuPdg; }
216  TVector3 & X3AtDetVolBoundary (void) { return fX3; }
217  TVector3 & P3AtDetVolBoundary (void) { return fP3; }
218  private:
219  double fStepSize;
220  int fNuPdg;
221  TVector3 fX3;
222  TVector3 fP3;
223  };
224 
225 };
226 
227 
228 //
229 // Concrete astrophysical flux drivers
230 //
231 
232 //............................................................................
233 // GENIE diffuse astrophysical neutrino flux driver
234 //
236 public :
239 
240  //
241  //
242 };
243 
244 //............................................................................
245 // GENIE concrete flux driver for astrophysical point neutrino sources
246 //
248 public :
251 
252  //
253  bool GenerateNext (void);
254 
255  //
256  void AddPointSource(string name, double ra, double dec, double rel_intensity);
257 
258 private:
259 
260  bool SelectSource(void);
261 
262  map<int, string> fPntSrcName; ///< point source name
263  map<int, double> fPntSrcRA; ///< right ascension
264  map<int, double> fPntSrcDec; ///< declination
265  map<int, double> fPntSrcRelI; ///< relative intensity
266  double fPntSrcTotI; ///< sum of all relative intensities
267 
268  unsigned int fSelSourceId;
269 };
270 //............................................................................
271 
272 } // flux namespace
273 } // genie namespace
274 
275 #endif // _GASTRO_FLUX_H_
276 
virtual const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GAstroFlux.h:142
void ForceMinEnergy(double emin)
Definition: GAstroFlux.cxx:150
TLorentzVector fgP4
(current) generated nu 4-momentum
Definition: GAstroFlux.h:178
virtual void Clear(Option_t *opt)
reset state variables based on opt
Definition: GAstroFlux.cxx:162
const XML_Char * name
Definition: expat.h:151
double fDetSize
(config) detector: size (detector should be enclosed in sphere of this radius)
Definition: GAstroFlux.h:188
double fMinEvCut
(config) user-defined minimum energy cut
Definition: GAstroFlux.h:183
bool fGenWeighted
(config) generate a weighted or unweighted flux?
Definition: GAstroFlux.h:184
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAstroFlux.cxx:54
void SetEnergyPowLawIdx(double n)
Definition: GAstroFlux.cxx:270
TRotation fRotTHz2User
(config) coord. system rotation: THZ -> Topocentric user-defined
Definition: GAstroFlux.h:191
double fMaxEvCut
(config) user-defined maximum energy cut
Definition: GAstroFlux.h:182
unsigned int nnue
Definition: runWimpSim.h:90
void SetUserCoordSystem(TRotation &rotation)
rotation Topocentric Horizontal -> User-defined Topocentric Coord System
Definition: GAstroFlux.cxx:293
const double kAstroDefMaxEv
Definition: GAstroFlux.h:119
void ForceMaxEnergy(double emax)
Definition: GAstroFlux.cxx:156
map< int, double > fRelNuPopulations
(config) relative neutrino populations
Definition: GAstroFlux.h:189
map< int, double > fPntSrcRelI
relative intensity
Definition: GAstroFlux.h:265
map< int, string > fPntSrcName
point source name
Definition: GAstroFlux.h:262
Loaders::FluxType flux
static constexpr Double_t km
Definition: Munits.h:148
const int kAstroNPhiBins
Definition: GAstroFlux.h:123
virtual bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GAstroFlux.cxx:59
A list of PDG codes.
Definition: PDGCodeList.h:33
bool SelectNuPdg(bool weighted, const map< int, double > &nupdgpdf, int &nupdg, double &wght)
Definition: GAstroFlux.cxx:375
virtual const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
Definition: GAstroFlux.h:137
const double emin
double fPntSrcTotI
sum of all relative intensities
Definition: GAstroFlux.h:266
double fDetGeoDepth
(config) detector: depth from surface
Definition: GAstroFlux.h:187
TRotation fRotGEF2THz
(config) coord. system rotation: GEF translated to detector centre -> THZ
Definition: GAstroFlux.h:190
NuGenerator * fNuGen
Definition: GAstroFlux.h:196
virtual double Weight(void)
returns the flux neutrino weight (if any)
Definition: GAstroFlux.h:141
Float_t E
Definition: plot.C:20
const double emax
double fDetGeoLongitude
(config) detector: geographic longitude
Definition: GAstroFlux.h:186
bool SelectEnergy(bool weighted, TH1D &log10epdf, double log10emin, double log10emax, double &log10e, double &wght)
Definition: GAstroFlux.cxx:423
virtual long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
Definition: GAstroFlux.h:145
double fDetGeoLatitude
(config) detector: geographic latitude
Definition: GAstroFlux.h:185
void SetDetectorPosition(double latitude, double longitude, double depth, double size)
Definition: GAstroFlux.cxx:174
PDGCodeList * fPdgCList
declared list of neutrino pdg-codes that can be thrown by current instance
Definition: GAstroFlux.h:176
int fgPdgC
(current) generated nu pdg-code
Definition: GAstroFlux.h:177
virtual const TLorentzVector & Position(void)
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
Definition: GAstroFlux.h:143
bool SelectOrigin(bool weighted, TH2D &opdf, double &phi, double &costheta, double &wght)
Definition: GAstroFlux.cxx:458
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAstroFlux.cxx:169
map< int, double > fPntSrcDec
declination
Definition: GAstroFlux.h:264
map< int, double > fPntSrcRA
right ascension
Definition: GAstroFlux.h:263
const int kAstroNCosThetaBins
Definition: GAstroFlux.h:122
const int kAstroNlog10EvBins
Definition: GAstroFlux.h:121
TLorentzVector fgX4
(current) generated nu 4-position
Definition: GAstroFlux.h:179
A base class for the concrete astrophysical neutrino flux drivers.
Definition: GAstroFlux.h:129
NuPropagator * fNuPropg
Definition: GAstroFlux.h:197
static const double GeV
Definition: Units.h:29
const double kAstroDefMinEv
Definition: GAstroFlux.h:120
virtual bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GAstroFlux.h:144
virtual int PdgCode(void)
returns the flux neutrino pdg code
Definition: GAstroFlux.h:140
void SetRelNuPopulations(double nnue=1, double nnumu=2, double nnutau=0, double nnuebar=1, double nnumubar=2, double nnutaubar=0)
Definition: GAstroFlux.cxx:220
double fgWeight
(current) generated nu weight
Definition: GAstroFlux.h:180
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
enum BeamMode string