DedxDistribution.cxx
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 /// \file DedxDistribution.cxx
4 /// \brief This is a helper class for ParticleIDAlg that provides a tidy structure
5 /// in which to hold the dE/dx histograms and provides access methods to the data
6 ///
7 ///
8 ///
9 /// \version $Id: DedxDistribution.cxx,v 1.0 2013-10-04 10:00:00 absmith Exp $
10 /// \author Alex Smith (smith@physics.umn.edu)
11 ////////////////////////////////////////////////////////////////////////
12 
15 #include "RecoBase/CellHit.h"
16 #include "RecoBase/WeightedHit.h"
17 #include "RecoBase/Cluster.h"
18 #include "RecoBase/Track.h"
19 #include "Utilities/func/MathUtil.h"
21 
22 #include "TGeoManager.h"
23 #include "TGeoNode.h"
24 
25 #include <iostream>
26 #include "TH1.h"
27 
28 namespace slid {
29 
30  /*********************************************************************************************/
31  ////////////////////////////////////////////////////////////////////////
32 
34  : fEmptyLongHist(0), fEmptyTransHist(0)
35  {
36  }
37 
38 
39  /*********************************************************************************************/
40  ////////////////////////////////////////////////////////////////////////
41 
43 
44  for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
45  for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
46  //
47  // delete longitudinal dE/dx histograms
48  //
49  for (int iPlane = 0; iPlane < kNumLongitudinalPlane; iPlane++) {
50  TH1F* h = fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane];
51  if(h != fEmptyLongHist) delete h;
52  }
53  //
54  // delete transverse dE/dx histograms
55  //
56 
57  for (int iPlane = 0; iPlane < kNumTransversePlane; iPlane++) {
58  TH1* h = fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iPlane];
59  if(h != fEmptyTransHist) delete h;
60  }
61  }
62  }
63 
64  delete fEmptyLongHist;
65  delete fEmptyTransHist;
66  }
67 
68  /*********************************************************************************************/
69  ////////////////////////////////////////////////////////////////////////
70 
72  std::string fileNameFclTag) {
73 
74  std::string libPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
75  std::string thistname = pset.get< std::string >(fileNameFclTag);
76  /// FD dE/dx input histogram files
77  std::string fileName = libPath + thistname;
78  LoadHistFromFile(fileName);
79  return;
80  }
81 
83 
84  TFile* fDistFile = new TFile(filePath.c_str(), "READ");
85  LoadDedxHistograms(fDistFile);
86  fDistFile->Close();
87  delete fDistFile;
88  }
89 
90  ////////////////////////////////////////////////////////////////////////
91  TH1F* DedxDistribution::TH1DToTH1F(TH1D* h, bool trans)
92  {
93  // A surprisingly large fraction (~1/2) of these histograms are empty. No
94  // need to have multiple seperate copies of identical empty histograms.
95  if(h->GetEntries() == 0){
96  // If we already have an appropriate histogram cached, that's what we'll
97  // return.
98  TH1F* ret = 0;
99  if(trans && fEmptyTransHist) ret = fEmptyTransHist;
100  if(!trans && fEmptyLongHist) ret = fEmptyLongHist;
101 
102  // Check the binning of the cached histogram in case we ever try to have
103  // variable binning within the longitudunal or transverse histogram sets.
104  if(ret &&
105  (h->GetNbinsX() != ret->GetNbinsX() ||
106  h->GetXaxis()->GetXmin() != ret->GetXaxis()->GetXmin() ||
107  h->GetXaxis()->GetXmax() != ret->GetXaxis()->GetXmax())){
108  std::cout << "Error in DedxDistribution::TH1DToTH1F. "
109  << "Cached empty " << (trans ? "transverse" : "longitudinal")
110  << " histogram binning does not match loaded histogram. "
111  << ret->GetNbinsX() << " bins from "
112  << ret->GetXaxis()->GetXmin() << " to "
113  << ret->GetXaxis()->GetXmax() << " vs "
114  << h->GetNbinsX() << " bins from "
115  << h->GetXaxis()->GetXmin() << " to "
116  << h->GetXaxis()->GetXmax() << std::endl;
117  abort();
118  }
119 
120  if(ret){
121  delete h;
122  return ret;
123  }
124  }
125 
126  // Make a histogram the same size, but with floats. Don't copy the title
127  // which is usually long and descriptive and takes up a surprising amount
128  // of storage.
129  TH1F* ret = new TH1F(h->GetName(),
130  "",
131  h->GetNbinsX(),
132  h->GetXaxis()->GetXmin(),
133  h->GetXaxis()->GetXmax());
134 
135  // Copy the bin contents across
136  for(int i = 0; i < h->GetNbinsX() + 2; ++i){
137  ret->SetBinContent(i, h->GetBinContent(i));
138  }
139 
140  // We also need to set the number of entries which is checked by some code
141  // that uses these histograms.
142  ret->SetEntries(h->GetEntries());
143  ret->SetDirectory(0);
144 
145  delete h;
146 
147  // If the histogram is empty and we got here, which means the cache is also
148  // empty, set the cache entry to this histogram.
149  if(ret->GetEntries() == 0){
150  if(trans)
152  else
154  }
155 
156  return ret;
157  }
158 
159  /*********************************************************************************************/
160  ////////////////////////////////////////////////////////////////////////
161 
162  bool DedxDistribution::LoadDedxHistograms(TFile* fDistFile) {
163 
164  //The detector is divided into 4 regions in XY, these regions extend in z.
165  //Region 1 is the top left quadrant in XY when looking at the front of the detector
166  //The training histograms are then divided into 11 energy regions from 0 to 5 GeV.
167  //In the training each shower is put into a bin based on its reco energy
168  //The first and last energy bin are 0-0.25 and 4.75-5 GeV. All bins inbetween are 0.5 GeV wide.
169  //For each region and each energy there is a de/dx histogram for each plane from 0 to 200
170  //the histograms being loaded below are for longitudinal dedx.
171  //These training histograms are also broken up by particle type (electron, muon, gamma, proton, pi0, pion, neutron)
172 
173  //NOTE: Jianming is exploring using the shower length as a training variable.
174  //In order to do this the variables below are initialized and then store the
175  //shower length at each energy in each region. This is done by finding the
176  //last plane in which there was energy deposited.
177  //These variables are not currently being used in the production version of the algorithm
178  for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
179  for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
180  fHtExpPlane[iXYRegion][iEnergyBin] = 0;
181  }
182  }
183  //
184  // Load histograms for longitudinal dE/dx
185  //
186  for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
187  for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
188  for (int iPlane = 0; iPlane < kNumLongitudinalPlane; iPlane++) {
189  char cht[200];
190  sprintf(cht, "ht_%d_%d_%d", iXYRegion, iEnergyBin, iPlane);
191  //convert TH1D to TH1F in order to save memory
192  fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane] = TH1DToTH1F((TH1D*)fDistFile->Get(cht), false);
193 
194  if (fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetEntries() > 0 && fHtExpPlane[iXYRegion][iEnergyBin] == 0) {
195  if (fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetBinContent(1) / (float) fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetEntries() > 0.997) {
196  fHtExpPlane[iXYRegion][iEnergyBin] = iPlane;
197  }
198  }
199  }
200  }
201  }
202 
203  //
204  // Load histograms for transverse dE/dx
205  //
206  for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
207  for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
208  for (int iPlane = 0; iPlane < kNumTransversePlane; iPlane++) {
209  char cht[200];
210  sprintf(cht, "httrans_%d_%d_%d", iXYRegion, iEnergyBin, iPlane);
211  //convert TH1D to TH1F in order to save memory.
212  fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iPlane] = TH1DToTH1F((TH1D*)fDistFile->Get(cht), true);
213  }
214  }
215  }
216  return true;
217  }
218 
219 
220 
221  /// \brief Transverse plane hit dE/dx.
222 
223  TH1F* DedxDistribution::GetLongitudinalPlaneHitDedxHist(int iXYRegion, int iEnergyBin, int iLongitudinalPlane) {
224  return fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iLongitudinalPlane];
225  }
226  /// \brief Transverse plane hit dE/dx.
227 
228  TH1F* DedxDistribution::GetTransversePlaneCellDedxHist(int iXYRegion, int iEnergyBin, int iTransversePlane) {
229  return fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iTransversePlane];
230  }
231 
232 
233 
234 
235 
236 }// End of namespace slid
void Initialize(const fhicl::ParameterSet &pset, std::string fileNameFclTag)
Load the histogram files.
void LoadHistFromFile(std::string filePath)
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
fileName
Definition: plotROC.py:78
unsigned fHtExpPlane[kNumXYRegion][kNumEnergyBin]
Need to look up what this in in JM&#39;s code and give it a more descriptive name (expected hits...
This is a helper class for ParticleIDAlg that provides a tidy structure in which to hold the dE/dx hi...
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
TH1F * GetLongitudinalPlaneHitDedxHist(int iXYRegion, int iEnergyBin, int iLongitudinalPlane)
Longitudinal plane hit dE/dx.
static const int kNumXYRegion
Number of XY regions into which detector is divided for dE/dx histograms.
cout<< t1-> GetEntries()
Definition: plottest35.C:29
TH1F * TH1DToTH1F(TH1D *h, bool trans)
T get(std::string const &key) const
Definition: ParameterSet.h:231
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool LoadDedxHistograms(TFile *fDistFile)
Load all of the histograms from the file into arrays of TH1F* objects. The details of handling the hi...
static const int kNumEnergyBin
Number of energy bins into which detector is divided for dE/dx histograms.
TH1F * fTransversePlaneCellDedx[kNumXYRegion][kNumEnergyBin][kNumTransversePlane]
Transverse plane hit dE/dx.
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
static const int kNumLongitudinalPlane
Number of longitudinal planes considered for dE/dx histograms.
TH1F * GetTransversePlaneCellDedxHist(int iXYRegion, int iEnergyBin, int iTransversePlane)
Transverse plane hit dE/dx ("planes" transverse to shower axis).
TH1F * fLongitudinalPlaneHitDedx[kNumXYRegion][kNumEnergyBin][kNumLongitudinalPlane]
Longitudinal plane hit dE/dx.