All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TCShowerElectronLikelihood_module.cc
Go to the documentation of this file.
1 // -------------------------------------------------
2 // makes shower profile templates
3 //
4 // Author: Rory Fitzpatrick (roryfitz@umich.edu)
5 // Created: 8/3/18
6 // -------------------------------------------------
7 
8 // Framework includes
9 #include "art/Framework/Core/EDAnalyzer.h"
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Principal/Handle.h"
13 #include "art/Framework/Services/Registry/ServiceHandle.h"
14 #include "art_root_io/TFileService.h"
15 #include "canvas/Persistency/Common/FindManyP.h"
16 #include "cetlib/pow.h"
17 #include "fhiclcpp/ParameterSet.h"
18 
19 #include "larcore/Geometry/Geometry.h"
20 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
21 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
22 #include "lardataobj/RecoBase/Hit.h"
23 #include "lardataobj/RecoBase/Shower.h"
24 #include "larreco/Calorimetry/CalorimetryAlg.h"
25 #include "nusimdata/SimulationBase/MCNeutrino.h"
26 #include "nusimdata/SimulationBase/MCParticle.h"
27 #include "nusimdata/SimulationBase/MCTruth.h"
28 
29 #include "TFile.h"
30 #include "TH1.h"
31 #include "TH3.h"
32 #include "TProfile.h"
33 #include "TProfile2D.h"
34 
35 namespace shower {
36 
37  class TCShowerElectronLikelihood : public art::EDAnalyzer {
38  public:
39  explicit TCShowerElectronLikelihood(fhicl::ParameterSet const& pset);
40 
41  private:
42  void beginJob() override;
43  void analyze(const art::Event& evt) override;
44 
45  void getShowerProfile(detinfo::DetectorClocksData const& clockdata,
46  detinfo::DetectorPropertiesData const& detProp,
47  std::vector<art::Ptr<recob::Hit>> showerhits,
48  TVector3 shwvtx,
49  TVector3 shwdir);
50  void findEnergyBin();
51  double getLongLikelihood();
52  double getTranLikelihood();
53 
54  void resetProfiles();
55 
56  std::string fTemplateFile;
57  std::string fROOTfile;
58 
59  TH3F* longTemplate;
60  TH3F* tranTemplate;
61  TH3F* tranTemplate_1;
62  TH3F* tranTemplate_2;
63  TH3F* tranTemplate_3;
64  TH3F* tranTemplate_4;
65  TH3F* tranTemplate_5;
66  TProfile2D* longTemplateProf2D;
67  TProfile2D* tranTemplateProf2D;
68  TProfile2D* tranTemplateProf2D_1;
69  TProfile2D* tranTemplateProf2D_2;
70  TProfile2D* tranTemplateProf2D_3;
71  TProfile2D* tranTemplateProf2D_4;
72  TProfile2D* tranTemplateProf2D_5;
73 
74  //TTree* fTree;
75  TH1F* energyDist;
76  TH1F* longLikelihoodHist;
77  TH1F* tranLikelihoodHist;
78  // just for printing purposes
79  TH1F* longProfHist;
80  TH1F* tranProfHist_1;
81  TH1F* tranProfHist_2;
82  TH1F* tranProfHist_3;
83  TH1F* tranProfHist_4;
84  TH1F* tranProfHist_5;
85 
86  TH1F* longProfile;
87  TH1F* tranProfile;
88  TH1F* tranProfile_1;
89  TH1F* tranProfile_2;
90  TH1F* tranProfile_3;
91  TH1F* tranProfile_4;
92  TH1F* tranProfile_5;
93 
94  int energyGuess;
95  double energyChi2;
96  int maxt;
97 
98  const int LBINS = 20;
99  const int LMIN = 0;
100  const int LMAX = 5;
101 
102  const int TBINS = 20;
103  const int TMIN = -5;
104  const int TMAX = 5;
105 
106  const double X0 = 14;
107 
108  std::string fHitModuleLabel;
109  std::string fShowerModuleLabel;
110  std::string fTemplateModuleLabel;
111  std::string fGenieGenModuleLabel;
112  std::string fDigitModuleLabel;
113 
114  calo::CalorimetryAlg fCalorimetryAlg;
115 
116  }; // class TCShowerElectronLikelihood
117 
118 } // shower
119 
120 // -------------------------------------------------
121 
122 shower::TCShowerElectronLikelihood::TCShowerElectronLikelihood(fhicl::ParameterSet const& pset)
123  : EDAnalyzer(pset)
124  , fHitModuleLabel(pset.get<std::string>("HitModuleLabel", "trajcluster"))
125  , fShowerModuleLabel(pset.get<std::string>("ShowerModuleLabel", "tcshower"))
126  , fGenieGenModuleLabel(pset.get<std::string>("GenieGenModuleLabel", "generator"))
127  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
128 {
129  fTemplateFile = pset.get<std::string>("TemplateFile");
130  cet::search_path sp("FW_SEARCH_PATH");
131  if (!sp.find_file(fTemplateFile, fROOTfile))
132  throw cet::exception("TCShowerElectronLikelihood")
133  << "cannot find the root template file: \n"
134  << fTemplateFile << "\n bail ungracefully.\n";
135 
136  std::unique_ptr<TFile> file{TFile::Open(fROOTfile.c_str())};
137 
138  longTemplate = file->Get<TH3F>("tcshowertemplate/fLongitudinal");
139  tranTemplate = file->Get<TH3F>("tcshowertemplate/fTransverse");
140  tranTemplate_1 = file->Get<TH3F>("tcshowertemplate/fTransverse_1");
141  tranTemplate_2 = file->Get<TH3F>("tcshowertemplate/fTransverse_2");
142  tranTemplate_3 = file->Get<TH3F>("tcshowertemplate/fTransverse_3");
143  tranTemplate_4 = file->Get<TH3F>("tcshowertemplate/fTransverse_4");
144  tranTemplate_5 = file->Get<TH3F>("tcshowertemplate/fTransverse_5");
145 
146  longTemplateProf2D = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoLong2D");
147  tranTemplateProf2D = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D");
148  tranTemplateProf2D_1 = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D_1");
149  tranTemplateProf2D_2 = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D_2");
150  tranTemplateProf2D_3 = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D_3");
151  tranTemplateProf2D_4 = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D_4");
152  tranTemplateProf2D_5 = file->Get<TProfile2D>("tcshowertemplate/fShowerProfileRecoTrans2D_5");
153 
154  longProfile = new TH1F("longProfile", "longitudinal shower profile;t;Q", LBINS, LMIN, LMAX);
155  tranProfile = new TH1F("tranProfile", "transverse shower profile;dist (cm);Q", TBINS, TMIN, TMAX);
156 
157  tranProfile_1 = new TH1F(
158  "tranProfile_1", "transverse shower profile [0 <= t < 1];dist (cm);Q", TBINS, TMIN, TMAX);
159  tranProfile_2 = new TH1F(
160  "tranProfile_2", "transverse shower profile [1 <= t < 2];dist (cm);Q", TBINS, TMIN, TMAX);
161  tranProfile_3 = new TH1F(
162  "tranProfile_3", "transverse shower profile [2 <= t < 3];dist (cm);Q", TBINS, TMIN, TMAX);
163  tranProfile_4 = new TH1F(
164  "tranProfile_4", "transverse shower profile [3 <= t < 4];dist (cm);Q", TBINS, TMIN, TMAX);
165  tranProfile_5 = new TH1F(
166  "tranProfile_5", "transverse shower profile [4 <= t < 5];dist (cm);Q", TBINS, TMIN, TMAX);
167 } // TCShowerElectronLikelihood
168 
169 // -------------------------------------------------
170 
171 void
172 shower::TCShowerElectronLikelihood::beginJob()
173 {
174 
175  art::ServiceHandle<art::TFileService const> tfs;
176 
177  energyDist = tfs->make<TH1F>("energyDist", "true energy - guess energy", 41, -20.5, 20.5);
178  longLikelihoodHist = tfs->make<TH1F>("longLikelihoodHist", "longitudinal likelihood", 20, 0, 2);
179  tranLikelihoodHist = tfs->make<TH1F>("tranLikelihoodHist", "transverse likelihood", 20, 0, 3);
180 
181  // just for printing purposes
182  longProfHist =
183  tfs->make<TH1F>("longProfHist", "longitudinal e- profile (reco);t;Q", LBINS, LMIN, LMAX);
184  tranProfHist_1 = tfs->make<TH1F>(
185  "tranProfHist", "transverse e- profile (reco) [0 <= t < 1];dist (cm);Q", TBINS, TMIN, TMAX);
186  tranProfHist_2 = tfs->make<TH1F>(
187  "tranProfHist", "transverse e- profile (reco) [1 <= t < 2];dist (cm);Q", TBINS, TMIN, TMAX);
188  tranProfHist_3 = tfs->make<TH1F>(
189  "tranProfHist", "transverse e- profile (reco) [2 <= t < 3];dist (cm);Q", TBINS, TMIN, TMAX);
190  tranProfHist_4 = tfs->make<TH1F>(
191  "tranProfHist", "transverse e- profile (reco) [3 <= t < 4];dist (cm);Q", TBINS, TMIN, TMAX);
192  tranProfHist_5 = tfs->make<TH1F>(
193  "tranProfHist", "transverse e- profile (reco) [4 <= t < 5];dist (cm);Q", TBINS, TMIN, TMAX);
194 
195 } // beginJob
196 
197 // -------------------------------------------------
198 
199 void
200 shower::TCShowerElectronLikelihood::analyze(const art::Event& evt)
201 {
202  resetProfiles();
203 
204  art::Handle<std::vector<recob::Hit>> hitListHandle;
205  std::vector<art::Ptr<recob::Hit>> hitlist;
206  if (evt.getByLabel(fHitModuleLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle);
207 
208  art::Handle<std::vector<recob::Shower>> showerListHandle;
209  std::vector<art::Ptr<recob::Shower>> showerlist;
210  if (evt.getByLabel(fShowerModuleLabel, showerListHandle))
211  art::fill_ptr_vector(showerlist, showerListHandle);
212 
213  art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
214  std::vector<art::Ptr<simb::MCTruth>> mclist;
215  if (evt.getByLabel(fGenieGenModuleLabel, mctruthListHandle))
216  art::fill_ptr_vector(mclist, mctruthListHandle);
217 
218  art::FindManyP<recob::Hit> shwfm(showerListHandle, evt, fShowerModuleLabel);
219 
220  if (showerlist.size()) {
221  auto const clock_data =
222  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
223  auto const det_prop =
224  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
225  std::vector<art::Ptr<recob::Hit>> showerhits = shwfm.at(0);
226  getShowerProfile(
227  clock_data, det_prop, showerhits, showerlist[0]->ShowerStart(), showerlist[0]->Direction());
228 
229  maxt = std::ceil((90 - showerlist[0]->ShowerStart().Z()) / X0);
230 
231  findEnergyBin();
232 
233  longLikelihoodHist->Fill(getLongLikelihood());
234  tranLikelihoodHist->Fill(getTranLikelihood());
235 
236  // check true shower energy
237  if (mclist.size()) {
238  art::Ptr<simb::MCTruth> mctruth = mclist[0];
239  if (mctruth->NeutrinoSet()) {
240  if (std::abs(mctruth->GetNeutrino().Nu().PdgCode()) == 12 &&
241  mctruth->GetNeutrino().CCNC() == 0) {
242  double elep = mctruth->GetNeutrino().Lepton().E();
243  energyDist->Fill(elep - energyGuess);
244  } // cc nue
245  } // neutrinoset
246  } // mclist
247  } // showerlist
248 
249  longProfHist = longProfile;
250  tranProfHist_1 = tranProfile_1;
251  tranProfHist_2 = tranProfile_2;
252  tranProfHist_3 = tranProfile_3;
253  tranProfHist_4 = tranProfile_4;
254  tranProfHist_5 = tranProfile_5;
255 } // analyze
256 
257 // -------------------------------------------------
258 
259 void
260 shower::TCShowerElectronLikelihood::resetProfiles()
261 {
262 
263  longProfile->Reset();
264  tranProfile->Reset();
265  tranProfile_1->Reset();
266  tranProfile_2->Reset();
267  tranProfile_3->Reset();
268  tranProfile_4->Reset();
269  tranProfile_5->Reset();
270 
271  energyGuess = -9999;
272  energyChi2 = -9999;
273  maxt = -9999;
274 } // resetProfiles
275 
276 // -------------------------------------------------
277 
278 void
279 shower::TCShowerElectronLikelihood::getShowerProfile(detinfo::DetectorClocksData const& clockdata,
280  detinfo::DetectorPropertiesData const& detProp,
281  std::vector<art::Ptr<recob::Hit>> showerhits,
282  TVector3 shwvtx,
283  TVector3 shwdir)
284 {
285  art::ServiceHandle<geo::Geometry const> geom;
286 
287  auto collectionPlane = geo::PlaneID(0, 0, 1);
288 
289  double shwVtxTime = detProp.ConvertXToTicks(shwvtx[0], collectionPlane);
290  double shwVtxWire = geom->WireCoordinate(shwvtx[1], shwvtx[2], collectionPlane);
291 
292  double shwTwoTime = detProp.ConvertXToTicks(shwvtx[0] + shwdir[0], collectionPlane);
293  double shwTwoWire =
294  geom->WireCoordinate(shwvtx[1] + shwdir[1], shwvtx[2] + shwdir[2], collectionPlane);
295 
296  for (size_t i = 0; i < showerhits.size(); ++i) {
297  if (showerhits[i]->WireID().Plane != collectionPlane.Plane) continue;
298 
299  double wirePitch = geom->WirePitch(showerhits[i]->WireID());
300  double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
301  tickToDist *= 1.e-3 * sampling_rate(clockdata); // 1e-3 is conversion of 1/us to 1/ns
302 
303  double xvtx = shwVtxTime * tickToDist;
304  double yvtx = shwVtxWire * wirePitch;
305 
306  double xtwo = shwTwoTime * tickToDist;
307  double ytwo = shwTwoWire * wirePitch;
308 
309  double xtwoorth = (ytwo - yvtx) + xvtx;
310  double ytwoorth = -(xtwo - xvtx) + yvtx;
311 
312  double xhit = showerhits[i]->PeakTime() * tickToDist;
313  double yhit = showerhits[i]->WireID().Wire * wirePitch;
314 
315  double ldist = std::abs((ytwoorth - yvtx) * xhit - (xtwoorth - xvtx) * yhit + xtwoorth * yvtx -
316  ytwoorth * xvtx) /
317  std::hypot(ytwoorth - yvtx, xtwoorth - xvtx);
318  double tdist = ((ytwo - yvtx) * xhit - (xtwo - xvtx) * yhit + xtwo * yvtx - ytwo * xvtx) /
319  std::hypot(ytwo - yvtx, xtwo - xvtx);
320 
321  double to3D = 1. / std::hypot(xvtx - xtwo,
322  yvtx - ytwo); // distance between two points in 3D space is one
323  ldist *= to3D;
324  tdist *= to3D;
325  double t = ldist / X0;
326 
327  double Q = showerhits[i]->Integral() *
328  fCalorimetryAlg.LifetimeCorrection(clockdata, detProp, showerhits[i]->PeakTime());
329 
330  longProfile->Fill(t, Q);
331  tranProfile->Fill(tdist, Q);
332 
333  if (t < 1)
334  tranProfile_1->Fill(tdist, Q);
335  else if (t < 2)
336  tranProfile_2->Fill(tdist, Q);
337  else if (t < 3)
338  tranProfile_3->Fill(tdist, Q);
339  else if (t < 4)
340  tranProfile_4->Fill(tdist, Q);
341  else if (t < 5)
342  tranProfile_5->Fill(tdist, Q);
343 
344  } // loop through showerhits
345 
346  return;
347 
348 } // getShowerProfile
349 
350 // -------------------------------------------------
351 
352 void
353 shower::TCShowerElectronLikelihood::findEnergyBin()
354 {
355 
356  if (longProfile->GetNbinsX() != longTemplate->GetNbinsX())
357  throw cet::exception("TCShowerElectronLikelihood")
358  << "Bin mismatch in longitudinal profile template \n";
359 
360  if (tranProfile->GetNbinsX() != tranTemplate->GetNbinsX())
361  throw cet::exception("TCShowerElectronLikelihood")
362  << "Bin mismatch in transverse profile template \n";
363 
364  double chi2min = 999999;
365  double bestbin = -1;
366 
367  int ebins = longTemplate->GetNbinsY();
368  int lbins = longTemplate->GetNbinsX();
369  int tbins = tranTemplate->GetNbinsX();
370 
371  for (int i = 0; i < ebins; ++i) {
372  double thischi2 = 0;
373 
374  TProfile* ltemp = longTemplateProf2D->ProfileX("_x", i + 1, i + 1);
375  TProfile* ttemp_1 = tranTemplateProf2D_1->ProfileX("_x_1", i + 1, i + 1);
376  TProfile* ttemp_2 = tranTemplateProf2D_2->ProfileX("_x_2", i + 1, i + 1);
377  TProfile* ttemp_3 = tranTemplateProf2D_3->ProfileX("_x_3", i + 1, i + 1);
378  TProfile* ttemp_4 = tranTemplateProf2D_4->ProfileX("_x_4", i + 1, i + 1);
379  TProfile* ttemp_5 = tranTemplateProf2D_5->ProfileX("_x_5", i + 1, i + 1);
380 
381  int nlbins = 0;
382  int ntbins = 0;
383 
384  for (int j = 0; j < lbins; ++j) {
385  double obs = longProfile->GetBinContent(j + 1);
386  double exp = ltemp->GetBinContent(j + 1);
387  if (obs != 0) {
388  thischi2 += cet::square(obs - exp) / exp;
389  ++nlbins;
390  }
391  } // loop through longitudinal bins
392 
393  for (int j = 0; j < tbins; ++j) {
394  double obs = tranProfile_1->GetBinContent(j + 1);
395  double exp = ttemp_1->GetBinContent(j + 1);
396  if (obs != 0) {
397  thischi2 += cet::square(obs - exp) / exp;
398  ++ntbins;
399  }
400 
401  obs = tranProfile_2->GetBinContent(j + 1);
402  exp = ttemp_2->GetBinContent(j + 1);
403  if (obs != 0) {
404  thischi2 += cet::square(obs - exp) / exp;
405  ++ntbins;
406  }
407 
408  obs = tranProfile_3->GetBinContent(j + 1);
409  exp = ttemp_3->GetBinContent(j + 1);
410  if (obs != 0) {
411  thischi2 += cet::square(obs - exp) / exp;
412  ++ntbins;
413  }
414 
415  obs = tranProfile_4->GetBinContent(j + 1);
416  exp = ttemp_4->GetBinContent(j + 1);
417  if (obs != 0) {
418  thischi2 += cet::square(obs - exp) / exp;
419  ++ntbins;
420  }
421 
422  obs = tranProfile_5->GetBinContent(j + 1);
423  exp = ttemp_5->GetBinContent(j + 1);
424  if (obs != 0) {
425  thischi2 += cet::square(obs - exp) / exp;
426  ++ntbins;
427  }
428  } // loop through longitudinal bins
429 
430  thischi2 /= (nlbins + ntbins);
431 
432  if (thischi2 < chi2min) {
433  chi2min = thischi2;
434  bestbin = i;
435  }
436 
437  } // loop through energy bins
438 
439  energyChi2 = chi2min;
440  energyGuess = bestbin + 1;
441 } // findEnergyBin
442 
443 // -------------------------------------------------
444 
445 double
446 shower::TCShowerElectronLikelihood::getLongLikelihood()
447 {
448  if (energyGuess < 0) return -9999.;
449  int energyBin = energyGuess;
450 
451  double longLikelihood = 0.;
452  int nbins = 0;
453 
454  for (int i = 0; i < LBINS; ++i) {
455  double qval = longProfile->GetBinContent(i + 1);
456  int qbin = longTemplate->GetZaxis()->FindBin(qval);
457  int binentries = longTemplate->GetBinContent(i + 1, energyBin, qbin);
458  int totentries = longTemplate->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
459  if (qval > 0) {
460  ++nbins;
461  double prob = (double)binentries / totentries * 100;
462  if (binentries > 0) longLikelihood += log(prob);
463  }
464  } // loop through
465 
466  longLikelihood /= nbins;
467 
468  std::cout << longLikelihood << std::endl;
469  return longLikelihood;
470 } // getLongLikelihood
471 
472 // -------------------------------------------------
473 
474 double
475 shower::TCShowerElectronLikelihood::getTranLikelihood()
476 {
477  if (energyGuess < 0) return -9999.;
478  int energyBin = energyGuess;
479 
480  double tranLikelihood_1 = 0;
481  double tranLikelihood_2 = 0;
482  double tranLikelihood_3 = 0;
483  double tranLikelihood_4 = 0;
484  double tranLikelihood_5 = 0;
485 
486  double qval;
487  int qbin, binentries, totentries;
488 
489  int nbins = 0;
490 
491  for (int i = 0; i < TBINS; ++i) {
492  qval = tranProfile_1->GetBinContent(i + 1);
493  qbin = tranTemplate_1->GetZaxis()->FindBin(qval);
494  binentries = tranTemplate_1->GetBinContent(i + 1, energyBin, qbin);
495  totentries = tranTemplate_1->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
496  if (qval > 0) {
497  ++nbins;
498  double prob = (double)binentries / totentries * 100;
499  if (binentries > 0) tranLikelihood_1 += log(prob);
500  }
501 
502  qval = tranProfile_2->GetBinContent(i + 1);
503  qbin = tranTemplate_2->GetZaxis()->FindBin(qval);
504  binentries = tranTemplate_2->GetBinContent(i + 1, energyBin, qbin);
505  totentries = tranTemplate_2->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
506  if (qval > 0) {
507  ++nbins;
508  double prob = (double)binentries / totentries * 100;
509  if (binentries > 0) tranLikelihood_2 += log(prob);
510  }
511 
512  qval = tranProfile_3->GetBinContent(i + 1);
513  qbin = tranTemplate_3->GetZaxis()->FindBin(qval);
514  binentries = tranTemplate_3->GetBinContent(i + 1, energyBin, qbin);
515  totentries = tranTemplate_3->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
516  if (qval > 0) {
517  ++nbins;
518  double prob = (double)binentries / totentries * 100;
519  if (binentries > 0) tranLikelihood_3 += log(prob);
520  }
521 
522  qval = tranProfile_4->GetBinContent(i + 1);
523  qbin = tranTemplate_4->GetZaxis()->FindBin(qval);
524  binentries = tranTemplate_4->GetBinContent(i + 1, energyBin, qbin);
525  totentries = tranTemplate_4->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
526  if (qval > 0) {
527  ++nbins;
528  double prob = (double)binentries / totentries * 100;
529  if (binentries > 0) tranLikelihood_4 += log(prob);
530  }
531 
532  qval = tranProfile_5->GetBinContent(i + 1);
533  qbin = tranTemplate_5->GetZaxis()->FindBin(qval);
534  binentries = tranTemplate_5->GetBinContent(i + 1, energyBin, qbin);
535  totentries = tranTemplate_5->Integral(i + 1, i + 1, energyBin, energyBin, 0, 100);
536  if (qval > 0) {
537  ++nbins;
538  double prob = (double)binentries / totentries * 100;
539  if (binentries > 0) tranLikelihood_5 += log(prob);
540  }
541 
542  } // loop through
543 
544  double const tranLikelihood =
545  (tranLikelihood_1 + tranLikelihood_2 + tranLikelihood_3 + tranLikelihood_4 + tranLikelihood_5) /
546  nbins;
547  std::cout << tranLikelihood << std::endl;
548  return tranLikelihood;
549 } // getTranLikelihood
550 
551 // -------------------------------------------------
552 
553 DEFINE_ART_MODULE(shower::TCShowerElectronLikelihood)
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
shower::TCShowerElectronLikelihood::getShowerProfile
void getShowerProfile(detinfo::DetectorClocksData const &clockdata, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> showerhits, TVector3 shwvtx, TVector3 shwdir)
Definition: TCShowerElectronLikelihood_module.cc:279
file
* file
Definition: file_to_url.sh:69
Hit.h
Declaration of signal hit object.
get
get()
Definition: PMTconfigurationExtraction_module.cc:168
detinfo::DetectorPropertiesData::Temperature
double Temperature() const
In kelvin.
Definition: DetectorPropertiesData.cc:81
shower
process_name shower
Definition: cheaterreco.fcl:51
shower::TCShowerElectronLikelihood::TCShowerElectronLikelihood
TCShowerElectronLikelihood(fhicl::ParameterSet const &pset)
Definition: TCShowerElectronLikelihood_module.cc:122
detinfo::DetectorPropertiesData::Efield
double Efield(unsigned int planegap=0) const
kV/cm
Definition: DetectorPropertiesData.cc:17
geo::fhicl::PlaneID
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
Definition: geo_types_fhicl.h:375
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
detinfo::DetectorPropertiesData::ConvertXToTicks
double ConvertXToTicks(double X, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:99
detinfo::DetectorPropertiesData::DriftVelocity
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Definition: DetectorPropertiesData.cc:23
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
shower::TCShowerElectronLikelihood::analyze
void analyze(const art::Event &evt) override
Definition: TCShowerElectronLikelihood_module.cc:200
calo::CalorimetryAlg::LifetimeCorrection
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:246
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
Generated by   doxygen 1.8.5