39 using namespace genie;
62 this->ResetSelection();
64 if(!fEnergySpectrum) {
67 if(!fSolidAngleAcceptance) {
70 if(fRelNuPopulations.size() == 0) {
79 double log10Emin = TMath::Log10(TMath::Max(
kAstroDefMinEv,fMinEvCut));
80 double log10Emax = TMath::Log10(TMath::Min(
kAstroDefMaxEv,fMaxEvCut));
82 double wght_species = 1.;
83 double wght_energy = 1.;
84 double wght_origin = 1.;
87 double log10E = -99999;
89 double costheta = -999999;
93 status = fNuGen->SelectNuPdg(
94 fGenWeighted, fRelNuPopulations, nupdg, wght_species);
99 status = fNuGen->SelectEnergy(
100 fGenWeighted, *fEnergySpectrum, log10Emin, log10Emax, log10E, wght_energy);
104 double Ev = TMath::Power(10.,log10E);
106 status = fNuGen->SelectOrigin(
107 fGenWeighted, *fSolidAngleAcceptance, phi, costheta, wght_origin);
117 status = fNuPropg->Go(phi, costheta, fDetCenter, fDetSize, nupdg, Ev);
122 int pnupdg = fNuPropg->NuPdgAtDetVolBoundary();
123 TVector3 & px3 = fNuPropg->X3AtDetVolBoundary();
124 TVector3 & pp3 = fNuPropg->P3AtDetVolBoundary();
130 px3 = fRotGEF2THz * px3;
131 pp3 = fRotGEF2THz * pp3;
134 px3 = fRotTHz2User * px3;
135 pp3 = fRotTHz2User * pp3;
140 fgWeight = wght_species * wght_energy * wght_origin;
145 fgP4.SetE(pp3.Mag());
166 LOG(
"Flux",
pERROR) <<
"No clear method implemented for option:"<<
opt;
171 fGenWeighted = gen_weighted;
175 double latitude,
double longitude,
double depth,
double size)
181 assert(longitude >= 0. && longitude < 2.*
kPi );
184 fDetGeoLatitude = latitude;
185 fDetGeoLongitude = longitude;
186 fDetGeoDepth = depth;
195 double radius = REarth - fDetGeoDepth;
197 double theta =
kPi/2. - fDetGeoLatitude;
198 double phi = fDetGeoLongitude;
200 double sintheta = TMath::Sin(theta);
201 double costheta = TMath::Cos(theta);
202 double sinphi = TMath::Sin(phi);
203 double cosphi = TMath::Cos(phi);
205 double xdc = radius * sintheta * cosphi;
206 double ydc = radius * sintheta * sinphi;
207 double zdc = radius * costheta;
209 fDetCenter.SetXYZ(xdc,ydc,zdc);
221 double nnue,
double nnumu,
double nnutau,
222 double nnuebar,
double nnumubar,
double nnutaubar)
224 fRelNuPopulations.clear();
228 fRelNuPopulations.insert(
229 map<int,double>::value_type(
kPdgNuE, nnue));
233 fRelNuPopulations.insert(
234 map<int,double>::value_type(
kPdgNuMu, nnumu));
238 fRelNuPopulations.insert(
239 map<int,double>::value_type(
kPdgNuTau, nnutau));
243 fRelNuPopulations.insert(
244 map<int,double>::value_type(
kPdgAntiNuE, nnuebar));
248 fRelNuPopulations.insert(
253 fRelNuPopulations.insert(
258 double tot = nnue + nnumu + nnutau + nnuebar + nnumubar + nnutaubar;
262 map<int,double>::iterator iter = fRelNuPopulations.begin();
263 for ( ; iter != fRelNuPopulations.end(); ++iter) {
264 double fraction = iter->second;
265 double norm_fraction = fraction/
tot;
266 fRelNuPopulations[iter->first] = norm_fraction;
272 if(fEnergySpectrum)
delete fEnergySpectrum;
279 fEnergySpectrum->SetDirectory(0);
282 double log10E = fEnergySpectrum->GetBinCenter(
i);
283 double Ev = TMath::Power(10., log10E);
284 double flux = TMath::Power(Ev, -1*n);
285 fEnergySpectrum->SetBinContent(
i,flux);
289 double max = fEnergySpectrum->GetMaximum();
290 fEnergySpectrum->Scale(1./max);
295 fRotTHz2User = rotation;
300 LOG(
"Flux",
pNOTICE) <<
"Initializing flux driver";
302 bool allow_dup =
false;
313 fDetGeoLatitude = -1.;
314 fDetGeoLongitude = -1.;
317 fDetCenter.SetXYZ(0,0,0);
321 fSolidAngleAcceptance = 0;
331 fRelNuPopulations.clear();
334 fRotGEF2THz.SetToIdentity();
335 fRotTHz2User.SetToIdentity();
342 this->ResetSelection();
350 fgP4.SetPxPyPzE (0.,0.,0.,0.);
351 fgX4.SetXYZT (0.,0.,0.,0.);
358 fRelNuPopulations.clear();
362 if(fEnergySpectrum)
delete fEnergySpectrum;
363 if(fSolidAngleAcceptance)
delete fSolidAngleAcceptance;
376 bool weighted,
const map<int,double> & nupdgpdf,
377 int & nupdg,
double & wght)
384 if(nupdgpdf.size() == 0) {
393 unsigned int nnu = nupdgpdf.size();
394 unsigned int inu = rnd->
RndFlux().Integer(nnu);
395 map<int,double>::const_iterator iter = nupdgpdf.begin();
404 double xrnd = rnd->
RndFlux().Uniform();
405 map<int,double>::const_iterator iter = nupdgpdf.begin();
406 for( ; iter != nupdgpdf.end(); ++iter) {
407 xsum += iter->second;
424 bool weighted, TH1D & log10Epdf,
double log10Emin,
double log10Emax,
425 double & log10E,
double & wght)
433 if(log10Emax <= log10Emin) {
441 log10E = log10Emin + (log10Emax-log10Emin) * rnd->
RndFlux().Rndm();
442 wght = log10Epdf.GetBinContent(log10Epdf.FindBin(log10E));
449 log10E = log10Epdf.GetRandom();
451 while(log10E < log10Emin || log10E > log10Emax);
459 bool weighted, TH2D & opdf,
460 double & phi,
double & costheta,
double & wght)
471 costheta = -1. + 2.*rnd->
RndFlux().Rndm();
472 wght = opdf.GetBinContent(opdf.FindBin(phi,costheta));
478 opdf.GetRandom2(phi,costheta);
486 double phi,
double costheta,
const TVector3 & detector_centre,
487 double detector_sz,
int nu_pdg,
double Ev)
495 double sintheta =
TMath::Sqrt(1-costheta*costheta);
496 double cosphi = TMath::Cos(phi);
497 double sinphi = TMath::Sin(phi);
499 double zs = REarth * costheta;
500 double ys = REarth * sintheta * cosphi;
501 double xs = REarth * sintheta * sinphi;
503 TVector3 start_position(xs,ys,zs);
504 fX3 = start_position - detector_centre;
509 TVector3 direction_unit_vec = -1. * fX3.Unit();
510 fP3 = Ev * direction_unit_vec;
516 LOG(
"Flux",
pWARN) <<
"|dist| = " << fX3.Mag();
517 LOG(
"Flux",
pWARN) <<
"|detsize| = " << detector_sz;
520 double currdist = fX3.Mag();
521 if(currdist <= detector_sz - 0.1)
break;
523 double stepsz = (currdist-detector_sz > fStepSize) ?
524 fStepSize : currdist-detector_sz;
525 if(stepsz <= 0.)
break;
536 fX3 += (stepsz * direction_unit_vec);
587 string name,
double ra,
double dec,
double rel_intensity)
590 (ra >= 0. && ra < 2.*
kPi) &&
591 (dec >= -
kPi/2. && dec <=
kPi/2.) &&
592 (rel_intensity > 0) &&
598 fPntSrcName.insert( map<int, string>::value_type(
id, name ) );
599 fPntSrcRA. insert( map<int, double>::value_type(
id, ra ) );
600 fPntSrcDec. insert( map<int, double>::value_type(
id, dec ) );
601 fPntSrcRelI.insert( map<int, double>::value_type(
id, rel_intensity) );
617 unsigned int srcid = 0;
625 srcid = rnd->
RndFlux().Integer(nsrc);
634 map<int,double>::const_iterator piter =
fPntSrcRelI.begin();
636 xsum += piter->second;
638 srcid = piter->first;
void ForceMinEnergy(double emin)
virtual void Clear(Option_t *opt)
reset state variables based on opt
bool fGenWeighted
(config) generate a weighted or unweighted flux?
THE MAIN GENIE PROJECT NAMESPACE
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
static RandomGen * Instance()
Access instance.
void SetEnergyPowLawIdx(double n)
vector< vector< double > > clear
void SetUserCoordSystem(TRotation &rotation)
rotation Topocentric Horizontal -> User-defined Topocentric Coord System
const double kAstroDefMaxEv
void ForceMaxEnergy(double emax)
map< int, double > fPntSrcRelI
relative intensity
A singleton holding random number generator classes. All random number generation in GENIE should tak...
map< int, string > fPntSrcName
point source name
void Initialize(Bool_t useTMVAStyle=kTRUE)
virtual bool GenerateNext(void)
generate the next flux neutrino (return false in err)
bool SelectNuPdg(bool weighted, const map< int, double > &nupdgpdf, int &nupdg, double &wght)
double fPntSrcTotI
sum of all relative intensities
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void ResetSelection(void)
bool SelectEnergy(bool weighted, TH1D &log10epdf, double log10emin, double log10emax, double &log10e, double &wght)
void SetDetectorPosition(double latitude, double longitude, double depth, double size)
void AddPointSource(string name, double ra, double dec, double rel_intensity)
unsigned int fSelSourceId
static const double kREarth
bool SelectOrigin(bool weighted, TH2D &opdf, double &phi, double &costheta, double &wght)
bool Go(double phi_start, double costheta_start, const TVector3 &detector_centre, double detector_sz, int nu_pdg, double Ev)
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Var Sqrt(const Var &v)
Use to take sqrt of a var.
map< int, double > fPntSrcDec
declination
map< int, double > fPntSrcRA
right ascension
assert(nhit_max >=nhit_nbins)
const int kAstroNlog10EvBins
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
A base class for the concrete astrophysical neutrino flux drivers.
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const double kAstroDefMinEv
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
void SetRelNuPopulations(double nnue=1, double nnumu=2, double nnutau=0, double nnuebar=1, double nnumubar=2, double nnutaubar=0)
double fgWeight
(current) generated nu weight