GSNovaModel.cxx
Go to the documentation of this file.
1 // C++ includes
2 #include <algorithm>
3 #include <exception>
4 
5 // ROOT includes
6 #include "TFile.h"
7 #include "TH2.h"
8 
9 // GENIE includes
10 #include "GENIE/Framework/Messenger/Messenger.h"
12 
13 using namespace genie;
14 using namespace genie::supernova;
15 
16 
17 //.......................................................................
19 {
20  this->ReadModel(filename);
21 }
22 
23 
24 //.......................................................................
25 TH2* ScaleHist(TH2* hist)
26 {
27  // The flux is normalized to a fluence through 1cm2 surface at 1kPc distance
28  const double kkpcPerCm = 3.086e21; /// kPc to cm
29  const double kLumiMag = 1e50; /// Luminosity magnitude
30  const double kSphereSolidAngle = 4 * TMath::Pi(); /// in str
31  const double kWindowArea = 1; /// in cm2
32 
33  const double norm = (kLumiMag / (kkpcPerCm * kkpcPerCm)) * kWindowArea / kSphereSolidAngle;
34 
35  hist->Scale(norm);
36 
37  return hist;
38 }
39 
40 
41 //.......................................................................
43 {
44  TFile f(filename.data(), "READ");
45  if (!f.IsOpen()) {
46  throw std::runtime_error("Error opening SN model file " + filename);
47  }
48 
49  std::map<int, std::string> histMap = {
50  { 12, "fluence_e" },
51  {-12, "fluence_ae"},
52  { 14, "fluence_x" },
53  {-14, "fluence_ax"}
54  };
55 
56  int histPDG;
58  for(auto& histMapEntry: histMap){
59  histPDG = histMapEntry.first;
60  histName = histMapEntry.second;
61 
62  TH2* hFlux = (TH2*)f.Get(histName.data());
63  if (!hFlux) {
64  throw std::runtime_error("Failed to read histogram \'"
65  + histName
66  + "\' from file \'"
67  + filename + "\'");
68  }
69 
70  this->AddGenerator(GFluxGenerator(histPDG, ScaleHist(hFlux)));
71  }
72 
73  this->SetProbScale(1.0);
74  f.Close();
75 }
76 
77 
78 //.......................................................................
80 {
81  return fPdgList;
82 }
83 
84 
85 //.......................................................................
87 {
88  // Add new generator to the model
89  LOG("GSNovaModel", pNOTICE) << "Add generator for pdg=" << generator.PDG()
90  << ", integral=" << generator.GetIntegral();
91 
92  //update generators list
93  fGeneratorsDead.push_back(generator);
94 
95  //update integral
96  fIntegral += generator.GetIntegral();
97 
98  //update pdg list
99  fPdgList.push_back(generator.PDG());
100 
101  //update limits
102  fLimsE.Update(generator.LimsE());
103  fLimsT.Update(generator.LimsT());
104 }
105 
106 
107 //.......................................................................
109 {
110  // Generate next particle and fill it to the argument
111  if (this->End()) return;
112 
113  // (Randomly) select from which generator we will take the input
114  auto itGenerator = this->SelectGenerator();
115  assert (itGenerator != fGeneratorsActive.end());
116 
117  GFluxGenerator& gen =* itGenerator;
118  LOG("GSNovaModel", pDEBUG) << "Using generator [" << gen.PDG() << "],"
119  << "t=" << gen.T()
120  << ",e=" << gen.E();
121 
122  // Fill the generated calues (T,E) to the part
123  gen.FillParticle(particle);
124 
125  //Generate a new particle in this generator
126  gen.GenerateNext();
127  if (gen.End()) {
128  LOG("GSNovaModel", pNOTICE) << "Finishing generator #" << gen.PDG();
129 
130  // Move generator to Dead list
131  fGeneratorsDead.splice(fGeneratorsDead.begin(), fGeneratorsActive, itGenerator);
132  LOG("GSNovaModel", pNOTICE) << fGeneratorsActive.size() << " left";
133  }
134 
135  LOG("GSNovaModel", pDEBUG) << "Done";
136 }
137 
138 
139 //.......................................................................
141 {
142  return fGeneratorsActive.empty();
143 }
144 
145 
146 //.......................................................................
148 {
149  // Move all generators to active state
151 
152  // Reload all generators
153  for (auto &gen: fGeneratorsActive) {
154  LOG("GSNovaModel", pDEBUG) << "Reloading generator " << gen.PDG();
155  gen.Reload();
156  gen.GenerateNext();
157  }
158 }
159 
160 
161 //.......................................................................
162 void GSNovaModel::SetEmin(double Emin)
163 {
164  for (auto& gen: fGeneratorsActive) gen.SetEmin(Emin);
165 }
166 
167 
168 //.......................................................................
170 {
171  LOG("GSNovaModel", pNOTICE) << "Set Probability scale:" << scale;
172 
173  for (auto& gen: fGeneratorsDead) gen.SetProbScale(scale);
174  for (auto& gen: fGeneratorsActive) gen.SetProbScale(scale);
175 
176  this->Reload();
177 }
178 
179 
180 //.......................................................................
181 bool compare(const GFluxGenerator& g1, const GFluxGenerator& g2)
182 {
183  return g1.T() < g2.T();
184 }
185 
186 
187 //.......................................................................
188 const std::list<GFluxGenerator> GSNovaModel::ActiveGenerators() {
189  return fGeneratorsActive;
190 }
191 
192 
193 //.......................................................................
194 std::list<GFluxGenerator>::iterator GSNovaModel::SelectGenerator()
195 {
196  // For timed distribution: select the earliest particle in time
197  LOG("GSNovaModel", pDEBUG) << "Select from " << fGeneratorsActive.size() << " active generators";
198 
199  return std::min_element(fGeneratorsActive.begin(), fGeneratorsActive.end(), compare);
200 }
bool End()
Check if there are no active generators left.
genie::PDGCodeList fPdgList
Definition: GSNovaModel.h:45
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
const Limits< double > & LimsT() const
std::list< GFluxGenerator > fGeneratorsDead
Definition: GSNovaModel.h:49
string filename
Definition: shutoffs.py:106
std::list< GFluxGenerator >::iterator SelectGenerator()
const Limits< double > & LimsE() const
A list of PDG codes.
Definition: PDGCodeList.h:33
const std::list< GFluxGenerator > ActiveGenerators()
Access to generators.
const double GetIntegral() const
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
std::list< GFluxGenerator > fGeneratorsActive
list of particles
Definition: GSNovaModel.h:48
const genie::PDGCodeList & PdgList()
Get model features.
Definition: GSNovaModel.cxx:79
Limits< double > fLimsE
Total flux integral.
Definition: GeneratorBase.h:22
void GenerateNext(GenParticle &)
Generate next particle.
Limits< double > fLimsT
Energy range.
Definition: GeneratorBase.h:23
void Update(const Limits< T > &other)
Definition: Limits.h:25
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Float_t norm
void Reload()
Reactivate all generators.
TH2 * ScaleHist(TH2 *hist)
Definition: GSNovaModel.cxx:25
bool compare(const GFluxGenerator &g1, const GFluxGenerator &g2)
void AddGenerator(GFluxGenerator)
Populate model with generators.
Definition: GSNovaModel.cxx:86
assert(nhit_max >=nhit_nbins)
void SetProbScale(double Scale)
#define pNOTICE
Definition: Messenger.h:62
void ReadModel(std::string filename)
Definition: GSNovaModel.cxx:42
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
GSNovaModel(std::string filename)
Definition: GSNovaModel.cxx:18
#define pDEBUG
Definition: Messenger.h:64