OscCalcSterile.cxx
Go to the documentation of this file.
1 #include "OscCalcSterile.h"
2 
3 #include <cassert>
4 #include <cstdlib>
5 #include <iostream>
6 #include <iomanip>
7 
8 namespace osc
9 {
11  : fPMNS_Sterile(0), fNFlavors(3), fDirty(true), fPrevE(0), fPrevAnti(0), fPrevFlavBefore(0)
12  {
14  }
15 
16  //---------------------------------------------------------------------------
17  std::vector<double> OscCalcSterile::GetState() const
18  {
19  std::vector<double> state;
20  state.push_back((double)fNFlavors);
21  state.push_back(fL);
22  state.push_back(fRho);
23  for(int i = 2; i <= fNFlavors; ++i) state.push_back(GetDm(i));
24  for(int j = 2; j <= fNFlavors; ++j) {
25  for(int i = 1; i < j; ++i) {
26  state.push_back(GetAngle(i, j));
27  if(i+1 != j) state.push_back(GetDelta(i, j));
28  }
29  }
30  return state;
31  }
32 
33  //---------------------------------------------------------------------------
34  void OscCalcSterile::SetState(std::vector<double> state)
35  {
36  int iState(0);
37  fDirty = true ;
38  SetNFlavors((int)state[iState++]);
39  SetL(state[iState++]);
40  SetRho(state[iState++]);
41  for (int i = 2; i <= fNFlavors; ++i) SetDm(i, state[iState++]);
42  for(int j = 2; j <= fNFlavors; ++j) {
43  for(int i = 1; i < j; ++i) {
44  SetAngle(i, j, state[iState++]);
45  if(i+1 != j) SetDelta(i, j, state[iState++]);
46  }
47  }
48  }
49 
50  //---------------------------------------------------------------------------
52  : OscCalcSterile()
53  {
54  std::vector<double> state = calc.GetState();
55  SetState(state);
56  }
57 
58  //---------------------------------------------------------------------------
60  {
61  TMD5* ret = new TMD5;
62  std::string txt = "PMNSSterile";
63  ret->Update((unsigned char*)txt.c_str(), txt.size());
64  std::vector<double> buf = GetState();
65  ret->Update((unsigned char*)&buf[0], sizeof(double)*buf.size());
66  ret->Final();
67  return ret;
68  }
69 
70  //---------------------------------------------------------------------------
72  {
73  delete fPMNS_Sterile;
74  }
75 
76  //---------------------------------------------------------------------------
78  {
79  return new OscCalcSterile(*this);
80  }
81 
82  //---------------------------------------------------------------------------
83  void OscCalcSterile::SetNFlavors(int nflavors)
84  {
85  fDirty = true;
86  delete fPMNS_Sterile;
87  fNFlavors = nflavors;
89  }
90 
91  //---------------------------------------------------------------------------
92  void OscCalcSterile::SetAngle(int i, int j, double th)
93  {
94  fDirty = true;
95  fPMNS_Sterile->SetAngle(i, j, th);
96  }
97 
98  //---------------------------------------------------------------------------
99  void OscCalcSterile::SetDelta(int i, int j, double delta)
100  {
101  fDirty = true;
102  fPMNS_Sterile->SetDelta(i, j, delta);
103  }
104 
105  //---------------------------------------------------------------------------
106  void OscCalcSterile::SetDm(int i, double dm)
107  {
108  fDirty = true;
109  fPMNS_Sterile->SetDm(i, dm);
110  }
111 
112  //---------------------------------------------------------------------------
113  void OscCalcSterile::SetDmsq21(const double& dmsq21)
114  {
115  std::cerr << "Must use SetDm!" << std::endl;
116  assert(false);
117  }
118 
119  //---------------------------------------------------------------------------
120  void OscCalcSterile::SetDmsq32(const double& dms32)
121  {
122  std::cerr << "Must use SetDm!" << std::endl;
123  assert(false);
124  }
125 
126  //---------------------------------------------------------------------------
127  void OscCalcSterile::SetTh12(const double& th12)
128  {
129  std::cerr << "Must use SetAngle!" << std::endl;
130  assert(false);
131  }
132 
133  //---------------------------------------------------------------------------
134  void OscCalcSterile::SetTh13(const double& th13)
135  {
136  std::cerr << "Must use SetAngle!" << std::endl;
137  assert(false);
138  }
139 
140  //---------------------------------------------------------------------------
141  void OscCalcSterile::SetTh23(const double& th23)
142  {
143  std::cerr << "Must use SetAngle!" << std::endl;
144  assert(false);
145  }
146 
147  //---------------------------------------------------------------------------
148  void OscCalcSterile::SetdCP(const double& dCP)
149  {
150  std::cerr << "Must use SetDelta!" << std::endl;
151  assert(false);
152  }
153 
154  //---------------------------------------------------------------------------
155  double OscCalcSterile::P(int flavBefore, int flavAfter, double E)
156  {
157  const int anti = (flavBefore > 0) ? +1 : -1;
158  //anti must be +/- 1 but flavAfter can be zero
159  assert(flavAfter/anti >= 0);
160  if (anti != fPrevAnti) fDirty = true;
161  fPrevAnti = anti;
162 
163  if (flavBefore != fPrevFlavBefore) fDirty = true;
164  fPrevFlavBefore = flavBefore;
165 
166  if (abs(fPrevE - E) > 1e-8) fDirty = true;
167  fPrevE = E;
168 
169  int i = -1, j = -1;
170  if(abs(flavBefore) == 12) i = 0;
171  if(abs(flavBefore) == 14) i = 1;
172  if(abs(flavBefore) == 16) i = 2;
173  if(abs(flavAfter) == 12) j = 0;
174  if(abs(flavAfter) == 14) j = 1;
175  if(abs(flavAfter) == 16) j = 2;
176  if(abs(flavAfter) == 0) j = 3;
177  assert(i >= 0 && j >= 0);
178 
179  if (fDirty) {
181  // Assume Z/A=0.5
182  const double Ne = fRho/2;
183  fPMNS_Sterile->PropMatter(fL, E, Ne, anti);
184  }
185 
186  // Return the active fraction
187  if (j == 3) return fPMNS_Sterile->P(0) + fPMNS_Sterile->P(1) + fPMNS_Sterile->P(2);
188  else return fPMNS_Sterile->P(j);
189  }
190 
191  //---------------------------------------------------------------------------
193  : OscCalcSterile()
194  {}
195 
196  //---------------------------------------------------------------------------
198  const OscCalcSterile& calc)
200  {
201  std::vector<double> state = calc.GetState();
202  SetState(state);
203  }
204 
205  //---------------------------------------------------------------------------
209  {
210  std::vector<double> state = calc.GetState();
211  SetState(state);
212  }
213 
214  //---------------------------------------------------------------------------
216  {
217  return new OscCalcSterileTrivial(*this);
218  }
219 
220  //---------------------------------------------------------------------------
222  int flavBefore, int flavAfter, double E)
223  {
224  return 1;
225  }
226 
227  //---------------------------------------------------------------------------
229  {
230  const OscCalcSterileTrivial* calc_trivial
231  = dynamic_cast<const OscCalcSterileTrivial*>(calc);
232  if (calc_trivial) return calc_trivial;
233  const OscCalcSterile* calc_sterile = dynamic_cast<const OscCalcSterile*>(calc);
234  if(calc_sterile) return calc_sterile;
235  else std::cout << "Input calculator was not of type OscCalcSterile." << std::endl;
236  return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
237  }
238 
239  //---------------------------------------------------------------------------
241  {
242  OscCalcSterileTrivial* calc_trivial
243  = dynamic_cast<OscCalcSterileTrivial*>(calc);
244  if (calc_trivial) return calc_trivial;
245  OscCalcSterile* calc_sterile = dynamic_cast<OscCalcSterile*>(calc);
246  if(calc_sterile) return calc_sterile;
247  else std::cout << "Input calculator was not of type OscCalcSterile." << std::endl;
248  return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
249  }
250 } // namespace
void SetState(std::vector< double > state)
double GetDelta(int i, int j) const
void SetNFlavors(int nflavors)
virtual void SetDelta(int i, int j, double delta)
virtual void SetAngle(int i, int j, double th)
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
virtual double P(int flv) const
virtual void SetTh12(const double &th12) override
virtual void PropMatter(double L, double E, double Ne, int anti=1)
virtual IOscCalcAdjustable * Copy() const override
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
void SetDelta(int i, int j, double delta)
virtual void SetDm(int i, double dm)
Adapt the PMNS_Sterile calculator to standard interface.
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
Version of OscCalcSterile that always returns probability of 1.
double th12
Definition: runWimpSim.h:98
osc::OscCalcDumb calc
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
virtual TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
virtual void SetdCP(const double &dCP) override
PMNS_Sterile * fPMNS_Sterile
std::vector< double > GetState() const
Float_t E
Definition: plot.C:20
double dCP
virtual void SetDmsq21(const double &dmsq21) override
const double j
Definition: BetheBloch.cxx:29
virtual void ResetToFlavour(int flv=1)
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
Oscillation probability calculators.
Definition: Calcs.h:5
const OscCalcSterile * DowncastToSterile(const IOscCalc *calc)
OStream cout
Definition: OStream.cxx:6
virtual void SetTh13(const double &th13) override
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetDm(int i, double dm)
virtual void SetL(double L) override
virtual void SetRho(double rho) override
assert(nhit_max >=nhit_nbins)
double GetDm(int i) const
virtual void SetDmsq32(const double &dmsq32) override
Float_t e
Definition: plot.C:35
virtual void SetTh23(const double &th23) override
double th13
Definition: runWimpSim.h:98
double GetAngle(int i, int j) const
virtual IOscCalcAdjustable * Copy() const override