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