CherenkovCalc.cxx
Go to the documentation of this file.
1 #include <string.h>
2 #include <signal.h>
3 #include <vector>
4 
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TH3.h"
8 #include "TString.h"
9 #include "TFile.h"
10 #include "TTreeReader.h"
11 #include "TTreeReaderValue.h"
12 #include "TMinuit.h"
13 #include "TRandom3.h"
14 
15 #include "consts.h"
16 #include "CherenkovCalc.h"
17 #include "Plot.h"
18 
19 // Remove for a quick and inaccurate fit
20 #define FULL
21 
22 const bool float_ndfx = false,
23  float_ndfy = false,
24 
25  float_fdfx = true,
26  float_fdfy = true,
27 
28  float_ec = false,
29 
30  float_attb = false,
31  float_att1 = false,
32  float_att2 = false,
33 
34  // So far most fits have preferred a value of 1, so leave it
35  // Based on Kurarray data, expect ~0.99 efficiency.
36  float_rngeff = false;
37 
38 static TRandom3 R(0); // 0 means random seed
39 
40 using std::vector;
41 
42 // We copy the pointers to the samples into these static variables so that
43 // MINUIT's objective function can get at them.
44 static vector<CherenkovCalc*> sfdSamples, sndSamples;
45 
46 static const unsigned int npar = 10;
47 
48 //#define USEPROTONWEIGHT
49 
50 #ifdef USEPROTONWEIGHT
51  // Weight true proton hits in the MC by this amount, to crudely adjust for
52  // cross section uncertainties.
53  const double PROTONWEIGHT = 1;
54 #endif
55 
56 // The minimum accepted reconstructed path through the cell.
57 const double MINPATHLENGTH = 1; // cm
58 
59 // For speed, don't accept MC events after this many, per sample.
60 static const unsigned int maxMChitsPERw = 1e7;
61 
62 // In order to prevent one sample from overpowering the rest (primarily this
63 // is a concern for FD cosmics), limit the number of data events per sample.
64 // A discussion of statistics philosophy should appear here.
65 static const unsigned int maxDATAhitsPERw = 1e6;
66 
67 // How many hits an MC bin has to have using the initial parameters for that
68 // bin to be used
69 static const int stat_req = 1000;
70 
71 // Warn if a bin's contents goes below this level as the fit proceeds
72 static const int stat_warn_level = 100;
73 
74 // If improvement is less than this after SCAN or SIMPLEX, or whatnot, don't
75 // trigger another round of trying to improve.
76 static const double epsilonll = 1e-4;
77 
78 // Project 'in' along XZ to 'out'. The user is responsible for making the
79 // number of X and Z bins match between the two histograms. All overflow bins
80 // are used. This is a replacement for ROOT's insane ProjectionX().
81 static void saneProjectionXZ(TH2D * out, TH3D * in)
82 {
83  double sumsum = 0;
84  for(int bz = 0; bz <= in->GetNbinsZ()+1; bz++){
85  for(int b = 0; b <= in->GetNbinsX()+1; b++){
86  double sum = 0;
87  for(int by = 0; by <= in->GetNbinsY()+1; by++)
88  sum += in->GetBinContent(b, by, bz);
89  out->SetBinContent(b, bz, sum);
90  sumsum += sum;
91  }
92  }
93  out->SetEntries(sumsum);
94 }
95 
96 CherenkovCalc::CherenkovCalc(const char * const dataFileName,
97  const char * const mcFileName,
98  const char * const dirName,
99  const char * const plot_name_,
100  const vector<double> & binX,
101  const vector<double> & binY,
102  const vector<double> & binZ,
103  const bool isFD_)
104 : isFD(isFD_), plot_name(plot_name_)
105 {
106  if(isFD) cellhalflength = 1549.4/2;
107  else cellhalflength = 399.28/2;
108 
110 
111  _dataFile = new TFile(dataFileName);
112  _mcFile = new TFile(mcFileName);
113 
114  if(_dataFile == NULL || _dataFile->IsZombie()){
115  fprintf(stderr, "Failed to open %s\n", dataFileName);
116  _exit(1);
117  }
118  if(_mcFile == NULL || _mcFile->IsZombie()){
119  fprintf(stderr, "Failed to open %s\n", mcFileName);
120  _exit(1);
121  }
122 
123  _dirName = dirName;
124 
125  if(binX.size() < 2 || binY.size() < 2 || binZ.size() < 2){
126  fprintf(stderr, "Must have at least two bin boundaries per dimension\n");
127  _exit(1);
128  }
129 
130  ndistbins = nx = binX.size()-1;
131  npebins = ny = binY.size()-1;
132  nwbins = nz = binZ.size()-1;
133 
134  distbins = binX;
135  pebins = binY;
136  wbins = binZ;
137 
138  // For smearing
139  extpebins = binY;
140  extpebins.insert(extpebins.begin(), pebins[0] - (pebins[1] - pebins[0]));
141  extpebins.push_back(pebins[ny] + pebins[ny] - pebins[ny-1]);
142 
143  const char * const title = ";Distance To Stop (cm);PE/Path (cm^{-1})";
144 
145  _mcX = new TH3D(Form("mcX_%s_%s", mcFileName, dirName),
146  title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
147  _mcY = new TH3D(Form("mcY_%s_%s", mcFileName, dirName),
148  title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
149  _dataX = new TH3D("dataX", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
150  _dataY = new TH3D("dataY", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
151 
152  _maskX = new TH3D("maskX", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
153  _maskY = new TH3D("maskY", title, nx, &binX[0], ny, &binY[0], nz, &binZ[0]);
154 
155  _mcXdist = new TH2D("mcXdist", "", nx, &binX[0], nz, &binZ[0]);
156  _mcYdist = new TH2D("mcYdist", "", nx, &binX[0], nz, &binZ[0]);
157  _dataXdist = new TH2D("dataXdist", "", nx, &binX[0], nz, &binZ[0]);
158  _dataYdist = new TH2D("dataYdist", "", nx, &binX[0], nz, &binZ[0]);
159 }
160 
162  const enum DET det)
163 {
164  _Init.Fx = det == ND? mcparams.NDFx: mcparams.FDFx;
165  _Init.Fy = det == ND? mcparams.NDFy: mcparams.FDFy;
166  _Init.attB = mcparams.attB;
167  _Init.inv_att1 = 1/mcparams.att1;
168  _Init.inv_att2 = 1/mcparams.att2;
169  _Init.Ys = mcparams.Ys;
170  _Init.Ec = mcparams.Ec;
171  _Init.rngeff = mcparams.rngeff;
172 }
173 
174 // Copied from Geo::getPigtail
175 double CherenkovCalc::getpigtail(const int cell, const bool isX)
176 {
177  int cellid = cell % kCellsPerModule;
178  // In vertical planes, high cell numbers have longer fibres.
179  // In horizontal planes it's the opposite, so flip it round.
180  if(!isX) cellid = kCellsPerModule-cellid-1;
181  // This should really never happen, but just to be safe...
182  if(cellid < 0 || cellid >= kCellsPerModule) return 100;
183  // Email from Tom Chase 2011-04-29
184  // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
185  // different happens.
186  return kPigtails[cellid];
187 }
188 
189 // Pure function that returns attenuation
190 static double _exact_atten(const double attB, const double inv_att1,
191  const double inv_att2, const double dist,
192  const double rngeff, const bool longway)
193 {
194  return (longway?rngeff:1)*
195  (attB * exp(-dist*inv_att1) + (1-attB) * exp(-dist*inv_att2));
196 }
197 
198 // Calculate the exact attenaution for a given distance 'dist' and the
199 // current parameters. Intended only for use to fill the cache.
200 double CherenkovCalc::exact_atten(const double dist, const bool longway)
201 {
203  dist, _Shift.rngeff, longway);
204 }
205 
206 // Because calculating attenuation exactly is expensive, make a table
207 // of attenuation to the nearest 'atten_step' (5cm as I write this).
209 {
211  atten_cache_long .resize(nattensteps);
212  for(int d = 0; d < nattensteps; d++){
214  atten_cache_long [d] = exact_atten(d*atten_step, true );
215  }
216 }
217 
218 // Make a TTreeReader given a file and the name of the desired tree.
219 // Be tolerant of name changes over time and try a few.
220 //
221 // Implemented with goto's for fun!
222 static TTreeReader * makereader(const char * const dir, TFile * file)
223 {
224  TTreeReader * reader = new TTreeReader(Form("cerenkov/%s", dir), file);
225  if(!reader->IsInvalid()) goto ok;
226 
227  reader->SetTree(Form("cerenkovselection/%s", dir), file);
228  if(!reader->IsInvalid()) goto ok;
229 
230  reader->SetTree(Form("cerenkovresponse/%s", dir), file);
231  if(!reader->IsInvalid()) goto ok;
232 
233  if(strcmp(dir, "cerenkovMuonTree")) goto fail;
234 
235  reader->SetTree("cerenkovresponse/cerenkovTree", file);
236  if(!reader->IsInvalid()) goto ok;
237 
238  goto fail;
239 
240  ok:
241 
242  if(reader->GetEntries(false) == 0){
243  fprintf(stderr, "Empty input from %s\n", file->GetName());
244  _exit(1);
245  }
246  return reader;
247 
248  fail:
249 
250  fprintf(stderr, "Couldn't get tree from %s, %s\n", file->GetName(), dir);
251  _exit(1);
252 }
253 
254 static inline double clamp(const double v, const double mi, const double ma)
255 {
256  return v < mi? mi: v > ma? ma: v;
257 }
258 
259 struct binw{
260  double prev, curr, next;
261 };
262 
263 #define SMEAR
264 
265 static binw smearbin(
266 #ifndef SMEAR
267  __attribute__((unused)) const vector<double> & bins,
268  __attribute__((unused)) const double val,
269  __attribute__((unused)) const int bin
270 #else
271  const vector<double> & bins, const double val, const int bin
272 #endif
273 )
274 {
275  binw ans;
276 
277  // This is the regular way of binning!
278 #ifndef SMEAR
279  ans.prev = ans.next = 0;
280  ans.curr = 1;
281  return ans;
282 #else
283  const double binlowedge = bins[bin];
284  const double binhigedge = bins[bin+1];
285  const double binfrac = clamp(( val - binlowedge)/
286  (binhigedge - binlowedge), -0.5, 1.5);
287  if(binfrac > 0.5){
288  ans.prev = 0;
289  ans.next = binfrac - 0.5;
290  ans.curr = 1 - ans.next;
291  }
292  else{
293  ans.prev = 0.5 - binfrac;
294  ans.curr = 1 - ans.prev;
295  ans.next = 0;
296  }
297 #endif
298 
299 #if 0
300  printf("%f %d:\t%f %f %f %f \n", val, bin, binfrac, ans.prev, ans.curr, ans.next);
301 #endif
302 
303  return ans;
304 }
305 
306 static void Add(TH3D * h, const int bx, const int by, const int bz,
307  const double w)
308 {
309  if(w == 0) return;
310  const int gbin = h->GetBin(bx, by, bz);
311  h->SetBinContent(gbin, h->GetBinContent(gbin) + w);
312 }
313 
314 // Get the zero-indexed bin number that 'pe' is in. Optimize by checking
315 // if it is the same as the previous answer. This works because we
316 // sorted by the initial PE.
317 static unsigned int getpebin(const double pe,
318  const std::vector<double> & bins)
319 {
320  static unsigned int prev = 0;
321  if(prev < bins.size()-1 && bins[prev] < pe && bins[prev+1] > pe)
322  return prev;
323 
324  // Don't try to optimize by restricting the search range. The pe value may
325  // have been modified so that it is out of range.
326  return prev = std::lower_bound(bins.begin(), bins.end(), pe)
327  - bins.begin() - 1;
328 }
329 
330 bool CherenkovCalc::fill_no_overflow(TH3D * h, const double x,
331  const double y, const double z,
332  const double weight)
333 {
334  const double xlo = distbins[0];
335  const double xhi = distbins[nx]; // sic, because nx is 1 less than the size
336  const double zlo = wbins[0];
337  const double zhi = wbins[nz];
338 
339  // Do not check y here; that's handled below
340  if(!(x < xhi && x >= xlo && z < zhi && z >= zlo)) return false;
341 
342  const int xbin = h->GetXaxis()->FindBin(x);
343  const int ybin = getpebin(y, extpebins);
344  const int zbin = h->GetZaxis()->FindBin(z);
345 
346  // extpebins is indexed like ROOT histograms with overflow: 0 to ny+1
347  const binw bw = smearbin(extpebins, y, ybin);
348 
349  if(ybin-1 >= 1) Add(h,xbin,ybin-1,zbin,bw.prev*weight);
350  if(ybin >= 1 && ybin <= (int)ny) Add(h,xbin,ybin ,zbin,bw.curr*weight);
351  if( ybin+1 <= (int)ny) Add(h,xbin,ybin+1,zbin,bw.next*weight);
352 
353  return true;
354 }
355 
356 bool CherenkovCalc::ok_w(const float w, __attribute__((unused)) const bool isX)
357 {
358  // From looking at the curves at
359  // https://nusoft.fnal.gov/nova/calibration/Attenuation_Calibration/SecondAna/Period3/Epoch3c
360  // Want to avoid the roll-off at the ends, even if we think we get
361  // some of it right.
362 
363 //#define RESTRICTFDW
364 
365  if(isFD){
366 #ifdef RESTRICTFDW
367  if(!isX) return w > fdwmin && w < fdwmax;
368  const double cen = 500;
369  return abs(w - cen) < 100;
370 #else
371  return w > fdwmin && w < fdwmax;
372 #endif
373  }
374  else return w > ndwmin && w < ndwmax;
375 }
376 
377 static bool goodfdplane(const int plane)
378 {
379 //#define BYDIBLOCK
380 
381 #ifdef BYDIBLOCK
382  const int diblock = 14;
383  return plane >= 64*(diblock-1) &&
384  plane < 64*diblock;
385 #else
386  return plane >= 64;
387 #endif
388 }
389 
390 static bool goodfdcell(__attribute__((unused)) const bool isX,
391  __attribute__((unused)) const int cell)
392 {
393  if(isX) return true;
394 
395 #define ALL
396 
397 #if defined ALL
398  return true;
399 #elif defined BOTTOMTHIRD
400  return cell < 32 * 4;
401 #elif defined MIDDLETHIRD
402  return cell >= 32 * 4 && cell < 32 * 8;
403 #elif defined TOPTHIRD
404  return cell >= 32 * 8 ;
405 #endif
406 }
407 
408 bool CherenkovCalc::accepthit(const bool isFD, const bool isX, const int plane,
409  const int cell,
410  const float pathlength, const float disttostop,
411  const float w, const float PE)
412 {
413  // remove first two blocks of FD, which have lower pseudocumene scintillator
414  if(isFD && !goodfdplane(plane)) return false;
415 
416  if(isFD && !goodfdcell(isX, cell)) return false;
417 
418  if(pathlength <= MINPATHLENGTH) return false;
419  if(disttostop > distbins[distbins.size()-1]) return false;
420  if(!ok_w(w, isX)) return false;
421 
422  // Data and MC disagree on where the threshold is, so remove the
423  // disagreeing part. Not sure if this is also true for the ND, but ND
424  // hits aren't mostly pasted up against the threshold, so it's not
425  // such an obvious problem.
426  if(isFD && PE < 30) return false;
427 
428  return true;
429 }
430 
432 {
433  TTreeReader * reader = makereader(_dirName, _dataFile);
434 
435  TTreeReaderValue<float> rvPathLength(*reader, "pathLength");
436  TTreeReaderValue<float> rvPE (*reader, "pe");
437  TTreeReaderValue<float> rvDistToStop(*reader, "distFromEnd");
438  TTreeReaderValue<char> rvIsX (*reader, "isX");
439  TTreeReaderValue<float> rvw (*reader, "w");
440  TTreeReaderValue<int> rvPlane (*reader, "plane");
441  TTreeReaderValue<int> rvCell (*reader, "cell");
442 
443  _dataX->Reset();
444  _dataY->Reset();
445  _dataXdist->Reset();
446  _dataYdist->Reset();
447 
448  printf("%9u %s %s\n", (unsigned)reader->GetEntries(false),
449  _dataFile->GetName(), _dirName);
450  fflush(stdout);
451 
452  unsigned int got[nwbins] = { 0 };
453  bool capped[nwbins] = { false };
454  unsigned int ncapped = 0;
455  while (reader->Next()) {
456  // remove first 2 blocks of FD, which have lower pseudocumene scintillator
457  if(!accepthit(isFD, *rvIsX, *rvPlane, *rvCell, *rvPathLength, *rvDistToStop, *rvw, *rvPE))
458  continue;
459 
460  // Zero-indexed bin number
461  const unsigned int wBin =
462  std::lower_bound(wbins.begin(), wbins.end(), *rvw) - wbins.begin() - 1;
463 
464  if(wBin >= wbins.size()) continue; // shouldn't be possible
465 
466  if(got[wBin] >= maxDATAhitsPERw){
467  if(!capped[wBin]){
468  capped[wBin] = true;
469  ncapped++;
470  }
471  if(nwbins == ncapped) break;
472  else continue;
473  }
474 
475  if(fill_no_overflow((*rvIsX)?_dataX:_dataY,
476  *rvDistToStop, (*rvPE)/(*rvPathLength), *rvw))
477  got[wBin]++;
478  }
479 
480  for(unsigned int i = 0; i < nwbins; i++)
481  printf("%9u data hits accepted in w bin %d%s\n", got[i],
482  i, capped[i]?" (capped)": Form(" (want %.1fx more)",
483  (float)maxDATAhitsPERw/got[i]));
484  printf("\n");
485 
486  delete reader;
487 
490 
491  _dataFile->Close();
492 }
493 
494 static TH3D * invmask(TH3D * mask)
495 {
496  TH3D * notmask = (TH3D*)mask->Clone("notX");
497  for(int i = 1; i <= mask->GetNbinsX(); i++)
498  for(int j = 1; j <= mask->GetNbinsY(); j++)
499  for(int k = 1; k <= mask->GetNbinsZ(); k++)
500  notmask->SetBinContent(i, j, k, !mask->GetBinContent(i, j, k));
501  return notmask;
502 }
503 
505 {
506  return invmask(_maskX);
507 }
508 
510 {
511  return invmask(_maskY);
512 }
513 
514 // Return true if any change was made
515 static bool fill_mask(TH3D * mask, TH3D * h)
516 {
517  printf("New mask for %s\n", h->GetName());
518  bool changed = false;
519  for(int k = 1; k <= h->GetNbinsZ(); k++){
520  for(int i = 1; i <= h->GetNbinsX(); i++){
521  for(int j = 1; j <= h->GetNbinsY(); j++){
522  const bool newval = h->GetBinContent(i, j, k) > stat_req;
523  const bool oldval = mask->GetBinContent(i, j, k);
524  printf("%c", newval?'.':'X');
525  if(oldval != newval){
526  changed = true;
527  mask->SetBinContent(i, j, k, newval);
528  }
529  }
530  puts("");
531  }
532  puts("");
533  }
534  return changed;
535 }
536 
537 
539 {
540  const bool changedx = fill_mask(_maskX, _mcX);
541  const bool changedy = fill_mask(_maskY, _mcY);
542  return changedx || changedy;
543 }
544 
545 static bool compare_dat(const dat & a, const dat & b)
546 {
547  return a.PE_PathLength < b.PE_PathLength;
548 }
549 
550 double CherenkovCalc::peWeight(const bool isX, const dat & d)
551 {
552  const double attenafter = atten_cache_short[d.distshort]
553  + atten_cache_long [d.distlong ];
554 
555  const double NphoShift = (isX?_Shift.Fx:_Shift.Fy)*
557 
558  return NphoShift * attenafter * d.InvTotPhoAtten;
559 }
560 
561 // Speeds up the whole program by 40% compared to naive use of TH2::Fill().
563  const std::vector<dat> & dats,
564  const bool isX)
565 {
566  double newhist[nx][ny][nz];
567  memset(newhist, 0, sizeof(double)*nx*ny*nz);
568 
569  for(const dat & d : dats){
570  const double weight = peWeight(isX, d);
571 
572  if(weight == 0) continue;
573 
574  // Apply a very rough threshold so that low energy hits drop out entirely
575  // Of course, this can only kill hits, not revive them.
576  const double threshold = 17.;
577  const double scaled_pe = d.PE * weight;
578  if(scaled_pe < threshold) continue;
579 
580  const double scaled_pe_pathLength = d.PE_PathLength * weight;
581 
582  // Range checked when we cached the hits: don't check again
583  const int xbin = d.DistBin, zbin = d.wBin;
584 
585  // Returns zero-indexed answer for extpebins, and subtract one to get
586  // zero-indexed answer for pebins.
587  const int ybin = getpebin(scaled_pe_pathLength, extpebins) - 1;
588 
589  if(!(ybin >= -1 && ybin < (int)ny+1 /* for overflow smear bins */)) continue;
590 
591  // Simple smearing over this bin and one adjacent bin so that events don't
592  // move abruptly and MINUIT can work better. Only needed for y (pe/path).
593  const binw bw = smearbin(extpebins, scaled_pe_pathLength, ybin+1 /* to index into extpebins */);
594 
595  #ifdef USEPROTONWEIGHT
596  if(d.proton && PROTONWEIGHT != 1.){
597  if(ybin-1>=0) newhist[xbin][ybin-1][zbin] += bw.prev*PROTONWEIGHT;
598  if(ybin >=0&&ybin <(int)ny)newhist[xbin][ybin ][zbin] += bw.curr*PROTONWEIGHT;
599  if( ybin+1<(int)ny)newhist[xbin][ybin+1][zbin] += bw.next*PROTONWEIGHT;
600  }
601  else
602  #endif
603  {
604  if(ybin-1>=0) newhist[xbin][ybin-1][zbin] += bw.prev;
605  if(ybin >=0&&ybin <(int)ny)newhist[xbin][ybin ][zbin] += bw.curr;
606  if( ybin+1<(int)ny)newhist[xbin][ybin+1][zbin] += bw.next;
607  }
608  }
609 
610  double tot = 0;
611  for(unsigned int i = 0; i < nx; i++){
612  for(unsigned int j = 0; j < ny; j++){
613  for(unsigned int k = 0; k < nz; k++){
614  tot += newhist[i][j][k];
615  hist->SetBinContent(i+1, j+1, k+1, newhist[i][j][k]);
616  }
617  }
618  }
619  hist->SetEntries(tot);
620 }
621 
623 {
624  TTreeReader * reader = makereader(_dirName, _mcFile);
625 
626  TTreeReaderValue<float> rvPathLength (*reader, "pathLength");
627  TTreeReaderValue<float> rvPE (*reader, "pe");
628  TTreeReaderValue<float> rvDistToStop (*reader, "distFromEnd");
629  TTreeReaderValue<char> rvIsX (*reader, "isX");
630  TTreeReaderValue<float> rvCerenkovPho(*reader, "cerenkovPho");
631  TTreeReaderValue<float> rvEBirks (*reader, "eBirks");
632  TTreeReaderValue<float> rvw (*reader, "w");
633  TTreeReaderValue<char> rvIsProton (*reader, "isProton");
634  TTreeReaderValue<int> rvPlane (*reader, "plane");
635  TTreeReaderValue<int> rvCell (*reader, "cell");
636 
637  _mcX->Reset();
638  _mcY->Reset();
639 
640  printf("%9u %s %s\n", (unsigned)reader->GetEntries(false),
641  _mcFile->GetName(), _dirName);
642 
643  unsigned int got[nwbins] = { 0 };
644  unsigned int gotall = 0, trueprotons = 0, truemuons = 0;
645  bool capped[nwbins] = { false };
646  unsigned int ncapped = 0;
647 
648  printf("Filling hit cache:\n");
649  fflush(stdout);
650  while (reader->Next()) {
651  // remove first two blocks of FD, which have lower pseudocumene scintillator
652  if(!accepthit(isFD, *rvIsX, *rvPlane, *rvCell, *rvPathLength, *rvDistToStop, *rvw, *rvPE)) continue;
653 
654  if(*rvEBirks <= 0) continue;
655 
656  dat d;
657  d.PE = *rvPE;
658  d.PE_PathLength = (*rvPE)/(*rvPathLength);
659  d.CherenkovPho = (*rvCerenkovPho)*_Init.Ec;
660 
661  // zero-indexed bins
662  d.DistBin = std::lower_bound(distbins.begin(), distbins.end(), *rvDistToStop)
663  - distbins.begin() - 1;
664  if(d.DistBin < 0 || d.DistBin >= (int)nx) continue;
665  d.wBin = std::lower_bound(wbins.begin(), wbins.end(), *rvw)
666  - wbins.begin() - 1;
667  if(d.wBin < 0 || d.wBin >= (int)nz) continue;
668  d.ScintPho = (*rvEBirks)*1000*_Init.Ys;
669  d.proton = *rvIsProton;
670 
671  const double distshort = cellhalflength - (*rvw) + getpigtail(*rvCell, *rvIsX);
672  const double distlong = 3*cellhalflength + (*rvw) + getpigtail(*rvCell, *rvIsX);
673 
674  const double attenbefore = _exact_atten(_Init.attB, _Init.inv_att1, _Init.inv_att2,
675  distshort, _Init.rngeff, false) +
677  distlong , _Init.rngeff, true );
678 
679  d.InvTotPhoAtten = 1/((d.ScintPho + d.CherenkovPho)*attenbefore);
680  if(d.InvTotPhoAtten <= 0) continue;
681 
682  d.distshort = clamp(0, distshort/atten_step + 0.5, nattensteps-1);
683  d.distlong = clamp(0, distlong /atten_step + 0.5, nattensteps-1);
684 
685  if(got[d.wBin] >= maxMChitsPERw){
686  if(!capped[d.wBin]){
687  capped[d.wBin] = true;
688  ncapped++;
689  }
690  if(nwbins == ncapped) break;
691  else continue;
692  }
693 
694  if(*rvIsX) datsx.push_back(d);
695  else datsy.push_back(d);
696 
697  if(fill_no_overflow((*rvIsX)?_mcX:_mcY, *rvDistToStop,
698  (*rvPE)/(*rvPathLength), *rvw,
699  #ifdef USEPROTONWEIGHT
700  (*rvIsProton)?PROTONWEIGHT:1)){
701  #else
702  1)){
703  #endif
704  got[d.wBin]++;
705  gotall++;
706  if(*rvIsProton) trueprotons++;
707  }
708 
709  if(gotall%1000000 == 0)
710  printf("Accepted %d hits and counting...\n", gotall);
711  }
712 
713  printf("Accepted %.0f%% true protons\n", 100.*trueprotons/gotall);
714  printf("Accepted %.0f%% true muons\n" , 100.*truemuons /gotall);
715 
716  delete reader;
717 
718  printf("Sorting hits\n");
719  // Optimize by sorting by PE/pathlength and then starting each bin search
720  // by guessing that the PE/pathlength bin is the same as the previous.
721  std::sort(datsx.begin(), datsx.end(), compare_dat);
722  std::sort(datsy.begin(), datsy.end(), compare_dat);
723  printf("Hit cache ready\n");
724 
725  for(unsigned int i = 0; i < nwbins; i++)
726  printf("%9u MC hits accepted in w bin %d%s\n", got[i],
727  i, capped[i]?" (capped)": Form(" (want %.1fx more)",
728  (float)maxMChitsPERw/got[i]));
729  printf("\n");
730 
732  saneProjectionXZ(_mcYdist, _mcY);
733 
734  MakeMasks();
735 }
736 
738  const bool changed_x, const bool changed_y)
739 {
740  _Shift.Fx = mcparams.Fx/_Init.Fx;
741  _Shift.Fy = mcparams.Fy/_Init.Fy;
742  _Shift.Ys = mcparams.Ys/_Init.Ys;
743  _Shift.Ec = mcparams.Ec/_Init.Ec;
744 
745  _Shift.attB = mcparams.attB;
746  _Shift.inv_att1 = mcparams.inv_att1;
747  _Shift.inv_att2 = mcparams.inv_att2;
748  _Shift.rngeff = mcparams.rngeff;
750 
752 
753  if(changed_x) FastMCFill(_mcX, datsx, true);
754  if(changed_y) FastMCFill(_mcY, datsy, false);
755 
756  // Do not need to fill projections again because they are unchanged.
757 }
758 
759 static bool complain_about_problems = true;
760 
761 static double
762 llike_hist(const TH3D * const hObs, const TH3D * const hMC,
763  const TH2D * const hObsDist, const TH2D * const hMCDist,
764  const TH3D * const mask)
765 {
766  double llike = 0;
767  int problems = 0;
768  for(int iZ = 1; iZ <= hMC->GetNbinsZ(); ++iZ) {
769  for(int iX = 1; iX <= hMC->GetNbinsX(); ++iX) {
770  const double data_at_dist = hObsDist->GetBinContent(iX, iZ);
771  const double mc_at_dist = hMCDist->GetBinContent(iX, iZ);
772  const double mc_scale = data_at_dist/mc_at_dist;
773 
774  for(int iY = 1; iY <= hMC->GetNbinsY(); ++iY) {
775  if(mask->GetBinContent(iX, iY, iZ) == 0) continue;
776 
777  const double data = hObs->GetBinContent(iX, iY, iZ);
778  const double rawmc = hMC->GetBinContent(iX, iY, iZ);
779 
780  if(rawmc < stat_warn_level){
781  problems++;
782  if(rawmc <= 0) continue;
783  }
784 
785  // normalize MC to match data
786  const double mc = rawmc * mc_scale;
787 
788  llike += mc - data + (data <= 0?0:data * log (data/mc));
789  }
790  }
791  }
792  if(complain_about_problems && problems)
793  printf("%d low-stat MC bin%s in %s!\n",
794  problems, problems == 1?"":"s", hMC->GetName());
795 
796  return llike;
797 }
798 
799 double CherenkovCalc::CalcChi2(const detparams & mcparams,
800  const bool changed_x, const bool changed_y)
801 {
802  FillMCHistos(mcparams, changed_x, changed_y);
805 }
806 
808  : _Init(mcparams)
809 {
810 }
811 
813 {
814  calc->SetInitialParams(_Init, FD);
815  calc->FillDataHistos();
816  _fdSamples.push_back(calc);
817 }
818 
820 {
821  calc->SetInitialParams(_Init, ND);
822  calc->FillDataHistos();
823  _ndSamples.push_back(calc);
824 }
825 
826 // Mechanism for the user to request plotting in the middle of the run
827 static CherenkovFitter * thefitter = NULL;
828 static bool plot = false;
829 
830 static void intermediate_plot(const int signal)
831 {
832  if(signal != SIGUSR1) return;
833  plot = true;
834 }
835 
836 // Keep track of the best values independently of MINUIT, because if
837 // MINUIT doesn't think a fit has converged, it may not leave the
838 // parameters at the best point seen so far (yes, this certainly
839 // happens).
840 static double lowest_llike = 1e300;
841 static double bestpars[npar] = { 0 };
842 
843 // The objective function for MINUIT.
844 static void objfcn(__attribute__((unused)) int & np,
845  __attribute__((unused)) double * gin,
846  double & llike,
847  double * par,
848  int flag)
849 {
850  // taking the absolute value is a nice way to prevent non-physical nonsense
851  // without using SET LIMITS. The fit will only hit negative values when it is
852  // going crazy, and this can give it a way to get back home.
853  detparams FD, ND;
854  static detparams old_FD, old_ND;
855  FD.Fx = fabs(par[0]);
856  FD.Fy = fabs(par[1]);
857  ND.Fx = fabs(par[2]);
858  ND.Fy = fabs(par[3]);
859 
860  FD.Ys = ND.Ys = fabs(par[4]);
861  FD.Ec = ND.Ec = fabs(par[5]);
862  FD.attB = ND.attB = fabs(par[6]);
863  FD.inv_att1 = ND.inv_att1 = 1/fabs(par[7]);
864  FD.inv_att2 = ND.inv_att2 = 1/fabs(par[8] + par[7]);
865  FD.rngeff = ND.rngeff = fabs(par[9]);
866 
867  double ndllike = 0, fdllike = 0;
868  static double old_ndllike = 0, old_fdllike = 0;
869 
870  const bool fdchanged = old_FD != FD;
871  const bool ndchanged = old_ND != ND;
872 
873  if(fdchanged){
874  const bool fdchanged_x = !old_FD.equalsx(FD);
875  const bool fdchanged_y = !old_FD.equalsy(FD);
876 
878  fdllike += calc->CalcChi2(FD, fdchanged_x, fdchanged_y);
879  old_fdllike = fdllike;
880  }
881  else{
882  fdllike = old_fdllike;
883  }
884 
885  if(ndchanged){
886  const bool ndchanged_x = !old_ND.equalsx(ND);
887  const bool ndchanged_y = !old_ND.equalsy(ND);
888 
890  ndllike += calc->CalcChi2(ND, ndchanged_x, ndchanged_y);
891  old_ndllike = ndllike;
892  }
893  else{
894  ndllike = old_ndllike;
895  }
896 
897  old_FD = FD;
898  old_ND = ND;
899 
900  llike = ndllike + fdllike;
901 
902  {
903  if(llike < lowest_llike){
904  lowest_llike = llike;
905  printf("lowest llike so far: %.9f%s\n", llike,
906  lowest_llike - llike < epsilonll?" (barely)":"");
907  lowest_llike = llike;
908  memcpy(bestpars, par, sizeof(double)*npar);
909  }
910  printf(".");
911  if(flag == 3) printf("\n");
912  fflush(stdout);
913  }
914  if(plot){
915  static int setn = 0;
916  thefitter->Plot(Form("intermediate%d", setn));
917  thefitter->Plot("intermediatelatest"); // inefficient, but ok...
918  setn++;
919  plot = false;
920  }
921 }
922 
923 // Get the zero-indexed value 'i'.
924 static double getpar(TMinuit & mn, const int i)
925 {
926  double val, err;
927  mn.GetParameter(i, val, err);
928  return val;
929 }
930 
931 // Get the zero-indexed error 'i'.
932 double geterr(TMinuit & mn, const int i)
933 {
934  double val, err;
935  mn.GetParameter(i, val, err);
936  return err;
937 }
938 
939 // Check if the value of the objective function at the current
940 // parameters is less than, equal to, or greater than the last time this
941 // function was called. If it greater, then reset the parameters to
942 // where they were before. This is useful because sometimes MIGRAD or
943 // SIMPLEX(?) get stuck and leave things worse than before.
944 void mincheck(TMinuit & mn, const bool reset = false)
945 {
946  static double oldllike = 0;
947 
948  if(oldllike == 0 || reset){
949  printf("\nllike = %8.1f\n", lowest_llike);
950  lowest_llike = 1e300;
951  memset(bestpars, 0, sizeof(double)*npar);
952  }
953  else if(oldllike < lowest_llike){
954  const int prevprint = mn.fISW[4];
955  mn.Command("SET PRINT -1");
956  if(prevprint >= 0)
957  printf("\nllike = %8.1f, %8.5f MORE than before (%8.1f). Resetting.\n",
958  lowest_llike, lowest_llike - oldllike, oldllike);
959  for(unsigned int i = 1; i <= npar; i++)
960  mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
961  mn.Command(Form("SET PRINT %d", prevprint));
962  return;
963  }
964  else if(lowest_llike == 0){
965  printf("\nllike = 0 which is nonsense. Resetting.\n");
966  for(unsigned int i = 1; i <= npar; i++)
967  mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
968  return;
969  }
970  else if(oldllike == lowest_llike){
971  printf("\nllike = %8.1f, which is the same as last time\n", lowest_llike);
972  }
973  else if(oldllike - lowest_llike > epsilonll){
974  printf("\nllike = %8.1f, which is %f less than last time\n",
975  lowest_llike, oldllike - lowest_llike);
976  }
977  else{
978  printf("\nllike = %8.1f, which is barely (%g) less than last time\n",
979  lowest_llike, oldllike - lowest_llike);
980  }
981 
982  if(lowest_llike < mn.fAmin){
983  const int prevprint = mn.fISW[4];
984  mn.Command("SET PRINT -1");
985  if(prevprint >= 0)
986  printf("\nNote: best = %f < fAmin = %f, using best pars instead\n",
987  lowest_llike, mn.fAmin);
988  for(unsigned int i = 1; i <= npar; i++)
989  mn.Command(Form("SET PAR %d %f\n", i, bestpars[i-1]));
990  mn.Command(Form("SET PRINT %d", prevprint));
991  }
992 
993  oldllike = mn.fAmin;
994  for(unsigned int i = 0; i < npar; i++) bestpars[i] = getpar(mn, i);
995 }
996 
997 void CherenkovFitter::Plot(const char * const tag)
998 {
999  for(CherenkovCalc * sample : _fdSamples){
1000  TwoDplots(*sample, Form("%s_%sfit", sample->plot_name, tag));
1001  Flatplots(*sample, Form("%s_3fitflat", sample->plot_name), tag);
1002  }
1003  for(CherenkovCalc * sample : _ndSamples){
1004  TwoDplots(*sample, Form("%s_%sfit", sample->plot_name, tag));
1005  Flatplots(*sample, Form("%s_3fitflat", sample->plot_name), tag);
1006  }
1007 }
1008 
1009 __attribute__((unused)) static void Console(TMinuit & mn)
1010 {
1011  std::string in;
1012  printf("MINIUT> ");
1013  while(getline(std::cin, in)){
1014  if(in == "quit" || in == "exit" || in == ".q") break;
1015  mn.Command(in.c_str());
1016  printf("MINIUT> ");
1017  }
1018 }
1019 
1020 static void setarand(TMinuit & mn, const int p, const double min,
1021  const double max)
1022 {
1023  const double v = min + (max-min)*R.Rndm();
1024  mn.Command(Form("SET PAR %d %f", p, v));
1025 }
1026 
1027 #ifndef FULL
1028 __attribute__((unused))
1029 #endif
1030 static void randomparameters(TMinuit & mn)
1031 {
1032  if(float_fdfx) setarand(mn, 1, 0.3, 0.6);
1033  if(float_fdfy) setarand(mn, 2, 0.3, 0.6);
1034  if(float_ndfx) setarand(mn, 3, 0.3, 0.6);
1035  if(float_ndfy) setarand(mn, 4, 0.3, 0.6);
1036 
1037  if(float_ec) setarand(mn, 6, 0.7, 0.9);
1038  if(float_attb) setarand(mn, 7, 0.1, 0.9);
1039  if(float_att1) setarand(mn, 8, 110, 400);
1040  if(float_att2) setarand(mn, 9, 110, 900);
1041  if(float_rngeff) setarand(mn, 10, 0.98, 1.0);
1042 }
1043 
1044 void CherenkovFitter::Fit(const params & mcparams)
1045 {
1046  sfdSamples = _fdSamples; // So the MINUIT FCN can get at them
1048 
1049  signal(SIGUSR1, intermediate_plot);
1050  thefitter = this; // so the signal handler can get at it. Oh boy.
1051 
1052  TMinuit mn(npar);
1053  mn.SetFCN(objfcn);
1054 
1055  mn.Command("SET ERR 0.5"); // Since we're using log likelihoods
1056 
1057  // Make SCAN plots smaller than default
1058  mn.Command("SET WIDTH 80");
1059  mn.Command("SET LINES 40");
1060 
1061  int ierr;
1062  mn.mnparm( 1 - 1, "FDFx", mcparams.FDFx, 0.01, 0.2, 2, ierr);
1063  mn.mnparm( 2 - 1, "FDFy", mcparams.FDFy, 0.01, 0.2, 2, ierr);
1064  mn.mnparm( 3 - 1, "NDFx", mcparams.NDFx, 0.01, 0, 2, ierr);
1065  mn.mnparm( 4 - 1, "NDFy", mcparams.NDFy, 0.01, 0, 2, ierr);
1066  mn.mnparm( 5 - 1, "Ys", mcparams.Ys , 10, 0, 0, ierr);
1067  mn.mnparm( 6 - 1, "Ec", mcparams.Ec , 0.1, 0, 0, ierr);
1068  mn.mnparm( 7 - 1, "attB", mcparams.attB, 0.01, 0, 1, ierr);
1069  mn.mnparm( 8 - 1, "att1", mcparams.att1, 10, 100, 2000, ierr);
1070 
1071  // Define this parameter to be the difference between the two attenuation
1072  // parameters. This is a handy way to keep them from becoming degenerate
1073  // and spoiling the fit
1074  mn.mnparm( 9 - 1, "Da12", mcparams.att2 - mcparams.att1, 10, 100, 2000, ierr);
1075  mn.mnparm(npar-1, "rngeff", mcparams.rngeff, 0.01, 0, 1.5, ierr);
1076  if(ierr > 0){
1077  fprintf(stderr, "Failed to set up parameters\n");
1078  _exit(1);
1079  }
1080 
1081  mn.Command("SET STRATEGY 2");
1082 
1083  // Reduce MINUIT's expectations of accuracy because the objective function is
1084  // not going to be very smooth since it involves events migrating over bin
1085  // boundaries. Experimentally, near a minimum, parameter changes of epsilon
1086  // result in FCN value changes of ~1E-6 of the FCN value.
1087  mn.Command("SET EPS 1E-6");
1088 
1089  mn.fGraphicsMode = false; // So SCAN prints a plot
1090 
1091  // Only three of {x, y, Ys} should be free. Which ones is arbitrary. It's
1092  // easiest to compare to previous results if we fix Ys.
1093  mn.Command("SET PAR 5 3151.04");
1094  mn.Command("FIX 5");
1095 
1096  if(!float_fdfx) mn.Command("FIX 1");
1097  if(!float_fdfy) mn.Command("FIX 2");
1098  if(!float_ndfx) mn.Command("FIX 3");
1099  if(!float_ndfy) mn.Command("FIX 4");
1100 
1101  if(!float_ec) mn.Command("FIX 6");
1102  if(!float_attb) mn.Command("FIX 7");
1103  if(!float_att1) mn.Command("FIX 8");
1104  if(!float_att2) mn.Command("FIX 9");
1105  if(!float_rngeff) mn.Command("FIX 10");
1106 
1107 #ifdef FULL
1108  mn.Command("CALL");
1109  for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
1110  for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();
1111  mn.Command("CALL");
1112 
1113  Plot("1ralsopre");
1114 
1115  // In case input parameters are crazy, do quick random search
1116  printf("WHAT WE ARE DOING: Trying some random parameters\n");
1117 
1118  // Start in random places instead of using any prior information
1119  #define NOHELP
1120 
1121  #ifndef NOHELP
1122  mincheck(mn, true);
1123  #endif
1124 
1125  const int prevprint = mn.fISW[4];
1126  mn.Command("SET PRINT -1");
1127  complain_about_problems = false;
1128  for(int i = 0; i < 1000; i++){
1129  randomparameters(mn);
1130  mn.Command("CALL");
1131  mincheck(mn,
1132  #ifdef NOHELP
1133  i == 0
1134  #else
1135  false
1136  #endif
1137  );
1138  }
1139  complain_about_problems = true;
1140  mn.Command(Form("SET PRINT %d", prevprint));
1141 
1142  for(int cornot = 0; cornot < 2; cornot++){
1143  if(float_ec){
1144  if(cornot == 0) mn.Command("FIX 6");
1145  if(cornot == 1) mn.Command("RELEASE 6");
1146  }
1147  if(cornot == 1 && !float_ec) break;
1148 
1149  printf("WHAT WE ARE DOING: Set masks at current value\n");
1150  mn.Command("CALL");
1151  for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
1152  for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();
1153  mn.Command("CALL");
1154 
1155  #if 0
1156  Console(mn);
1157  #endif
1158 
1159  mincheck(mn, true);
1160 
1161  printf("WHAT WE ARE DOING: Do a little SIMPLEX to get in a nice spot\n");
1162  mn.Command("SIMPLEX 100");
1163  mincheck(mn);
1164 
1165  printf("WHAT WE ARE DOING: See if MIGRAD helps\n");
1166  mn.Command("MIGRAD");
1167  mincheck(mn);
1168 
1169  for(int i = 0; i < 2; i++){
1170  printf("WHAT WE ARE DOING: Reset masks at current value and minimize again\n");
1171  mn.Command("CALL");
1172  bool changed = false;
1173  for (CherenkovCalc* calc : sfdSamples) changed |= calc->MakeMasks();
1174  for (CherenkovCalc* calc : sndSamples) changed |= calc->MakeMasks();
1175  mn.Command("CALL");
1176 
1177  if(!changed){
1178  printf("No change in masks, so no need to minimize again\n");
1179  break;
1180  }
1181 
1182  mincheck(mn, true);
1183 
1184  mn.Command("SIMPLEX");
1185  mincheck(mn);
1186  mn.Command("SEEK 300 8");
1187  mincheck(mn);
1188  mn.Command("MIGRAD");
1189  mincheck(mn);
1190  mn.Command("IMPROVE");
1191  mincheck(mn);
1192  }
1193 
1194  while(true){
1195  {
1196  const double oldmin = mn.fAmin;
1197 
1198  printf("WHAT WE ARE DOING: Scan all to make sure things are sensible\n");
1199  const double r = 0.02, lo = 1-r, hi = 1+r;
1200  if(float_fdfx) mn.Command(Form("SCAN 1 81 %f %f",getpar(mn,0)*lo,getpar(mn,0)*hi));
1201  if(float_fdfy) mn.Command(Form("SCAN 2 81 %f %f",getpar(mn,1)*lo,getpar(mn,1)*hi));
1202  if(float_ndfx) mn.Command(Form("SCAN 3 81 %f %f",getpar(mn,2)*lo,getpar(mn,2)*hi));
1203  if(float_ndfy) mn.Command(Form("SCAN 4 81 %f %f",getpar(mn,3)*lo,getpar(mn,3)*hi));
1204 
1205  if(float_ec) mn.Command(Form("SCAN 6 81 %f %f",getpar(mn,5)*lo,getpar(mn,5)*hi));
1206  if(float_attb) mn.Command(Form("SCAN 7 81 %f %f",getpar(mn,6)*lo,getpar(mn,6)*hi));
1207  if(float_att1) mn.Command(Form("SCAN 8 81 %f %f",getpar(mn,7)*lo,getpar(mn,7)*hi));
1208  if(float_att2) mn.Command(Form("SCAN 9 81 %f %f",getpar(mn,8)*lo,getpar(mn,8)*hi));
1209  if(float_rngeff)mn.Command(Form("SCAN 10 81 %f %f",getpar(mn,9)*lo,getpar(mn,9)*hi));
1210 
1211  const double newmin = mn.fAmin;
1212  const bool significantreduction = oldmin - newmin > epsilonll;
1213  if(newmin < oldmin)
1214  printf("SCANs reduced llike by %g%s\n", oldmin - newmin,
1215  significantreduction?"":" -- not significant");
1216  if(!significantreduction) break;
1217  }
1218 
1219  {
1220  printf("WHAT WE ARE DOING: found new min, so SIMPLEX, SEEK, MIGRAD, IMPROVE\n");
1221  const double oldmin = mn.fAmin;
1222  mn.Command("SIMPLEX");
1223  mincheck(mn);
1224  mn.Command("SEEK 300 8");
1225  mincheck(mn);
1226  mn.Command("MIGRAD");
1227  mincheck(mn);
1228  mn.Command("IMPROVE");
1229  mincheck(mn);
1230  const double newmin = mn.fAmin;
1231  const bool significantreduction = oldmin - newmin > epsilonll;
1232  if(newmin < oldmin)
1233  printf("SIMPLEX, etc. reduced llike by %g%s\n", oldmin - newmin,
1234  significantreduction?"":" -- not significant");
1235  if(!significantreduction) break;
1236  }
1237  }
1238  }
1239 
1240  printf("WHAT WE ARE DOING: Call once more at min to make plots right\n");
1241  mn.Command("CALL");
1242 #else
1243  mn.Command("SET PAR 1 6.41795e-01");
1244  mn.Command("SET PAR 2 5.53340e-01");
1245  mn.Command("SET PAR 3 5.09120e-01");
1246  mn.Command("SET PAR 4 4.94346e-01");
1247  mn.Command("SET PAR 5 3.15104e+03");
1248  mn.Command("SET PAR 6 8.10680e-01");
1249  mn.Command("SET PAR 7 4.46867e-01");
1250  mn.Command("SET PAR 8 2.03410e+02");
1251  mn.Command("SET PAR 9 5.48778e+02");
1252  mn.Command("SET PAR 10 1.00000e+00 ");
1253 
1254  const double besta1 = getpar(mn, 8-1), bestda = getpar(mn, 9-1);
1255 
1256  printf("WHAT WE ARE DOING: Reset masks at current value\n");
1257  mn.Command("CALL");
1258  for (CherenkovCalc* calc : sfdSamples) calc->MakeMasks();
1259  for (CherenkovCalc* calc : sndSamples) calc->MakeMasks();
1260 
1261  // Helps?
1262  mincheck(mn, true);
1263  mn.Command("SIMPLEX 30");
1264  mincheck(mn);
1265 #if 0
1266  mn.Command("MIGRAD");
1267  mincheck(mn);
1268 #endif
1269 
1270  mn.Command("FIX 8 9");
1271 
1272  mn.Command("FIX 5"); // Fix Ys as in the regular fit
1273  mn.Command("FIX 10"); // not messing with this at the moment
1274  mn.Command("SET ERR 5"); // just to cut off SIMPLEX earlier without
1275  // changing the number of calls
1276 
1277  const int nside = 28;
1278  const double a1step = 0.003, dastep = 0.005/4;
1279  printf("contcenter2: %f %f %d %f %f\n",
1280  besta1, bestda, nside, a1step, dastep);
1281  for(int i = nside; i >= -nside; i--){
1282  for(int j = -nside; j <= nside; j++){
1283  const double a1rat = 1 + i*a1step;
1284  const double darat = 1 + j*dastep;
1285  const double a1 = besta1*a1rat;
1286  const double da = bestda*darat;
1287  mn.Command("SET PAR 1 6.41795e-01");
1288  mn.Command("SET PAR 2 5.53340e-01");
1289  mn.Command("SET PAR 3 5.09120e-01");
1290  mn.Command("SET PAR 4 4.94346e-01");
1291  mn.Command("SET PAR 5 3.15104e+03");
1292  mn.Command("SET PAR 6 8.10680e-01");
1293  mn.Command("SET PAR 7 4.46867e-01");
1294 
1295  mn.Command(Form("SET PAR 8 %f", a1));
1296  mn.Command(Form("SET PAR 9 %f", da));
1297 
1298  mn.Command("FIX 1 2 3 4 5 6 7 8 9");
1299 
1300  // Have to do this every time, because sometimes MINUIT
1301  // decides 0 isn't good enough, and switches to 1, and it
1302  // is a permanent switch!
1303  mn.Command("SET STRATEGY 0");
1304  mn.Command("CALL");
1305 #if 0
1306  printf("cont: %10f %10f %10f\n", a1, da, mn.fAmin);
1307 #endif
1308  mincheck(mn, true);
1309 #if 1
1310  mn.Command("RELEASE 7");
1311  mn.Command("SIMPLEX 10");
1312  mn.Command("MIGRAD 100");
1313  mincheck(mn);
1314 #endif
1315 #if 1
1316  printf("cont2: %10f %10f %10f\n", a1, da, mn.fAmin);
1317 #endif
1318 #if 0
1319  mn.Command("RELEASE 1 2 3 4");
1320  mn.Command("SIMPLEX 10");
1321  //mn.Command("MIGRAD 50");
1322  mincheck(mn);
1323 
1324  mn.Command("RELEASE 6");
1325  mn.Command("SIMPLEX 10");
1326  mn.Command("MIGRAD 50");
1327  mincheck(mn);
1328 
1329  printf("cont3: %10f %10f %10f\n", a1, da, lowest_llike);
1330 #endif
1331  }
1332  }
1333 
1334  Console(mn);
1335 
1336  #if 0
1337  mincheck(mn, true);
1338  mn.Command("SIMPLEX");
1339 
1340  printf("WHAT WE ARE DOING: Scan all to make sure things are sensible\n");
1341  mn.Command(Form("SCAN 1 21 %f %f", getpar(mn, 0)*0.98, getpar(mn, 0)*1.02));
1342  mn.Command(Form("SCAN 2 21 %f %f", getpar(mn, 1)*0.98, getpar(mn, 1)*1.02));
1343  mn.Command(Form("SCAN 3 21 %f %f", getpar(mn, 2)*0.98, getpar(mn, 2)*1.02));
1344  mn.Command(Form("SCAN 4 21 %f %f", getpar(mn, 3)*0.98, getpar(mn, 3)*1.02));
1345  mn.Command(Form("SCAN 6 21 %f %f", getpar(mn, 5)*0.98, getpar(mn, 5)*1.02));
1346  mn.Command(Form("SCAN 7 21 %f %f", getpar(mn, 6)*0.98, getpar(mn, 6)*1.02));
1347  mn.Command(Form("SCAN 8 21 %f %f", getpar(mn, 7)*0.98, getpar(mn, 7)*1.02));
1348  mn.Command(Form("SCAN 9 21 %f %f", getpar(mn, 8)*0.98, getpar(mn, 8)*1.02));
1349  mn.Command(Form("SCAN 10 21 %f %f", getpar(mn, 9)*0.98, getpar(mn, 9)*1.02));
1350 
1351  printf("WHAT WE ARE DOING: Call once more at min to make plots right\n");
1352  mn.Command("CALL");
1353  #endif
1354 #endif
1355 
1356  printf("\n\nFD Fx: %.4f Fy: %.4f\n"
1357  "ND Fx: %.4f Fy: %.4f\n"
1358  "Ys: %.2f (by definition)\n"
1359  "Ec: %.4f\n"
1360  "attB: %6.4f att1: %6.1f att2: %6.1f\n",
1361  getpar(mn, 0), getpar(mn, 1), getpar(mn, 2), getpar(mn, 3),
1362  getpar(mn, 4), getpar(mn, 5), getpar(mn, 6), getpar(mn, 7),
1363  getpar(mn, 7) + getpar(mn, 8));
1364 
1365  Plot("2post");
1366 }
static void setarand(TMinuit &mn, const int p, const double min, const double max)
double FDFx
Definition: CherenkovCalc.h:12
std::vector< double > pebins
static bool fill_mask(TH3D *mask, TH3D *h)
detparams _Shift
const double MINPATHLENGTH
static void objfcn(__attribute__((unused)) int &np, __attribute__((unused)) double *gin, double &llike, double *par, int flag)
double Ys
Definition: CherenkovCalc.h:21
TSpline3 lo("lo", xlo, ylo, 12,"0")
static TRandom3 R(0)
const double fdwmax
Definition: consts.h:2
double geterr(TMinuit &mn, const int i)
diblock
print "ROW IS " print row
Definition: geo2elec.py:31
void FastMCFill(TH3D *hist, const std::vector< dat > &dats, const bool isX)
void mincheck(TMinuit &mn, const bool reset=false)
int distshort
Definition: CherenkovCalc.h:83
std::vector< CherenkovCalc * > _ndSamples
const bool float_ec
#define NOHELP
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
def fail(msg)
double inv_att1
Definition: CherenkovCalc.h:21
unsigned int ny
double Ec
Definition: CherenkovCalc.h:21
static binw smearbin(const vector< double > &bins, const double val, const int bin)
int distlong
Definition: CherenkovCalc.h:84
unsigned int nz
const Var weight
double FDFy
Definition: CherenkovCalc.h:12
void Plot(const char *const tag)
static TH3D * invmask(TH3D *mask)
double inv_att2
Definition: CherenkovCalc.h:21
void AddNDSample(CherenkovCalc *calc)
static bool compare_dat(const dat &a, const dat &b)
void SetInitialParams(const params &mcparams, const enum DET det)
static double bestpars[npar]
const char * p
Definition: xmltok.h:285
double NDFx
Definition: CherenkovCalc.h:12
static bool goodfdplane(const int plane)
Int_t par
Definition: SimpleIterate.C:24
double att1
Definition: CherenkovCalc.h:12
static double llike_hist(const TH3D *const hObs, const TH3D *const hMC, const TH2D *const hObsDist, const TH2D *const hMCDist, const TH3D *const mask)
static vector< CherenkovCalc * > sfdSamples
void AddFDSample(CherenkovCalc *calc)
void abs(TH1 *hist)
CherenkovCalc(const char *const dataFileName, const char *const mcFileName, const char *const dirName, const char *const plot_name, const std::vector< double > &binX, const std::vector< double > &binY, const std::vector< double > &binZ, const bool isFD_)
double rngeff
Definition: CherenkovCalc.h:21
static const int stat_req
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 bool float_ndfy
const double kPigtails[kCellsPerModule]
static const unsigned int npar
const char * _dirName
double dist
Definition: runWimpSim.h:113
const XML_Char const XML_Char * data
Definition: expat.h:268
bool proton
Definition: CherenkovCalc.h:82
static double clamp(const double v, const double mi, const double ma)
double prev
std::vector< double > atten_cache_long
unsigned int nx
TSpline3 hi("hi", xhi, yhi, 18,"0")
double peWeight(const bool isX, const dat &d)
double Ys
Definition: CherenkovCalc.h:12
const double atten_step
void prev()
Definition: show_event.C:91
static const unsigned int maxDATAhitsPERw
void init_atten_cache()
static const int kCellsPerModule
const double ndwmin
Definition: consts.h:3
static const unsigned int maxMChitsPERw
double Fx
Definition: CherenkovCalc.h:21
double exact_atten(const double dist, const bool longway)
std::vector< double > wbins
static bool goodfdcell(__attribute__((unused)) const bool isX, __attribute__((unused)) const int cell)
TH1F * a1
Definition: f2_nu.C:476
static void saneProjectionXZ(TH2D *out, TH3D *in)
static CherenkovFitter * thefitter
const double a
const bool float_fdfy
const double ndwmax
Definition: consts.h:3
int DistBin
Definition: CherenkovCalc.h:80
const double maxpigtail
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
CherenkovFitter(const params &mcparams)
printf("%d Experimental points found\n", nlines)
double rngeff
Definition: CherenkovCalc.h:12
const bool float_att1
bool accepthit(const bool isFD, const bool isX, const int plane, const int cell, const float pathlength, const float disttostop, const float w, const float PE)
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
const Binning bins
static vector< CherenkovCalc * > sndSamples
#define SMEAR
DET
Definition: CherenkovCalc.h:8
double NDFy
Definition: CherenkovCalc.h:12
float bin[41]
Definition: plottest35.C:14
z
Definition: test.py:28
static double _exact_atten(const double attB, const double inv_att1, const double inv_att2, const double dist, const double rngeff, const bool longway)
const bool float_fdfx
const bool float_rngeff
double getpigtail(const int cell, const bool isX)
double Ec
Definition: CherenkovCalc.h:12
bool fill_no_overflow(TH3D *h, const double x, const double y, const double z, const double weight=1)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
std::vector< dat > datsx
__attribute__((unused)) static void Console(TMinuit &mn)
static double lowest_llike
double curr
bool equalsx(const detparams &b)
Definition: CherenkovCalc.h:47
std::vector< CherenkovCalc * > _fdSamples
void fill_hit_cache_and_init_things()
Float_t iX
ifstream in
Definition: comparison.C:7
unsigned int npebins
std::string dirName
Definition: PlotSpectra.h:47
double next
TDirectory * dir
Definition: macro.C:5
float InvTotPhoAtten
Definition: CherenkovCalc.h:79
static const int stat_warn_level
static bool complain_about_problems
static unsigned int getpebin(const double pe, const std::vector< double > &bins)
void Flatplots(CherenkovCalc &calc, const char *const title1, const char *const title2)
Definition: Plot.cxx:173
void TwoDplots(CherenkovCalc &calc, const char *const title)
Definition: Plot.cxx:27
std::vector< double > atten_cache_short
const hit & b
Definition: hits.cxx:21
detparams _Init
const bool float_ndfx
TFile * file
Definition: cellShifts.C:17
double att2
Definition: CherenkovCalc.h:12
TRandom3 r(0)
float ScintPho
Definition: CherenkovCalc.h:78
float PE
Definition: CherenkovCalc.h:76
static double getpar(TMinuit &mn, const int i)
const bool float_attb
static void Add(TH3D *h, const int bx, const int by, const int bz, const double w)
std::vector< double > distbins
TFile * _dataFile
std::vector< double > extpebins
bool ok_w(const float w, const bool isX)
const double fdwmin
Definition: consts.h:2
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
static const double epsilonll
static void randomparameters(TMinuit &mn)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
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
void next()
Definition: show_event.C:84
void FillMCHistos(const detparams &mcparams, const bool changed_x, const bool changed_y)
double cellhalflength
unsigned int ndistbins
const bool float_att2
double CalcChi2(const detparams &mcparams, const bool changed_x, const bool changed_y)
void Fit(const params &mcparams)
static TTreeReader * makereader(const char *const dir, TFile *file)
double Fy
Definition: CherenkovCalc.h:21
static bool plot
unsigned int nwbins
static void intermediate_plot(const int signal)
enum BeamMode string