GAstroFlux.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Mar 27, 2010 - CA
14  This class was first added in 2.7.1.
15  @ Feb 22, 2011 - JD
16  Implemented dummy versions of the new GFluxI::Clear and GFluxI::Index as
17  these methods needed for pre-generation of flux interaction probabilities
18  in GMCJDriver.
19 
20 */
21 //____________________________________________________________________________
22 
23 #include <cassert>
24 
25 #include <TH1D.h>
26 #include <TH2D.h>
27 #include <TMath.h>
28 
30 #include "Tools/Flux/GAstroFlux.h"
38 
39 using namespace genie;
40 using namespace genie::flux;
41 using namespace genie::constants;
42 
43 //____________________________________________________________________________
45 {
46  this->Initialize();
47 }
48 //___________________________________________________________________________
50 {
51  this->CleanUp();
52 }
53 //___________________________________________________________________________
55 {
56  return TMath::Min(kAstroDefMaxEv, fMaxEvCut);
57 }
58 //___________________________________________________________________________
60 {
61  // Reset previously generated neutrino code / 4-p / 4-x
62  this->ResetSelection();
63 
64  if(!fEnergySpectrum) {
65  return false;
66  }
67  if(!fSolidAngleAcceptance) {
68  return false;
69  }
70  if(fRelNuPopulations.size() == 0) {
71  return false;
72  }
73 
74  //
75  // Generate neutrino energy & starting position at the Geocentric
76  // coordinate system
77  //
78 
79  double log10Emin = TMath::Log10(TMath::Max(kAstroDefMinEv,fMinEvCut));
80  double log10Emax = TMath::Log10(TMath::Min(kAstroDefMaxEv,fMaxEvCut));
81 
82  double wght_species = 1.;
83  double wght_energy = 1.;
84  double wght_origin = 1.;
85 
86  int nupdg = 0;
87  double log10E = -99999;
88  double phi = -999999;
89  double costheta = -999999;
90 
91  bool status = true;
92 
93  status = fNuGen->SelectNuPdg(
94  fGenWeighted, fRelNuPopulations, nupdg, wght_species);
95  if(!status) {
96  return false;
97  }
98 
99  status = fNuGen->SelectEnergy(
100  fGenWeighted, *fEnergySpectrum, log10Emin, log10Emax, log10E, wght_energy);
101  if(!status) {
102  return false;
103  }
104  double Ev = TMath::Power(10.,log10E);
105 
106  status = fNuGen->SelectOrigin(
107  fGenWeighted, *fSolidAngleAcceptance, phi, costheta, wght_origin);
108  if(!status) {
109  return false;
110  }
111 
112  //
113  // Propagate through the Earth: Get position, 4-momentum and neutrino
114  // pdg code at the boundary of the detector volume
115  //
116 
117  status = fNuPropg->Go(phi, costheta, fDetCenter, fDetSize, nupdg, Ev);
118  if(!status) {
119  return false;
120  }
121 
122  int pnupdg = fNuPropg->NuPdgAtDetVolBoundary();
123  TVector3 & px3 = fNuPropg->X3AtDetVolBoundary();
124  TVector3 & pp3 = fNuPropg->P3AtDetVolBoundary();
125 
126  //
127  // Rotate vectors:
128 
129  // GEF translated to detector centre -> THZ
130  px3 = fRotGEF2THz * px3;
131  pp3 = fRotGEF2THz * pp3;
132 
133  // THZ -> Topocentric user-defined detetor system
134  px3 = fRotTHz2User * px3;
135  pp3 = fRotTHz2User * pp3;
136 
137  //
138  // Set position, momentum, pdg code and weight variables reported back
139  //
140  fgWeight = wght_species * wght_energy * wght_origin;
141  fgPdgC = pnupdg;
142  fgX4.SetVect(px3*(units::m/units::km));
143  fgX4.SetT(0.);
144  fgP4.SetVect(pp3);
145  fgP4.SetE(pp3.Mag());
146 
147  return true;
148 }
149 //___________________________________________________________________________
151 {
152  emin = TMath::Max(0., emin/units::GeV);
153  fMinEvCut = emin;
154 }
155 //___________________________________________________________________________
157 {
158  emax = TMath::Max(0., emax/units::GeV);
159  fMaxEvCut = emax;
160 }
161 //___________________________________________________________________________
162 void GAstroFlux::Clear(Option_t * opt)
163 {
164 // Dummy clear method needed to conform to GFluxI interface
165 //
166  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
167 }
168 //___________________________________________________________________________
169 void GAstroFlux::GenerateWeighted(bool gen_weighted)
170 {
171  fGenWeighted = gen_weighted;
172 }
173 //___________________________________________________________________________
175  double latitude, double longitude, double depth, double size)
176 {
177  depth = TMath::Max(0., depth/units::km);
178  size = TMath::Max(0., size /units::km);
179 
180  assert(latitude >= -kPi/2. && latitude <= kPi/2.);
181  assert(longitude >= 0. && longitude < 2.*kPi );
182 
183  // set inputs
184  fDetGeoLatitude = latitude;
185  fDetGeoLongitude = longitude;
186  fDetGeoDepth = depth;
187  fDetSize = size;
188 
189  //
190  // Compute detector/topocentric coordinate system center in the
191  // geocentric coordinate system.
192  //
193 
194  double REarth = constants::kREarth/units::km;
195  double radius = REarth - fDetGeoDepth;
196 
197  double theta = kPi/2. - fDetGeoLatitude;
198  double phi = fDetGeoLongitude;
199 
200  double sintheta = TMath::Sin(theta);
201  double costheta = TMath::Cos(theta);
202  double sinphi = TMath::Sin(phi);
203  double cosphi = TMath::Cos(phi);
204 
205  double xdc = radius * sintheta * cosphi;
206  double ydc = radius * sintheta * sinphi;
207  double zdc = radius * costheta;
208 
209  fDetCenter.SetXYZ(xdc,ydc,zdc);
210 
211  //
212  // Coordinate System Rotation:
213  // GEF translated to detector centre -> THZ
214  //
215  // ...
216  // ... TODO
217  // ...
218 }
219 //___________________________________________________________________________
221  double nnue, double nnumu, double nnutau,
222  double nnuebar, double nnumubar, double nnutaubar)
223 {
224  fRelNuPopulations.clear();
225  fPdgCList->clear();
226 
227  if(nnue>0.) {
228  fRelNuPopulations.insert(
229  map<int,double>::value_type(kPdgNuE, nnue));
230  fPdgCList->push_back(kPdgNuE);
231  }
232  if(nnumu>0.) {
233  fRelNuPopulations.insert(
234  map<int,double>::value_type(kPdgNuMu, nnumu));
235  fPdgCList->push_back(kPdgNuMu);
236  }
237  if(nnutau>0.) {
238  fRelNuPopulations.insert(
239  map<int,double>::value_type(kPdgNuTau, nnutau));
240  fPdgCList->push_back(kPdgNuTau);
241  }
242  if(nnuebar>0.) {
243  fRelNuPopulations.insert(
244  map<int,double>::value_type(kPdgAntiNuE, nnuebar));
245  fPdgCList->push_back(kPdgAntiNuE);
246  }
247  if(nnumubar>0.) {
248  fRelNuPopulations.insert(
249  map<int,double>::value_type(kPdgAntiNuMu, nnumubar));
250  fPdgCList->push_back(kPdgAntiNuMu);
251  }
252  if(nnutaubar>0.) {
253  fRelNuPopulations.insert(
254  map<int,double>::value_type(kPdgAntiNuTau, nnutaubar));
255  fPdgCList->push_back(kPdgAntiNuTau);
256  }
257 
258  double tot = nnue + nnumu + nnutau + nnuebar + nnumubar + nnutaubar;
259  assert(tot>0.);
260 
261  // normalize to 1.
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;
267  }
268 }
269 //___________________________________________________________________________
271 {
272  if(fEnergySpectrum) delete fEnergySpectrum;
273 
274  double log10Emin = TMath::Log10(kAstroDefMinEv);
275  double log10Emax = TMath::Log10(kAstroDefMaxEv);
276 
277  fEnergySpectrum =
278  new TH1D("fEnergySpectrum","",kAstroNlog10EvBins,log10Emin,log10Emax);
279  fEnergySpectrum->SetDirectory(0);
280 
281  for(int i=1; i<=kAstroNlog10EvBins; i++) {
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);
286  }
287 
288  // normalize
289  double max = fEnergySpectrum->GetMaximum();
290  fEnergySpectrum->Scale(1./max);
291 }
292 //___________________________________________________________________________
293 void GAstroFlux::SetUserCoordSystem(TRotation & rotation)
294 {
295  fRotTHz2User = rotation;
296 }
297 //___________________________________________________________________________
299 {
300  LOG("Flux", pNOTICE) << "Initializing flux driver";
301 
302  bool allow_dup = false;
303  fPdgCList = new PDGCodeList(allow_dup);
304 
305  // Default: No min/max energy cut
306  this->ForceMinEnergy(kAstroDefMinEv);
307  this->ForceMaxEnergy(kAstroDefMaxEv);
308 
309  // Generate weighted or un-weighted flux?
310  fGenWeighted = true;
311 
312  // Detector position & size
313  fDetGeoLatitude = -1.;
314  fDetGeoLongitude = -1.;
315  fDetGeoDepth = -1.;
316  fDetSize = -1.;
317  fDetCenter.SetXYZ(0,0,0); // in the geocentric coord system
318 
319  // Normalized 2-D histogram (phi,costheta): detector solid angle
320  // as seen from different positions across the face of the Earth
321  fSolidAngleAcceptance = 0;
322 
323  // Neutrino energy spectrum
324  // To be set via SetEnergyPowLawIdx()
325  // Can be trivially modified to accomodate different spectra
326  fEnergySpectrum = 0;
327 
328  // Relative neutrino populations
329  // To be set via SetRelNuPopulations()
330  // Default nue:numu:nutau:nuebar:numubar:nutaubar = 1:2:0:1:2:0
331  fRelNuPopulations.clear();
332 
333  // Rotations
334  fRotGEF2THz.SetToIdentity();
335  fRotTHz2User.SetToIdentity();
336 
337  // Utility objects for generating and propagating neutrinos
338  fNuGen = new NuGenerator();
339  fNuPropg = new NuPropagator(1.0*units::km);
340 
341  // Reset `current' selected flux neutrino
342  this->ResetSelection();
343 }
344 //___________________________________________________________________________
346 {
347 // initializing running neutrino pdg-code, 4-position, 4-momentum
348 
349  fgPdgC = 0;
350  fgP4.SetPxPyPzE (0.,0.,0.,0.);
351  fgX4.SetXYZT (0.,0.,0.,0.);
352 }
353 //___________________________________________________________________________
355 {
356  LOG("Flux", pNOTICE) << "Cleaning up...";
357 
358  fRelNuPopulations.clear();
359  fPdgCList->clear();
360 
361  delete fPdgCList;
362  if(fEnergySpectrum) delete fEnergySpectrum;
363  if(fSolidAngleAcceptance) delete fSolidAngleAcceptance;
364 
365  delete fNuGen;
366  delete fNuPropg;
367 }
368 //___________________________________________________________________________
369 
370 //
371 // GAstroFlux utility class method definitions
372 // ..........................................................................
373 //
374 //___________________________________________________________________________
376  bool weighted, const map<int,double> & nupdgpdf,
377  int & nupdg, double & wght)
378 {
379 // select neutrino species based on relative neutrino species populations
380 //
381  nupdg = 0;
382  wght = 0;
383 
384  if(nupdgpdf.size() == 0) {
385  return false;
386  }
387 
389 
390  // Generate weighted flux:
391  //
392  if(weighted) {
393  unsigned int nnu = nupdgpdf.size();
394  unsigned int inu = rnd->RndFlux().Integer(nnu);
395  map<int,double>::const_iterator iter = nupdgpdf.begin();
396  advance(iter,inu);
397  nupdg = iter->first;
398  wght = iter->second;
399  }
400  // Generate un-weighted flux:
401  //
402  else {
403  double xsum = 0.;
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;
408  if(xrnd < xsum) {
409  nupdg = iter->first;
410  break;
411  }
412  }
413  wght = 1.;
414  }
415 
416  if(nupdg==0) {
417  return false;
418  }
419 
420  return true;
421 }
422 //___________________________________________________________________________
424  bool weighted, TH1D & log10Epdf, double log10Emin, double log10Emax,
425  double & log10E, double & wght)
426 {
427 // select neutrino energy
428 //
429 
430  log10E = -9999999;
431  wght = 0;
432 
433  if(log10Emax <= log10Emin) {
434  return false;
435  }
436 
437  // Generate weighted flux:
438  //
439  if(weighted) {
441  log10E = log10Emin + (log10Emax-log10Emin) * rnd->RndFlux().Rndm();
442  wght = log10Epdf.GetBinContent(log10Epdf.FindBin(log10E));
443  }
444 
445  // Generate un-weighted flux:
446  //
447  else {
448  do {
449  log10E = log10Epdf.GetRandom();
450  }
451  while(log10E < log10Emin || log10E > log10Emax);
452  wght = 1.;
453  }
454 
455  return true;
456 }
457 //___________________________________________________________________________
459  bool weighted, TH2D & opdf,
460  double & phi, double & costheta, double & wght)
461 {
462  wght = 0;
463  costheta = -999999;
464  phi = -999999;
465 
466  // Generate weighted flux:
467  //
468  if(weighted) {
470  phi = 2.*kPi * rnd->RndFlux().Rndm();
471  costheta = -1. + 2.*rnd->RndFlux().Rndm();
472  wght = opdf.GetBinContent(opdf.FindBin(phi,costheta));
473  }
474 
475  // Generate un-weighted flux:
476  //
477  else {
478  opdf.GetRandom2(phi,costheta);
479  wght = 1.;
480  }
481 
482  return true;
483 }
484 //___________________________________________________________________________
486  double phi, double costheta, const TVector3 & detector_centre,
487  double detector_sz, int nu_pdg, double Ev)
488 {
489  // initialize neutrino code
490  fNuPdg = nu_pdg;
491 
492  //
493  // initialize neutrino position vector
494  //
495  double sintheta = TMath::Sqrt(1-costheta*costheta);
496  double cosphi = TMath::Cos(phi);
497  double sinphi = TMath::Sin(phi);
498  double REarth = constants::kREarth/units::km;
499  double zs = REarth * costheta;
500  double ys = REarth * sintheta * cosphi;
501  double xs = REarth * sintheta * sinphi;
502 
503  TVector3 start_position(xs,ys,zs);
504  fX3 = start_position - detector_centre;
505 
506  //
507  // initialize neutrino momentum 4-vector
508  //
509  TVector3 direction_unit_vec = -1. * fX3.Unit();
510  fP3 = Ev * direction_unit_vec;
511 
512  //
513  // step through the Earth
514  //
515 
516  LOG("Flux", pWARN) << "|dist| = " << fX3.Mag();
517  LOG("Flux", pWARN) << "|detsize| = " << detector_sz;
518 
519  while(1) {
520  double currdist = fX3.Mag();
521  if(currdist <= detector_sz - 0.1) break;
522 
523  double stepsz = (currdist-detector_sz > fStepSize) ?
524  fStepSize : currdist-detector_sz;
525  if(stepsz <= 0.) break;
526 
527 // LOG("Flux", pWARN) << "Stepping by... |dr| = " << stepsz;
528 
529  //
530  // check earth density at current position, calculate interaction
531  // probability over next step size, decide whether it interacts and
532  // what happens if it does...
533  //
534  // ... todo ...
535 
536  fX3 += (stepsz * direction_unit_vec);
537  }
538 
539  return true;
540 }
541 //___________________________________________________________________________
542 
543 //
544 // GDiffuseAstroFlux concrete flux driver
545 // ..........................................................................
546 //
547 //___________________________________________________________________________
549 GAstroFlux()
550 {
551 
552 }
553 //___________________________________________________________________________
555 {
556 
557 }
558 //___________________________________________________________________________
559 
560 //
561 // GPointSourceAstroFlux concrete flux driver
562 // ..........................................................................
563 //
564 //___________________________________________________________________________
566 GAstroFlux()
567 {
568  fPntSrcName.clear();
569  fPntSrcRA. clear();
570  fPntSrcDec. clear();
571  fPntSrcRelI.clear();
572 
573  fPntSrcTotI = 0;
574 }
575 //___________________________________________________________________________
577 {
578 
579 }
580 //___________________________________________________________________________
582 {
583  return true;
584 }
585 //___________________________________________________________________________
587  string name, double ra, double dec, double rel_intensity)
588 {
589  bool accept =
590  (ra >= 0. && ra < 2.*kPi) &&
591  (dec >= -kPi/2. && dec <= kPi/2.) &&
592  (rel_intensity > 0) &&
593  (name.size() > 0);
594 
595  if(accept) {
596  int id = fPntSrcName.size();
597 
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) );
602 
603  fPntSrcTotI += rel_intensity;
604  }
605 }
606 //___________________________________________________________________________
608 {
609  if(fPntSrcRelI.size() == 0) {
610  return false;
611  }
612 
613  if(fPntSrcTotI <= 0.) {
614  return false;
615  }
616 
617  unsigned int srcid = 0;
618  double wght = 0;
619 
620  // Generate weighted flux:
621  //
622  if(fGenWeighted) {
624  unsigned int nsrc = fPntSrcName.size();
625  srcid = rnd->RndFlux().Integer(nsrc);
626  wght = fPntSrcRelI[srcid] / fPntSrcTotI;
627  }
628  // Generate un-weighted flux:
629  //
630  else {
632  double xsum = 0.;
633  double xrnd = fPntSrcTotI * rnd->RndFlux().Uniform();
634  map<int,double>::const_iterator piter = fPntSrcRelI.begin();
635  for( ; piter != fPntSrcRelI.end(); ++piter) {
636  xsum += piter->second;
637  if(xrnd < xsum) {
638  srcid = piter->first;
639  break;
640  }
641  }
642  wght = 1.;
643  }
644 
645  //
646  fSelSourceId = srcid;
647 
648  //
649  fgWeight *= wght;
650 
651  return true;
652 }
653 //___________________________________________________________________________
654 
const double kPi
void ForceMinEnergy(double emin)
Definition: GAstroFlux.cxx:150
virtual void Clear(Option_t *opt)
reset state variables based on opt
Definition: GAstroFlux.cxx:162
const XML_Char * name
Definition: expat.h:151
static const double m
Definition: Units.h:79
Basic constants.
const int kPdgNuE
Definition: PDGCodes.h:28
int status
Definition: fabricate.py:1613
bool fGenWeighted
(config) generate a weighted or unweighted flux?
Definition: GAstroFlux.h:184
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAstroFlux.cxx:54
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const int kPdgAntiNuE
Definition: PDGCodes.h:29
void SetEnergyPowLawIdx(double n)
Definition: GAstroFlux.cxx:270
vector< vector< double > > clear
const int kPdgNuMu
Definition: PDGCodes.h:30
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 > fPntSrcRelI
relative intensity
Definition: GAstroFlux.h:265
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
map< int, string > fPntSrcName
point source name
Definition: GAstroFlux.h:262
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
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
const double emin
GENIE flux drivers.
static const double km
Definition: Units.h:72
double fPntSrcTotI
sum of all relative intensities
Definition: GAstroFlux.h:266
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const double emax
bool SelectEnergy(bool weighted, TH1D &log10epdf, double log10emin, double log10emax, double &log10e, double &wght)
Definition: GAstroFlux.cxx:423
Double_t radius
void SetDetectorPosition(double latitude, double longitude, double depth, double size)
Definition: GAstroFlux.cxx:174
void AddPointSource(string name, double ra, double dec, double rel_intensity)
Definition: GAstroFlux.cxx:586
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:61
const int kPdgNuTau
Definition: PDGCodes.h:32
static const double kREarth
Definition: Constants.h:111
bool SelectOrigin(bool weighted, TH2D &opdf, double &phi, double &costheta, double &wght)
Definition: GAstroFlux.cxx:458
bool Go(double phi_start, double costheta_start, const TVector3 &detector_centre, double detector_sz, int nu_pdg, double Ev)
Definition: GAstroFlux.cxx:485
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAstroFlux.cxx:169
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
map< int, double > fPntSrcDec
declination
Definition: GAstroFlux.h:264
map< int, double > fPntSrcRA
right ascension
Definition: GAstroFlux.h:263
assert(nhit_max >=nhit_nbins)
const int kAstroNlog10EvBins
Definition: GAstroFlux.h:121
#define pNOTICE
Definition: Messenger.h:62
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
A base class for the concrete astrophysical neutrino flux drivers.
Definition: GAstroFlux.h:129
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GAstroFlux.cxx:581
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static const double GeV
Definition: Units.h:29
const double kAstroDefMinEv
Definition: GAstroFlux.h:120
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)
Definition: GAstroFlux.cxx:220
double fgWeight
(current) generated nu weight
Definition: GAstroFlux.h:180