CherenkovCalc.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <string.h>
4 #include "TH3.h"
5 
6 class TFile;
7 
8 enum DET { ND, FD };
9 
10 // Parameters for the whole fit
11 struct params{
12  double FDFx, FDFy,
13  NDFx, NDFy,
14  Ys, Ec,
15  attB, att1, att2,
16  rngeff;
17 };
18 
19 // Parameters just for one detector
20 struct detparams{
21  double Fx, Fy,
22  Ys, Ec,
23  attB,
24  // For efficiency, store the inverses
25  inv_att1, inv_att2,
26  rngeff;
27 
28  // Does this happen by default? C++ is so inconsistent about zero-inits...
30  {
31  memset(this, 0, sizeof(detparams));
32  }
33 
34  bool operator !=(const detparams b)
35  {
36  return b.Fx != this->Fx ||
37  b.Fy != this->Fy ||
38  b.Ys != this->Ys ||
39  b.Ec != this->Ec ||
40  b.attB != this->attB ||
41  b.inv_att1 != this->inv_att1 ||
42  b.inv_att2 != this->inv_att2 ||
43  b.rngeff != this->rngeff;
44  }
45 
46  // True if this == b except for Fy, which we don't care about
47  bool equalsx(const detparams & b)
48  {
49  return b.Fx == this->Fx &&
50  // don't check Fy
51  b.Ys == this->Ys &&
52  b.Ec == this->Ec &&
53  b.attB == this->attB &&
54  b.inv_att1 == this->inv_att1 &&
55  b.inv_att2 == this->inv_att2 &&
56  b.rngeff == this->rngeff;
57  }
58 
59  // True if this == b except for Fx, which we don't care about
60  bool equalsy(const detparams & b)
61  {
62  return // don't check Fx
63  b.Fy == this->Fy &&
64  b.Ys == this->Ys &&
65  b.Ec == this->Ec &&
66  b.attB == this->attB &&
67  b.inv_att1 == this->inv_att1 &&
68  b.inv_att2 == this->inv_att2 &&
69  b.rngeff == this->rngeff;
70  }
71 };
72 
73 // Used for caching selected data from the input files.
74 struct dat{
75  float PE_PathLength; // PE divided by path length in cm
76  float PE; // Raw PE
77  float CherenkovPho; // Number of Cherenkov photons
78  float ScintPho; // Number of scintillation photons
79  float InvTotPhoAtten;// 1/((ScintPho + CherenkovPho)*original attenuation)
80  int DistBin; // Which bin in distance to put this hit
81  int wBin; // Which bin in w to put this hit
82  bool proton; // Whether this is truly a proton hit
83  int distshort; // distance to APD the short way, NOTE: in units of atten_step
84  int distlong; // ... the long way
85 };
86 
88 {
89  public:
90  CherenkovCalc(const char * const dataFileName, const char * const mcFileName,
91  const char * const dirName, const char * const plot_name,
92  const std::vector<double> & binX,
93  const std::vector<double> & binY,
94  const std::vector<double> & binZ, const bool isFD_);
95 
96  // Stores the initial parameters affecting one detector
97  void SetInitialParams(const params & mcparams, const enum DET det);
98 
99  void FillDataHistos();
100  void FillMCHistos(const detparams & mcparams,
101  const bool changed_x, const bool changed_y);
102 
103  double peWeight(const bool isX, const dat & d);
104  double CalcChi2(const detparams & mcparams,
105  const bool changed_x, const bool changed_y);
106 
107  TH3D* dataX() {return _dataX;}
108  TH3D* dataY() {return _dataY;}
109  TH3D* mcX() {return _mcX;}
110  TH3D* mcY() {return _mcY;}
111  TH3D* maskX() {return _maskX;}
112  TH3D* maskY() {return _maskY;}
113 
114  // Return the logical NOT of the mask
115  TH3D* invmaskX();
116  TH3D* invmaskY();
117 
118  // Fill masks at the current values of the parameters, and return whether they
119  // were changed.
120  bool MakeMasks();
121 
122  bool isFD;
123 
124  const char * plot_name;
125 
126  private:
127 
128  // Run the first time we go to fill MC histograms. Caches all the hits.
129  // Fills the projection histograms used for the likelihood. Sets masks for
130  // low-stat bins.
131  void fill_hit_cache_and_init_things();
132 
133  // Effectively round all w values to the nearest 5cm for attenuation
134  // purposes. That should be plenty accurate. And since w's aren't varied
135  // in the fit, there's no problem quantizing them.
136  const double atten_step = 3.;
137 
140 
141  static const int kCellsPerModule = 32;
142  const double kPigtails[kCellsPerModule] = {
143  34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
144  64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
145  95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
146  125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
147  };
148 
149  const double maxpigtail = kPigtails[kCellsPerModule-1];
150 
151  std::vector<double> atten_cache_short, atten_cache_long;
152 
153  double getpigtail(const int cell, const bool isX);
154 
155  // dist in cm
156  double exact_atten(const double dist, const bool longway);
157 
158  void init_atten_cache();
159 
160  // Returns true if hits from this 'w' values are thought to be likely modeled
161  // well in MC and not part of the poorly-modeled roll-off at the ends.
162  bool ok_w(const float w, const bool isX);
163 
164  // Returns true if we should accept this hit
165  bool accepthit(const bool isFD, const bool isX, const int plane, const int cell,
166  const float pathlength, const float disttostop,
167  const float w, const float PE);
168 
169  // Optimized filling of MC histograms, avoiding TH2::Fill().
170  void FastMCFill(TH3D * hist,
171  const std::vector<dat> & dats,
172  const bool isX);
173 
174  // Does h->Fill() if in range. Returns true iff it is. 'h' must be one of
175  // the member TH3Ds.
176  bool fill_no_overflow(TH3D * h, const double x, const double y,
177  const double z, const double weight = 1);
178 
181 
182  TFile* _dataFile;
183  TFile* _mcFile;
184  const char * _dirName; // Directory inside the input root files
185 
186  // The top three are synonyms for the bottom three and are the sizes
187  // of the vectors just below, minus one. Of course we can derive these
188  // at any time from the histograms, but this keeps the code clean.
189  unsigned int nx, ny, nz;
190  unsigned int ndistbins, npebins, nwbins;
191 
192  // Bin boundaries for each direction. Size is nbins+1.
193  std::vector<double> distbins, pebins, wbins;
194 
195  // For smearing
196  std::vector<double> extpebins;
197 
198  TH3D* _mcX;
199  TH3D* _mcY;
200  TH3D* _dataX;
201  TH3D* _dataY;
202 
203  // Stores for each MC bin whether there is sufficient statistics (as
204  // determined by the initial parameters) to use it.
205  TH3D* _maskX;
206  TH3D* _maskY;
207 
208  // Projections of the data and MC distributions onto distance-from-end axis
209  TH2D* _mcXdist;
210  TH2D* _mcYdist;
211  TH2D* _dataXdist;
212  TH2D* _dataYdist;
213 
214  std::vector<dat> datsx;
215  std::vector<dat> datsy;
216 };
217 
219 {
220  public:
221  CherenkovFitter(const params & mcparams);
222 
223  void AddFDSample(CherenkovCalc* calc);
224  void AddNDSample(CherenkovCalc* calc);
225 
226  void Fit(const params & mcparams);
227 
228  // write out PDFs, some of which get 'tag' in their file name
229  void Plot(const char * const tag);
230 
231  private:
232 
234 
235  std::vector<CherenkovCalc*> _fdSamples;
236  std::vector<CherenkovCalc*> _ndSamples;
237 };
double FDFx
Definition: CherenkovCalc.h:12
detparams _Shift
double Ys
Definition: CherenkovCalc.h:21
TH3D * maskY()
struct Plot Plot
Definition: DrawUtils.h:15
int distshort
Definition: CherenkovCalc.h:83
std::vector< CherenkovCalc * > _ndSamples
double inv_att1
Definition: CherenkovCalc.h:21
void CalcChi2(const Spectrum &sAll, const Spectrum &sData, FILE *textOFS)
double Ec
Definition: CherenkovCalc.h:21
int distlong
Definition: CherenkovCalc.h:84
unsigned int nz
const Var weight
double FDFy
Definition: CherenkovCalc.h:12
const Binning distbins
double inv_att2
Definition: CherenkovCalc.h:21
fVtxDx Fit("f")
double NDFx
Definition: CherenkovCalc.h:12
double att1
Definition: CherenkovCalc.h:12
double rngeff
Definition: CherenkovCalc.h:21
std::vector< dat > datsy
double attB
Definition: CherenkovCalc.h:12
float CherenkovPho
Definition: CherenkovCalc.h:77
int wBin
Definition: CherenkovCalc.h:81
osc::OscCalcDumb calc
const char * _dirName
double dist
Definition: runWimpSim.h:113
bool proton
Definition: CherenkovCalc.h:82
double Ys
Definition: CherenkovCalc.h:12
double Fx
Definition: CherenkovCalc.h:21
std::vector< double > wbins
int DistBin
Definition: CherenkovCalc.h:80
double rngeff
Definition: CherenkovCalc.h:12
Float_t d
Definition: plot.C:236
DET
Definition: CherenkovCalc.h:8
double NDFy
Definition: CherenkovCalc.h:12
z
Definition: test.py:28
double Ec
Definition: CherenkovCalc.h:12
std::vector< dat > datsx
bool operator!=(View lhs, View rhs)
Definition: BaseProducts.h:155
const char * plot_name
TH3D * dataX()
bool equalsx(const detparams &b)
Definition: CherenkovCalc.h:47
std::vector< CherenkovCalc * > _fdSamples
std::string dirName
Definition: PlotSpectra.h:47
TH3D * maskX()
float InvTotPhoAtten
Definition: CherenkovCalc.h:79
std::vector< double > atten_cache_short
const hit & b
Definition: hits.cxx:21
detparams _Init
double att2
Definition: CherenkovCalc.h:12
float ScintPho
Definition: CherenkovCalc.h:78
float PE
Definition: CherenkovCalc.h:76
TFile * _dataFile
std::vector< double > extpebins
double attB
Definition: CherenkovCalc.h:21
bool equalsy(const detparams &b)
Definition: CherenkovCalc.h:60
Float_t w
Definition: plot.C:20
float PE_PathLength
Definition: CherenkovCalc.h:75
TH3D * dataY()
double cellhalflength
double Fy
Definition: CherenkovCalc.h:21
unsigned int nwbins