AttenFit_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file AttenFit_module.cc
3 // \brief Fit attenuation behaviour from AttenHists files
4 // \author Christopher Backhouse, Dan Pershey - bckhouse@caltech.edu, pershey@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 // C/C++ includes
8 #include <cmath>
9 #include <iostream>
10 
11 // ROOT includes
12 #include "TFile.h"
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TF1.h"
16 #include "TString.h"
17 #include "TGraph.h"
18 #include "TCanvas.h"
19 #include "TLegend.h"
20 #include "TProfile.h"
21 #include "TFitResult.h"
22 #include "TFitResultPtr.h"
23 #include "TStyle.h"
24 #include "TSystem.h"
25 
26 // Framework includes
32 
33 // NOvASoft includes
36 #include "Geometry/Geometry.h"
41 #include "RecoBase/CellHit.h"
42 #include "RecoBase/Track.h"
44 #include "Utilities/func/MathUtil.h"
45 #include "NovaDAQConventions/DAQConventions.h"
47 
48 namespace calib
49 {
50  // Parameters for fits, should be moved to .fcl file
51  int kStatsLimit = 1;
52  //const double kTrueShortAtten = 289.5;
53 
54  class AttenFit: public art::EDAnalyzer
55  {
56  public:
57  explicit AttenFit(const fhicl::ParameterSet& pset);
58  ~AttenFit();
59 
60  void reconfigure(const fhicl::ParameterSet& pset);
61  virtual void beginRun(const art::Run& run);
62  virtual void analyze(const art::Event& evt);
63  virtual void endSubRun(const art::SubRun& sr);
64  virtual void endRun(const art::Run& run);
65  virtual void endJob();
66 
67  protected:
68  int MinBlock() const;
69  int MaxBlock() const;
71  void ConsolidateViews();
72 
73  caldp::AttenHists GetChannelHists(bool isData, int plane, int cell);
74  caldp::AttenProfiles GetChannelProfiles(bool isData, int plane, int cell);
75 
77  TH1* MedianProfile(TH2* f);
78  TH1* MeanProfile(TH2* f);
79  TH1* TruncatedMeanProfile(TH2* h);
80  void ApplyThresholdCorrection(int view, int cell, TH1* prof);
81  TF1* expFit(TH1* prof, TString opt, geo::OfflineChan chan, calib::AttenCurve*& res);
82  TF1* rolloffFit(TH1* prof, TString opt, geo::OfflineChan chan, calib::AttenCurve*& res);
83  TGraph* lowessFit(TH1* prof, int plane, calib::AttenCurve*& res);
84  void fit_channel_prof(TH1* atten, int plane, int cell, calib::AttenCurve*& res);
85  void DrawDetectorEdges(double range);
86  void DrawDetectorEdges(double rangedown, double rangeup);
87 
88  double FitQuality(TH1* prof, TF1* fit) const;
89 
90  /// Ensure every cell has a CSV entry even if it couldn't be fit
91  void WriteDummyCSVRow(FILE* f, int plane, int cell) const;
92 
104 
105  bool fIsMC;
106 
108 
109  /// This is where \ref GetChannelHist stores its histograms. Map is from block
111  std::map<int, caldp::AttenProfilesMap> fChannelMapProf;
112  };
113 
114  //......................................................................
116  : EDAnalyzer(pset)
117  {
118  reconfigure(pset);
119 
120  fThreshCorrMap = 0;
121 
122  gStyle->SetOptStat(0);
123  }
124 
125  //......................................................................
127  {
128  delete fThreshCorrMap;
129  }
130 
131  //......................................................................
133  {
134  fAttenHistsMapLabel = pset.get<std::string>("AttenHistsMapLabel");
135  fThresholdCorrMapFile = pset.get<std::string>("ThresholdCorrMapFile");
136  fPlotsDirectory = pset.get<std::string>("PlotsDirectory");
137  fProfileMode = pset.get<bool>("ProfileMode");
138  fConsolidateViewsMode = pset.get<bool>("ConsolidateViewsMode");
139  fIncludeRolloffFit = pset.get<bool>("IncludeRolloffFit");
140  fIncludeLowessFit = pset.get<bool>("IncludeLowessFit");
141  fPlaneMin = pset.get<int>("PlaneMin");
142  fPlaneMax = pset.get<int>("PlaneMax");
143  fLookInSubRuns = pset.get<bool>("LookInSubRuns");
144  }
145 
146  //......................................................................
148  {
150  fInitEdgeX = geom->DetHalfWidth();
151  fInitEdgeY = geom->DetHalfHeight();
152  if(geom->DetId() == novadaq::cnv::kNEARDET){
153  fInitEdgeX-=50;
154  fInitEdgeY-=50;
155  fInitEdgeMuXMin=-150;
156  fInitEdgeMuXMax=50;
157  }
158 
160  fIsMC = (rh->GetDataType() == nova::dbi::kMCOnly);
162  }
163 
164  //......................................................................
166  {
167  // If we get given any events, make sure they match what we guessed for the
168  // run.
169  assert(fIsMC != evt.isRealData());
170  }
171 
172  //......................................................................
174  {
175  if(!fLookInSubRuns) return;
176 
177  if(fProfileMode){
178  for(int block = MinBlock(); block <= MaxBlock(); ++block){
180  sr.getByLabel(fAttenHistsMapLabel, TString::Format("block%d", block).Data(), profsmap);
181  if(profsmap.failedToGet()) continue;
182 
183  // Accumulate across runs
184  if(fChannelMapProf.find(block) == fChannelMapProf.end())
185  fChannelMapProf.insert(std::make_pair(block, caldp::AttenProfilesMap(profsmap->MinW(), profsmap->MaxW())));
186  auto it = fChannelMapProf.find(block);
187  it->second += *profsmap;
188 
189  // No one else will ever want this. Save the memory
190  sr.removeCachedProduct(profsmap);
191  }
192  }
193  else{
195  sr.getByLabel(fAttenHistsMapLabel, histsmap);
196 
197  // Accumulate across runs
198  fChannelMap += *histsmap;
199  }
200  }
201 
202  //......................................................................
204  {
205  // Profiles could be stored in subruns or runs, depending whether summing
206  // has been done.
207  for(int block = MinBlock(); block <= MaxBlock(); ++block){
209  run.getByLabel(fAttenHistsMapLabel, TString::Format("block%d", block).Data(), profsmap);
210 
211  if(profsmap.failedToGet()) continue;
212 
213  // Accumulate across runs
214  if(fChannelMapProf.find(block) == fChannelMapProf.end())
215  fChannelMapProf.insert(std::make_pair(block, caldp::AttenProfilesMap(profsmap->MinW(), profsmap->MaxW())));
216  auto it = fChannelMapProf.find(block);
217  it->second += *profsmap;
218 
219  // No one else will ever want this. Save the memory
220  run.removeCachedProduct(profsmap);
221  }
222  }
223 
224  //......................................................................
226  {
228 
230 
231  TFile* outfile = new TFile("Atten_Fit_Results.root","RECREATE");
232 
233  unsigned int plane_max = fPlaneMax;
234  bool isND=geom->DetId() == novadaq::cnv::kNEARDET;
235  if(fPlaneMax < 0) plane_max = geom->NPlanes()-1;
236  if(fConsolidateViewsMode) plane_max = 1;
237  if(fConsolidateViewsMode&&isND) plane_max = 3;
238 
239  const TString mcStr = fIsMC ? "_mc" : "";
240  FILE* fConstsCSV = fopen("calib_atten_consts"+mcStr+".csv","w");
241  FILE* fPointsCSV = fopen("calib_atten_points"+mcStr+".csv","w");
242 
243  fprintf(fConstsCSV, "#coeff_exp,atten_length,background,edge_low,edge_high,coeff_low,coeff_high,chisq,channel,tv\n");
244  fprintf(fPointsCSV, "#w,factor,channel,tv\n");
245 
246  for (unsigned int plane = fPlaneMin; plane <= plane_max; plane++){
247  mf::LogVerbatim("AttenFit") << "Fitting plane "
248  << plane << "/" << plane_max;
249  const int block = plane/32;
250  int view = GetView(plane,geom);
251  if(fConsolidateViewsMode) view = plane;
252  auto it = fChannelMapProf.find(block);
253  std::cout<<"it == fChannelMapProf.end()"<<(it == fChannelMapProf.end())<<std::endl;
254  if(it == fChannelMapProf.end()) continue;
255  caldp::AttenProfilesMap& profmap = it->second;
256 
257  TDirectory* dir = 0;
258  const unsigned int cell_max = geom->Plane(plane)->Ncells();
259  std::cout<<"CellMax: "<<cell_max<<std::endl;
260  for(unsigned int cell = 0; cell < cell_max; ++cell){
261  mf::LogVerbatim("AttenFit") << " Fitting cell "
262  << cell << "/" << cell_max;
263 
264  TH1* prof = 0;
265  if(fProfileMode){
266  if(!profmap.HasProfiles(plane, cell)){
267  mf::LogWarning("AttenFit") << "Nothing to fit for plane "
268  << plane << ", cell " << cell;
269  WriteDummyCSVRow(fConstsCSV, plane, cell);
270  continue;
271  }
273 
274  prof = v.WPE_corr_xy.ToTH1F();
275  // if(plane==2&&cell==33){
276 
277  // std::cout<<"Integral "<<prof->Integral()<<std::endl;
278  // std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl;
279  // std::cout<<"GetEffectiveEntries "<<prof->GetEffectiveEntries()<<std::endl;
280 
281  // std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl;
282 
283  // std::cout<<"TotalEntries z "<<v.WPE_corr_z.TotalEntries()<<std::endl;
284 
285  // std::cout<<"TotalEntries "<<v.WPE_corr.TotalEntries()<<std::endl;
286 
287  // assert(0);
288  // }
289 
290  // TODO: figure out where these bad axis limits come from
291 
292  if(prof->GetXaxis()->GetXmin() == 0 ||
293  prof->GetXaxis()->GetXmax() == 0 ||
295  // std::cout<<"Stats Fail "<<v.WPE_corr_xy.TotalEntries()<<" "<<kStatsLimit<<std::endl;
296 
297  // std::cout<<"Integral "<<prof->Integral()<<std::endl;
298  // std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl;
299  // std::cout<<"GetEffectiveEntries "<<prof->GetEffectiveEntries()<<std::endl;
300 
301  // std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl;
302 
303  // std::cout<<"TotalEntries z "<<v.WPE_corr_z.TotalEntries()<<std::endl;
304 
305  // std::cout<<"TotalEntries "<<v.WPE_corr.TotalEntries()<<std::endl;
306 
307  prof = v.WPE_corr_z.ToTH1F();
308 
309  if(prof->GetXaxis()->GetXmin() == 0 ||
310  prof->GetXaxis()->GetXmax() == 0 ||
312  std::cout<<"Stats Double Fail "<<v.WPE_corr_z.TotalEntries()<<" "<<kStatsLimit<<std::endl;
313  WriteDummyCSVRow(fConstsCSV, plane, cell);
314  //continue;
315  }
316  }
317  }
318  else{
321  TH2* h_WPE = GetBestHistogram(v,u);
322  if(!h_WPE){
323  WriteDummyCSVRow(fConstsCSV, plane, cell);
324  continue;
325  }
326  prof = MedianProfile(h_WPE);
327  }
328 
329  // Never evaluated for NDOS
330  if(geom->DetId() == novadaq::cnv::kFARDET ||
331  geom->DetId() == novadaq::cnv::kNEARDET){
332  ApplyThresholdCorrection(view, cell, prof);
333  }
334 
336  fit_channel_prof(prof, plane, cell, res);
337  outfile->cd();
338  if(!dir) dir = outfile->mkdir(Form("plane_%03d", plane));
339  dir->cd();
340  prof->Write(Form("Atten_Fit_%03d_%03d",plane,cell));
341  // Encoding to indicate it's the backup averaged constants
343  res->chan = geo::OfflineChan(1000+view, cell);
344  res->WriteToCSVs(fConstsCSV, fPointsCSV, fIsMC);
345 
346  delete res;
347  delete prof;
348  }
349  }
350  fclose(fConstsCSV);
351  fclose(fPointsCSV);
352  }
353 
354  //......................................................................
355  int AttenFit::MinBlock() const
356  {
357  return fPlaneMin/32;
358  }
359 
360  //......................................................................
361  int AttenFit::MaxBlock() const
362  {
363  if(fPlaneMax < 0) return 27;
364  return fPlaneMax/32;
365  }
366  //......................................................................
367  int AttenFit::GetView(int plane, art::ServiceHandle<geo::Geometry> geom){//get the view for the current combination of view and plane (x,y,mux, or muy)
368  int view=0;
369  if(geom->DetId() == novadaq::cnv::kNEARDET){
370  int firstPlane=geom->FirstPlaneInMuonCatcher();
371  int isMuonCatcher=firstPlane <= plane;
372  view = geom->Plane(plane)->View()+(2*isMuonCatcher);
373  }
374  else{
375  view = geom->Plane(plane)->View();
376  }
377  return view;
378  }
379  //......................................................................
381  {
383 
384  // Has consolidation already been done by an earlier step?
385  bool isConsolidated = true;
386  bool isND=geom->DetId() == novadaq::cnv::kNEARDET;
387  if(fChannelMapProf.find(0) != fChannelMapProf.end()){
388  const std::vector<geo::OfflineChan> chans = fChannelMapProf.find(0)->second.Channels();
389  for(geo::OfflineChan chan: chans){
390  if((chan.Plane() > 3 && isND) || (chan.Plane() > 1 && isND==false)){
391  isConsolidated = false;
392  break;
393  }
394  }
395  }
396 
397  if(isConsolidated){
398  mf::LogInfo("AttenFit") << "Input file appears to already be consolidated by view. Skipping that step.";
399  return;
400  }
401 
402  // Need to do it for ourselves
403  mf::LogInfo("AttenFit") << "Input file is not consolidated by view. Consolidating";
404 
405  // All x-planes will go into plane 0 of this, y to 1.
407 
408  for(int block = MinBlock(); block <= MaxBlock(); ++block){
409  auto it = fChannelMapProf.find(block);
410  if(it == fChannelMapProf.end()) continue;
411  caldp::AttenProfilesMap& profmap = it->second;
412  const std::vector<geo::OfflineChan> chans = profmap.Channels();
413 
414  for(geo::OfflineChan chan: chans){
415  const caldp::AttenProfiles& profs = profmap.GetProfiles(chan);
416  if(!sum) sum = new caldp::AttenProfilesMap(profmap.MinW(),
417  profmap.MaxW());
418 
419  const int view = GetView(chan.Plane(),geom);
420 
421  sum->GetProfiles(view, chan.Cell()) += profs;
422  } // end for chan
423 
424  profmap.Clear();
425  } // end for block
426 
427  fChannelMapProf.insert(std::make_pair(0, *sum));
428  delete sum;
429  }
430 
431  //......................................................................
432  double AttenFit::FitQuality(TH1* prof, TF1* fit) const
433  {
434  // chisq lets bad fits to low stats data through, and cuts out good(ish)
435  // fits to high stats data. Try this metric, mean fractional deviation of
436  // the data from the fit (in quadrature).
437 
438  double xmin, xmax;
439  fit->GetRange(xmin, xmax);
440  int tot = 0;
441  double diff = 0;
442 //std::cout << "xmin = " << xmin << " , xmax = " << xmax << std::endl;
443  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
444  const double x = prof->GetXaxis()->GetBinCenter(i);
445 //std::cout << "i = " << i << " , x = " << x << " ,xmin = " << xmin << " , xmax = " << xmax << " , Fit value = " << fit->Eval(x) << " , Bin content = " << prof->GetBinContent(i) << ", tot = " << tot << " and diff = " << diff << std::endl;
446  if(x < xmin || x > xmax) {
447 // std::cout << "i = " << i << " , x = " << x << " ,xmin = " << xmin << " , xmax = " << xmax << " , Fit value = " << fit->Eval(x) << " , Bin content = " << prof->GetBinContent(i) << ", tot = " << tot << " and diff = " << diff << std::endl;
448  continue;
449  }
450  ++tot;
451  diff += util::sqr((fit->Eval(x)-prof->GetBinContent(i))/fit->Eval(x));
452  }
453 
454  const double ret = sqrt(diff/tot);
455  if(std::isnan(ret) || std::isinf(ret)) {
456 // std::cout << "xmin = " << xmin << " , xmax = " << xmax << " , tot = " << tot << " , diff = " << diff << " and ret = " << ret << std::endl;
457  return 2;
458  }
459  return ret;
460  }
461 
462  //......................................................................
464  {
465  if(!v.WPE) return 0;
466  if(u.WPE_corr_xy.TotalEntries() > kStatsLimit) return v.WPE_corr_xy;
467  if(u.WPE_corr_z.TotalEntries() > kStatsLimit) return v.WPE_corr_z;
468  if(u.WPE_corr.TotalEntries() > kStatsLimit) return v.WPE_corr;
469  return 0;
470  }
471 
472  //......................................................................
474  {
475  TProfile* p = h->ProfileX();
476 
477  TAxis* xax = h->GetXaxis();
478  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
479  ret->GetXaxis()->SetTitle(xax->GetTitle());
480  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
481  ret->SetBinContent(ix, p->GetBinContent(ix));
482  ret->SetBinError(ix, p->GetBinError(ix));
483  }
484  return ret;
485  }
486 
487  //......................................................................
489  {
491  TAxis* xax = h->GetXaxis();
492  TAxis* yax = h->GetYaxis();
493  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
494  ret->GetXaxis()->SetTitle(xax->GetTitle());
495  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
496  double tot = 0;
497  // Make sure to include overflow bins
498  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
499  tot += h->GetBinContent(ix, iy);
500  }
501  double seen = 0;
502  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
503  seen += h->GetBinContent(ix, iy);
504  if(seen > tot/2){
505  const double z = h->GetBinContent(ix, iy);
506  // Interpolate
507  const double frac = 1-(seen-tot/2)/z;
508  ret->SetBinContent(ix, yax->GetBinLowEdge(iy)+frac*yax->GetBinWidth(iy));
509  // Intro of http://www.jstor.org/stable/2286545
510  ret->SetBinError(ix, sqrt(tot)*yax->GetBinWidth(iy)/(2*z));
511  break;
512  }
513  } // end for iy
514  } // end for ix
515  return ret;
516  }
517 
518  //......................................................................
520  {
521  TAxis* xax = h->GetXaxis();
522  TAxis* yax = h->GetYaxis();
523  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
524  ret->GetXaxis()->SetTitle(xax->GetTitle());
525  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
526  double tot = 0;
527  // Make sure to include overflow bins
528  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
529  tot += h->GetBinContent(ix, iy);
530  }
531 
532  double avg = 0;
533  double rms = 0;
534  double seen = 0;
535  double inc = 0;
536  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
537  const double z = h->GetBinContent(ix, iy);
538  if(!z) continue;
539 
540  seen += z;
541 
542  // Interpolation
543  double weight_lo = (seen-(.25*tot))/z;
544  double weight_hi = 1-(seen-(.75*tot))/z;
545  if(weight_lo > 1) weight_lo = 1;
546  if(weight_hi > 1) weight_hi = 1;
547  double weight = weight_lo+weight_hi-1;
548  if(weight < 0) weight = 0;
549 
550  inc += z*weight;
551  const double c = yax->GetBinCenter(iy);
552  avg += z*weight*c;
553  rms += z*weight*util::sqr(c);
554  } // end for iy
555 
556  // In principle this should be problematic because of floating point errors
557  assert(inc == tot/2);
558 
559  if(inc){
560  avg /= inc;
561  rms /= inc;
562  rms -= util::sqr(avg);
563  assert(!std::isinf(avg) && !std::isnan(avg));
564  assert(!std::isinf(rms) && !std::isnan(rms));
565  assert(rms >= 0);
566  ret->SetBinContent(ix, avg);
567  ret->SetBinError(ix, sqrt(rms)/sqrt(tot/2));
568  }
569  else{
570  ret->SetBinContent(ix, 0);
571  ret->SetBinError(ix, 1e10);
572  }
573  }
574 
575  return ret;
576  }
577 
578  //......................................................................
579  void AttenFit::ApplyThresholdCorrection(int view, int cell, TH1* prof)
580  {
581 //std::cout << "////////////////////////////// Apply Threshold correction /////////////////////////////" << std::endl;
582  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
583  const double w = prof->GetXaxis()->GetBinCenter(i);
584  const double y = prof->GetBinContent(i);
585  const double corr = fThreshCorrMap->ThresholdCorrAt(view, cell, w);
586 //std::cout << "i = " << i << " , Bin center = " << w << " , Bincontent = " << y << " , threshold corr = " << corr << " , view = " << view << " , cell " << cell << " , new Bin content " << y/corr << " , old bin error = "<< prof->GetBinError(i) << " and new bin error " << prof->GetBinError(i)/corr << std::endl;
587  prof->SetBinContent(i, y/corr);
588  prof->SetBinError(i, prof->GetBinError(i)/corr);
589  }
590 //std::cout << "////////////////////////////// Done Threshold correction /////////////////////////////" << std::endl;
591  }
592 
593  //......................................................................
594  TF1* AttenFit::expFit(TH1* prof, TString opt, geo::OfflineChan chan,
595  calib::AttenCurve*& res)
596  {
598 
599  int view = GetView(chan.Plane(),geom);
600  double range = geom->Plane(chan.Plane())->Cell(chan.Cell())->HalfL();
602  view = (chan.Plane());
603  range = (view == geo::kX) ? geom->DetHalfWidth() : geom->DetHalfHeight();
604  }
605  res->cell_length = 2*range;
606  if(geom->DetId() == novadaq::cnv::kNEARDET){
607  range -= 50;
608  }
609  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
610  res->cell_length =range+50+fInitEdgeMuXMax;
611  }
612  res->coeff_low = 0;
613  res->coeff_high = 0;
614  res->edge_low = 0;
615  res->edge_high = 0;
616 
617  struct ExpFitAdaptor
618  {
619  ExpFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
620  double operator()(double* xs, double* ps)
621  {
622  fRes->coeff_exp = ps[0];
623  fRes->atten_length = ps[1];
624  fRes->background = ps[2];
625  return fRes->MeanPEPerCmAt(xs[0]);
626  }
627  calib::AttenCurve* fRes;
628  };
629  ExpFitAdaptor* adapt = new ExpFitAdaptor(res);
630 
631  TString fitStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
632  TString fit_name = fitStrArr[view];
633  double rangeup=range;
634  if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
635  rangeup=fInitEdgeMuXMax;
636  }
637 
638  TF1* fit = new TF1(fit_name, adapt, -range, +rangeup, 3);
639  // Empirically, the fits end up close to here
640  fit->SetParameter(0, 11);
641  if(geom->DetId() == novadaq::cnv::kNEARDET){
642  fit->SetParameter(0, 43);
643  fit->SetParLimits(0, 0, 1000);
644  fit->SetParLimits(1, 0, 1000);
645  fit->SetParLimits(2, 0, 1000);
646  }
647  if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
648  fit->SetParameter(0, 70);
649  }
650  fit->SetParLimits(0, 0, 1000);
651  fit->SetParLimits(1, 0, 1000);
652  fit->SetParLimits(2, 0, 1000);
653 
654  fit->SetParameter(1, 400);//kTrueShortAtten);
655  fit->SetParameter(2, 7.5);
656  fit->SetParName(0, "norm");
657  fit->SetParName(1, "atten");
658  fit->SetParName(2, "bkgd");
659 
660  fit->SetRange(-range, +rangeup);
661 
662  res->chisq = 0; // Need to count as calibrated during fit
663  TFitResultPtr fr = prof->Fit(fit, "RS"+opt);
664  //std::cout<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell()<<std::endl;
665  //std::cout<<fit->GetParameter(0)<<" "<<fit->GetParameter(1)<<" "<<fit->GetParameter(2)<<std::endl;
666 
667  if(FitQuality(prof, fit)>0.2 && view==2){
668  //std::cout<<"Default "<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell()<<std::endl;
669  //70 400 7.5
670  fit->SetParameter(0, 43);
671  fit->SetParameter(1, 452);//kTrueShortAtten);
672  fit->SetParameter(2, 7.6);
673  res->atten_length = 73;
674  res->coeff_exp = 452;
675  res->background = 7.6;
676  res->chisq = FitQuality(prof, fit);
677  //res->chisq = 0;
678 
679  return fit;
680  }
681 
682  res->atten_length = fit->GetParameter(0);
683  res->coeff_exp = fit->GetParameter(1);
684  res->background = fit->GetParameter(2);
685  res->chisq = FitQuality(prof, fit);
686 
687 
688  // const int dof = prof->FindBin(+range)-prof->FindBin(-range)-3;
689  // res->chisq =(fr == 0) ? fr->Chi2()/dof : -1;
690 
691  return fit;
692  }
693 
694  //......................................................................
695  TF1* AttenFit::rolloffFit(TH1* prof, TString opt, geo::OfflineChan chan,
696  calib::AttenCurve*& res)
697  {
698  struct RolloffFitAdaptor
699  {
700  RolloffFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
701  double operator()(double* xs, double* ps)
702  {
703  fRes->coeff_low = ps[0];
704  fRes->coeff_high = ps[1];
705  fRes->edge_low = ps[2];
706  fRes->edge_high = ps[3];
707  return fRes->MeanPEPerCmAt(xs[0]);
708  }
709  calib::AttenCurve* fRes;
710  };
711  RolloffFitAdaptor* adapt = new RolloffFitAdaptor(res);
712 
714  int view = GetView(chan.Plane(),geom);
715  if(fConsolidateViewsMode) view = (chan.Plane());
716 
717  const double edge = (view == geo::kX) ? fInitEdgeX : fInitEdgeY;
718  const double range = (view == geo::kX) ? geom->DetHalfWidth() : geom->DetHalfHeight();
719  double edgeup=edge;
720  double rangeup=range;
721  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
722  edgeup=fInitEdgeMuXMax;
723  rangeup=fInitEdgeMuXMax+50;
724  }
725 
726  TF1* fit = new TF1(view == geo::kX ? "rollFitX" : "rollFitY",
727  adapt, -range, +rangeup, 4);
728  fit->SetParName(0, "coeffL");
729  fit->SetParName(1, "coeffR");
730  fit->SetParName(2, "edgeL");
731  fit->SetParName(3, "edgeR");
732  fit->SetParameter(0, 0);
733  fit->SetParameter(1, 0);
734  fit->SetParameter(2, -edge);
735  fit->SetParameter(3, edgeup);
736  fit->SetRange(-range, +rangeup);
737 
738  if(fIncludeRolloffFit&&view!=2){
739  res->chisq = 0; // Need to count as calibrated during fit
740  TFitResultPtr fr = prof->Fit(fit, "RS"+opt);
741 
742  res->coeff_low = fit->GetParameter(0);
743  res->coeff_high = fit->GetParameter(1);
744  res->edge_low = fit->GetParameter(2);
745  res->edge_high = fit->GetParameter(3);
746 
747  res->chisq = FitQuality(prof, fit);
748 
749  // const int dof = prof->FindBin(+range)-prof->FindBin(-range)-7;
750  // res->chisq = (fr == 0) ? fr->Chi2()/dof : -1;
751  }
752 
753  return fit;
754  }
755 
756  //......................................................................
757  TGraph* AttenFit::lowessFit(TH1* prof, int plane, calib::AttenCurve*& res)
758  {
760  int view = GetView(plane, geom);
761  TGraph* ret = new TGraph;
762  double x0 = (view == geo::kX) ? -geom->DetHalfHeight() : -geom->DetHalfWidth();
763  double x1 = (view == geo::kX) ? +geom->DetHalfHeight() : +geom->DetHalfWidth();
764  const int kNumControlPts = 20;
765  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
766  x0 = -geom->DetHalfHeight();
767  x1 = fInitEdgeMuXMax+50;
768  }
769  // Range of measurements that influence each control point
770  const double kRange = 1.5*(x1-x0)/kNumControlPts;
771 
772  for(int n = 0; n < kNumControlPts; ++n){
773  const double x = x0+double(n)/(kNumControlPts-1)*(x1-x0);
774 
775  std::vector<double> xs, ys, ws;
776 
777  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
778  const double xi = prof->GetBinCenter(i);
779  if(xi < x0 || xi > x1) continue;
780  if(fabs(xi-x) > kRange) continue;
781 
782  xs.push_back(xi);
783  ys.push_back(prof->GetBinContent(i));
784  // Apparently "tricube" weights are traditional
785  ws.push_back(util::cube(1-util::cube(fabs(xi-x)/kRange)));
786  }
787 
788  double m, c;
789 
790  util::LinFit(xs, ys, ws, m, c);
791 
792  const double y = m*x+c;
793 
794  ret->SetPoint(n, x, y);
795 
796  res->AddInterpPoint(x, y);
797  }
798  ret->SetLineWidth(2);
799  ret->SetMarkerColor(kRed);
800  ret->SetLineColor(kRed);
801  ret->SetMarkerStyle(kFullCircle);
802  return ret;
803  }
804 
805  //......................................................................
807  {
808  for(int sign = -1; sign <= +1; sign += 2){
809  TGraph* g = new TGraph;
810  g->SetPoint(0, sign*range, 0);
811  g->SetPoint(1, sign*range, 1e4);
812  g->SetLineStyle(7);
813  g->Draw("l same");
814  }
815  }
816 
817  //......................................................................
818  void AttenFit::DrawDetectorEdges(double rangedown, double rangeup)
819  {
820  TGraph* g = new TGraph;
821  g->SetPoint(0, rangedown, 0);
822  g->SetPoint(1, rangedown, 1e4);
823  g->SetLineStyle(7);
824  g->Draw("l same");
825  TGraph* g2 = new TGraph;
826  g2->SetPoint(0, rangeup, 0);
827  g2->SetPoint(1, rangeup, 1e4);
828  g2->SetLineStyle(7);
829  g2->Draw("l same");
830 
831 
832 
833  }
834 
835  //......................................................................
836  void AttenFit::fit_channel_prof(TH1* atten, int plane, int cell, calib::AttenCurve*& res)
837  {
838 
840  int view = GetView(plane,geom);
841  if(fConsolidateViewsMode) view = (plane);
842 
843 
844  const TString dataMCStr = fIsMC ? "_mc" : "_data";
845  const TString viewStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
846  const TString viewStr =viewStrArr[view];
847 
848  res->initialized = true;
849 
850  double cell_length = (view == geo::kX) ? 2*geom->DetHalfHeight() : 2*geom->DetHalfWidth();
851  res->cell_length = cell_length;
852  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
853  res->cell_length =(cell_length/2)+50+fInitEdgeMuXMax;
854  }
855 
856  TString suffix = dataMCStr+viewStr;
857  if(cell >= 0)
858  suffix += TString::Format("_%03d_%03d", plane, cell);
859 
860  TString quiet;
861  if(cell >= 0) quiet = "Q";
862 
863  // TString title = "NDOS cosmic ";
864  TString title;
865  if (geom->DetId() == novadaq::cnv::kFARDET) title = "FD cosmic ";
866  if (geom->DetId() == novadaq::cnv::kNEARDET) title = "ND cosmic ";
867  if (geom->DetId() == novadaq::cnv::kNDOS) title = "NDOS cosmic ";
868  if (fIsMC) title += "Monte Carlo"; else title += "data";
869 
870  TString viewStrArrVertHor[4]={"vertical","horizontal","vertical","horizontal"};
871  TString viewStrVertHor =viewStrArrVertHor[view];
872 
873  title += TString::Format(" - plane %d (%s), cell %d",
874  plane,
875  viewStrVertHor.Data(),
876  cell);
877 
878  new TCanvas;
879 
880  TF1* fit = expFit(atten, quiet, geo::OfflineChan(plane, cell), res);
881 
882  atten->Draw();
883  atten->SetTitle(title+";Distance from center (cm);Mean PE / cm");
884 
885  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
887  }
888  else{
889  DrawDetectorEdges(cell_length/2);
890  }
891 
892  const double attenErr = fit->GetParError(0); /* unused */
893 
894  atten->GetYaxis()->SetRangeUser(0, 120);
895 
896  TLegend* leg = new TLegend(.3, .15, .7, .35);
897  leg->SetFillStyle(0);
898 
899  leg->AddEntry((TObject*)0, TString::Format("Attenuation length = %d#pm%.1lf cm", int(res->atten_length), attenErr), "");
900  leg->AddEntry((TObject*)0, TString::Format("Background = %.1lf", res->background), "");
901  leg->Draw();
902 
903  std::string plotsDirectory = fPlotsDirectory;
904  if(!fPlotsDirectory.empty()){
905  plotsDirectory = TString::Format("%s/%03d", fPlotsDirectory.c_str(), plane);
906  gSystem->mkdir(plotsDirectory.c_str());
907  gPad->Print(plotsDirectory+"/atten"+suffix+".eps");
908  }
909 
910  new TCanvas;
911  TH1D* profCorr = (TH1D*)atten->Clone();
912  profCorr->GetYaxis()->SetTitle("Ratio to exponential fit");
913  fit->SetRange(-1.1*cell_length, +1.1*cell_length);
914  profCorr->Divide(fit);
915  profCorr->Draw();
916 
917  TF1* rollFit = rolloffFit(atten, quiet, geo::OfflineChan(plane, cell), res);
918 
919  if(res->coeff_exp == 0)
920  mf::LogWarning("AttenFit") << "Fit failed: plane = " << plane << ", cell = " << cell;
921 
922  profCorr->GetYaxis()->SetRangeUser(.4, 1.1);
923 
924  leg = new TLegend(.3, .15, .7, .35);
925  leg->SetFillStyle(0);
926 
927  leg->AddEntry((TObject*)0, TString::Format("Left edge = %.1lf cm", res->edge_low), "");
928  leg->AddEntry((TObject*)0, TString::Format("Right edge = %.1lf cm", res->edge_high), "");
929  leg->Draw();
930 
931  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
933  }
934  else{
935  DrawDetectorEdges(cell_length/2);
936  }
937 
938  if(!plotsDirectory.empty()){
939  gPad->Print(plotsDirectory+"/rolloff"+suffix+".eps");
940  }
941 
942 
943  new TCanvas;
944  TH1* profFullCorr = (TH1*)atten->Clone();
945  profFullCorr->GetYaxis()->SetTitle("Ratio to full fit");
946  rollFit->SetRange(-800, +800);
947  if(geom->DetId() == novadaq::cnv::kNEARDET){
948  rollFit->SetRange(-250, +250);
949  }
950  profFullCorr->Divide(rollFit);
951  profFullCorr->GetYaxis()->SetRangeUser(0.9, 1.1);
952  profFullCorr->Draw();
953 
954  TGraph* one = new TGraph;
955  one->SetPoint(0, -1e4, 1);
956  one->SetPoint(1, +1e4, 1);
957  one->SetLineStyle(7);
958  one->Draw("l same");
959 
960  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
962  }
963  else{
964  DrawDetectorEdges(cell_length/2);
965  }
966 
967  if(fIncludeLowessFit){
968  TGraph* lowess = lowessFit(profFullCorr, plane, res);
969  lowess->Draw("lp same");
970  }
971 
972  res->initialized = true;
973 
974  if(!plotsDirectory.empty()){
975  gPad->Print(plotsDirectory+"/residual"+suffix+".eps");
976  }
977 
978 
979  new TCanvas;
980  TH1* plotProf = (TH1*)atten->Clone();
981  plotProf->Draw();
982 
983  struct FullFitAdaptor
984  {
985  FullFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
986  double operator()(double* xs, double* ps)
987  {
988  return fRes->MeanPEPerCmAt(xs[0]);
989  }
990  calib::AttenCurve* fRes;
991  } adapt(res);
992 
993  const double range = cell_length/2;
994  double rangeup=range;
995  if(view==2&& geom->DetId() == novadaq::cnv::kNEARDET){
996  rangeup=fInitEdgeMuXMax+50;
997  }
998  TF1* fullFit = new TF1("fullFit", adapt, -range, +rangeup, 0);
999  fullFit->SetLineColor(kBlue);
1000  fullFit->SetRange(-range, +rangeup);
1001  fullFit->Draw("same");
1002 
1003  plotProf->SetLineColor(kBlack);
1004  plotProf->GetYaxis()->SetTitle("Mean PE / cm");
1005  plotProf->GetYaxis()->SetRangeUser(0, 150);
1006 
1007  // HACKs for making a pretty plot for blessing
1008  plotProf->GetXaxis()->CenterTitle();
1009  plotProf->GetYaxis()->CenterTitle();
1010  plotProf->GetYaxis()->SetTitleSize(plotProf->GetXaxis()->GetTitleSize());
1011  plotProf->GetYaxis()->SetTitleOffset(plotProf->GetXaxis()->GetTitleOffset());
1012 
1013  if(geom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1015  }
1016  else{
1017  DrawDetectorEdges(cell_length/2);
1018  }
1019 
1020  if(!plotsDirectory.empty()){
1021  gPad->Print(plotsDirectory+"/totfit"+suffix+".eps");
1022  gPad->Print(plotsDirectory+"/totfit"+suffix+".png");
1023  }
1024 
1025  new TCanvas;
1026  TH1* resid = new TH1F("", ";W (cm);Ratio to total fit", atten->GetNbinsX(), atten->GetXaxis()->GetXmin(), atten->GetXaxis()->GetXmax());
1027 
1028  for(int n = 0; n < atten->GetNbinsX()+2; ++n){
1029  resid->SetBinContent(n, atten->GetBinContent(n)/res->MeanPEPerCmAt(atten->GetBinCenter(n)));
1030  }
1031  resid->GetYaxis()->SetRangeUser(.9, 1.1);
1032  resid->Draw();
1033  one->Draw("l same");
1034 
1035  if(!plotsDirectory.empty()){
1036  gPad->Print(plotsDirectory+"/residual_final"+suffix+".eps");
1037  }
1038 
1039  delete fit;
1040  delete rollFit;
1041  delete profCorr;
1042  delete profFullCorr;
1043  delete one;
1044  delete leg;
1045  delete plotProf;
1046  delete fullFit;
1047  delete resid;
1048  }
1049 
1050  //......................................................................
1051  void AttenFit::WriteDummyCSVRow(FILE* f, int plane, int cell) const
1052  {
1053  const int vldTime = 1; // 1970
1054 
1055  geo::OfflineChan chan(plane, cell);
1056 
1057  if(fConsolidateViewsMode && !fIsMC){
1058  assert(plane == 0 || plane == 1 || plane == 2 || plane == 3);
1059  chan = geo::OfflineChan(1000+plane, cell);
1060  }
1061 
1062  fprintf(f, "-1,-1,-1,-1,-1,-1,-1,1000,%d,%d\n",
1063  chan.ToDBValidityChan(), vldTime);
1064  }
1065 
1067 
1068 } // end namespace
1069 ////////////////////////////////////////////////////////////////////////
caldp::AttenHists GetChannelHists(bool isData, int plane, int cell)
double MinW() const
void fit_channel_prof(TH1 *atten, int plane, int cell, calib::AttenCurve *&res)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
int TotalEntries() const
bool HasProfiles(const geo::OfflineChan &c) const
std::map< std::string, double > xmax
int MinBlock() const
TH1F * ToTH1F() const
TH2 * rh
Definition: drawXsec.C:5
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void DrawDetectorEdges(double range)
virtual void beginRun(const art::Run &run)
const Var weight
virtual void endJob()
double ThresholdCorrAt(int view, int cell, double w) const
LiteProfile WPE_corr_xy
Definition: AttenProfiles.h:85
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void reconfigure(const fhicl::ParameterSet &pset)
Float_t x1[n_points_granero]
Definition: compare.C:5
TF1 * expFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
static AttenCurve * Uninitialized(int det, geo::OfflineChan chan)
Return a new AttenCurve objects with fields uninitialized.
Definition: AttenCurve.cxx:25
TH2 * GetBestHistogram(caldp::AttenHists v, caldp::AttenProfiles u)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
TGraph * lowessFit(TH1 *prof, int plane, calib::AttenCurve *&res)
AttenProfiles for many channels.
Definition: AttenProfiles.h:89
virtual void endSubRun(const art::SubRun &sr)
TH1 * MedianProfile(TH2 *f)
double corr(TGraph *g, double thresh)
std::vector< geo::OfflineChan > Channels() const
List of channels with profiles on.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
TH1 * TruncatedMeanProfile(TH2 *h)
const PlaneGeo * Plane(unsigned int i) const
std::string fThresholdCorrMapFile
bool isRealData() const
Definition: Event.h:83
TH1 * MeanProfile(TH2 *f)
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Definition: Run.h:31
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
AttenFit(const fhicl::ParameterSet &pset)
fclose(fg1)
int ToDBValidityChan() const
Definition: OfflineChan.cxx:19
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
CDPStorage service.
block
print "ROW IS " print row
Definition: elec2geo.py:31
std::map< int, caldp::AttenProfilesMap > fChannelMapProf
void WriteDummyCSVRow(FILE *f, int plane, int cell) const
Ensure every cell has a CSV entry even if it couldn&#39;t be fit.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
AttenHists & GetHists(const geo::OfflineChan &c)
Definition: AttenHists.h:81
Near Detector in the NuMI cavern.
virtual void analyze(const art::Event &evt)
void ApplyThresholdCorrection(int view, int cell, TH1 *prof)
Histograms used by attenuation calibration.
Definition: AttenHists.h:27
caf::StandardRecord * sr
std::string fAttenHistsMapLabel
int MaxBlock() const
double frac(double x)
Fractional part.
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
z
Definition: test.py:28
AttenHists for many channels.
Definition: AttenHists.h:77
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
unsigned short Plane() const
Definition: OfflineChan.h:31
Profiles used by attenuation calibration.
Definition: AttenProfiles.h:61
geo::OfflineChan chan
Definition: AttenCurve.h:53
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
unsigned short Cell() const
Definition: OfflineChan.h:32
caldp::AttenHistsMap fChannelMap
This is where GetChannelHist stores its histograms. Map is from block.
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
LiteProfile WPE_corr_z
Definition: AttenProfiles.h:85
A (plane, cell) pair.
Definition: OfflineChan.h:17
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
LiteProfile WPE_corr
Definition: AttenProfiles.h:85
void geom(int which=0)
Definition: geom.C:163
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
virtual void endRun(const art::Run &run)
assert(nhit_max >=nhit_nbins)
int GetView(int plane, art::ServiceHandle< geo::Geometry > geom)
float MeanPEPerCmAt(double w) const
Mean response of this channel at this distance from detector centre.
Definition: AttenCurve.cxx:37
caldp::AttenProfiles GetChannelProfiles(bool isData, int plane, int cell)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
bool removeCachedProduct(Handle< PROD > &h) const
Definition: DataViewImpl.h:551
auto one()
Definition: PMNS.cxx:49
unsigned int NPlanes() const
nova::dbi::DataType GetDataType()
Definition: RunHistory.h:452
Double_t sum
Definition: plot.C:31
std::string fPlotsDirectory
Float_t w
Definition: plot.C:20
TF1 * rolloffFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
AttenProfiles & GetProfiles(const geo::OfflineChan &c)
double FitQuality(TH1 *prof, TF1 *fit) const
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
void WriteToCSVs(FILE *fConsts, FILE *fPoints, bool mc) const
Definition: AttenCurve.cxx:129
double MaxW() const
FILE * outfile
Definition: dump_event.C:13
def sign(x)
Definition: canMan.py:197
ThresholdCorrMap * fThreshCorrMap
TH2F * WPE_corr_xy
Definition: AttenHists.h:67
void AddInterpPoint(float w, float factor)
Definition: AttenCurve.cxx:122
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.