AttenProfiles.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 /// \brief Histograms used by attenuation calibration
3 /// \author Christopher Backhouse - bckhouse@caltech.edu
4 /// \date Apr 2013
5 ///////////////////////////////////////////////////////////////////////////
8 
9 namespace caldp
10 {
11  //----------------------------------------------------------------------------
12  LiteProfile::LiteProfile(float w0, float w1)
13  : fMinW(w0)
14  , fMaxW(w1)
15  {
16  }
17 
18  //----------------------------------------------------------------------------
19  void LiteProfile::Fill(float x, float y)
20  {
21  const unsigned int bin = ((x-fMinW)*kNumBins)/(fMaxW-fMinW);
22  if(x < fMinW || bin > kNumBins) return; // Drop overflows
23 
24  assert(y >= 0);
25 
26  // only reserve the memory necessary for the
27  // the vectors when needed, not at construction
28  // assume if the fCount vector is not at the size
29  // of kNumBins, then none of them are
30  if(fCount.size() < kNumBins){
31  fCount.resize(kNumBins, 0 );
32  fSum .resize(kNumBins, 0.);
33  fSqSum.resize(kNumBins, 0.);
34  }
35 
36  ++fCount[bin];
37  fSum[bin] += y;
38  fSqSum[bin] += y*y;
39  }
40 
41  //----------------------------------------------------------------------------
43  {
44  assert(rhs.fMinW == fMinW);
45  assert(rhs.fMaxW == fMaxW);
46 
47  // only reserve the memory necessary for the
48  // the vectors when needed, not at construction
49  // assume if the fCount vector is not at the size
50  // of kNumBins, then none of them are
51  if(fCount.size() < kNumBins){
52  fCount.resize(kNumBins, 0 );
53  fSum .resize(kNumBins, 0.);
54  fSqSum.resize(kNumBins, 0.);
55  }
56 
57  // only add the profiles if they have the same
58  // number of entried for each vector
59  if(fCount.size() == rhs.fCount.size() &&
60  fSum.size() == rhs.fSum.size() &&
61  fSqSum.size() == rhs.fSqSum.size() ){
62 
63  for(size_t i = 0; i < fCount.size(); ++i){
64  LOG_DEBUG("AttenProfiles") << i << " " << rhs.fCount[i] << " " << rhs.fSum[i] << " " << rhs.fSqSum[i];
65  fCount[i] += rhs.fCount[i];
66  fSum[i] += rhs.fSum[i];
67  fSqSum[i] += rhs.fSqSum[i];
68  }
69  }
70 
71  return *this;
72  }
73 
74  //----------------------------------------------------------------------------
75  TH1F* LiteProfile::ToTH1F() const
76  {
77  TH1F* h = new TH1F("", "", kNumBins, fMinW, fMaxW);
78 
79  for(size_t i = 0; i < kNumBins; ++i){
80  const int N = fCount[i];
81  if(N < 1) continue;
82  h->SetEntries(h->GetEntries()+N);
83 
84  h->SetBinContent(i+1, fSum[i]/N);
85 
86  const float arg = fSqSum[i]-fSum[i]*fSum[i]/N;
87  if(arg >= 0){
88  h->SetBinError(i+1, std::sqrt(arg)/N);
89  }
90  else{
91  h->SetBinError(i+1, fSum[i]/N);
92  }
93  }
94 
95  return h;
96  }
97 
98  //----------------------------------------------------------------------------
100  {
101  int ret = 0;
102  for(int n: fSum) ret += n;
103  return ret;
104  }
105 
106  //----------------------------------------------------------------------------
108  {
109  int ret =0;
110  for(size_t i = 0; i < kNumBins; ++i){
111  ret+=fCount[i];
112  }
113 
114  return ret;
115  }
116 
117  //----------------------------------------------------------------------------
118  std::vector<int> const& LiteProfile::Count() const
119  {
120  return fCount;
121  }
122 
123  //----------------------------------------------------------------------------
124  std::vector<float> const& LiteProfile::Sum() const
125  {
126  return fSum;
127  }
128 
129  //----------------------------------------------------------------------------
130  std::vector<float> const& LiteProfile::SqSum() const
131  {
132  return fSqSum;
133  }
134 
135  //----------------------------------------------------------------------------
136  AttenProfiles::AttenProfiles(float minW, float maxW)
137  : WPE(minW, maxW)
138  , WPE_corr(minW, maxW)
139  , WPE_corr_xy(minW, maxW)
140  , WPE_corr_z(minW, maxW)
141  , WPE_corr_traj(minW, maxW)
142  , WPE_corr_xy_truth(minW, maxW)
143  {
144  }
145 
146  //----------------------------------------------------------------------------
148  {
149  WPE += rhs.WPE;
150  WPE_corr += rhs.WPE_corr;
151  WPE_corr_xy += rhs.WPE_corr_xy;
152  WPE_corr_z += rhs.WPE_corr_z;
155 
156  return *this;
157  }
158 
159  //----------------------------------------------------------------------------
160  AttenProfilesMap::AttenProfilesMap(float w0, float w1) : fMinW(w0), fMaxW(w1) {}
161 
162  //----------------------------------------------------------------------------
164  {
165  return fMap.find(c) != fMap.end();
166  }
167 
168  //----------------------------------------------------------------------------
170  {
171  return HasProfiles(geo::OfflineChan(plane, cell));
172  }
173 
174  //----------------------------------------------------------------------------
176  {
177  auto it = fMap.find(c);
178  if(it != fMap.end()) return it->second;
179 
181  return fMap.find(c)->second;
182  }
183 
184  //----------------------------------------------------------------------------
186  {
187  return GetProfiles(geo::OfflineChan(plane, cell));
188  }
189 
190  //----------------------------------------------------------------------------
191  /// List of channels with profiles on
192  std::vector<geo::OfflineChan> AttenProfilesMap::Channels() const
193  {
194  std::vector<geo::OfflineChan> ret;
195  ret.reserve(fMap.size());
196  for(auto& it: fMap) ret.push_back(it.first);
197  return ret;
198  }
199 
200  //----------------------------------------------------------------------------
201  /// Intended for writing out/drawing all the histograms
202  std::map<int, std::map<int, std::vector<TH1F*>>> AttenProfilesMap::GetAllProfilesByPlaneAndCell()
203  {
204  std::map<int, std::map<int, std::vector<TH1F*>>> ret;
205  for(auto& it: fMap){
206  const int plane = it.first.Plane();
207  const int cell = it.first.Cell();
208  const AttenProfiles& hs = it.second;
209  std::vector<TH1F*>& v = ret[plane][cell];
210  /* if(hs.WPE.TotalHits()){
211  v.push_back(hs.WPE.ToTH1F());
212  v.back()->SetName(TString::Format("h_WPE_%03d_%03d", plane, cell));
213  v.back()->SetTitle(TString::Format("Plane %d, Cell %d;W;PE", plane, cell));
214  }
215 
216  if(hs.WPE_corr.TotalHits()){
217  v.push_back(hs.WPE_corr.ToTH1F());
218  v.back()->SetName(TString::Format("h_WPE_corr_avg_%03d_%03d", plane, cell));
219  v.back()->SetTitle(TString::Format("Avg Corr, Plane %d, Cell %d;W;PE/cm", plane, cell));
220  }*/
221 
222  if(hs.WPE_corr_xy.TotalHits()){
223  v.push_back(hs.WPE_corr_xy.ToTH1F());
224  v.back()->SetName(TString::Format("h_WPE_corr_xy_%03d_%03d", plane, cell));
225  v.back()->SetTitle(TString::Format("XY, Plane %d, Cell %d;W;PE/cm", plane, cell));
226  }
227 
228  if(hs.WPE_corr_z.TotalHits()){
229  v.push_back(hs.WPE_corr_z.ToTH1F());
230  v.back()->SetName(TString::Format("h_WPE_corr_z_%03d_%03d", plane, cell));
231  v.back()->SetTitle(TString::Format("Z, Plane %d, Cell %d;W;PE/cm", plane, cell));
232  }
233 
234  if(hs.WPE_corr_traj.TotalHits()){
235  v.push_back(hs.WPE_corr_traj.ToTH1F());
236  v.back()->SetName(TString::Format("h_WPE_corr_traj_%03d_%03d", plane, cell));
237  v.back()->SetTitle(TString::Format("Traj, Plane %d, Cell %d;W;PE/cm", plane, cell));
238  }
239 
240  /*
241  if(hs.WPE_corr_xy_truth.TotalHits()){
242  v.push_back(hs.WPE_corr_xy_truth.ToTH1F());
243  v.back()->SetName(TString::Format("h_WPE_corr_xy_truth_%03d_%03d", plane, cell));
244  v.back()->SetTitle(TString::Format("XY, Plane %d, Cell %d;W;MeV/cm", plane, cell));
245  }*/
246  }
247  return ret;
248  }
249 
250  //----------------------------------------------------------------------------
252  {
253  for(auto & it: rhs.fMap){
254  GetProfiles(it.first) += it.second;
255  }
256  return *this;
257  }
258 
259  //----------------------------------------------------------------------------
261  {
262  fMap.clear();
263  }
264 
265 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int TotalEntries() const
bool HasProfiles(const geo::OfflineChan &c) const
TH1F * ToTH1F() const
set< int >::iterator it
LiteProfile WPE_corr_xy
Definition: AttenProfiles.h:85
std::vector< int > fCount
Definition: AttenProfiles.h:55
T sqrt(T number)
Definition: d0nt_math.hpp:156
void Fill(float x, float y)
std::vector< int > const & Count() const
AttenProfiles for many channels.
Definition: AttenProfiles.h:89
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
AttenProfilesMap & operator+=(const AttenProfilesMap &rhs)
int TotalHits() const
LiteProfile WPE_corr_traj
Definition: AttenProfiles.h:85
LiteProfile WPE_corr_xy_truth
Definition: AttenProfiles.h:85
std::vector< float > const & Sum() const
LiteProfile & operator+=(LiteProfile const &rhs)
std::map< geo::OfflineChan, AttenProfiles > fMap
float bin[41]
Definition: plottest35.C:14
static const size_t kNumBins
Definition: AttenProfiles.h:53
AttenProfiles & operator+=(AttenProfiles const &rhs)
std::vector< float > fSum
Definition: AttenProfiles.h:56
Profiles used by attenuation calibration.
Definition: AttenProfiles.h:61
Histograms used by attenuation calibration.
Lightweight TProfile equivalent.
Definition: AttenProfiles.h:24
std::vector< float > fSqSum
Definition: AttenProfiles.h:57
std::map< int, std::map< int, std::vector< TH1F * > > > GetAllProfilesByPlaneAndCell()
Intended for writing out/drawing all the histograms.
LiteProfile WPE_corr_z
Definition: AttenProfiles.h:85
A (plane, cell) pair.
Definition: OfflineChan.h:17
LiteProfile WPE_corr
Definition: AttenProfiles.h:85
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
AttenProfiles & GetProfiles(const geo::OfflineChan &c)
std::vector< float > const & SqSum() const