AttenuationFit_plugin.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file AttenuationFit_plugin.cc
3 // \brief Fit attenuation behaviour from AttenHists files
4 // \author Brian Rebel - brebel@fnal.gov
5 // Based on AttenFit_module.cc by Christopher Backhouse, Dan Pershey
6 // - bckhouse@caltech.edu, pershey@caltech.edu
7 ////////////////////////////////////////////////////////////////////////
8 
9 // C/C++ includes
10 #include <cmath>
11 #include <iostream>
12 #include <sstream>
13 
14 // ROOT includes
15 #include "TFile.h"
16 #include "TH1.h"
17 #include "TH2.h"
18 #include "TF1.h"
19 #include "TString.h"
20 #include "TGraph.h"
21 #include "TCanvas.h"
22 #include "TLegend.h"
23 #include "TProfile.h"
24 #include "TFitResult.h"
25 #include "TFitResultPtr.h"
26 #include "TStyle.h"
27 #include "TSystem.h"
28 
29 // Framework includes
36 
37 // NOvASoft includes
41 #include "Geometry/Geometry.h"
46 #include "RecoBase/CellHit.h"
47 #include "RecoBase/Track.h"
49 #include "Utilities/func/MathUtil.h"
50 #include "NovaDAQConventions/DAQConventions.h"
54 
55 namespace calib
56 {
57  // Parameters for fits, should be moved to .fcl file
58  int kStatsLimit = 1;
59  //const double kTrueShortAtten = 289.5;
60 
62  {
63  public:
64  explicit AttenuationFit(fhicl::ParameterSet const& pset);
66 
67  void reconfigure(fhicl::ParameterSet const& pset);
68  void clear();
69  void writeResults(art::Results & r);
70  void readResults(art::Results const& run);
71  void beginRun(art::Run const& run);
72 
73  protected:
74  int MinBlock() const;
75  int MaxBlock() const;
76  int GetView(int plane);
77  //int GetFiberBrightness(int plane, int cell);
78  void ConsolidateViews();
79 
81  int plane,
82  int cell);
84  int plane,
85  int cell);
86 
89  TH1* MedianProfile(TH2* f);
90  TH1* MeanProfile(TH2* f);
91  TH1* TruncatedMeanProfile(TH2* h);
93  int cell,
94  TH1* prof);
95  void ApplyThresholdCorrection(int view,
96  int cell,
97  TH1* prof,
98  int fiberbrightness);
99  TF1* expFit(TH1* prof,
100  TString opt,
102  calib::AttenCurve*& res);
103  TF1* rolloffFit(TH1* prof,
104  TString opt,
105  geo::OfflineChan chan,
106  calib::AttenCurve*& res);
107  TGraph* lowessFit(TH1* prof,
108  int plane,
109  calib::AttenCurve*& res);
110  void fit_channel_prof(TH1* atten,
111  int plane,
112  int cell,
113  calib::AttenCurve*& res,
114  int fiberbrightness = 0);
115  void DrawDetectorEdges(double range);
116  void DrawDetectorEdges(double rangedown,
117  double rangeup);
118 
119  double FitQuality(TH1* prof,
120  TF1* fit) const;
121 
122  /// Ensure every cell has a CSV entry even if it couldn't be fit
123  void WriteDummyCSVRow(FILE* f,
124  int plane,
125  int cell) const;
126 
127  double fInitEdgeX;
128  double fInitEdgeY;
138  bool fHybridMode;//For MC when we want the final results separated by fb, but no fb in the correction factor
146  bool fIsMC;
147 
149 
150  /// This is where \ref GetChannelHist stores its histograms. Map is from block
152  std::map<int, caldp::AttenProfilesMap> fChannelMapProf;
153 
157  };
158 
159  //......................................................................
161  {
162  reconfigure(pset);
163 
164  fThreshCorrMap = 0;
165 
166  gStyle->SetOptStat(0);
167  }
168 
169  //......................................................................
171  {
172  delete fThreshCorrMap;
173  }
174 
175  //......................................................................
177  {
178  fAttenHistsMapLabel = pset.get<std::string>("AttenHistsMapLabel" );
179  fAttenHistsMapSuffix = pset.get<std::string>("AttenHistsMapSuffix", "block" );
180  fThresholdCorrMapFile = pset.get<std::string>("ThresholdCorrMapFile" );
181  fPlotsDirectory = pset.get<std::string>("PlotsDirectory" );
182  fProfileType = pset.get<std::string>("ProfileType", "XY" );
183  fIsMC = pset.get<bool >("IsMC", true );
184  fFiberBrightnessMode = pset.get<bool >("FiberBrightnessMode", false );
185  fMuCFiberBrightnessMode = pset.get<bool >("MuCFiberBrightnessMode", false);
186  fHybridMode = pset.get<bool >("HybridMode", false );
187  fProfileMode = pset.get<bool >("ProfileMode" );
188  fConsolidateViewsMode = pset.get<bool >("ConsolidateViewsMode" );
189  fMuCConsolidateViewsMode = pset.get<bool >("MuCConsolidateViewsMode", true);
190  fIncludeRolloffFit = pset.get<bool >("IncludeRolloffFit" );
191  fIncludeLowessFit = pset.get<bool >("IncludeLowessFit" );
192  fPlaneMin = pset.get<int >("PlaneMin" );
193  fPlaneMax = pset.get<int >("PlaneMax" );
194 
195  if(fHybridMode) fFiberBrightnessMode = false;
196 
197  LOG_VERBATIM("AttenuationFit") << "PlaneMin: "<< fPlaneMin;
198  LOG_VERBATIM("AttenuationFit") << "PlaneMax: "<< fPlaneMax;
199 
200  LOG_VERBATIM("AttenuationFit") << "07/2/2019 updates";
201  if (fIsMC) LOG_VERBATIM("AttenuationFit") << "File is MC";
202  else LOG_VERBATIM("AttenuationFit") << "File is data";
203 
204  if(fFiberBrightnessMode) LOG_VERBATIM("AttenuationFit") << "FiberBrightnessMode is ON";
205  else if (fHybridMode) LOG_VERBATIM("AttenuationFit") << "HybridMode is ON";
206  else LOG_VERBATIM("AttenuationFit") << "FiberBrightnessMode and HybridMode are OFF";
207 
208  LOG_VERBATIM("AttenuationFit") << "Profile Type = " << fProfileType;
209 
210  LOG_VERBATIM("AttenuationFit") << "ThresholdCorrMapFile = " << fThresholdCorrMapFile;
211  }
212 
213  //......................................................................
215  {
216  std::cout<<"beginRun"<<std::endl;
217  // fInitEdgeX = fGeom->DetHalfWidth();
218  // fInitEdgeY = fGeom->DetHalfHeight();
222  fInitEdgeX-=50;
223  fInitEdgeY-=50;
224  fInitEdgeMuXMin=-150;
225  fInitEdgeMuXMax=50;
226  }
227 
228  if(!fThreshCorrMap)
230  if(!fThreshCorrMap)
231  LOG_VERBATIM("AttenuationFit") << "Can't find fThreshCorrMap in " << fThresholdCorrMapFile << "!";
232  }
233 
234  //......................................................................
236  {
237  return;
238  }
239 
240  //......................................................................
242  {
243  std::cout<<"writeResults"<<std::endl;
244 
245  LOG_VERBATIM("AttenuationFit") << "Writing results";
246 
248 
249  TFile* outfile = new TFile("Atten_Fit_Results.root","RECREATE");
250 
251  //Look at a single cell from a couple of FBs and views
252  TCanvas * c_XView_fb0 = new TCanvas;
253  TCanvas * c_XView_fb8 = new TCanvas;
254  TCanvas * c_YView_fb0 = new TCanvas;
255  TCanvas * c_YView_fb8 = new TCanvas;
256  TCanvas * c_XView_ratio_fb0 = new TCanvas;
257  TCanvas * c_XView_ratio_fb8 = new TCanvas;
258  TCanvas * c_YView_ratio_fb0 = new TCanvas;
259  TCanvas * c_YView_ratio_fb8 = new TCanvas;
260  TH1F* profclone = 0;
261  TLegend * l_PECorrections = new TLegend(0.7, 0.1, 0.9, 0.3);
262  l_PECorrections->SetFillStyle(0);
263  l_PECorrections->SetBorderSize(0);
264 
265  unsigned int plane_max = fPlaneMax;
266  bool isND = fGeom->DetId() == novadaq::cnv::kNEARDET;
267  bool isTB = fGeom->DetId() == novadaq::cnv::kTESTBEAM;
268 
269  unsigned int fiberbrightness_max = 0;
270  if(fIsMC && (fFiberBrightnessMode || fHybridMode)) fiberbrightness_max = 11;
271  LOG_VERBATIM("AttenuationFit") << "Number of fiber brightness bins : " << fiberbrightness_max + 1;
272 
273  // Need more info on stats and quality of fit of MC
274  // For each view, make a 2d fb vs cell plot
275  // with z = number of events (aka TotalEntries)
276  // Same for FitQuality, aka chi sq
277  // Can we make this into a map of the detector??
278  TH2F * h_nEvents_x = new TH2F("h_nEvents_x", "Number of events in XView cells;Cell;FB", 384, 0, 384, fiberbrightness_max+1, 0, fiberbrightness_max+1);
279  TH2F * h_nEvents_y = new TH2F("h_nEvents_y", "Number of events in YView cells;Cell;FB", 384, 0, 384, fiberbrightness_max+1, 0, fiberbrightness_max+1);
280  TH2F * h_chiSq_x = new TH2F("h_chiSq_x", "#chi^{2} of fits in XView cells;Cell;FB", 384, 0, 384, fiberbrightness_max+1, 0, fiberbrightness_max+1);
281  TH2F * h_chiSq_y = new TH2F("h_chiSq_y", "#chi^{2} of fits in YView cells;Cell;FB", 384, 0, 384, fiberbrightness_max+1, 0, fiberbrightness_max+1);
282 
283  if(fPlaneMax < 0) {
284  plane_max = fGeom->NPlanes()-1;
285  if(fConsolidateViewsMode) plane_max = 1;
286  if(fConsolidateViewsMode && isND) plane_max = 3;
288  }
289 
290  const TString mcStr = fIsMC ? "_mc" : "";
291  FILE* fConstsCSV[fiberbrightness_max+1];
292  FILE* fPointsCSV[fiberbrightness_max+1];
293 
294  for(unsigned int i=0; i<=fiberbrightness_max; i++)
295  {
296  TString fbStr = "";
297  if(fIsMC && (fFiberBrightnessMode || fHybridMode)) fbStr.Form("_fb%i", i);
298  //LOG_VERBATIM("AttenuationFit") << "Making CSV files for " << fbStr;
299  fConstsCSV[i] = fopen("calib_atten_consts"+mcStr+fbStr+".csv","w");
300  fPointsCSV[i] = fopen("calib_atten_points"+mcStr+fbStr+".csv","w");
301 
302  fprintf(fConstsCSV[i], "#coeff_exp,atten_length,background,edge_low,edge_high,coeff_low,coeff_high,chisq,channel,tv\n");
303  fprintf(fPointsCSV[i], "#w,factor,channel,tv\n");
304  }
305 
306  //LOG_VERBATIM("AttenuationFit") << "Fitting plane " << fPlaneMin << " to plane " << plane_max;
307 
308  for (unsigned int fiberbrightness = 0; fiberbrightness <= fiberbrightness_max; ++fiberbrightness) {
309  TString fbStr;
310  fbStr.Form("FB%i", fiberbrightness);
311 
312  for (unsigned int plane = fPlaneMin; plane <= plane_max; ++plane){
313 
315  LOG_VERBATIM("AttenuationFit")
316  << "Fitting fb "
317  << fiberbrightness << "/"
318  << fiberbrightness_max;
319  }
320 
321  const int block = plane/32;
322  int view = GetView(plane);
323 
325  LOG_VERBATIM("AttenuationFit") << "Consolidated by views";
326  view = plane;
327  }
328 
330  view = plane - fGeom->FirstPlaneInMuonCatcher() + 2;
331  }
332 
333  LOG_VERBATIM("AttenuationFit") << "View = " << view;
334 
335  int vartofind;
336  if(fIsMC && (fFiberBrightnessMode || fHybridMode)) vartofind = fiberbrightness;
337  else vartofind = block;
338  if((!fMuCFiberBrightnessMode || !fIsMC) && view > geo::kY) vartofind = 100;
339 
340  auto it = fChannelMapProf.find(vartofind);
341 
342  if(fIsMC && (fFiberBrightnessMode || fHybridMode) && it == fChannelMapProf.end()){
343  LOG_VERBATIM("AttenuationFit")
344  << "fiberbrightness " << vartofind
345  << " it == fChannelMapProf.end()";
346  continue;
347  } else if(it == fChannelMapProf.end()){
348  LOG_VERBATIM("AttenuationFit")
349  << "block " << vartofind
350  << " it == fChannelMapProf.end()";
351  continue;
352  }
353 
354  caldp::AttenProfilesMap& profmap = it->second;
355 
357  LOG_DEBUG("AttenuationFit")
358  << "profmap for fiberbrightness "
359  << vartofind << " has "
360  << profmap.Channels().size() << " channels";
361 
362  if(profmap.Channels().size() < 1){
363  LOG_WARNING("AttenuationFit")
364  << "no channels for " << vartofind
365  << " don't try to fit any";
366  continue;
367  }
368  } else {
369  LOG_DEBUG("AttenuationFit")
370  << "profmap for block "
371  << vartofind << " has "
372  << profmap.Channels().size() << " channels";
373 
374  if(profmap.Channels().size() < 1){
375  LOG_WARNING("AttenuationFit")
376  << "no channels for " << vartofind
377  << " don't try to fit any";
378  continue;
379  }
380  }//end if(fibrebrightness){}else{}
381 
382  TDirectory* dir = 0;
383  unsigned int cell_max = fGeom->Plane(plane)->Ncells();
384  // offset by 1 if only consolidating muon catcher
386  cell_max = fGeom->Plane(plane+1)->Ncells();
387 
388 
389  LOG_DEBUG("AttenuationFit")
390  << "CellMax: "
391  << cell_max;
392  std::cout << "CellMax: " << cell_max << std::endl;
393 
394  for(unsigned int cell = 0; cell < cell_max; ++cell){
395  LOG_DEBUG("AttenuationFit")
396  << " Fitting cell "
397  << cell << "/" << cell_max;
398 
399  //Make plots of raw pe, threshold/shadowing corrected pe, pecorr, for 4 example cells (x/y, fb0/fb8)
400  bool fb0_xview_example = false;
401  bool fb8_xview_example = false;
402  bool fb0_yview_example = false;
403  bool fb8_yview_example = false;
404  bool fb0_nd_xview_example = false;
405  bool fb8_nd_xview_example = false;
406  bool fb0_nd_yview_example = false;
407  bool fb8_nd_yview_example = false;
408 
409  int data_fiberbrightness = 0;
410 
411  std::vector<float> PESum;
412  std::vector<int> PECount;
413  std::vector<float> MeVSum;
414  std::vector<int> MeVCount;
415 
416  if(!fIsMC && (fFiberBrightnessMode || fHybridMode)) {
417  data_fiberbrightness = fFBhandle->getBrightnessIndex(plane,cell);
418  }
419 
420  TH1* prof = 0;
421  if(fProfileMode){
422  LOG_VERBATIM("AttenuationFit") << "ProfileMode is ON";
423  int pln = plane;
426  if(!profmap.HasProfiles(pln, cell)){
427  LOG_DEBUG("AttenuationFit")
428  << "Nothing to fit for plane "
429  << plane << ", cell " << cell;
430  WriteDummyCSVRow(fConstsCSV[fiberbrightness], plane, cell);
431  continue;
432  }
433  caldp::AttenProfiles v = profmap.GetProfiles(pln, cell);
434 
435  if(fIsMC) {
436  PESum = v.WPE_corr_xy.Sum();
437  PECount = v.WPE_corr_xy.Count();
438  MeVSum = v.WPE_corr_xy_truth.Sum();
439  MeVCount = v.WPE_corr_xy_truth.Count();
440  }
441 
442  if(fProfileType.compare("Trajectory")==0) prof = v.WPE_corr_traj.ToTH1F();
443  else if(fProfileType.compare("XY")==0){
444  prof = v.WPE_corr_xy.ToTH1F();
445 
446  // LOG_VERBATIM("AttenuationFit") << "Got prof = v.WPE_corr_xy.ToTH1F()";
447 
448  // TODO: figure out where these bad axis limits come from
449 
450  //Could be an issue, seems like FD code...
451  if(fIsMC && !isND && !isTB){
452  int binx = h_nEvents_x->GetBin(h_nEvents_x->GetXaxis()->FindBin(cell), h_nEvents_x->GetYaxis()->FindBin(fiberbrightness));
453  int biny = h_nEvents_y->GetBin(h_nEvents_y->GetXaxis()->FindBin(cell), h_nEvents_y->GetYaxis()->FindBin(fiberbrightness));
454  if(view == 0)h_nEvents_x->SetBinContent(binx, v.WPE_corr_xy.TotalEntries());
455  else h_nEvents_y->SetBinContent(biny, v.WPE_corr_xy.TotalEntries());
456  }
457 
458  // std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl; //Remove
459  // std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl; //Remove
460 
461  if(prof->GetXaxis()->GetXmin() == 0 ||
462  prof->GetXaxis()->GetXmax() == 0 ||
464  LOG_VERBATIM("AttenuationFit") << "Stats Fail "<<v.WPE_corr_xy.TotalEntries()<<" < "<<kStatsLimit<<std::endl;
465 
466  // std::cout<<"Integral "<<prof->Integral()<<std::endl;
467  // std::cout<<"GetEntries "<<prof->GetEntries()<<std::endl;
468  // std::cout<<"GetEffectiveEntries "<<prof->GetEffectiveEntries()<<std::endl;
469 
470  // std::cout<<"TotalEntries xy "<<v.WPE_corr_xy.TotalEntries()<<std::endl;
471 
472  // std::cout<<"TotalEntries z "<<v.WPE_corr_z.TotalEntries()<<std::endl;
473 
474  // std::cout<<"TotalEntries "<<v.WPE_corr.TotalEntries()<<std::endl;
475 
476  prof = v.WPE_corr_z.ToTH1F();
477 
478  if(prof->GetXaxis()->GetXmin() == 0 ||
479  prof->GetXaxis()->GetXmax() == 0 ||
481  LOG_VERBATIM("AttenuationFit")<<"Stats Double Fail "<<v.WPE_corr_z.TotalEntries()<<" "<<kStatsLimit<<std::endl;
482  WriteDummyCSVRow(fConstsCSV[fiberbrightness], plane, cell);
483  //continue;
484 
485  }
486  }
487  LOG_VERBATIM("AttenuationFit") << "Finished if XY";
488  }//end of trajectory based condition
489  }//end of ProfileMode condition
490  else{
491  int pln = plane;
495  caldp::AttenProfiles u = profmap.GetProfiles(pln, cell);
496  TH2* h_WPE = GetBestHistogram(v,u);
497  if(!h_WPE){
498  WriteDummyCSVRow(fConstsCSV[fiberbrightness], plane, cell);
499  continue;
500  }
501  prof = MedianProfile(h_WPE);
502  }
503 
504  // Never evaluated for NDOS
505  if(fGeom->DetId() == novadaq::cnv::kFARDET ||
508 
509  //DON'T COMMENT OUT THESE THREE LINES
510  //If you do, you'll get a segfault. I have no idea why...
511  if(prof){
512  LOG_VERBATIM("AttenuationFit") << "Drawing prof";
513  new TCanvas;
514  prof->Draw();
515  //gPad->Print("prof.png");
516  } else {
517  LOG_VERBATIM("AttenuationFit") << "Failed if (prof)";
518  }
519 
520  //Make plots of raw pe, threshold/shadowing corrected pe, pecorr, for 4 example cells (x/y, fb0/fb8)
521  /* if(!fIsMC){
522  if (plane == 451 && cell == 226) fb0_xview_example = true;
523  if (plane == 451 && cell == 190) fb8_xview_example = true;
524  if (plane == 452 && cell == 280) fb0_yview_example = true;
525  if (plane == 452 && cell == 309) fb8_yview_example = true;
526 
527  if (isND && plane == 87 && cell == 50) fb0_xview_example = true;
528  if (isND && plane == 95 && cell == 50) fb8_xview_example = true;
529  if (isND && plane == 90 && cell == 50) fb0_yview_example = true;
530  if (isND && plane == 108 && cell == 50) fb8_yview_example = true;
531  } else {
532  if (!isND && plane == 0 && fiberbrightness == 0 && cell == 226) fb0_xview_example = true;
533  if (!isND && plane == 0 && fiberbrightness == 8 && cell == 190) fb8_xview_example = true;
534  if (!isND && plane == 1 && fiberbrightness == 0 && cell == 280) fb0_yview_example = true;
535  if (!isND && plane == 1 && fiberbrightness == 8 && cell == 309) fb8_yview_example = true;
536 
537  if (isND && plane == 0 && fiberbrightness == 0 && cell == 50) fb0_nd_xview_example = true;
538  if (isND && plane == 0 && fiberbrightness == 8 && cell == 50) fb8_nd_xview_example = true;
539  if (isND && plane == 1 && fiberbrightness == 0 && cell == 50) fb0_nd_yview_example = true;
540  if (isND && plane == 1 && fiberbrightness == 8 && cell == 50) fb8_nd_yview_example = true;
541  }*/
542  if(fb0_xview_example || fb8_xview_example || fb0_yview_example || fb8_yview_example ||
543  fb0_nd_xview_example || fb8_nd_xview_example || fb0_nd_yview_example || fb8_nd_yview_example) {
544  profclone = (TH1F*)prof->Clone();
545  profclone->SetLineColor(kRed);
546  profclone->SetMaximum(100);
547  if(fb0_xview_example){
548  c_XView_fb0->cd();
549  profclone->SetTitle("XView, Cell 226, FB0;W (cm);PE/cm");
550  } else if(fb8_xview_example){
551  c_XView_fb8->cd();
552  profclone->SetTitle("XView, Cell 190, FB8;W (cm);PE/cm");
553  } else if(fb0_yview_example){
554  c_YView_fb0->cd();
555  profclone->SetTitle("YView, Cell 280, FB0;W (cm);PE/cm");
556  } else if(fb8_yview_example){
557  c_YView_fb8->cd();
558  profclone->SetTitle("YView, Cell 309, FB8;W (cm);PE/cm");
559  } else if(fb0_nd_xview_example){
560  c_XView_fb0->cd();
561  profclone->SetTitle("XView, Cell 50, FB0;W (cm);PE/cm");
562  } else if(fb8_nd_xview_example){
563  c_XView_fb8->cd();
564  profclone->SetTitle("XView, Cell 50, FB8;W (cm);PE/cm");
565  } else if(fb0_nd_yview_example){
566  c_YView_fb0->cd();
567  profclone->SetTitle("YView, Cell 50, FB0;W (cm);PE/cm");
568  } else if(fb8_nd_yview_example){
569  c_YView_fb8->cd();
570  profclone->SetTitle("YView, Cell 50, FB8;W (cm);PE/cm");
571  }
572  profclone->Draw();
573  l_PECorrections->AddEntry(profclone, "Raw PE", "l");
574  }
575 
576  if((view <= geo::kY && fFiberBrightnessMode) || (view > geo::kY && fMuCFiberBrightnessMode)){
577  //LOG_VERBATIM("AttenuationFit") << "About to apply threshold corrections (with fibre brightness)";
578  if(fIsMC) ApplyThresholdCorrection(view, cell, prof, fiberbrightness);
579  else ApplyThresholdCorrection(view, cell, prof, data_fiberbrightness);
580  } else {
581  //LOG_VERBATIM("AttenuationFit") << "About to apply threshold corrections";
582  ApplyThresholdCorrection(view, cell, prof);
583  }
584 
585  if(fb0_xview_example || fb8_xview_example || fb0_yview_example || fb8_yview_example ||
586  fb0_nd_xview_example || fb8_nd_xview_example || fb0_nd_yview_example || fb8_nd_yview_example) {
587  if(fb0_xview_example || fb0_nd_xview_example) c_XView_fb0->cd();
588  else if(fb8_xview_example || fb8_nd_xview_example) c_XView_fb8->cd();
589  else if(fb0_yview_example || fb0_nd_yview_example) c_YView_fb0->cd();
590  else if(fb8_yview_example || fb8_nd_yview_example) c_YView_fb8->cd();
591  prof->SetLineColor(kBlue);
592  prof->Draw("same");
593  l_PECorrections->AddEntry(prof, "Threshold corrected PE", "l");
594 
595  if(fb0_xview_example || fb0_nd_xview_example) c_XView_ratio_fb0->cd();
596  else if(fb8_xview_example || fb8_nd_xview_example) c_XView_ratio_fb8->cd();
597  else if(fb0_yview_example || fb0_nd_yview_example) c_YView_ratio_fb0->cd();
598  else if(fb8_yview_example || fb8_nd_yview_example) c_YView_ratio_fb8->cd();
599  TH1F* profratio = profclone;
600  profratio->Divide(prof);
601  profratio->SetMaximum(1.6);
602  profratio->SetMinimum(1.1);
603  profratio->Draw();
604  profratio->SetTitle("Ratio before / after threshold corrections");
605  if(fb0_xview_example || fb0_nd_xview_example){
606  profratio->SetTitle("XView FB0: Ratio before / after threshold corrections");
607  c_XView_ratio_fb0->SaveAs("plots/PEThresholdCorrRatio_xview_fb0.pdf");
608  } else if(fb8_xview_example || fb8_nd_xview_example){
609  profratio->SetTitle("XView FB8: Ratio before / after threshold corrections");
610  c_XView_ratio_fb8->SaveAs("plots/PEThresholdCorrRatio_xview_fb8.pdf");
611  } else if(fb0_yview_example || fb0_nd_yview_example){
612  profratio->SetTitle("YView FB0: Ratio before / after threshold corrections");
613  c_YView_ratio_fb0->SaveAs("plots/PEThresholdCorrRatio_yview_fb0.pdf");
614  } else if(fb8_yview_example || fb8_nd_yview_example){
615  profratio->SetTitle("YView FB8: Ratio before / after threshold corrections");
616  c_YView_ratio_fb8->SaveAs("plots/PEThresholdCorrRatio_yview_fb8.pdf");
617  }
618  }
619 
620  }//end if(TB or ND or FD)
621 
622  // LOG_VERBATIM("AttenuationFit") << "About to declare res";
623  calib::AttenCurve* res;
624  //int nviews = 2;
625  //if(fGeom->DetId() == novadaq::cnv::kNEARDET) nviews = 4;
626  //if((fFiberBrightnessMode || fHybridMode) && fIsMC) res = calib::AttenCurve::Uninitialized(fGeom->DetId(), geo::OfflineChan((fiberbrightness*nviews)+plane, cell));
627  //else
629  LOG_VERBATIM("AttenuationFit") << "Plane: " << plane << "\t Cell: " << cell;
631  fit_channel_prof(prof, plane, cell, res, fiberbrightness);
632  int bin = h_chiSq_x->GetBin(h_chiSq_x->GetXaxis()->FindBin(cell), h_chiSq_x->GetYaxis()->FindBin(fiberbrightness));
633  if(view == 0) h_chiSq_x->SetBinContent(bin, res->chisq);
634  else if(view == 1) h_chiSq_y->SetBinContent(bin, res->chisq);
635 
636  } else if (!fIsMC && (fFiberBrightnessMode || fHybridMode)) fit_channel_prof(prof, plane, cell, res, data_fiberbrightness);
637  else {
638  //LOG_VERBATIM("AttenuationFit") << "About to fit channel";
639  fit_channel_prof(prof, plane, cell, res);
640  }
641  //LOG_VERBATIM("AttenuationFit") << "Go to outfile";
642  outfile->cd();
643  //LOG_VERBATIM("AttenuationFit") << "about to save prof in outfile";
645  if(!dir) dir = outfile->mkdir(Form("fb%i_plane_%03d", fiberbrightness, plane));
646  dir->cd();
647  prof->Write(Form("Atten_Fit_fb%i_%03d_%03d",fiberbrightness,plane,cell));
648  } else {
649  if(!dir) dir = outfile->mkdir(Form("plane_%03d", plane));
650  dir->cd();
651  prof->Write(Form("Atten_Fit_%03d_%03d",plane,cell));
652  }
653  // Encoding to indicate it's the backup averaged constants
655  res->chan = geo::OfflineChan(1000+view, cell);
656  if(fMuCConsolidateViewsMode && !fIsMC && view > geo::kY)
657  res->chan = geo::OfflineChan(1000+view, cell);
658 
659  res->WriteToCSVs(fConstsCSV[fiberbrightness], fPointsCSV[fiberbrightness], fIsMC);
660 
661  if(fb0_xview_example || fb8_xview_example || fb0_yview_example || fb8_yview_example ||
662  fb0_nd_xview_example || fb8_nd_xview_example || fb0_nd_yview_example || fb8_nd_yview_example) {
664  TH1F* pecorrprof = (TH1F*)prof->Clone();
665  for(int profentry = 0; profentry < pecorrprof->GetNbinsX(); ++profentry){
666  double W = pecorrprof->GetBinCenter(profentry);
667  float attenfactor = calibhandle->GetAttenScale(res, -fInitEdgeX, fInitEdgeX, W);//pecorrprof->GetBinCenter(profentry));
668  pecorrprof->SetBinContent(profentry, pecorrprof->GetBinContent(profentry)*attenfactor);
669  }
670  pecorrprof->SetLineColor(kBlack);
671  l_PECorrections->AddEntry(pecorrprof, "PECorr", "l");
672 
673  if(fb0_xview_example || fb0_nd_xview_example) c_XView_fb0->cd();
674  if(fb8_xview_example || fb8_nd_xview_example) c_XView_fb8->cd();
675  if(fb0_yview_example || fb0_nd_yview_example) c_YView_fb0->cd();
676  if(fb8_yview_example || fb8_nd_yview_example) c_YView_fb8->cd();
677 
678  pecorrprof->Draw("same");
679  l_PECorrections->Draw();
680 
681  if(fb0_xview_example) c_XView_fb0->SaveAs("plots/PEpercm_xview_cell226.pdf");
682  if(fb8_xview_example) c_XView_fb8->SaveAs("plots/PEpercm_xview_cell190.pdf");
683  if(fb0_yview_example) c_YView_fb0->SaveAs("plots/PEpercm_yview_cell280.pdf");
684  if(fb8_yview_example) c_YView_fb8->SaveAs("plots/PEpercm_yview_cell309.pdf");
685  if(fb0_nd_xview_example) c_XView_fb0->SaveAs("plots/PEpercm_xview_cell50.pdf");
686  if(fb8_nd_xview_example) c_XView_fb8->SaveAs("plots/PEpercm_xview_cell50.pdf");
687  if(fb0_nd_yview_example) c_YView_fb0->SaveAs("plots/PEpercm_yview_cell50.pdf");
688  if(fb8_nd_yview_example) c_YView_fb8->SaveAs("plots/PEpercm_yview_cell50.pdf");
689 
690  l_PECorrections->Clear();
691  fb0_xview_example = false;
692  fb8_xview_example = false;
693  fb0_yview_example = false;
694  fb8_yview_example = false;
695 
696  fb0_nd_xview_example = false;
697  fb8_nd_xview_example = false;
698  fb0_nd_yview_example = false;
699  fb8_nd_yview_example = false;
700  }
701 
702  /*for(unsigned int profentry = 0; profentry < PESum.size(); ++profentry) {
703  int N_PE = PECount.at(profentry);
704  if(N_PE < 1) N_PE = 1;
705  double W = -fInitEdgeX + ((2.*fInitEdgeX) / profentry);
706  float attenfactor = calibhandle->GetAttenScale(res, -fInitEdgeX, fInitEdgeX, W);
707  float PECorrSum = PESum.at(profentry) * attenfactor;
708  int N_MeV = MeVCount.at(profentry);
709  if(N_MeV < 1) N_MeV = 1;
710  }*/
711 
712  delete res;
713  delete prof;
714  } // end for cell
715  } // end for plane
716 
717  } // end for fiberbrightness
718 
719  if(fIsMC && !isTB){
720  TDirectory* statsdir = outfile->mkdir("Stats");
721  statsdir->cd();
722  h_nEvents_x->Write("nEvents_xview");
723  h_nEvents_y->Write("nEvents_yview");
724  h_chiSq_x->Write("chisq_xview");
725  h_chiSq_y->Write("chisq_yview");
726 
727  //map it out on the detector
728  int nplanes = 896;
729  int ncells = 384;
730  if(isND){
731  nplanes = 214;
732  ncells = 96;
733  }
734 
735  TH2F * h_nEvents_map = new TH2F("h_nEvents_map", "Number of events in each plane/cell;Plane;Cell", nplanes, 0, nplanes, ncells, 0, ncells);
736  TH2F * h_chiSq_map = new TH2F("h_chiSq_map", "#chi^{2} in each plane/cell;Plane;Cell", nplanes, 0, nplanes, ncells, 0, ncells);
737  for(int plane = 0; plane < nplanes; ++plane){
738  for(int cell = 0; cell < ncells; ++cell){
740  if(fb > (int)fiberbrightness_max) continue;//useful when testing with fewer than 9 FBs, to speed things up
741  int bintofill = h_nEvents_map->GetBin(h_nEvents_map->GetXaxis()->FindBin(plane), h_nEvents_map->GetYaxis()->FindBin(cell));
742  int bintoread = h_nEvents_x->GetBin(h_nEvents_x->GetXaxis()->FindBin(cell), h_nEvents_x->GetYaxis()->FindBin(fb));
743  if(GetView(plane) == 0) {
744  h_nEvents_map->SetBinContent(bintofill, h_nEvents_x->GetBinContent(bintoread));
745  h_chiSq_map->SetBinContent(bintofill, h_chiSq_x->GetBinContent(bintoread));
746  if(h_chiSq_x->GetBinContent(bintoread) == 0 || h_chiSq_x->GetBinContent(bintoread) > 0.2)
747  LOG_VERBATIM("AttenuationFit") << "Plane " << plane << " Cell " << cell << " FB " << fb << " chisq " << h_chiSq_x->GetBinContent(bintoread);
748  } else if(GetView(plane) == 1){
749  h_nEvents_map->SetBinContent(bintofill, h_nEvents_y->GetBinContent(bintoread));
750  h_chiSq_map->SetBinContent(bintofill, h_chiSq_y->GetBinContent(bintoread));
751  if(h_chiSq_y->GetBinContent(bintoread) == 0 || h_chiSq_y->GetBinContent(bintoread) > 0.2)
752  LOG_VERBATIM("AttenuationFit") << "Plane " << plane << " Cell " << cell << " FB " << fb << " chisq " << h_chiSq_y->GetBinContent(bintoread);
753  }
754  }
755  }
756 
757  h_nEvents_map->Write("nEvents_map");
758  h_chiSq_map->Write("chisq_map");
759  }
760 
761  for (unsigned int fiberbrightness = 0; fiberbrightness <= fiberbrightness_max; ++fiberbrightness){
762  fclose(fConstsCSV[fiberbrightness]);
763  fclose(fPointsCSV[fiberbrightness]);
764  }
765 
766  return;
767  }
768 
769  //......................................................................
771  {
772  std::cout<<"readResults"<<std::endl;
773 
774  int min, max;
776  min = 0;
777  max = 11;
778  } else {
779  min = MinBlock();
780  max = MaxBlock();
781  }
782 
783  std::map<int, std::string> iLabels;
784  for(int i = min; i <= max; i++){
785  std::ostringstream iLabel;
786  iLabel << fAttenHistsMapSuffix << i;
787  iLabels.insert(std::make_pair(i, iLabel.str()));
788  }
790  iLabels.insert(std::make_pair(100, "nofbmuC"));
791 
792 
793  for(auto &iLabelMap: iLabels){
794  int i = iLabelMap.first;
795  std::string iLabel = iLabelMap.second;
797 
798  //Commented out because fIsMC isn't set yet
799  //std::ostringstream iLabel;
800  //if (fIsMC && (fFiberBrightnessMode || fHybridMode)) iLabel << "fiberbrightness" << i;
801  //else iLabel << "block" << i;
802  //art::InputTag tag(fAttenHistsMapLabel, iLabel.str());
803 
805 
806  r.getByLabel(tag, profsmap);
807 
808  if(profsmap.failedToGet()){
809  LOG_WARNING("AttenuationFit")
810  << "unable to find profile map in "
811  << tag.label() << " "
812  << tag.instance();
813  continue;
814  } else {
815  LOG_DEBUG("AttenuationFit")
816  << "retrieving "
817  << tag.label() << " "
818  << tag.instance() << "; "
819  << profsmap->Channels().size() << " channels";
820  }
821 
822  // Accumulate across results
823  auto it = fChannelMapProf.find(i);
824  if(it == fChannelMapProf.end()){
825  fChannelMapProf.insert(std::make_pair(i, caldp::AttenProfilesMap(profsmap->MinW(), profsmap->MaxW())));
826  it = fChannelMapProf.find(i);
827  }
828 
829  it->second += *profsmap;
830 
831  }
832  }
833 
834  //......................................................................
836  {
837  return fPlaneMin/32;
838  }
839 
840  //......................................................................
842  {
843  if(fPlaneMax < 0) return 27;
844  return fPlaneMax/32;
845  }
846 
847  //......................................................................
849 
850  //get the view for the current combination of view and plane (x,y,mux, or muy)
851  int view=0;
853  int isMuonCatcher = (int)fGeom->FirstPlaneInMuonCatcher() <= plane;
854  view = fGeom->Plane(plane)->View() + (2 * isMuonCatcher);
855  if(isMuonCatcher) LOG_VERBATIM("AttenuationFit") << "Is muon catcher, view " << view;
856  }
857  else{
858  view = fGeom->Plane(plane)->View();
859  }
860  return view;
861  }
862 
863  //......................................................................
865  {
866  std::cout<<"ConsolidateViews"<<std::endl;
867 
868  // Has consolidation already been done by an earlier step?
869  bool isConsolidated = true;
870  bool isND = fGeom->DetId() == novadaq::cnv::kNEARDET;
871 
872  if(fChannelMapProf.find(0) != fChannelMapProf.end()){
873  for(auto chan : fChannelMapProf.find(0)->second.Channels() ){
874  if((chan.Plane() > 3 && isND) || (chan.Plane() > 1 && isND==false)){
875  isConsolidated = false;
876  break;
877  }
878  }
879  }
880 
881  if(isConsolidated){
882  LOG_INFO("AttenuationFit")
883  << "Input file appears to already be consolidated by view. Skipping that step.";
884  return;
885  }
886 
887  // Need to do it for ourselves
888  LOG_INFO("AttenuationFit")
889  << "Input file is not consolidated by view. Consolidating";
890 
891  //art::ServiceHandle<photrans::FiberBrightness> fbhandle;
892 
893  if(!fFiberBrightnessMode) LOG_VERBATIM("AttenuationFit") << "Consolidating by view";
894  if(fFiberBrightnessMode || fHybridMode) LOG_VERBATIM("AttenuationFit") << "Consolidating by view+fiber brightness";
895 
896  // if(!fiberbrightnessmode) All x-planes will go into plane 0 of this, y to 1.
897  // if(fiberbrightnessmode) Same, but with 9 fibre brightnesses
898  caldp::AttenProfilesMap sum(fChannelMapProf.begin()->second.MinW(),
899  fChannelMapProf.begin()->second.MaxW());
900 
901  //If not using fiberbrightness or looking at data, then the ChannelMapProf is sorted by block
902  //If it's mc and we're using fiberbrightness, it's sorted by fiberbrightness
903  int varMin = MinBlock();
904  int varMax = MaxBlock();
906  varMin = 0;
907  varMax = 11;
908  }
909 
910  for(int var = varMin; var <= varMax; ++var){
911  auto it = fChannelMapProf.find(var);
912  if(it == fChannelMapProf.end()) continue;
913 
914  caldp::AttenProfilesMap& profmap = it->second;
915 
916  LOG_DEBUG("AttenuationFit")
917  << "consolidating var "
918  << var << " with "
919  << profmap.Channels().size() << " channels";
920 
921  for(auto chan : profmap.Channels() ){
922  const caldp::AttenProfiles& profs = profmap.GetProfiles(chan);
923 
924  const int view = GetView(chan.Plane());
925  int index = view;
926 
928  int fiberbrightness;
929  if(!fIsMC) fiberbrightness = fFBhandle->getBrightnessIndex(chan.Plane(), chan.Cell());
930  else fiberbrightness = var;
931  index = view + fiberbrightness * (2 + isND*2);
932  //LOG_VERBATIM("AttenuationFit")
933  //<< "View: " << view
934  //<< "\t FB: " << fiberbrightness
935  //<< "\t Index: " << index;
936  }
937 
938  sum.GetProfiles(index, chan.Cell()) += profs;
939  } // end for chan
940 
941  profmap.Clear();
942  } // end for var (either block or fiberbrightness)
943 
944  /*for(int block = MinBlock(); block <= MaxBlock(); ++block){
945  auto it = fChannelMapProf.find(block);
946  if(it == fChannelMapProf.end()) continue;
947 
948  caldp::AttenProfilesMap& profmap = it->second;
949 
950  LOG_VERBATIM("AttenuationFit")//was DEBUG
951  << "consolidating block "
952  << block << " with "
953  << profmap.Channels().size() << " channels";
954 
955  for(auto chan : profmap.Channels() ){
956  const caldp::AttenProfiles& profs = profmap.GetProfiles(chan);
957 
958  const int view = GetView(chan.Plane());
959  int index = view;
960 
961  if(fFiberBrightnessMode || fHybridMode){
962  const int fiberbrightness = fFBhandle->getBrightnessIndex(chan.Plane(), chan.Cell());
963  index = view*9 + fiberbrightness;
964  }
965 
966  sum.GetProfiles(index, chan.Cell()) += profs;
967  } // end for chan
968 
969  profmap.Clear();
970  } // end for block*/
971 
972  if(sum.Channels().size() > 0){
973 
974  fChannelMapProf.clear();
975  fChannelMapProf.insert(std::make_pair(0, sum));
976 
977  auto it = fChannelMapProf.begin();
978  LOG_VERBATIM("AttenuationFit")
979  << "after consolidation, there are "
980  << it->second.Channels().size()
981  << " channels in the profile map";
982  }
983  else
984  throw cet::exception("AttenuationFit")
985  << "unable to consolidate blocks, summed profile is empty";
986 
987  return;
988  }
989 
990  //......................................................................
991  double AttenuationFit::FitQuality(TH1* prof, TF1* fit) const
992  {
993  std::cout<<"FitQuality"<<std::endl;
994 
995  // chisq lets bad fits to low stats data through, and cuts out good(ish)
996  // fits to high stats data. Try this metric, mean fractional deviation of
997  // the data from the fit (in quadrature).
998 
999  double xmin = 1.e9;
1000  double xmax = -1.e9;
1001  fit->GetRange(xmin, xmax);
1002 
1003  int tot = 0;
1004  double diff = 0;
1005 
1006  //std::cout << "xmin = " << xmin << " , xmax = " << xmax << std::endl;
1007  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
1008  const double x = prof->GetXaxis()->GetBinCenter(i);
1009  if(x < xmin || x > xmax) {
1010  continue;
1011  }
1012  ++tot;
1013  diff += util::sqr((fit->Eval(x)-prof->GetBinContent(i))/fit->Eval(x));
1014  }
1015 
1016  const double ret = sqrt(diff/tot);
1017  if(std::isnan(ret) || std::isinf(ret)) {
1018  return 2;
1019  }
1020  return ret;
1021  }
1022 
1023  //......................................................................
1025  {
1026  std::cout<<"GetBestHist"<<std::endl;
1027 
1028  if(!v.WPE) return 0;
1029  if(fProfileType.compare("Trajectory")==0){
1031  }
1032  else if(fProfileType.compare("XY")==0){
1033  if(u.WPE_corr_xy.TotalEntries() > kStatsLimit) return v.WPE_corr_xy;
1034  if(u.WPE_corr_z.TotalEntries() > kStatsLimit) return v.WPE_corr_z;
1035  if(u.WPE_corr.TotalEntries() > kStatsLimit) return v.WPE_corr;
1036  }
1037  return 0;
1038  }
1039 
1040  //......................................................................
1042  {
1043  TProfile* p = h->ProfileX();
1044 
1045  TAxis* xax = h->GetXaxis();
1046  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
1047  ret->GetXaxis()->SetTitle(xax->GetTitle());
1048  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
1049  ret->SetBinContent(ix, p->GetBinContent(ix));
1050  ret->SetBinError(ix, p->GetBinError(ix));
1051  }
1052  return ret;
1053  }
1054 
1055  //......................................................................
1057  {
1058  TAxis* xax = h->GetXaxis();
1059  TAxis* yax = h->GetYaxis();
1060  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
1061  ret->GetXaxis()->SetTitle(xax->GetTitle());
1062  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
1063  double tot = 0;
1064  // Make sure to include overflow bins
1065  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
1066  tot += h->GetBinContent(ix, iy);
1067  }
1068  double seen = 0;
1069  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
1070  seen += h->GetBinContent(ix, iy);
1071  if(seen > tot/2){
1072  const double z = h->GetBinContent(ix, iy);
1073  // Interpolate
1074  const double frac = 1-(seen-tot/2)/z;
1075  ret->SetBinContent(ix, yax->GetBinLowEdge(iy)+frac*yax->GetBinWidth(iy));
1076  // Intro of http://www.jstor.org/stable/2286545
1077  ret->SetBinError(ix, sqrt(tot)*yax->GetBinWidth(iy)/(2*z));
1078  break;
1079  }
1080  } // end for iy
1081  } // end for ix
1082  return ret;
1083  }
1084 
1085  //......................................................................
1087  {
1088  TAxis* xax = h->GetXaxis();
1089  TAxis* yax = h->GetYaxis();
1090  TH1* ret = new TH1D("", "", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
1091  ret->GetXaxis()->SetTitle(xax->GetTitle());
1092  for(int ix = 0; ix < xax->GetNbins()+2; ++ix){
1093  double tot = 0;
1094  // Make sure to include overflow bins
1095  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
1096  tot += h->GetBinContent(ix, iy);
1097  }
1098 
1099  double avg = 0;
1100  double rms = 0;
1101  double seen = 0;
1102  double inc = 0;
1103  for(int iy = 0; iy < yax->GetNbins()+2; ++iy){
1104  const double z = h->GetBinContent(ix, iy);
1105  if(!z) continue;
1106 
1107  seen += z;
1108 
1109  // Interpolation
1110  double weight_lo = (seen-(.25*tot))/z;
1111  double weight_hi = 1-(seen-(.75*tot))/z;
1112  if(weight_lo > 1) weight_lo = 1;
1113  if(weight_hi > 1) weight_hi = 1;
1114  double weight = weight_lo+weight_hi-1;
1115  if(weight < 0) weight = 0;
1116 
1117  inc += z*weight;
1118  const double c = yax->GetBinCenter(iy);
1119  avg += z*weight*c;
1120  rms += z*weight*util::sqr(c);
1121  } // end for iy
1122 
1123  // In principle this should be problematic because of floating point errors
1124  assert(inc == tot/2);
1125 
1126  if(inc){
1127  avg /= inc;
1128  rms /= inc;
1129  rms -= util::sqr(avg);
1130  assert(!std::isinf(avg) && !std::isnan(avg));
1131  assert(!std::isinf(rms) && !std::isnan(rms));
1132  assert(rms >= 0);
1133  ret->SetBinContent(ix, avg);
1134  ret->SetBinError(ix, sqrt(rms)/sqrt(tot/2));
1135  }
1136  else{
1137  ret->SetBinContent(ix, 0);
1138  ret->SetBinError(ix, 1e10);
1139  }
1140  }
1141 
1142  return ret;
1143  }
1144 
1145  //......................................................................
1147  {
1148  //std::cout << "////////////////////////////// Apply Threshold correction /////////////////////////////" << std::endl;
1149  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
1150  const double w = prof->GetXaxis()->GetBinCenter(i);
1151  const double y = prof->GetBinContent(i);
1152  const double corr = fThreshCorrMap->ThresholdCorrAt(view, cell, w);
1153  //LOG_VERBATIM("AttenuationFit") << "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;
1154  prof->SetBinContent(i, y/corr);
1155  prof->SetBinError(i, prof->GetBinError(i)/corr);
1156  }
1157  //std::cout << "////////////////////////////// Done Threshold correction /////////////////////////////" << std::endl;
1158  }
1159 
1160  //......................................................................
1161  void AttenuationFit::ApplyThresholdCorrection(int view, int cell, TH1* prof, int fiberbrightness)
1162  {
1163  //std::cout << "////////////////////////////// Apply Threshold correction /////////////////////////////" << std::endl;
1164  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
1165  const double w = prof->GetXaxis()->GetBinCenter(i);
1166  const double y = prof->GetBinContent(i);
1167  const double corr = fThreshCorrMap->ThresholdCorrAt(view, cell, w, fiberbrightness);
1168  //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;
1169  prof->SetBinContent(i, y/corr);
1170  prof->SetBinError(i, prof->GetBinError(i)/corr);
1171  }
1172  //std::cout << "////////////////////////////// Done Threshold correction /////////////////////////////" << std::endl;
1173  }
1174 
1175  //......................................................................
1176  TF1* AttenuationFit::expFit(TH1* prof,
1177  TString opt,
1179  calib::AttenCurve*& res)
1180  {
1181  std::cout<<"expFit"<<std::endl;
1182 
1183  int view = GetView(chan.Plane());
1184  double range = 0.;
1185  if(fConsolidateViewsMode /*&& !(fFiberBrightnessMode||fHybridMode)*/){
1186  view = (chan.Plane());
1187  // range = (view == geo::kX) ? fGeom->DetHalfWidth() : fGeom->DetHalfHeight();
1188  range = (view == geo::kX) ? fGeom->DetHalfHeight() : fGeom->DetHalfWidth();
1189  }
1191  view = chan.Plane()-fGeom->FirstPlaneInMuonCatcher()+2;
1192  range = (view == geo::kX) ? fGeom->DetHalfHeight() : fGeom->DetHalfWidth();
1193  }
1194  else{
1195  range = fGeom->Plane(chan.Plane())->Cell(chan.Cell())->HalfL();
1196  }
1197 
1198  res->cell_length = 2*range;
1200  range -= 50;
1201  }
1202  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1203  res->cell_length =range+50+fInitEdgeMuXMax;
1204  }
1205  res->coeff_low = 0;
1206  res->coeff_high = 0;
1207  res->edge_low = 0;
1208  res->edge_high = 0;
1209 
1210  struct ExpFitAdaptor
1211  {
1212  ExpFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
1213  double operator()(double* xs, double* ps)
1214  {
1215  fRes->coeff_exp = ps[0];
1216  fRes->atten_length = ps[1];
1217  fRes->background = ps[2];
1218  return fRes->MeanPEPerCmAt(xs[0]);
1219  }
1220  calib::AttenCurve* fRes;
1221  };
1222  ExpFitAdaptor* adapt = new ExpFitAdaptor(res);
1223 
1224  TString fitStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
1225  TString fit_name = fitStrArr[view];
1226  double rangeup=range;
1227  std::cout << " View : " << view << std::endl;
1228  if(view == 2 && fGeom->DetId() == novadaq::cnv::kNEARDET){
1229  std::cout << "****************" << std::endl;
1230  rangeup=fInitEdgeMuXMax;
1231  }
1232 
1233  TF1* fit = new TF1(fit_name, adapt, -range, +rangeup, 3);
1234  // Empirically, the fits end up close to here
1235  fit->SetParameter(0, 11);
1237  fit->SetParameter(0, 43);
1238  fit->SetParLimits(0, 0, 1000);
1239  fit->SetParLimits(1, 0, 1000);
1240  fit->SetParLimits(2, 0, 1000);
1241  }
1242  if(view == 2 && fGeom->DetId() == novadaq::cnv::kNEARDET){
1243  fit->SetParameter(0, 70);
1244  }
1246  fit->SetParameter(0, 48);
1247  }
1248  fit->SetParLimits(0, 0, 1000);
1249  fit->SetParLimits(1, 0, 1000);
1250  fit->SetParLimits(2, 0, 1000);
1251 
1252  fit->SetParameter(1, 400);//kTrueShortAtten);
1253  fit->SetParameter(2, 7.5);
1254  fit->SetParName(0, "norm");
1255  fit->SetParName(1, "atten");
1256  fit->SetParName(2, "bkgd");
1257 
1258  fit->SetRange(-range, +rangeup);
1259  res->chisq = 0; // Need to count as calibrated during fit
1260  TFitResultPtr fr = prof->Fit(fit, "RS"+opt);
1261  LOG_VERBATIM("AttenuationFit")<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell();
1262  LOG_VERBATIM("AttenuationFit")<<fit->GetParameter(0)<<" "<<fit->GetParameter(1)<<" "<<fit->GetParameter(2);
1263 
1264  if(FitQuality(prof, fit)>0.2 && view==2){
1265  //std::cout<<"Default "<<FitQuality(prof, fit)<<" "<<chan.Plane()<<" "<<chan.Cell()<<std::endl;
1266  //70 400 7.5
1267  fit->SetParameter(0, 43);
1268  fit->SetParameter(1, 452);//kTrueShortAtten);
1269  fit->SetParameter(2, 7.6);
1270  res->atten_length = 73;
1271  res->coeff_exp = 452;
1272  res->background = 7.6;
1273  res->chisq = FitQuality(prof, fit);
1274  //res->chisq = 0;
1275 
1276  return fit;
1277  }
1278 
1279  res->atten_length = fit->GetParameter(0);
1280  res->coeff_exp = fit->GetParameter(1);
1281  res->background = fit->GetParameter(2);
1282  res->chisq = FitQuality(prof, fit);
1283 
1284 
1285  // const int dof = prof->FindBin(+range)-prof->FindBin(-range)-3;
1286  // res->chisq =(fr == 0) ? fr->Chi2()/dof : -1;
1287 
1288  return fit;
1289  }
1290 
1291  //......................................................................
1293  TString opt,
1295  calib::AttenCurve*& res)
1296  {
1297  struct RolloffFitAdaptor
1298  {
1299  RolloffFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
1300  double operator()(double* xs, double* ps)
1301  {
1302  fRes->coeff_low = ps[0];
1303  fRes->coeff_high = ps[1];
1304  fRes->edge_low = ps[2];
1305  fRes->edge_high = ps[3];
1306  return fRes->MeanPEPerCmAt(xs[0]);
1307  }
1308  calib::AttenCurve* fRes;
1309  };
1310  RolloffFitAdaptor* adapt = new RolloffFitAdaptor(res);
1311 
1312  int view = GetView(chan.Plane());
1313  if(fConsolidateViewsMode/*&&!(fFiberBrightnessMode||fHybridMode)*/) view = (chan.Plane());
1315 
1316  const double edge = (view == geo::kX) ? fInitEdgeX : fInitEdgeY;
1317  // const double range = (view == geo::kX) ? fGeom->DetHalfWidth() : fGeom->DetHalfHeight();
1318  const double range = (view == geo::kX) ? fGeom->DetHalfHeight() : fGeom->DetHalfWidth();
1319  double edgeup=edge;
1320  double rangeup=range;
1321  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1322  edgeup=fInitEdgeMuXMax;
1323  rangeup=fInitEdgeMuXMax+50;
1324  }
1325 
1326  TF1* fit = new TF1(view == geo::kX ? "rollFitX" : "rollFitY",
1327  adapt, -range, +rangeup, 4);
1328  fit->SetParName(0, "coeffL");
1329  fit->SetParName(1, "coeffR");
1330  fit->SetParName(2, "edgeL");
1331  fit->SetParName(3, "edgeR");
1332  fit->SetParameter(0, 0);
1333  fit->SetParameter(1, 0);
1334  fit->SetParameter(2, -edge);
1335  fit->SetParameter(3, edgeup);
1336  fit->SetRange(-range, +rangeup);
1337 
1338  if(fIncludeRolloffFit&&view!=2){
1339  res->chisq = 0; // Need to count as calibrated during fit
1340  TFitResultPtr fr = prof->Fit(fit, "RS"+opt);
1341 
1342  res->coeff_low = fit->GetParameter(0);
1343  res->coeff_high = fit->GetParameter(1);
1344  res->edge_low = fit->GetParameter(2);
1345  res->edge_high = fit->GetParameter(3);
1346 
1347  res->chisq = FitQuality(prof, fit);
1348 
1349  // const int dof = prof->FindBin(+range)-prof->FindBin(-range)-7;
1350  // res->chisq = (fr == 0) ? fr->Chi2()/dof : -1;
1351  }
1352 
1353  return fit;
1354  }
1355 
1356  //......................................................................
1357  TGraph* AttenuationFit::lowessFit(TH1* prof, int plane, calib::AttenCurve*& res)
1358  {
1359  int view = GetView(plane);
1360  if(fConsolidateViewsMode) view = plane;
1362  TGraph* ret = new TGraph;
1363  double x0 = (view == geo::kX) ? -fGeom->DetHalfHeight() : -fGeom->DetHalfWidth();
1364  double x1 = (view == geo::kX) ? +fGeom->DetHalfHeight() : +fGeom->DetHalfWidth();
1365  const int kNumControlPts = 20;
1366  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1367  x0 = -fGeom->DetHalfHeight();
1368  x1 = fInitEdgeMuXMax+50;
1369  }
1370  // Range of measurements that influence each control point
1371  const double kRange = 1.5*(x1-x0)/kNumControlPts;
1372 
1373  for(int n = 0; n < kNumControlPts; ++n){
1374  const double x = x0+double(n)/(kNumControlPts-1)*(x1-x0);
1375 
1376  std::vector<double> xs, ys, ws;
1377 
1378  for(int i = 0; i < prof->GetNbinsX()+2; ++i){
1379  const double xi = prof->GetBinCenter(i);
1380  if(xi < x0 || xi > x1) continue;
1381  if(fabs(xi-x) > kRange) continue;
1382 
1383  xs.push_back(xi);
1384  ys.push_back(prof->GetBinContent(i));
1385  // Apparently "tricube" weights are traditional
1386  ws.push_back(util::cube(1-util::cube(fabs(xi-x)/kRange)));
1387  }
1388 
1389  double m, c;
1390 
1391  util::LinFit(xs, ys, ws, m, c);
1392 
1393  const double y = m*x+c;
1394 
1395  ret->SetPoint(n, x, y);
1396 
1397  res->AddInterpPoint(x, y);
1398  }
1399  ret->SetLineWidth(2);
1400  ret->SetMarkerColor(kRed);
1401  ret->SetLineColor(kRed);
1402  ret->SetMarkerStyle(kFullCircle);
1403  return ret;
1404  }
1405 
1406  //......................................................................
1408  {
1409  for(int sign = -1; sign <= +1; sign += 2){
1410  TGraph* g = new TGraph;
1411  g->SetPoint(0, sign*range, 0);
1412  g->SetPoint(1, sign*range, 1e4);
1413  g->SetLineStyle(7);
1414  g->Draw("l same");
1415  }
1416  }
1417 
1418  //......................................................................
1419  void AttenuationFit::DrawDetectorEdges(double rangedown, double rangeup)
1420  {
1421  TGraph* g = new TGraph;
1422  g->SetPoint(0, rangedown, 0);
1423  g->SetPoint(1, rangedown, 1e4);
1424  g->SetLineStyle(7);
1425  g->Draw("l same");
1426  TGraph* g2 = new TGraph;
1427  g2->SetPoint(0, rangeup, 0);
1428  g2->SetPoint(1, rangeup, 1e4);
1429  g2->SetLineStyle(7);
1430  g2->Draw("l same");
1431 
1432  return;
1433  }
1434 
1435  //......................................................................
1437  int plane,
1438  int cell,
1439  calib::AttenCurve*& res,
1440  int fiberbrightness)
1441  {
1442 
1443  std::cout<<"fit_channel_prof"<<std::endl;
1444 
1445  //art::ServiceHandle<photrans::FiberBrightness> fbhandle;
1446  //int fiberbrightness = fbhandle->getBrightnessIndex(plane, cell);
1447 
1448  int view = GetView(plane);
1449  if(fConsolidateViewsMode) view = plane;
1451 
1452  if(fFiberBrightnessMode)LOG_VERBATIM("AttenuationFit") << "\t View: " << view << ", fiber brightness: " << fiberbrightness << ", cell: " << cell;
1453  else LOG_VERBATIM("AttenuationFit") << "\t View: " << view << ", cell: " << cell;
1454 
1455  const TString dataMCStr = fIsMC ? "_mc" : "_data";
1456  const TString viewStrArr[4]={"fitX","fitY","fitMuX","fitMuY"};
1457  const TString viewStr =viewStrArr[view];
1458 
1459  LOG_VERBATIM("AttenuationFit") << "\t data? " << dataMCStr << " viewStr? " << viewStr;
1460 
1461  res->initialized = true;
1462 
1463  double cell_length = (view == geo::kX) ? 2 * fGeom->DetHalfHeight() : 2 * fGeom->DetHalfWidth();
1464  res->cell_length = cell_length;
1465  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1466  res->cell_length =(cell_length/2)+50+fInitEdgeMuXMax;
1467  }
1468 
1469  TString suffix = dataMCStr+viewStr;
1470  if(cell >= 0){
1471  if(fFiberBrightnessMode) suffix += TString::Format("_fb%i", fiberbrightness);
1472  suffix += TString::Format("_%03d_%03d", plane, cell);
1473  }
1474  TString quiet;
1475  if(cell >= 0) quiet = "Q";
1476 
1477  // TString title = "NDOS cosmic ";
1478  TString title;
1479  if (fGeom->DetId() == novadaq::cnv::kFARDET) title = "FD cosmic ";
1480  if (fGeom->DetId() == novadaq::cnv::kNEARDET) title = "ND cosmic ";
1481  if (fGeom->DetId() == novadaq::cnv::kTESTBEAM) title = "TB cosmic ";
1482  if (fGeom->DetId() == novadaq::cnv::kNDOS) title = "NDOS cosmic ";
1483  if (fIsMC) title += "Monte Carlo"; else title += "data";
1484 
1485  TString viewStrArrVertHor[4]={"vertical","horizontal","vertical","horizontal"};
1486  TString viewStrVertHor =viewStrArrVertHor[view];
1487 
1489  title += TString::Format(" - %s view, cell %d, FB %d",
1490  viewStrVertHor.Data(),
1491  cell,
1492  fiberbrightness);
1493  } else {
1494  title += TString::Format(" - plane %d (%s), cell %d",
1495  plane,
1496  viewStrVertHor.Data(),
1497  cell);
1498  }
1499 
1500  LOG_VERBATIM("AttenuationFit") << "\t title: " << title;
1501 
1502  new TCanvas;
1503  //TCanvas * c1 = new TCanvas("c1");
1504  //c1->cd();
1505 
1506  TF1* fit = expFit(atten, quiet, geo::OfflineChan(plane, cell), res);
1507 
1508  atten->Draw();
1509  atten->SetTitle(title+";Distance from center (cm);Mean PE / cm");
1510 
1511  LOG_VERBATIM("AttenuationFit") << "\t atten drawn";
1512 
1513  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1515  }
1516  else{
1517  DrawDetectorEdges(cell_length/2);
1518  }
1519 
1520  const double attenErr = fit->GetParError(0); /* unused */
1521 
1522  atten->GetYaxis()->SetRangeUser(0, 120);
1523 
1524  TLegend* leg = new TLegend(.3, .15, .7, .35);
1525  leg->SetFillStyle(0);
1526 
1527  leg->AddEntry((TObject*)0, TString::Format("Attenuation length = %d#pm%.1lf cm", int(res->atten_length), attenErr), "");
1528  leg->AddEntry((TObject*)0, TString::Format("Background = %.1lf", res->background), "");
1529  leg->Draw();
1530 
1531  std::string plotsDirectory = fPlotsDirectory;
1532  //plotsDirectory = TString::Format("%s/fb%i", fPlotsDirectory.c_str(), fiberbrightness);
1533  //gSystem->mkdir(plotsDirectory.c_str());
1534  plotsDirectory = TString::Format("%s/%03d", fPlotsDirectory.c_str(), plane);
1535  gSystem->mkdir(plotsDirectory.c_str());
1536 
1537  if(!fPlotsDirectory.empty()){
1538  gPad->Print(plotsDirectory+"/atten"+suffix+".eps");
1539  }
1540 
1541  new TCanvas;
1542  TH1D* profCorr = (TH1D*)atten->Clone();
1543  profCorr->GetYaxis()->SetTitle("Ratio to exponential fit");
1544  fit->SetRange(-1.1*cell_length, +1.1*cell_length);
1545  profCorr->Divide(fit);
1546  profCorr->Draw();
1547 
1548  TF1* rollFit = rolloffFit(atten, quiet, geo::OfflineChan(plane, cell), res);
1549 
1550  if(res->coeff_exp == 0)
1551  mf::LogWarning("AttenuationFit") << "Fit failed: plane = " << plane << ", cell = " << cell;
1552 
1553  profCorr->GetYaxis()->SetRangeUser(.4, 1.1);
1554 
1555  leg = new TLegend(.3, .15, .7, .35);
1556  leg->SetFillStyle(0);
1557 
1558  leg->AddEntry((TObject*)0, TString::Format("Left edge = %.1lf cm", res->edge_low), "");
1559  leg->AddEntry((TObject*)0, TString::Format("Right edge = %.1lf cm", res->edge_high), "");
1560  leg->Draw();
1561 
1562  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1564  }
1565  else{
1566  DrawDetectorEdges(cell_length/2);
1567  }
1568 
1569  if(!plotsDirectory.empty()){
1570  gPad->Print(plotsDirectory+"/rolloff"+suffix+".eps");
1571  }
1572 
1573 
1574  new TCanvas;
1575  TH1* profFullCorr = (TH1*)atten->Clone();
1576  profFullCorr->GetYaxis()->SetTitle("Ratio to full fit");
1577  rollFit->SetRange(-800, +800);
1579  rollFit->SetRange(-250, +250);
1580  }
1582  rollFit->SetRange(-150, +150);
1583  }
1584  profFullCorr->Divide(rollFit);
1585  profFullCorr->GetYaxis()->SetRangeUser(0.9, 1.1);
1586  profFullCorr->Draw();
1587 
1588  TGraph* one = new TGraph;
1589  one->SetPoint(0, -1e4, 1);
1590  one->SetPoint(1, +1e4, 1);
1591  one->SetLineStyle(7);
1592  one->Draw("l same");
1593 
1594  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1596  }
1597  else{
1598  DrawDetectorEdges(cell_length/2);
1599  }
1600 
1601  if(fIncludeLowessFit){
1602  TGraph* lowess = lowessFit(profFullCorr, plane, res);
1603  lowess->Draw("lp same");
1604  }
1605 
1606  res->initialized = true;
1607 
1608  if(!plotsDirectory.empty()){
1609  gPad->Print(plotsDirectory+"/residual"+suffix+".eps");
1610  }
1611 
1612 
1613  new TCanvas;
1614  TH1* plotProf = (TH1*)atten->Clone();
1615  plotProf->Draw();
1616 
1617  struct FullFitAdaptor
1618  {
1619  FullFitAdaptor(calib::AttenCurve* res) : fRes(res) {}
1620  double operator()(double* xs, double* ps)
1621  {
1622  return fRes->MeanPEPerCmAt(xs[0]);
1623  }
1624  calib::AttenCurve* fRes;
1625  } adapt(res);
1626 
1627  const double range = cell_length/2;
1628  double rangeup=range;
1629  if(view == 2 && fGeom->DetId() == novadaq::cnv::kNEARDET){
1630  rangeup = fInitEdgeMuXMax+50;
1631  }
1632  TF1* fullFit = new TF1("fullFit", adapt, -range, +rangeup, 0);
1633  fullFit->SetLineColor(kBlue);
1634  fullFit->SetRange(-range, +rangeup);
1635  fullFit->Draw("same");
1636 
1637  plotProf->SetLineColor(kBlack);
1638  plotProf->GetYaxis()->SetTitle("Mean PE / cm");
1639  plotProf->GetYaxis()->SetRangeUser(0, 150);
1640 
1641  // HACKs for making a pretty plot for blessing
1642  plotProf->GetXaxis()->CenterTitle();
1643  plotProf->GetYaxis()->CenterTitle();
1644  plotProf->GetYaxis()->SetTitleSize(plotProf->GetXaxis()->GetTitleSize());
1645  plotProf->GetYaxis()->SetTitleOffset(plotProf->GetXaxis()->GetTitleOffset());
1646 
1647  if(fGeom->DetId() == novadaq::cnv::kNEARDET&&view==2){
1649  }
1650  else{
1651  DrawDetectorEdges(cell_length/2);
1652  }
1653 
1654  if(!plotsDirectory.empty()){
1655  gPad->Print(plotsDirectory+"/totfit"+suffix+".eps");
1656  gPad->Print(plotsDirectory+"/totfit"+suffix+".png");
1657  }
1658 
1659  new TCanvas;
1660  TH1* resid = new TH1F("", ";W (cm);Ratio to total fit", atten->GetNbinsX(), atten->GetXaxis()->GetXmin(), atten->GetXaxis()->GetXmax());
1661 
1662  for(int n = 0; n < atten->GetNbinsX()+2; ++n){
1663  resid->SetBinContent(n, atten->GetBinContent(n)/res->MeanPEPerCmAt(atten->GetBinCenter(n)));
1664  }
1665  resid->GetYaxis()->SetRangeUser(.9, 1.1);
1666  resid->Draw();
1667  one->Draw("l same");
1668 
1669  if(!plotsDirectory.empty()){
1670  gPad->Print(plotsDirectory+"/residual_final"+suffix+".eps");
1671  }
1672 
1673  delete fit;
1674  delete rollFit;
1675  delete profCorr;
1676  delete profFullCorr;
1677  delete one;
1678  delete leg;
1679  delete plotProf;
1680  delete fullFit;
1681  delete resid;
1682  }
1683 
1684  //......................................................................
1685  void AttenuationFit::WriteDummyCSVRow(FILE* f, int plane, int cell) const
1686  {
1687  std::cout<<"WriteDummyCSVrow"<<std::endl;
1688 
1689  const int vldTime = 1; // 1970
1690 
1691  geo::OfflineChan chan(plane, cell);
1692 
1693  if(fConsolidateViewsMode && !fIsMC){
1694  assert(plane == 0 || plane == 1 || plane == 2 || plane == 3);
1695  chan = geo::OfflineChan(1000+plane, cell);
1696  }
1698  int pln_muc = plane - fGeom->FirstPlaneInMuonCatcher() + 2;
1699  chan = geo::OfflineChan(1000+pln_muc, cell);
1700  }
1701 
1702  fprintf(f, "-1,-1,-1,-1,-1,-1,-1,1000,%d,%d\n",
1703  chan.ToDBValidityChan(), vldTime);
1704  }
1705 
1707 
1708 } // end namespace
1709 ////////////////////////////////////////////////////////////////////////
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double MinW() const
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
int TotalEntries() const
void readResults(art::Results const &run)
enum BeamMode kRed
bool HasProfiles(const geo::OfflineChan &c) const
std::map< std::string, double > xmax
TH1F * ToTH1F() const
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TGraph * lowessFit(TH1 *prof, int plane, calib::AttenCurve *&res)
const Var weight
double ThresholdCorrAt(int view, int cell, double w) const
LiteProfile WPE_corr_xy
Definition: AttenProfiles.h:85
Float_t x1[n_points_granero]
Definition: compare.C:5
art::ServiceHandle< geo::Geometry > fGeom
static AttenCurve * Uninitialized(int det, geo::OfflineChan chan)
Return a new AttenCurve objects with fields uninitialized.
Definition: AttenCurve.cxx:25
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
std::vector< int > const & Count() const
void writeResults(art::Results &r)
AttenProfiles for many channels.
Definition: AttenProfiles.h:89
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
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
const PlaneGeo * Plane(unsigned int i) const
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
TH2 * GetBestHistogram(caldp::AttenHists v, caldp::AttenProfiles u)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
std::map< int, caldp::AttenProfilesMap > fChannelMapProf
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
art::ServiceHandle< nova::dbi::RunHistoryService > fRunHist
LiteProfile WPE_corr_traj
Definition: AttenProfiles.h:85
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
fclose(fg1)
LiteProfile WPE_corr_xy_truth
Definition: AttenProfiles.h:85
caldp::AttenHistsMap fChannelMap
This is where GetChannelHist stores its histograms. Map is from block.
TF1 * rolloffFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
std::vector< float > const & Sum() const
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
void beginRun(art::Run const &run)
ThresholdCorrMap * fThreshCorrMap
art::ServiceHandle< photrans::FiberBrightness > fFBhandle
double GetAttenScale(rb::CellHit const &cellhit, double w)
for PE->PECorr conversion
Prototype Near Detector on the surface at FNAL.
#define DEFINE_ART_RESULTS_PLUGIN(klass)
void reconfigure(fhicl::ParameterSet const &pset)
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
std::string const & instance() const noexcept
Definition: InputTag.h:60
void WriteDummyCSVRow(FILE *f, int plane, int cell) const
Ensure every cell has a CSV entry even if it couldn&#39;t be fit.
double FitQuality(TH1 *prof, TF1 *fit) const
AttenHists & GetHists(const geo::OfflineChan &c)
Definition: AttenHists.h:81
Near Detector in the NuMI cavern.
Histograms used by attenuation calibration.
Definition: AttenHists.h:27
double frac(double x)
Fractional part.
double DetHalfHeight() const
caldp::AttenProfiles GetChannelProfiles(bool isData, int plane, int cell)
caldp::AttenHists GetChannelHists(bool isData, int plane, int cell)
float bin[41]
Definition: plottest35.C:14
TF1 * expFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
z
Definition: test.py:28
AttenHists for many channels.
Definition: AttenHists.h:77
Definition: run.py:1
#define LOG_WARNING(category)
OStream cout
Definition: OStream.cxx:6
unsigned short Plane() const
Definition: OfflineChan.h:31
int nplanes
Definition: geom.C:145
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Profiles used by attenuation calibration.
Definition: AttenProfiles.h:61
geo::OfflineChan chan
Definition: AttenCurve.h:53
unsigned short Cell() const
Definition: OfflineChan.h:32
void ApplyThresholdCorrection(int view, int cell, TH1 *prof)
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
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
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
int ncells
Definition: geom.C:124
float MeanPEPerCmAt(double w) const
Mean response of this channel at this distance from detector centre.
Definition: AttenCurve.cxx:37
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
#define LOG_VERBATIM(category)
auto one()
Definition: PMNS.cxx:49
unsigned int NPlanes() const
#define LOG_INFO(stream)
Definition: Messenger.h:144
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue
Float_t w
Definition: plot.C:20
TH2F * WPE_corr_traj
Definition: AttenHists.h:67
#define W(x)
AttenProfiles & GetProfiles(const geo::OfflineChan &c)
std::string const & label() const noexcept
Definition: InputTag.h:55
AttenuationFit(fhicl::ParameterSet const &pset)
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
void DrawDetectorEdges(double range)
def sign(x)
Definition: canMan.py:197
void fit_channel_prof(TH1 *atten, int plane, int cell, calib::AttenCurve *&res, int fiberbrightness=0)
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.
enum BeamMode string