ThresholdCorrMap.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ThresholdCorrMap.cxx
3 // \brief Parameterization for use in threshold and shadowing corrections
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
9 #include "TFile.h"
10 #include "TTree.h"
11 
12 #include <cassert>
13 #include <cmath>
14 
15 #include "cetlib/search_path.h"
16 #include "cetlib_except/exception.h"
17 
18 namespace calib
19 {
20  //......................................................................
22  std::string const& fileName,
23  bool FiberBrightnessMode)
24  {
25  std::cout << "ThresholdCorrMap fileName: " << fileName << std::endl;
26 
27  /* std::string fname;
28  std::string filePath("Calibration/");
29  filePath += fileName;
30  cet::search_path sp("FW_SEARCH_PATH");
31  sp.find_file(filePath.c_str(), fname);
32 
33  std::cout << "sp.find_file(" << filePath.c_str() << ", " << fname << std::endl;
34  */
35 
36  switch(det){
38  if(fileName.find("nd") == std::string::npos)//was filePath
39  throw cet::exception("ThresholdCorrMap")
40  << "threshold correction map file " << fileName //was filePath
41  << " does not appear to be valid for current detector"
42  << det;
43  break;
45  if(fileName.find("fd") == std::string::npos)//was filePath
46  throw cet::exception("ThresholdCorrMap")
47  << "threshold correction map file " << fileName //was filePath
48  << " does not appear to be valid for current detector"
49  << det;
50  break;
52  if(fileName.find("tb") == std::string::npos)//was filePath
53  throw cet::exception("ThresholdCorrMap")
54  << "threshold correction map file " << fileName //was filePath
55  << " does not appear to be valid for current detector"
56  << det;
57  break;
58 
59  default:
60  assert(0 && "Unhandled detector, tb test");
61  }
62  assert(!fileName.empty());//was fname
63 
64  TFile corrFile(fileName.c_str());//was fname
65  assert(!corrFile.IsZombie());
66 
67  TTree* tr = (TTree*)corrFile.Get("tree");
68  int view;
69  int fiberbrightness=0;
70  tr->SetBranchAddress("view", &view);
71  if(FiberBrightnessMode) tr->SetBranchAddress("fiberbrightness", &fiberbrightness);
72  Coeffs cs;
73  double cutw, cutwmax;
74  tr->SetBranchAddress("Na",&fNa);
75  tr->SetBranchAddress("Nb",&fNb);
76  tr->SetBranchAddress("Nflag",&fNflag);
77  tr->SetBranchAddress("wpow", &cs.wpow);
78  tr->SetBranchAddress("cellpow", &cs.cellpow);
79  tr->SetBranchAddress("modpow", &cs.modpow);
80  tr->SetBranchAddress("coeff", &cs.coeff);
81  tr->SetBranchAddress("dcellMax", &fCellMax);
82  tr->SetBranchAddress("dcutW", &cutw);
83  tr->SetBranchAddress("dcutWMax", &cutwmax);
84 
85  for(int n = 0; n < tr->GetEntries(); ++n){
86  tr->GetEntry(n);
87  fCoeffs[view][fiberbrightness].push_back(cs);
88  fCutW[view] = cutw;
89  fCutWMax[view] = cutwmax;
90  }
91 
92  corrFile.Close();
93  }
94 
95  //......................................................................
97  {
98  }
99 
100  namespace{
101  /// Raise a double to a (small) integer
102  inline double ipow(double a, int b)
103  {
104  double ret = 1;
105  for(int i = 0; i < b; ++i) ret *= a;
106  return ret;
107  }
108  }
109 
110  //......................................................................
111  double ThresholdCorrMap::ThresholdCorrAt(int view, int cell, double w) const
112  {
113  return ThresholdCorrAt(view,cell,w,0);
114  }
115 
116  //......................................................................
117  double ThresholdCorrMap::ThresholdCorrAt(int view, int cell, double w, int fiberbrightness) const
118  {
119 
120  double kMaxW = fCutWMax[view];
121  double kCutW = fCutW[view];
122  if(w < -kMaxW) w = -kMaxW;
123  if(w > +kCutW) w = +kCutW;
124 
125  const int mod = cell%32; // Position within module
126  double ret = 0;
127 
128  if(fCoeffs[view][fiberbrightness].size()){
129  for(int a = 0; a < fNa; ++a){
130  for(int b = 0; b < fNb; ++b){
131  ret += fCoeffs[view][fiberbrightness][fNb*a+b].coeff*pow(w/kMaxW, a)*pow(cell/fCellMax, b);
132  }
133  }
134  ret+=fCoeffs[view][fiberbrightness][fNa*fNb+mod].coeff;
135  }
136  std::cout << "ThresholdCorrMapAt(" << view << "," << cell << "," << w << "," << fiberbrightness << ") returns " << ret << std::endl;
137 
138  return ret;
139 
140  }
141 }
fileName
Definition: plotROC.py:78
ThresholdCorrMap(novadaq::cnv::DetId det, std::string const &fileName, bool FiberBrightness=false)
double ThresholdCorrAt(int view, int cell, double w) const
constexpr T pow(T x)
Definition: pow.h:75
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Far Detector at Ash River, MN.
CDPStorage service.
const double a
std::vector< Coeffs > fCoeffs[4][12]
Near Detector in the NuMI cavern.
T ipow(T x, unsigned int n)
More efficient exponentiation function than pow(x,n) for small n.
Definition: MathUtil.h:29
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
Float_t w
Definition: plot.C:20