FlashMatchAna_module.cc
Go to the documentation of this file.
1 // -*- mode: c++; c-basic-offset: 2; -*-
2 // This analyzer writes out a TTree for studying the matching
3 // between flashes and events
4 //
5 
6 #ifndef FlashMatchAna_H
7 #define FlashMatchAna_H 1
8 
9 // ROOT includes
10 #include "TH1.h"
11 #include "TEfficiency.h"
12 #include "TTree.h"
13 
14 // C++ includes
15 #include <map>
16 #include <vector>
17 #include <iostream>
18 #include <cstring>
19 #include <sstream>
20 #include "math.h"
21 #include <climits>
22 
23 // LArSoft includes
34 
35 // ART includes.
38 #include "fhiclcpp/ParameterSet.h"
43 #include "art_root_io/TFileService.h"
44 #include "art_root_io/TFileDirectory.h"
47 #include "canvas/Persistency/Common/FindManyP.h"
48 
49 namespace opdet {
50 
52  public:
53 
54  // Standard constructor and destructor for an ART module.
56  virtual ~FlashMatchAna();
57 
58  // This method is called once, at the start of the job. In this
59  // example, it will define the histogram we'll write.
60  void beginJob();
61 
62  // The analyzer routine, called once per event.
63  void analyze (const art::Event&);
64 
65  private:
66 
67  // The stuff below is the part you'll most likely have to change to
68  // go from this custom example to your own task.
69 
70  // The parameters we'll read from the .fcl file.
71  std::string fOpFlashModuleLabel; // Input tag for OpFlash collection
72  std::string fOpHitModuleLabel; // Input tag for OpHit collection
73  std::string fSignalLabel; // Input tag for the signal generator label
74  std::string fGeantLabel; // Input tag for GEANT
75 
76  TTree * fFlashMatchTree;
79 
80 
81  TEfficiency * fRecoEfficiencyVsE;
82  TEfficiency * fRecoEfficiencyVsX;
83  TEfficiency * fRecoEfficiencyVsXandE;
84  TEfficiency * fLargestEfficiencyVsE;
85  TEfficiency * fLargestEfficiencyVsX;
87  TEfficiency * fSelectedEfficiencyVsE;
88  TEfficiency * fSelectedEfficiencyVsX;
90 
91  // Parameters from the fhicl
92  int fNBinsE;
93  float fLowE;
94  float fHighE;
95  int fNBinsX;
96  float fLowX;
97  float fHighX;
98  float fDistanceCut;
99 
100  Int_t fEventID;
101 
102  Float_t fTrueX;
103  Float_t fTrueY;
104  Float_t fTrueZ;
105  Float_t fTrueT;
106  Float_t fDetectedT;
107  Float_t fTrueE;
108  Int_t fTruePDG;
109  Float_t fRecoX;
110 
111  std::vector< Float_t > fTruePxallpart;
112  std::vector< Float_t > fTruePyallpart;
113  std::vector< Float_t > fTruePzallpart;
114  std::vector< Float_t > fTrueEallpart;
115  std::vector< Int_t > fTrueAllPDG;
116 
117  Int_t fNFlashes;
118  std::vector< Int_t > fFlashIDVector;
119  std::vector< Float_t > fYCenterVector;
120  std::vector< Float_t > fZCenterVector;
121  std::vector< Float_t > fYWidthVector;
122  std::vector< Float_t > fZWidthVector;
123  std::vector< Float_t > fTimeVector;
124  std::vector< Float_t > fRecoXVector;
125  std::vector< Float_t > fTimeWidthVector;
126  std::vector< Float_t > fTimeDiffVector;
127  std::vector< Int_t > fFlashFrameVector;
128  std::vector< Bool_t > fInBeamFrameVector;
129  std::vector< Int_t > fOnBeamTimeVector;
130  std::vector< Float_t > fTotalPEVector;
131  Float_t fSumPE;
132  std::vector< Float_t > fPurityVector;
133  std::vector< Float_t > fDistanceVector;
134  Int_t fNOpDets;
135  std::vector<Int_t> fNHitOpDetVector;
136 
137  Int_t fFlashID;
138  Float_t fYCenter;
139  Float_t fZCenter;
140  Float_t fYWidth;
141  Float_t fZWidth;
142  Float_t fTime;
143  Float_t fTimeWidth;
144  Float_t fTimeDiff;
145  //Int_t fFlashFrame; // unused
146  //Bool_t fInBeamFrame; // unused
147  //Int_t fOnBeamTime; // unused
148  Float_t fTotalPE;
149  Float_t fPurity;
150  Float_t fDistance;
151  Int_t fNHitOpDets;
152  std::vector< Float_t > fPEsPerOpDetVector;
153 
154  // For counting waveforms
155 
157  float fBaseline;
158  float fPE;
159  TTree * fCountTree;
163 
164  };
165 
166 }
167 
168 #endif // FlashMatchAna_H
169 
170 namespace opdet {
171 
172  //-----------------------------------------------------------------------
173  // Constructor
175  : EDAnalyzer(pset)
176  {
177 
178  // Indicate that the Input Module comes from .fcl
179  fOpFlashModuleLabel = pset.get<std::string>("OpFlashModuleLabel");
180  fOpHitModuleLabel = pset.get<std::string>("OpHitModuleLabel");
181  fSignalLabel = pset.get<std::string>("SignalLabel");
182  fGeantLabel = pset.get<std::string>("GeantLabel");
183  fNBinsE = pset.get<int>("NBinsE");
184  fLowE = pset.get<float>("LowE");
185  fHighE = pset.get<float>("HighE");
186  fNBinsX = pset.get<int>("NBinsX");
187  fLowX = pset.get<float>("LowX");
188  fHighX = pset.get<float>("HighX");
189  fDistanceCut = pset.get<float>("DistanceCut");
190 
191  fOpDetWaveformLabel = pset.get<std::string>("OpDetWaveformLabel","");
192  fBaseline = pset.get<float>("Baseline", 1500.);
193  fPE = pset.get<float>("PE", 18.);
194 
196 
197  fFlashMatchTree = tfs->make<TTree>("FlashMatchTree","FlashMatchTree");
198  fFlashMatchTree->Branch("EventID", &fEventID, "EventID/I");
199  fFlashMatchTree->Branch("TrueX", &fTrueX, "TrueX/F");
200  fFlashMatchTree->Branch("TrueY", &fTrueY, "TrueY/F");
201  fFlashMatchTree->Branch("TrueZ", &fTrueZ, "TrueZ/F");
202  fFlashMatchTree->Branch("TrueT", &fTrueT, "TrueT/F");
203  fFlashMatchTree->Branch("DetectedT", &fDetectedT, "DetectedT/F");
204  fFlashMatchTree->Branch("TrueE", &fTrueE, "TrueE/F");
205  fFlashMatchTree->Branch("TruePDG", &fTruePDG, "TruePDG/I");
206  fFlashMatchTree->Branch("NFlashes", &fNFlashes, "NFlashes/I");
207  fFlashMatchTree->Branch("FlashIDVector", &fFlashIDVector);
208  fFlashMatchTree->Branch("YCenterVector", &fYCenterVector);
209  fFlashMatchTree->Branch("ZCenterVector", &fZCenterVector);
210  fFlashMatchTree->Branch("YWidthVector", &fYWidthVector);
211  fFlashMatchTree->Branch("ZWidthVector", &fZWidthVector);
212  fFlashMatchTree->Branch("TimeVector", &fTimeVector);
213  fFlashMatchTree->Branch("TimeWidthVector", &fTimeWidthVector);
214  fFlashMatchTree->Branch("TimeDiffVector", &fTimeDiffVector);
215  fFlashMatchTree->Branch("TotalPEVector", &fTotalPEVector);
216  fFlashMatchTree->Branch("SumPE", &fSumPE, "SumPE/F");
217  fFlashMatchTree->Branch("NOpDets", &fNOpDets, "NOpDets/I");
218  fFlashMatchTree->Branch("NHitOpDetVector", &fNHitOpDetVector);
219  fFlashMatchTree->Branch("Purity", &fPurityVector);
220  fFlashMatchTree->Branch("Distance", &fDistanceVector);
221  fFlashMatchTree->Branch("RecoXVector", &fRecoXVector);
222 
223 
224  fLargestFlashTree = tfs->make<TTree>("LargestFlashTree","LargestFlashTree");
225  fLargestFlashTree->Branch("EventID", &fEventID, "EventID/I");
226  fLargestFlashTree->Branch("TrueX", &fTrueX, "TrueX/F");
227  fLargestFlashTree->Branch("TrueY", &fTrueY, "TrueY/F");
228  fLargestFlashTree->Branch("TrueZ", &fTrueZ, "TrueZ/F");
229  fLargestFlashTree->Branch("TrueT", &fTrueT, "TrueT/F");
230  fLargestFlashTree->Branch("DetectedT", &fDetectedT, "DetectedT/F");
231  fLargestFlashTree->Branch("TrueE", &fTrueE, "TrueE/F");
232  fLargestFlashTree->Branch("TruePDG", &fTruePDG, "TruePDG/I");
233  fLargestFlashTree->Branch("NFlashes", &fNFlashes, "NFlashes/I");
234  fLargestFlashTree->Branch("FlashID", &fFlashID, "FlashID/I");
235  fLargestFlashTree->Branch("YCenter", &fYCenter, "YCenter/F");
236  fLargestFlashTree->Branch("ZCenter", &fZCenter, "ZCenter/F");
237  fLargestFlashTree->Branch("YWidth", &fYWidth, "YWidth/F");
238  fLargestFlashTree->Branch("ZWidth", &fZWidth, "ZWidth/F");
239  fLargestFlashTree->Branch("Time", &fTime, "Time/F");
240  fLargestFlashTree->Branch("TimeWidth", &fTimeWidth, "TimeWidth/F");
241  fLargestFlashTree->Branch("TimeDiff", &fTimeDiff, "TimeDiff/F");
242  fLargestFlashTree->Branch("TotalPE", &fTotalPE, "TotalPE/F");
243  fLargestFlashTree->Branch("NOpDets", &fNOpDets, "NOpDets/I");
244  fLargestFlashTree->Branch("NHitOpDets", &fNHitOpDets,"NHitOpDets/I");
245  fLargestFlashTree->Branch("PEsPerOpDetVector", &fPEsPerOpDetVector);
246  fLargestFlashTree->Branch("Purity", &fPurity, "Purity/F");
247  fLargestFlashTree->Branch("Distance", &fDistance, "Distance/F");
248  fLargestFlashTree->Branch("RecoX", &fRecoX, "RecoX/F");
249 
250 
251  fSelectedFlashTree = tfs->make<TTree>("SelectedFlashTree","SelectedFlashTree");
252  fSelectedFlashTree->Branch("EventID", &fEventID, "EventID/I");
253  fSelectedFlashTree->Branch("TrueX", &fTrueX, "TrueX/F");
254  fSelectedFlashTree->Branch("TrueY", &fTrueY, "TrueY/F");
255  fSelectedFlashTree->Branch("TrueZ", &fTrueZ, "TrueZ/F");
256  fSelectedFlashTree->Branch("TrueT", &fTrueT, "TrueT/F");
257  fSelectedFlashTree->Branch("DetectedT", &fDetectedT, "DetectedT/F");
258  fSelectedFlashTree->Branch("TrueE", &fTrueE, "TrueE/F");
259  fSelectedFlashTree->Branch("TruePDG", &fTruePDG, "TruePDG/I");
260  fSelectedFlashTree->Branch("NFlashes", &fNFlashes, "NFlashes/I");
261  fSelectedFlashTree->Branch("FlashID", &fFlashID, "FlashID/I");
262  fSelectedFlashTree->Branch("YCenter", &fYCenter, "YCenter/F");
263  fSelectedFlashTree->Branch("ZCenter", &fZCenter, "ZCenter/F");
264  fSelectedFlashTree->Branch("YWidth", &fYWidth, "YWidth/F");
265  fSelectedFlashTree->Branch("ZWidth", &fZWidth, "ZWidth/F");
266  fSelectedFlashTree->Branch("Time", &fTime, "Time/F");
267  fSelectedFlashTree->Branch("TimeWidth", &fTimeWidth, "TimeWidth/F");
268  fSelectedFlashTree->Branch("TimeDiff", &fTimeDiff, "TimeDiff/F");
269  fSelectedFlashTree->Branch("TotalPE", &fTotalPE, "TotalPE/F");
270  fSelectedFlashTree->Branch("NOpDets", &fNOpDets, "NOpDets/I");
271  fSelectedFlashTree->Branch("NHitOpDets", &fNHitOpDets,"NHitOpDets/I");
272  fSelectedFlashTree->Branch("PEsPerOpDetVector", &fPEsPerOpDetVector);
273  fSelectedFlashTree->Branch("Purity", &fPurity, "Purity/F");
274  fSelectedFlashTree->Branch("Distance", &fDistance, "Distance/F");
275  fSelectedFlashTree->Branch("RecoX", &fRecoX, "RecoX/F");
276  fSelectedFlashTree->Branch("TruePxallpart", &fTruePxallpart);
277  fSelectedFlashTree->Branch("TruePyallpart", &fTruePyallpart);
278  fSelectedFlashTree->Branch("TruePzallpart", &fTruePzallpart);
279  fSelectedFlashTree->Branch("TrueEallpart", &fTrueEallpart);
280  fSelectedFlashTree->Branch("TrueAllPDG", &fTrueAllPDG);
281 
282  if (!fOpDetWaveformLabel.empty()) {
283  fCountTree = tfs->make<TTree>("CountWaveforms","CountWaveforms");
284  fCountTree->Branch("EventID", &fEventID, "EventID/I");
285  fCountTree->Branch("nwaveforms1pe", &fnwaveforms1pe, "nwaveforms1pe/I");
286  fCountTree->Branch("nwaveforms2pe", &fnwaveforms2pe, "nwaveforms2pe/I");
287  fCountTree->Branch("nwaveforms3pe", &fnwaveforms3pe, "nwaveforms3pe/I");
288  }
289 
290 
291 
292  fRecoEfficiencyVsE = tfs->make<TEfficiency>("recoEfficiencyVsE", ";Energy (GeV);Efficiency", fNBinsE, fLowE, fHighE);
293  fRecoEfficiencyVsX = tfs->make<TEfficiency>("recoEfficiencyVsX", ";Position (cm);Efficiency", fNBinsX, fLowX, fHighX);
294  fRecoEfficiencyVsXandE = tfs->make<TEfficiency>("recoEfficiencyVsXandE", ";Position (cm);Energy (GeV);Efficiency", fNBinsX, fLowX, fHighX, fNBinsE, fLowE, fHighE);
295  fLargestEfficiencyVsE = tfs->make<TEfficiency>("largestEfficiencyVsE", ";Energy (GeV);Efficiency", fNBinsE, fLowE, fHighE);
296  fLargestEfficiencyVsX = tfs->make<TEfficiency>("largestEfficiencyVsX", ";Position (cm);Efficiency", fNBinsX, fLowX, fHighX);
297  fLargestEfficiencyVsXandE = tfs->make<TEfficiency>("largestEfficiencyVsXandE", ";Position (cm);Energy (GeV);Efficiency", fNBinsX, fLowX, fHighX, fNBinsE, fLowE, fHighE);
298  fSelectedEfficiencyVsE = tfs->make<TEfficiency>("selectedEfficiencyVsE", ";Energy (GeV);Efficiency", fNBinsE, fLowE, fHighE);
299  fSelectedEfficiencyVsX = tfs->make<TEfficiency>("selectedEfficiencyVsX", ";Position (cm);Efficiency", fNBinsX, fLowX, fHighX);
300  fSelectedEfficiencyVsXandE = tfs->make<TEfficiency>("selectedEfficiencyVsXandE", ";Position (cm);Energy (GeV);Efficiency", fNBinsX, fLowX, fHighX, fNBinsE, fLowE, fHighE);
301  }
302 
303  //-----------------------------------------------------------------------
304  // Destructor
306  {}
307 
308  //-----------------------------------------------------------------------
310  {}
311 
312  //-----------------------------------------------------------------------
314  {
315  // Get the required services
320  //pbt->Rebuild(evt);
321 
322 
323  // Record the event ID
324  fEventID = evt.id().event();
325 
326 
327 
328 
329  ///////////////////////////////////////////////////
330  // Count waveforms if a waveform label was given //
331  ///////////////////////////////////////////////////
332 
333  if (!fOpDetWaveformLabel.empty()) {
334  fnwaveforms1pe = 0;
335  fnwaveforms2pe = 0;
336  fnwaveforms3pe = 0;
337  auto wfHandle = evt.getHandle< std::vector< raw::OpDetWaveform > >(fOpDetWaveformLabel);
338  if (wfHandle) {
339  fnwaveforms1pe = wfHandle->size();
340 
341  for (auto wf: *wfHandle) {
342  auto it = max_element(std::begin(wf), std::end(wf));
343  double peak = *it - fBaseline;
344  if ( peak > (1.5*fPE)) {
345  ++fnwaveforms2pe;
346 
347  if ( peak > (2.5*fPE) )
348  ++fnwaveforms3pe;
349  }
350  }
351  fCountTree->Fill();
352  }
353  }
354 
355 
356 
357 
358  //////////////////////////////////////
359  // Access all the Flash Information //
360  //////////////////////////////////////
361 
362  // Get flashes from event
363  std::vector<art::Ptr<recob::OpFlash> > flashlist;
364  auto FlashHandle = evt.getHandle< std::vector< recob::OpFlash > >(fOpFlashModuleLabel);
365  if (FlashHandle) {
366  art::fill_ptr_vector(flashlist, FlashHandle);
367  std::sort(flashlist.begin(), flashlist.end(), recob::OpFlashPtrSortByPE);
368  }
369  else {
370  mf::LogWarning("FlashMatchAna") << "Cannot load any flashes. Failing";
371  return;
372  }
373 
374  std::vector<art::Ptr<recob::OpHit> > hitlist;
375  auto HitHandle = evt.getHandle< std::vector< recob::OpHit > >(fOpHitModuleLabel);
376  if (HitHandle) {
377  art::fill_ptr_vector(hitlist, HitHandle);
378  }
379 
380  // Get total PE in all hits
381  fSumPE = 0;
382  for (auto hit: hitlist)
383  fSumPE += hit->PE();
384 
385  // Get assosciations between flashes and hits
386  //art::FindManyP< recob::OpHit > Assns(flashlist, evt, fOpFlashModuleLabel);
387 
388 
389 
390  //////////////////////////////////////
391  // Access all the truth information //
392  //////////////////////////////////////
393 
394  std::set<int> signal_trackids;
395  geo::PlaneID planeid;
396 
397  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
398  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
399 
400  try {
401  auto MClistHandle = evt.getValidHandle<std::vector<simb::MCTruth> >(fSignalLabel);
402 
403  art::Ptr<simb::MCTruth> mctruth(MClistHandle, 0);
404  if (mctruth->NParticles() == 0) {
405  mf::LogError("FlashMatchAna") << "No MCTruth Particles";
406  }
407 
408  // Get all the track ids associated with the signal event.
409  art::FindManyP<simb::MCParticle> SignalGeantAssns(MClistHandle,evt,fGeantLabel);
410  for ( size_t i = 0; i < SignalGeantAssns.size(); i++) {
411  auto parts = SignalGeantAssns.at(i);
412  for (auto part = parts.begin(); part != parts.end(); part++) {
413  signal_trackids.emplace((*part)->TrackId());
414  }
415  }
416 
417  // Get just the neutrino, entry 0 from the list, and record its properties
418  const simb::MCParticle& part(mctruth->GetParticle(0));
419  fTrueX = part.Vx();
420  fTrueY = part.Vy();
421  fTrueZ = part.Vz();
422  fTrueT = part.T()*1000; // ns -> us
423  fTrueE = part.E();
424  fTruePDG = part.PdgCode();
425 
426  // Get all the paricle including neutrino, and record its properties
427  unsigned int const nParticles = mctruth->NParticles();
428  for (unsigned int i = 0; i < nParticles; ++i) {
429  simb::MCParticle const& particle = mctruth->GetParticle(i);
430  fTruePxallpart .emplace_back(particle.Px());
431  fTruePyallpart .emplace_back(particle.Py());
432  fTruePzallpart .emplace_back(particle.Pz());
433  fTrueEallpart .emplace_back(particle.E());
434  fTrueAllPDG .emplace_back(particle.PdgCode());
435  }
436 
437  // Get the PlaneID which describes the location of the true vertex
438  int plane = 0;
439  double loc[] = {part.Vx(), part.Vy(), part.Vz()};
440  geo::TPCID tpc = geom->FindTPCAtPosition(loc);
441  if (! geom->HasTPC(tpc) ) {
442  mf::LogInfo("FlashMatchAna") << "No valid TPC for " << tpc;
443  return;
444  }
445  geo::PlaneID tempid(tpc, plane);
446  planeid = tempid;
447 
448  // Convert true X to would-be charge arrival time, and convert from ticks to us, add to MC time
449  double deltaTicks = detProp.ConvertXToTicks(part.Vx(), planeid);
450  double deltaT = clockData.TPCTick2Time(deltaTicks);
451  fDetectedT = fTrueT + deltaT;
452  }
453  catch (art::Exception const& err)
454  {
455  // If the error isn't that a product wasn't found, throw it back up.
456  if ( err.categoryCode() != art::errors::ProductNotFound ) throw;
457 
458  // Otherwise, this event just doesn't have signal. Fill default values
459  fTrueX = 0;
460  fTrueY = 0;
461  fTrueZ = 0;
462  fTrueT = 0;
463  fTrueE = 0;
464  fTruePDG = 0;
465  fDetectedT = 0;
466  }
467 
468  // Get the maximum possible time difference by getting number of ticks corresponding to
469  // one full drift distance, and converting to time.
470  double maxT = clockData.TPCTick2Time(detProp.NumberTimeSamples());
471 
472 
473 
474  /////////////////////////
475  // Analyze the flashes //
476  /////////////////////////
477 
478 
479  // Set up some flags to fill as we loop
480  // through flashes. These will control
481  // filling of efficiency plots after the loop.
482  bool AnyReconstructed = false;
483  bool LargestFound = false;
484  bool LargestRight = false;
485  bool SelectedFound = false;
486  bool SelectedRight = false;
487 
488 
489  // For every OpFlash in the vector
490  fNOpDets = geom->NOpDets();
491  fNFlashes = flashlist.size();
492  for(unsigned int i = 0; i < flashlist.size(); ++i)
493  {
494  // Get OpFlash and associated hits
495  recob::OpFlash TheFlash = *flashlist[i];
496  art::Ptr<recob::OpFlash> FlashP = flashlist[i];
497  std::vector< art::Ptr<recob::OpHit> > hitFromFlash = pbt->OpFlashToOpHits_Ps(FlashP);
498  std::vector< art::Ptr<recob::OpHit> > matchedHits = pbt->OpFlashToOpHits_Ps(FlashP);
499  //std::vector< art::Ptr<recob::OpHit> > matchedHits = Assns.at(i);
500 
501  // Calculate the flash purity
502  double purity = pbt->OpHitCollectionPurity(signal_trackids, matchedHits);
503 
504  // Calcuate relative detection time
505  double timeDiff = fDetectedT - TheFlash.Time();
506 
507  if (!planeid) {
508  // planeid isn't valid
509  fRecoX = 0;
510  }
511  else {
512  double ticks = clockData.Time2Tick(timeDiff);
513  fRecoX = detProp.ConvertTicksToX(ticks, planeid);
514  }
515 
516  // Check if this is a possible flash (w/in 1 drift window)
517  // Otherwise, skip it
518  if (timeDiff < -10 || timeDiff > maxT)
519  continue;
520 
521  // Put flash info into variables
522  fFlashID = i;
523  fYCenter = TheFlash.YCenter();
524  fZCenter = TheFlash.ZCenter();
525  fYWidth = TheFlash.YWidth();
526  fZWidth = TheFlash.ZWidth();
527  fTime = TheFlash.Time();
528  fTimeWidth = TheFlash.TimeWidth();
529  fTimeDiff = timeDiff;
530  fTotalPE = TheFlash.TotalPE();
531  fPurity = purity;
532 
533  // Calculate distance from MC truth vertex in the Y-Z plane
534  fDistance = sqrt( pow(fTrueY-fYCenter,2) + pow(fTrueZ-fZCenter,2) );
535 
536 
537 
538  // Loop through all the opdets with hits in this flash
539  fPEsPerOpDetVector.clear();
540  for(unsigned int iOD = 0; iOD < geom->NOpDets(); ++iOD){
541  fPEsPerOpDetVector.emplace_back(0);
542  }
543  for(unsigned int iC=0; iC < geom->NOpChannels(); ++iC)
544  {
545  unsigned int iOD = geom->OpDetFromOpChannel(iC);
546  fPEsPerOpDetVector[iOD] += TheFlash.PE(iC);
547  }
548 
549  fNHitOpDets = 0;
550  for(unsigned int iOD = 0; iOD < geom->NOpDets(); ++iOD){
551  if (fPEsPerOpDetVector[iOD] > 0) ++fNHitOpDets;
552  }
553  fNHitOpDetVector.emplace_back(fNHitOpDets);
554 
555 
556  // Add flash info to the tree of all possible flashes
557  fFlashIDVector .emplace_back(fFlashID);
558  fYCenterVector .emplace_back(fYCenter);
559  fZCenterVector .emplace_back(fZCenter);
560  fYWidthVector .emplace_back(fYWidth);
561  fZWidthVector .emplace_back(fZWidth);
562  fTimeVector .emplace_back(fTime);
563  fTimeWidthVector .emplace_back(fTimeWidth);
564  fTimeDiffVector .emplace_back(fTimeDiff);
565  fTotalPEVector .emplace_back(fTotalPE);
566  fPurityVector .emplace_back(fPurity);
567  fDistanceVector .emplace_back(fDistance);
568  fRecoXVector .emplace_back(fRecoX);
569 
570 
571  // Did we reconstruct any flashes with signal in them?
572  if (fPurity > 0) AnyReconstructed = true;
573 
574 
575  // First == Largest, so if this is the first flash it is also the largest.
576  // So, fill the LargestFlash tree and the LargestFlash efficiency plots
577  if (!LargestFound) {
578 
579  // Write out the info into the tree for the largest flash
580  fLargestFlashTree->Fill();
581 
582  // Record that we found the largest flash
583  // and if we got it right
584  LargestFound = true;
585  if (fPurity > 0) LargestRight = true;
586  }
587 
588 
589  // The first time we get into here we have the largest flash that is
590  // within the distance cut. So, fill the SelectedFlash tree and the
591  // selected flash efficiency plots
592  if (!SelectedFound && fDistance < fDistanceCut) {
593 
594  // Write out the info for the selected flash
595  fSelectedFlashTree->Fill();
596 
597  // Record that we found the selected flash
598  // and if we got it right
599  SelectedFound = true;
600  if (fPurity > 0) SelectedRight = true;
601  }
602 
603  }
604 
605  // Fill these TEfficiencies once for every event
606  // but use the booleans to decide if it was
607  // "selected" or not.
608  fRecoEfficiencyVsE->Fill(AnyReconstructed, fTrueE);
609  fRecoEfficiencyVsX->Fill(AnyReconstructed, fTrueX);
610  fRecoEfficiencyVsXandE->Fill(AnyReconstructed, fTrueX, fTrueE);
611 
612  fLargestEfficiencyVsE->Fill(LargestRight, fTrueE);
613  fLargestEfficiencyVsX->Fill(LargestRight, fTrueX);
614  fLargestEfficiencyVsXandE->Fill(LargestRight, fTrueX, fTrueE);
615 
616  fSelectedEfficiencyVsE->Fill(SelectedRight, fTrueE);
617  fSelectedEfficiencyVsX->Fill(SelectedRight, fTrueX);
618  fSelectedEfficiencyVsXandE->Fill(SelectedRight, fTrueX, fTrueE);
619 
620 
621 
622 
623  ///////////////////////////////////////////////
624  // Write out the FlashMatchTree and clean up //
625  ///////////////////////////////////////////////
626 
627  fFlashMatchTree->Fill();
628  fFlashIDVector .clear();
629  fYCenterVector .clear();
630  fZCenterVector .clear();
631  fYWidthVector .clear();
632  fZWidthVector .clear();
633  fTimeVector .clear();
634  fTimeWidthVector .clear();
635  fTimeDiffVector .clear();
636  fTotalPEVector .clear();
637  fNHitOpDetVector .clear();
638  fPurityVector .clear();
639  fDistanceVector .clear();
640  fRecoXVector .clear();
641  fTruePxallpart .clear();
642  fTruePyallpart .clear();
643  fTruePzallpart .clear();
644  fTrueEallpart .clear();
645  fTrueAllPDG .clear();
646  }
647 } // namespace opdet
648 
649 namespace opdet {
651 }
double E(const int i=0) const
Definition: MCParticle.h:233
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
std::vector< Int_t > fFlashFrameVector
int PdgCode() const
Definition: MCParticle.h:212
std::vector< Int_t > fTrueAllPDG
std::vector< Float_t > fTrueEallpart
double Py(const int i=0) const
Definition: MCParticle.h:231
std::vector< Float_t > fPurityVector
const std::vector< art::Ptr< recob::OpHit > > OpFlashToOpHits_Ps(art::Ptr< recob::OpFlash > &flash_P)
double TimeWidth() const
Definition: OpFlash.h:107
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< Float_t > fDistanceVector
TEfficiency * fLargestEfficiencyVsXandE
constexpr T pow(T x)
Definition: pow.h:72
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
double Px(const int i=0) const
Definition: MCParticle.h:230
std::vector< Int_t > fFlashIDVector
double PE(unsigned int i) const
Definition: OpFlash.h:110
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
std::vector< Float_t > fTruePzallpart
const double OpHitCollectionPurity(std::set< int > const &tkIds, std::vector< art::Ptr< recob::OpHit > > const &opHits_Ps)
unsigned int NOpChannels() const
Number of electronics channels for all the optical detectors.
tick ticks
Alias for common language habits.
Definition: electronics.h:78
art framework interface to geometry description
TEfficiency * fLargestEfficiencyVsE
unsigned int OpDetFromOpChannel(int opChannel) const
Convert unique channel to detector number.
std::vector< Float_t > fTimeDiffVector
TEfficiency * fSelectedEfficiencyVsX
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
double ZCenter() const
Definition: OpFlash.h:117
double Time() const
Definition: OpFlash.h:106
void analyze(const art::Event &)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
FlashMatchAna(const fhicl::ParameterSet &)
TEfficiency * fLargestEfficiencyVsX
std::vector< Int_t > fOnBeamTimeVector
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::vector< Float_t > fZCenterVector
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
struct recob::OpFlashPtrSortByPE_t OpFlashPtrSortByPE
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
TEfficiency * fRecoEfficiencyVsE
TEfficiency * fSelectedEfficiencyVsXandE
void err(const char *fmt,...)
Definition: message.cpp:226
unsigned int NOpDets() const
Number of OpDets in the whole detector.
Detector simulation of raw signals on wires.
std::vector< Float_t > fZWidthVector
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
double YWidth() const
Definition: OpFlash.h:116
TEfficiency * fRecoEfficiencyVsXandE
std::vector< Float_t > fRecoXVector
std::vector< Float_t > fYWidthVector
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Pz(const int i=0) const
Definition: MCParticle.h:232
TEfficiency * fSelectedEfficiencyVsE
EventNumber_t event() const
Definition: EventID.h:116
std::vector< Float_t > fTimeVector
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
std::vector< Float_t > fTruePxallpart
double TotalPE() const
Definition: OpFlash.cxx:68
TCEvent evt
Definition: DataStructs.cxx:7
std::vector< Float_t > fTruePyallpart
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< Bool_t > fInBeamFrameVector
std::vector< Float_t > fYCenterVector
double YCenter() const
Definition: OpFlash.h:115
EventID id() const
Definition: Event.cc:34
std::vector< Float_t > fPEsPerOpDetVector
std::vector< Float_t > fTotalPEVector
double ZWidth() const
Definition: OpFlash.h:118
TEfficiency * fRecoEfficiencyVsX
std::vector< Float_t > fTimeWidthVector
std::vector< Int_t > fNHitOpDetVector