GFluxGenerator.cxx
Go to the documentation of this file.
1 // ROOT includes
2 #include "TH2.h"
3 
4 // GENIE includes
5 #include "GENIE/Framework/Messenger/Messenger.h"
6 #include "GENIE/Framework/Conventions/Units.h"
7 
8 // NOvASoft includes
10 
11 using namespace genie::supernova;
12 
13 
14 //.......................................................................
15 GFluxGenerator::GFluxGenerator(int pdg, TH2* hist) : fPdg(pdg), fT(0), fE(0)
16 {
17  // Set the limits
18  fLimsT.Set(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
19  fLimsE.Set(hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
20 
21  // Set generator integral
22  fIntegral = hist->Integral("width");
23 
24  // Init the time sampler
25  TH1D* histT = hist->ProjectionX(Form("HistT_%d", pdg));
26  fTSampler.Init(histT);
27  int NBins = hist->GetNbinsX();
28 
29  // Init the energy sampler (per each time bin)
30  fESamplers.resize(NBins);
31  for (int nbin = 0; nbin < NBins; ++nbin) {
32  TH1D* h = hist->ProjectionY(Form("HistE_%d_%d", pdg, nbin), nbin, nbin);
33  fESamplers[nbin].Init(h);
34  }
35 }
36 
37 
38 //.......................................................................
40 {
41  fE = 0;
42  do {
43  // Generate next event time
45  if (this->End()) return;
46 
47  // Get current time bin
48  int bin = fTSampler.GetBin();
49  LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
50  << "] T=" << fT
51  << " (bin=" << bin <<")";
52 
53  // Generate event energy using current bin distribution
54  fE = fESamplers[bin].GenerateNext();
55  LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
56  << "] E=" << fE
57  << " vs Emin=" << fEmin;
58 
59  } while (fE <= fEmin);
60 
62 
63  LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
64  << "]: generated T=" << fT
65  << "; E=" << fE;
66 }
67 
68 
69 //.......................................................................
71 {
72  particle.SetPDG(fPdg);
73  particle.SetE(fE);
74  particle.SetT(fT);
75 }
76 
77 
78 //.......................................................................
80 {
81  fTSampler.Reload();
82 }
83 
84 
85 //.......................................................................
87 {
88  fTSampler.SetProbScale(scale);
89 }
90 
91 
92 //.......................................................................
93 void GFluxGenerator::SetEmin(double Emin)
94 {
95  fEmin = Emin;
96 }
97 
98 
99 //.......................................................................
101 {
102  return fTSampler.End();
103 }
104 
105 
106 //.......................................................................
108 {
109  return fPdg;
110 }
111 
112 
113 //.......................................................................
114 double GFluxGenerator::T() const
115 {
116  return fT;
117 }
118 
119 
120 //.......................................................................
121 double GFluxGenerator::E() const
122 {
123  return fE;
124 }
125 
126 
127 //.......................................................................
129 {
130  return fTSampler;
131 }
132 
133 
134 //.......................................................................
136 {
137  return fESamplers[i];
138 }
const OrderedSampler & GetSamplerT() const
const RandomSampler & GetSamplerE(size_t i) const
void Set(double val1, double val2)
Definition: Limits.h:19
void SetProbScale(double scale)
Definition: Sampler.cxx:121
static const double MeV
Definition: Units.h:130
GFluxGenerator(int pdg, TH2 *hist)
Double_t scale
Definition: plot.C:25
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Limits< double > fLimsE
Total flux integral.
Definition: GeneratorBase.h:22
float bin[41]
Definition: plottest35.C:14
Limits< double > fLimsT
Energy range.
Definition: GeneratorBase.h:23
std::vector< RandomSampler > fESamplers
void FillParticle(GenParticle &particle)
void Init(TH1 *hist, bool integral=false)
Definition: Sampler.cxx:52
#define pDEBUG
Definition: Messenger.h:64