BPFdEdxHistoMaker_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BPFdEdxHistoMaker_module.cc
3 // \brief Module to create dE/dx vs. beta*gamma plots from BPF tracks
4 // \version
5 // \author $Author: mbaird42 $
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <cmath>
10 #include <iostream>
11 #include <vector>
12 
13 // ROOT includes
14 #include "TH2F.h"
15 
16 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
28 
29 // NOvASoft includes
30 #include "RecoBase/Cluster.h"
31 #include "RecoBase/Prong.h"
32 #include "RecoBase/Track.h"
33 #include "RecoBase/FitSum.h"
34 
35 #include "Utilities/AssociationUtil.h"
36 #include "Geometry/Geometry.h"
37 #include "Calibrator/Calibrator.h"
38 
40 
41 
42 
43 // BPFdEdxHistoMaker header
44 namespace bpfit {
46  public:
47  explicit BPFdEdxHistoMaker(fhicl::ParameterSet const& pset);
49 
50  void analyze(const art::Event& evt);
51  void reconfigure(const fhicl::ParameterSet& pset);
52  void beginJob();
53  void endJob();
54 
55  private:
60 
61  double fdxTOL; ///< cut off value for dx, if dx is less than this, don't compute dE/dx
62  bool fSmooth; ///< Smooth the likelihood histos at the end?
63  unsigned int fNsmooth; ///< number of bins on either side of the current bin to smooth over
64  bool fNormalize; ///< Normalize the likelihood histos at the end?
65 
66  TH2F *fdEdxVSbgMu[2][2]; ///< Computed dE/dx for cell hits (vs. beta*gamma) for the muon assupmtion.
67  TH2F *fdEdxVSbgPi[2][2]; ///< Computed dE/dx for cell hits (vs. beta*gamma) for the pion assupmtion.
68  TH2F *fdEdxVSbgPr[2][2]; ///< Computed dE/dx for cell hits (vs. beta*gamma) for the proton assupmtion.
69 
70  bpfit::dEdxCalculator fdEdxCalc; ///< helper class for computing dEdx values
71 
72  };
73 }
74 
75 
76 
77 namespace bpfit
78 {
79  //.......................................................................
81  : EDAnalyzer(pset),
82  fSlicerToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("SlicerLabel"))),
83  fdEdxCalc()
84  {
85  this->reconfigure(pset);
86  }
87 
88  //......................................................................
90  {
91  //======================================================================
92  // Clean up any memory allocated by your module... ...if you dare...
93  //======================================================================
94  }
95 
96  //......................................................................
98  {
99  fProngLabel = pset.get<std::string> ("ProngLabel");
100  fProngInstance = pset.get<std::string> ("ProngInstance");
101  fTrackLabel = pset.get<std::string> ("TrackLabel");
102  fdxTOL = pset.get<double> ("dxTOL");
103  fSmooth = pset.get<bool> ("Smooth");
104  fNsmooth = pset.get<unsigned int>("Nsmooth");
105  fNormalize = pset.get<bool> ("Normalize");
106  }
107 
108  //......................................................................
110  {
112 
113  //
114  // Book the NTuples and histograms...
115  //
116 
117  // We are dividing the dE/dx info from cell hits into multiple histograms
118  // based on whether or not the track is steep (large dx) and whether the
119  // cell hit is near to (W <= 0.0) or far from (W > 0.0) the readout. This
120  // is to reflect the fact that the dE/dx calculations might be slightly
121  // different depending on these conditions.
122 
123  // A note on naming of the histos below:
124  //
125  // * {Mu,Pi,Pr} obviously refers to the particle assumption
126  // * W{0,1} refers to whether the track W values was <= 0 (W0) or > 0 (W1)
127  // * X{0,1} refers to whether the dx value was <= 8 cm (X0) or > 8 cm (X1)
128 
129  for(unsigned int w = 0; w < 2; ++w) {
130  for(unsigned int x = 0; x < 2; ++x) {
131 
132  char histoname[30];
133 
134  sprintf(histoname, "fdEdxVSbgMu_W%.1u_X%.1u", w, x);
135  fdEdxVSbgMu[w][x] = tfs->make<TH2F>(histoname,
136  "dE/dx vs. #beta#gamma (muon track assumption);ln(#beta#gamma);dE/dx [GeV/cm]",
137  160,-2.0,6.0,400,0.0,0.04);
138 
139  sprintf(histoname, "fdEdxVSbgPi_W%.1u_X%.1u", w, x);
140  fdEdxVSbgPi[w][x] = tfs->make<TH2F>(histoname,
141  "dE/dx vs. #beta#gamma (pion track assumption);ln(#beta#gamma);dE/dx [GeV/cm]",
142  160,-2.0,6.0,400,0.0,0.04);
143 
144  sprintf(histoname, "fdEdxVSbgPr_W%.1u_X%.1u", w, x);
145  fdEdxVSbgPr[w][x] = tfs->make<TH2F>(histoname,
146  "dE/dx vs. #beta#gamma (proton track assumption);ln(#beta#gamma);dE/dx [GeV/cm]",
147  160,-2.0,6.0,400,0.0,0.04);
148 
149 
150  } // end loop over x
151  } // end loop over w
152 
153  }
154 
155  //......................................................................
157  {
158 
159  //
160  // Apply simple smoothing if requested.
161  //
162  // The smoothing is simple: set each bin to the average of that bin
163  // with the bin immedately above and below. Do NOT include the
164  // overflow bins in the smoothing.
165  //
166 
167  if(fSmooth) {
168 
169  unsigned int NXbins = 0.0;
170  unsigned int NYbins = 0.0;
171  double ave = 0.0;
172 
173  for(unsigned int w = 0; w < 2; ++w) {
174  for(unsigned int x = 0; x < 2; ++x) {
175 
176  // muon histos
177  NXbins = fdEdxVSbgMu[w][x]->GetNbinsX();
178  NYbins = fdEdxVSbgMu[w][x]->GetNbinsY();
179  for(unsigned int i = 1; i <= NXbins; ++i) {
180  for(unsigned int j = 1+fNsmooth; j+fNsmooth <= NYbins; ++j) {
181  ave = 0.0;
182  for(unsigned int s = j-fNsmooth; s <= j+fNsmooth; ++s)
183  ave += fdEdxVSbgMu[w][x]->GetBinContent(i,s);
184  ave = ave/(2*fNsmooth+1);
185  fdEdxVSbgMu[w][x]->SetBinContent(i,j,ave);
186  }
187  }
188 
189  // pion histos
190  NXbins = fdEdxVSbgPi[w][x]->GetNbinsX();
191  NYbins = fdEdxVSbgPi[w][x]->GetNbinsY();
192  for(unsigned int i = 1; i <= NXbins; ++i) {
193  for(unsigned int j = 1+fNsmooth; j+fNsmooth <= NYbins; ++j) {
194  ave = 0.0;
195  for(unsigned int s = j-fNsmooth; s <= j+fNsmooth; ++s)
196  ave += fdEdxVSbgPi[w][x]->GetBinContent(i,s);
197  ave = ave/(2*fNsmooth+1);
198  fdEdxVSbgPi[w][x]->SetBinContent(i,j,ave);
199  }
200  }
201 
202  // proton histos
203  NXbins = fdEdxVSbgPr[w][x]->GetNbinsX();
204  NYbins = fdEdxVSbgPr[w][x]->GetNbinsY();
205  for(unsigned int i = 1; i <= NXbins; ++i) {
206  for(unsigned int j = 1+fNsmooth; j+fNsmooth <= NYbins; ++j) {
207  ave = 0.0;
208  for(unsigned int s = j-fNsmooth; s <= j+fNsmooth; ++s)
209  ave += fdEdxVSbgPr[w][x]->GetBinContent(i,s);
210  ave = ave/(2*fNsmooth+1);
211  fdEdxVSbgPr[w][x]->SetBinContent(i,j,ave);
212  }
213  }
214 
215  } // end loop over x
216  } // end loop over w
217 
218  } // end smoothing
219 
220 
221 
222  //
223  // Normalize the likelihood histos if requested.
224  //
225  // The normalization is so that the area of each vertical stip
226  // (one bin wide) is 1.0. Since each bin represents a probability,
227  // we will include the overflow bins in the normalization.
228  //
229 
230  if(fNormalize) {
231 
232  unsigned int NXbins = 0.0;
233  unsigned int NYbins = 0.0;
234 
235  for(unsigned int w = 0; w < 2; ++w) {
236  for(unsigned int x = 0; x < 2; ++x) {
237 
238  // muon histos
239  NXbins = fdEdxVSbgMu[w][x]->GetNbinsX();
240  NYbins = fdEdxVSbgMu[w][x]->GetNbinsY();
241  for(unsigned int i = 0; i <= NXbins+1; ++i) {
242  float integral = 0.0;
243  for(unsigned int j = 0; j <= NYbins+1; ++j)
244  integral += fdEdxVSbgMu[w][x]->GetBinContent(i,j);
245  if(integral == 0.0) continue;
246  for(unsigned int j = 0; j <= NYbins+1; ++j)
247  fdEdxVSbgMu[w][x]->SetBinContent(i,j,fdEdxVSbgMu[w][x]->GetBinContent(i,j)/integral);
248  }
249 
250  // pion histos
251  NXbins = fdEdxVSbgPi[w][x]->GetNbinsX();
252  NYbins = fdEdxVSbgPi[w][x]->GetNbinsY();
253  for(unsigned int i = 0; i <= NXbins+1; ++i) {
254  float integral = 0.0;
255  for(unsigned int j = 0; j <= NYbins+1; ++j)
256  integral += fdEdxVSbgPi[w][x]->GetBinContent(i,j);
257  if(integral == 0.0) continue;
258  for(unsigned int j = 0; j <= NYbins+1; ++j)
259  fdEdxVSbgPi[w][x]->SetBinContent(i,j,fdEdxVSbgPi[w][x]->GetBinContent(i,j)/integral);
260  }
261 
262  // proton histos
263  NXbins = fdEdxVSbgPr[w][x]->GetNbinsX();
264  NYbins = fdEdxVSbgPr[w][x]->GetNbinsY();
265  for(unsigned int i = 0; i <= NXbins+1; ++i) {
266  float integral = 0.0;
267  for(unsigned int j = 0; j <= NYbins+1; ++j)
268  integral += fdEdxVSbgPr[w][x]->GetBinContent(i,j);
269  if(integral == 0.0) continue;
270  for(unsigned int j = 0; j <= NYbins+1; ++j)
271  fdEdxVSbgPr[w][x]->SetBinContent(i,j,fdEdxVSbgPr[w][x]->GetBinContent(i,j)/integral);
272  }
273 
274  } // end loop over x
275  } // end loop over w
276 
277  } // end normalization
278 
279  }
280 
281  //......................................................................
283  {
284 
285  // get the slices
287  evt.getByToken(fSlicerToken, slices);
288 
289  // create the FindMany object for getting prongs
291 
292  // make a geometry service handle
294 
295  // make a calibrator service handle
297 
298 
299 
300  // loop over all slices and get the fuzzyk 3D prongs
301  for(unsigned int islice = 0; islice < slices->size(); ++islice) {
302 
303  // As usual, skip the noise slice...
304  if((*slices)[islice].IsNoise()) continue;
305 
306  // NOTE: Currently, the fundamental object of analysis is the slice.
307  // If in the future we change things so that the fundamental object
308  // of analysis becomes the vertex, then we will have to get prongs
309  // associated with the vertex instead.
310 
311  // get the 3D prongs associated with this slice
312  std::vector< art::Ptr< rb::Prong > > prongs;
313  if(prongFMP.isValid()) prongs = prongFMP.at(islice);
314 
315  // create the FindMany object for getting tracks
316  art::FindManyP<rb::Track> trackFMP(prongs, evt, fTrackLabel);
317 
318  // loop over all prongs and get the tracks
319  for(unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
320 
321  // get the tracks associated with this prong
322  std::vector< art::Ptr< rb::Track > > tracks;
323  if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
324 
325  // create the FindMany object for getting FitSum objects
326  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
327 
328  // Loop over all tracks...
329  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
330 
331  // get the FitSum object associated with this track
332  std::vector< art::Ptr< rb::FitSum > > fitsums;
333  if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
334 
335  // There can be only one!!!
336  // (something went really wrong with BPF if there's not one and
337  // only one of these per track.)
338  if(fitsums.size() != 1) abort();
339 
340  // Note to the user:
341  // It is intended that this program only runs over a sample of
342  // muons that have been tracked only under the muon assumption. I
343  // haven't forced that asumption here.
344 
345  //
346  // Perform the dE/dx calculations with the helper class.
347  //
348  fdEdxCalc.computeDEDX(tracks[itrack],fitsums[0]->PDG());
349 
350  std::vector<double> dE;
351  std::vector<double> dx;
352  std::vector<double> W;
353  std::vector<double> BG;
354  std::vector<double> s; // this is a throw away variable...
355 
356  fdEdxCalc.getResults(dE,dx,W,BG,s,fdxTOL);
357 
358  for(unsigned int a = 0; a < dE.size(); ++a) {
359 
360  unsigned int w = 0;
361  unsigned int x = 0;
362 
363  if(W[a] > 0.0) w = 1;
364  if(dx[a] > 8.0) x = 1;
365 
366  if(fitsums[0]->PDG() == 13)
367  fdEdxVSbgMu[w][x]->Fill(log(BG[a]),dE[a]/dx[a]);
368  else if(fitsums[0]->PDG() == 211)
369  fdEdxVSbgPi[w][x]->Fill(log(BG[a]),dE[a]/dx[a]);
370  else if(fitsums[0]->PDG() == 2212)
371  fdEdxVSbgPr[w][x]->Fill(log(BG[a]),dE[a]/dx[a]);
372  else {
373  std::cout << "\n\nERROR: PDG code " << fitsums[0]->PDG()
374  << " is not supported! Aborting...\n\n";
375  abort();
376  }
377  }
378 
379  } // end loop over tracks (itrack)
380 
381  } // end loop over prongs (iprong)
382 
383  } // end loop over slices (islice)
384 
385  return;
386 
387  }
388 
389 } // end namespace bpfit
390 ////////////////////////////////////////////////////////////////////////
391 
392 
393 
395 
396 namespace bpfit
397 {
399 }
void analyze(const art::Event &evt)
TH2F * fdEdxVSbgPr[2][2]
Computed dE/dx for cell hits (vs. beta*gamma) for the proton assupmtion.
void reconfigure(const fhicl::ParameterSet &pset)
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
void computeDEDX(art::Ptr< rb::Track > &track, int pdg=13)
Compute dE/dx for all cells along this track and the fitsum object that went with it...
DEFINE_ART_MODULE(TestTMapFile)
double dE
Calculates dE/dx in cells passed through by a track.
double fdxTOL
cut off value for dx, if dx is less than this, don&#39;t compute dE/dx
const XML_Char * s
Definition: expat.h:262
BPFdEdxHistoMaker(fhicl::ParameterSet const &pset)
double dx[NP][NC]
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
bool fNormalize
Normalize the likelihood histos at the end?
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Perform a "2 point" Hough transform on a collection of hits.
bpfit::dEdxCalculator fdEdxCalc
helper class for computing dEdx values
OStream cout
Definition: OStream.cxx:6
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
bool fSmooth
Smooth the likelihood histos at the end?
T * make(ARGS...args) const
void geom(int which=0)
Definition: geom.C:163
TH2F * fdEdxVSbgPi[2][2]
Computed dE/dx for cell hits (vs. beta*gamma) for the pion assupmtion.
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Float_t w
Definition: plot.C:20
#define W(x)
ProductToken< T > consumes(InputTag const &)
void getResults(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &Wtemp, std::vector< double > &bgtemp, std::vector< double > &stemp, double dxTol)
Return the results of the dE/dx calculation.
Encapsulate the geometry of one entire detector (near, far, ndos)
TH2F * fdEdxVSbgMu[2][2]
Computed dE/dx for cell hits (vs. beta*gamma) for the muon assupmtion.
unsigned int fNsmooth
number of bins on either side of the current bin to smooth over
enum BeamMode string