OscCalcSterile.cxx
Go to the documentation of this file.
2 
3 #include <cassert>
4 
5 namespace osc
6 {
8  : fPMNS_Sterile(0), fNFlavors(3), fPrevE(0), fPrevAnti(0), fPrevFlavBefore(0)
9  {
11  fDirty = true;
12  }
13 
14  //---------------------------------------------------------------------------
15  std::vector<double> OscCalcSterile::GetState() const
16  {
17  std::vector<double> state;
18  state.push_back((double)fNFlavors);
19  state.push_back(fL);
20  state.push_back(fRho);
21  for(int i = 2; i <= fNFlavors; ++i) state.push_back(GetDm(i));
22  for(int j = 2; j <= fNFlavors; ++j) {
23  for(int i = 1; i < j; ++i) {
24  state.push_back(GetAngle(i, j));
25  if(i+1 != j) state.push_back(GetDelta(i, j));
26  }
27  }
28  return state;
29  }
30 
31  //---------------------------------------------------------------------------
32  void OscCalcSterile::SetState(std::vector<double> state)
33  {
34  int iState(0);
35  fDirty = true ;
36  SetNFlavors((int)state[iState++]);
37  SetL(state[iState++]);
38  SetRho(state[iState++]);
39  for (int i = 2; i <= fNFlavors; ++i) SetDm(i, state[iState++]);
40  for(int j = 2; j <= fNFlavors; ++j) {
41  for(int i = 1; i < j; ++i) {
42  SetAngle(i, j, state[iState++]);
43  if(i+1 != j) SetDelta(i, j, state[iState++]);
44  }
45  }
46  }
47 
48  //---------------------------------------------------------------------------
50  : OscCalcSterile()
51  {
52  std::vector<double> state = calc.GetState();
53  SetState(state);
54  }
55 
56  //---------------------------------------------------------------------------
58  {
59  TMD5* ret = new TMD5;
60  std::string txt = "PMNSSterile";
61  ret->Update((unsigned char*)txt.c_str(), txt.size());
62  std::vector<double> buf = GetState();
63  ret->Update((unsigned char*)&buf[0], sizeof(double)*buf.size());
64  ret->Final();
65  return ret;
66  }
67 
68  //---------------------------------------------------------------------------
70  {
71  delete fPMNS_Sterile;
72  }
73 
74  //---------------------------------------------------------------------------
76  {
77  return new OscCalcSterile(*this);
78  }
79 
80  //---------------------------------------------------------------------------
81  void OscCalcSterile::SetNFlavors(int nflavors)
82  {
83  fDirty = true;
84  delete fPMNS_Sterile;
85  fNFlavors = nflavors;
87  }
88 
89  //---------------------------------------------------------------------------
90  void OscCalcSterile::SetAngle(int i, int j, double th)
91  {
92  fDirty = true;
93  fPMNS_Sterile->SetAngle(i, j, th);
94  }
95 
96  //---------------------------------------------------------------------------
97  void OscCalcSterile::SetDelta(int i, int j, double delta)
98  {
99  fDirty = true;
100  fPMNS_Sterile->SetDelta(i, j, delta);
101  }
102 
103  //---------------------------------------------------------------------------
104  void OscCalcSterile::SetDm(int i, double dm)
105  {
106  fDirty = true;
107  fPMNS_Sterile->SetDm(i, dm);
108  }
109 
110  //---------------------------------------------------------------------------
111  double OscCalcSterile::P(int flavBefore, int flavAfter, double E)
112  {
113  const int anti = (flavBefore > 0) ? +1 : -1;
114  //anti must be +/- 1 but flavAfter can be zero
115  assert(flavAfter/anti >= 0);
116  if (anti != fPrevAnti) fDirty = true;
117  fPrevAnti = anti;
118 
119  if (flavBefore != fPrevFlavBefore) fDirty = true;
120  fPrevFlavBefore = flavBefore;
121 
122  if (abs(fPrevE - E) > 1e-8) fDirty = true;
123  fPrevE = E;
124 
125  int i = -1, j = -1;
126  if(abs(flavBefore) == 12) i = 0;
127  if(abs(flavBefore) == 14) i = 1;
128  if(abs(flavBefore) == 16) i = 2;
129  if(abs(flavAfter) == 12) j = 0;
130  if(abs(flavAfter) == 14) j = 1;
131  if(abs(flavAfter) == 16) j = 2;
132  if(abs(flavAfter) == 0) j = 3;
133  assert(i >= 0 && j >= 0);
134 
135  if (fDirty) {
137  // Assume Z/A=0.5
138  const double Ne = fRho/2;
139  fPMNS_Sterile->PropMatter(fL, E, Ne, anti);
140  }
141 
142  // Return the active fraction
143  if (j == 3) return fPMNS_Sterile->P(0) + fPMNS_Sterile->P(1) + fPMNS_Sterile->P(2);
144  else return fPMNS_Sterile->P(j);
145  }
146 
147 } // namespace
void SetState(std::vector< double > state)
virtual void SetDm(int i, double dm) override
void SetNFlavors(int nflavors)
virtual void SetDelta(int i, int j, double delta)
virtual void SetAngle(int i, int j, double th)
double delta
Definition: runWimpSim.h:98
virtual double P(int flv) const
virtual double GetDm(int i) const override
virtual void PropMatter(double L, double E, double Ne, int anti=1)
virtual void SetDm(int i, double dm)
Adapt the PMNS_Sterile calculator to standard interface.
void abs(TH1 *hist)
osc::OscCalcDumb calc
virtual TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
virtual void SetRho(double rho) override
PMNS_Sterile * fPMNS_Sterile
virtual double GetDelta(int i, int j) const override
std::vector< double > GetState() const
Float_t E
Definition: plot.C:20
virtual void SetAngle(int i, int j, double th) 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
assert(nhit_max >=nhit_nbins)
virtual void SetDelta(int i, int j, double delta) override
virtual double GetAngle(int i, int j) const override
Float_t e
Definition: plot.C:35
virtual void SetL(double L) override
virtual IOscCalcAdjustable * Copy() const override
enum BeamMode string