BPFCVNAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BPFCVNAna_module.cc
3 // \brief Analyzer module for CVN prongs with BPF goodness
4 // \author psihas@fnal.gov essmith@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 
7 // C/C++ includes
8 #include <iostream>
9 #include <string>
10 #include <sstream>
11 #include <math.h>
12 
13 
14 // ROOT includes
15 #include "TTree.h"
16 #include "TFile.h"
17 
18 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
30 
31 // NOvASoft includes
32 #include "ReMId/ReMId.h"
33 #include "RecoBase/CellHit.h"
34 #include "RecoBase/Cluster.h"
35 #include "RecoBase/FitSum.h"
36 #include "RecoBase/PID.h"
37 #include "RecoBase/Prong.h"
38 #include "RecoBase/Shower.h"
39 #include "RecoBase/FilterList.h"
40 #include "Utilities/AssociationUtil.h"
41 #include "ShowerLID/ShowerLID.h"
45 
46 #include "Calibrator/Calibrator.h"
47 #include "MCCheater/BackTracker.h"
48 
51 #include "CVN/func/TrainingData.h"
53 #include "CVN/func/ProngType.h"
54 #include "CVN/func/AssignLabels.h"
55 #include "CVN/func/Result.h"
56 
57 #include "NumuEnergy/NumuE.h"
58 #include "NumuEnergy/NumuEAlg.h"
59 
60 // -------------------------------------------------------------
61 // TO DO
62 //
63 // (fill out my to do list...)
64 //
65 // -------------------------------------------------------------
66 
67 namespace bpf {
68 
69  class BPFCVNAna : public art::EDAnalyzer {
70  public:
71  explicit BPFCVNAna(fhicl::ParameterSet const& pset);
72  ~BPFCVNAna();
73  void analyze(const art::Event& evt);
74  void reconfigure(const fhicl::ParameterSet& pset);
75  void beginJob();
76  void endJob();
77  void beginRun(art::Run const&);
78  static bool CompareE( const art::Ptr<rb::Prong>& a, const art::Ptr<rb::Prong>& b);
79 
80  private:
82 
96  std::vector<std::string> fPreselectionLabels;
98 
100 
101  TTree *fCVNBPFslice;
103 
104  unsigned int Nprongs;
105  unsigned int Nprongs2D;
106 
107  //TFile* fFile;
108 
109  static bool byScore (std::pair<int,float> i,std::pair<int,float> j){
110  return ( (i.second > j.second) );
111  }
112  // Will use to order pairs of < PDG, score >
113  static bool scoreCompare(const art::Ptr<rb::PID>& pidA, const art::Ptr<rb::PID>& pidB) {
114  return pidA->Value() > pidB->Value();
115  }
116 
117  };
118 
119 
120  //.......................................................................
122  : EDAnalyzer(pset)
123  , fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel")))
124 // , fFile(nullptr)
125  {
126  this->reconfigure(pset);
127 
128  }
129 
130  //......................................................................
132  {
133  }
134 
135  //......................................................................
137  {
138  fPixelMapInput = pset.get<std::string>("PixelMapInput");
139  fVertexLabel = pset.get<std::string>("VertexLabel");
140  fProngModLabel = pset.get<std::string>("ProngModLabel");
141  fProng3DLabel = pset.get<std::string>("Prong3DLabel");
142  fProng2DLabel = pset.get<std::string>("Prong2DLabel");
143  fCVNLabel = pset.get<std::string>("CVNLabel");
144  fCVNSliceLabel = pset.get<std::string>("CVNSliceLabel");
145  fFitSumLabel = pset.get<std::string>("FitSumLabel");
146  fBPFPidLabel = pset.get<std::string>("BPFPidLabel");
147  fTrackLabel = pset.get<std::string>("TrackLabel");
148  fEnergyModuleLabel = pset.get<std::string>("EnergyModuleLabel");
149  fShowerLabel = pset.get<std::string>("ShowerLabel");
150  fShowerLIDLabel = pset.get<std::string>("ShowerLIDLabel");
151  fPreselectionLabels = pset.get<std::vector<std::string>>("PreselectionLabels");
152  fCosRejLabel = pset.get<std::string>("CosRejLabel");
153  }
154 
155  //......................................................................
157  {
159 
160  fCVNBPFslice = tfs->make<TTree>("CVNBPFnu","CVNBPFnu");
161  fCVNBPFslice->Branch("fSlice", "bpf::SliceSummary", &fInfo);
162 
163  }
164 
165  //......................................................................
167  {
168 
169  }
170 
171  //......................................................................
172 
174  {
175  }
176 
178  {
179  return a->TotalGeV() > b->TotalGeV();
180  }
181 
182  //......................................................................
184  {
185 
187 
189  evt.getByToken(fClusterToken, sliceHandle);
190 
191  const std::vector<rb::Cluster> & clusters = *sliceHandle;
193  for ( unsigned int i=0; i<sliceHandle->size(); ++i )
194  slices.push_back(art::Ptr<rb::Cluster>(sliceHandle, i));
195 
196  art::FindOneP<cvn::Result> foSliceCVN(sliceHandle, evt, fCVNSliceLabel);
197  if ( !foSliceCVN.isValid() ) return;
198 
199  // Get Numu Energy objects
200  art::FindOneP<numue::NumuE> foNumuE(sliceHandle, evt, fEnergyModuleLabel);
201  if ( !foNumuE.isValid() ) return;
202 
203  // Get MCNeutrino Idx vector matched to the slices
204  std::vector<cheat::NeutrinoWithIndex> nus = bt->allMCTruth();
205  std::vector<std::vector<cheat::NeutrinoEffPur>> slTruth = bt->SlicesToMCTruthsTable(slices);
206  std::vector<int> sliceNuId = bt->SliceToOrderedNuIds(nus, slTruth, cheat::EnergyMetric, cheat::EnergyMetric);
207  art::FindManyP<rb::Track> fmSliceKalman(sliceHandle, evt, "kalmantrackmerge");
208  for ( size_t iSlice = 0; iSlice < clusters.size(); ++iSlice ) {
209 
210  fInfo = new SliceSummary;
211 
212  const rb::Cluster& slice = clusters[iSlice];
213 
214  if ( rb::IsFiltered(evt,sliceHandle,iSlice,fPreselectionLabels) ) continue;
215  if ( sliceNuId[iSlice] == -1 ) continue;
216  if ( slice.IsNoise() ) continue;
217 
218  // Fill slice info
219  fInfo->slice = iSlice;
220  fInfo->slice_totalGeV = slice.TotalGeV();
221  fInfo->slice_nhits = slice.NCell(geo::kXorY);
222  truth = slTruth[iSlice][sliceNuId[iSlice]].neutrinoInt;
223  if ( !truth->NeutrinoSet() ) continue;
224 
225  // Fill true nu info
226  simb::MCNeutrino thisNu = truth->GetNeutrino();
227  fInfo->trueE_nu = thisNu.Nu().E();
228  fInfo->trueE_lep = thisNu.Lepton().E();
229  fInfo->trueE_had = thisNu.Nu().E() - thisNu.Lepton().E();
230  fInfo->truePDG_lep = thisNu.Lepton().PdgCode();
231  fInfo->true_CCNC = thisNu.CCNC();
232 
233  // Fill CVN ID
234  art::Ptr<cvn::Result> cvnSliceResult = foSliceCVN.at(iSlice);
235  if (!cvnSliceResult) continue;
236 
237  // from CAFMaker::FillCVNResultVars
238  fInfo->cvnID_NUMU = cvnSliceResult->fOutput[0] + cvnSliceResult->fOutput[1]
239  + cvnSliceResult->fOutput[2] + cvnSliceResult->fOutput[3];
240  fInfo->cvnID_NUE = cvnSliceResult->fOutput[4] + cvnSliceResult->fOutput[5]
241  + cvnSliceResult->fOutput[6] + cvnSliceResult->fOutput[7];
242 
243  // Get 3D Prongs
244  std::vector<art::Ptr<rb::Prong>> prongs3D;
245  art::FindManyP<rb::Prong> fmProng3D(sliceHandle, evt,
247  if ( !fmProng3D.isValid() ) continue;
248  if ( fmProng3D.isValid() ) prongs3D = fmProng3D.at(iSlice);
249 
250  // Get prongs by energy
251  std::sort(prongs3D.begin(),prongs3D.end(), CompareE);
252  Nprongs = std::min((int)prongs3D.size(), 20);
253 
254  // Find CVN result and BPF track
255  art::FindManyP<rb::PID> fmCVNID(prongs3D, evt, fCVNLabel);
256  art::FindOneP<cvn::Result> foCVN(prongs3D, evt, fCVNLabel);
257  art::FindManyP<rb::Track> fmTrack(prongs3D, evt, fTrackLabel);
258  art::FindOneP<rb::Shower> foSh(prongs3D, evt, fShowerLabel);
259 
260  // Make container to contain all hits
261  rb::Cluster hitsMPP; //hits for muon, proton, pions
262  rb::Cluster hits3DProngsOther;
263 
264  for ( unsigned int iProng = 0; iProng < Nprongs; ++iProng ){
265 
266  if (fmCVNID.isValid()) {
267  std::vector< art::Ptr<rb::PID> > cvnPIDs = fmCVNID.at(iProng);
268 
269  // Get cvn results by PID value
270  std::sort(cvnPIDs.begin(), cvnPIDs.end(), scoreCompare);
271  for( unsigned int iPID = 0; iPID < cvnPIDs.size(); ++iPID ){
272  fInfo->prng_cvnscore[iProng][iPID] = cvnPIDs[iPID]->Value();
273  fInfo->prng_cvnpPDG[iProng][iPID] = cvnPIDs[iPID]->Pdg();
274  }
275 
276  } // end if(fmCVNID.isValid())
277 
278  // Get 2D and 3D truth and max purity
279  cvn::ProngType ptype3D, ptypeX, ptypeY;
280  double purity3D, purityX, purityY, recE;
281  unsigned int ncellX, ncellY;
282  bool isprimary;
283 
284  ptype3D = ProngClassify(*prongs3D[iProng], &ptype3D, &ptypeX, &ptypeY, &isprimary,
285  &purity3D, &purityX, &purityY, &recE, &ncellX, &ncellY);
286 
287  // Get particle contrubutions by purity
288  std::vector< std::pair<int, double> > purPDG = cvn::GetProngPurityByPDG(*prongs3D[iProng]);
289  std::vector< std::pair<int, double> > enePDG = cvn::GetProngEnergyByPDG(*prongs3D[iProng]);
290  std::vector< std::pair<int, double> > effPDG = cvn::GetProngEfficiencyByPDG(*prongs3D[iProng], slice);
291 
292  if(purPDG.size()==0) continue;
293  for( int iPur = 0; iPur < std::min((int)purPDG.size(),20); ++iPur ){
294  fInfo->prng_truePDG[iProng][iPur] = (int)purPDG[iPur].first;
295  fInfo->prng_truepur[iProng][iPur] = (double)purPDG[iPur].second;
296  fInfo->prng_trueeff[iProng][iPur] = (double)effPDG[iPur].second;
297  fInfo->prng_trueE[iProng][iPur] = (double)enePDG[iPur].second;
298  }
299 
300  // account for hits that are in prongs that are muon, proton, or pion (by truth)
301  // TODO - this, but for prongs that are muon, proton, or pion by CVN...?
302 
303  if ((int)purPDG[0].first == 13 || (int)purPDG[0].first == 2212 || (int)purPDG[0].first == 211){
304 
305  // loop over all hits in prong
306  for(unsigned int h = 0; h < prongs3D[iProng]->NCell(); ++h){
307 
308  const rb::CellHit *pronghit = prongs3D[iProng]->Cell(h).get();
309 
310  //check to see if this hit is already in the list
311  bool onlist = false;
312  for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
313  const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
314  if((*pronghit) == (*prongHitInList)) {
315  onlist = true;
316  break;
317  }
318  } // end loop over hits in list (hpr)
319 
320  // Add the hit if it's not already on the MPP list
321  if(!onlist) hitsMPP.Add(prongs3D[iProng]->Cell(h));
322  } // end loop over hits
323  }
324 
325  // hits that are in 3D prongs that are *not* muon, proton, or pion by truth
326  else {
327 
328  // loop over all hits in prong
329  for(unsigned int h = 0; h < prongs3D[iProng]->NCell(); ++h){
330 
331  const rb::CellHit *pronghit = prongs3D[iProng]->Cell(h).get();
332 
333  // check to see if this hit is already in the MPP list
334  bool onlist = false;
335  for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
336  const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
337  if((*pronghit) == (*prongHitInList)) {
338  onlist = true;
339  break;
340  }
341  } // end loop over hits in MPP list (hpr)
342 
343  // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
344  if (!onlist){
345  for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
346  const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
347  if((*pronghit) == (*prongHitInList)) {
348  onlist = true;
349  break;
350  }
351  }
352  }
353 
354  // Add the hit if its not already on the "other" list
355  if(!onlist) hits3DProngsOther.Add(prongs3D[iProng]->Cell(h));
356  } // end loop over hits
357  } //end else for non-{muon, proton, pion}
358 
359  // get LID info
360  if ( foSh.isValid() ){
361  cet::maybe_ref<art::Ptr<rb::Shower> const> roSh(foSh.at(iProng));
362  art::Ptr<rb::Shower> shower = roSh.ref();
363 
364  std::vector<art::Ptr<rb::Shower>> showers;
365  showers.push_back(shower);
366 
367  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
368  if ( foShLID.isValid() ){
369  cet::maybe_ref<art::Ptr<slid::ShowerLID> const> roShLID(foShLID.at(0));
370  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
371  if ( roShLID.ref() ) {
372  fInfo->shwLID_ePID[iProng] = showerLID->Value();
373  fInfo->shwLID_shwE[iProng] = showerLID->ShowerEnergy();
374 
375  // Get particle contrubutions by purity
376  std::vector< std::pair<int,double> > purPDGs=cvn::GetProngPurityByPDG(*showers[0]);
377  std::vector< std::pair<int,double> > enePDGs=cvn::GetProngEnergyByPDG(*showers[0]);
378  std::vector< std::pair<int,double> > effPDGs=cvn::GetProngEfficiencyByPDG(*showers[0], slice);
379 
380  if(purPDGs.size()==0) continue;
381  for( int iPur = 0; iPur < std::min((int)purPDGs.size(),20); ++iPur ){
382  fInfo->shwLID_truePDG[iProng][iPur] = (int)purPDGs[iPur].first;
383  fInfo->shwLID_truepur[iProng][iPur] = (double)purPDGs[iPur].second;
384  fInfo->shwLID_trueeff[iProng][iPur] = (double)effPDGs[iPur].second;
385  }
386  }
387  }
388  } // end shower LID
389 
390  art::Ptr<rb::Prong> prong = prongs3D[iProng];
391  // Get the tracks associated with this prong
392  std::vector< art::Ptr< rb::Track > > tracks;
393  if (!fmTrack.isValid()) continue;
394  tracks = fmTrack.at(iProng);
395 
396  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
397  // loop over BPF tracks
398  art::FindManyP<bpfit::BPFPId> fmBPFTrackBPFPid(tracks, evt, fBPFPidLabel);
399  for ( unsigned int itrack = 0; itrack < tracks.size(); ++itrack ){
400  std::vector< art::Ptr< rb::FitSum > > fitsums;
401  if ( fitsumFMP.isValid() ) fitsums = fitsumFMP.at(itrack);
402  if ( fitsums.size() != 1 ) continue;
403 
404  art::Ptr< rb::Track > track = tracks.at(itrack);
405 
406  double phi = atan2(track->Dir().Y(), track->Dir().X())*180.0/M_PI;
407  if(phi<0.0) phi += 360.0;
408 
409  fInfo->trkBPF_PDG[iProng][itrack] = fitsums[0]->PDG();
410  fInfo->trkBPF_E[iProng][itrack] = fitsums[0]->FourMom().E();
411  fInfo->trkBPF_len[iProng][itrack] = track->TotalLength();
412  fInfo->trkBPF_theta[iProng][itrack] = acos(track->Dir().Z())*180.0/M_PI;
413  fInfo->trkBPF_phi[iProng][itrack] = phi;
414 
415  std::vector< std::pair<int, double> > purPDG = cvn::GetProngPurityByPDG(*(tracks[itrack]));
416  std::vector< std::pair<int, double> > effPDG = cvn::GetProngEfficiencyByPDG(*(tracks[itrack]), slice);
417 
418  // loop over particle contributors to this track
419  if(purPDG.size() == 0) continue;
420  for(int iMatch = 0; iMatch < std::min((int)purPDG.size(),20); ++iMatch) {
421 
422  fInfo->trkBPF_truePDG[iProng][itrack][iMatch] = (int)purPDG[iMatch].first;
423  fInfo->trkBPF_truepur[iProng][itrack][iMatch] = (double)purPDG[iMatch].second;
424  fInfo->trkBPF_trueeff[iProng][itrack][iMatch] = (double)effPDG[iMatch].second;
425 
426  } // end loop over matches (iMatch)
427 
428  // get the BPFpid value for this track...
429  if(fmBPFTrackBPFPid.isValid()) {
430  std::vector<art::Ptr<bpfit::BPFPId>> bpfpids = fmBPFTrackBPFPid.at(itrack);
431  if(bpfpids.size() < 1) continue;
432  // NOTE: There should only be one PID associated with each track. PIDs for different particle assumptions
433  // are expected to be associated with the appropriate BPF track
434  fInfo->trkBPF_pID[iProng][itrack] = bpfpids[0]->Value();
435  }
436  }//bpf tracks
437 
438  // Fill other prong information
439  fInfo->prng_trueIDx[iProng] = ptypeX;
440  fInfo->prng_trueIDy[iProng] = ptypeY;
441  fInfo->prng_trueID[iProng] = ptype3D;
442  fInfo->prng_nhitsX[iProng] = ncellX;
443  fInfo->prng_nhitsY[iProng] = ncellY;
444  fInfo->prng_nhits3D[iProng] = ncellX+ncellY;
445  fInfo->prng_purX[iProng] = purityX;
446  fInfo->prng_purY[iProng] = purityY;
447  fInfo->prng_pur3D[iProng] = purity3D;
448  fInfo->prng_length[iProng] = prong->TotalLength();
449  fInfo->prng_totalGeV[iProng] = prong->TotalGeV();
450  fInfo->prng_missingpln[iProng] = prong->NMissingPlanes(geo::kXorY);
451 
452  }// over prongs
453 
454  // Get nu, muon, hadronic energies as calculated by SA numu analysis
455  const art::Ptr<numue::NumuE> sliceEnergy = foNumuE.at(iSlice);
456  if ( sliceEnergy ){
457  fInfo->numuE_nu = sliceEnergy->E();
458  fInfo->numuE_lep = sliceEnergy->RecoMuonE();
459  fInfo->numuE_had = sliceEnergy->RecoTrkCCHadE();
460  }
461 
462  // Fill info for Kalman tracks
463  std::vector< art::Ptr<rb::Track> > sliceKalmanTracks;
464  if ( fmSliceKalman.isValid() ){
465 
466  sliceKalmanTracks = fmSliceKalman.at(iSlice);
467  art::FindOneP<remid::ReMId> foKalmanRemid(sliceKalmanTracks, evt, "remid");
468  // loop over all tracks
469  for(int iTrk = 0; iTrk < std::min((int)sliceKalmanTracks.size(),20); ++iTrk) {
470  art::Ptr< rb::Track > track = sliceKalmanTracks.at(iTrk);
471 
472  std::vector< std::pair<int, double> > purPDG = cvn::GetProngPurityByPDG(*track);
473  std::vector< std::pair<int, double> > effPDG = cvn::GetProngEfficiencyByPDG(*track, slice);
474 
475  double phi = atan2(track->Dir().Y(), track->Dir().X())*180.0/M_PI;
476  if(phi<0.0) phi += 360.0;
477 
478  fInfo->trkKal_len[iTrk] = track->TotalLength();
479  fInfo->trkKal_theta[iTrk] = acos(track->Dir().Z())*180.0/M_PI;
480  fInfo->trkKal_phi[iTrk] = phi;
481 
482  // loop over particle contributors to this track
483  if(purPDG.size() == 0) continue; // don't need this line
484  for(int iMatch = 0; iMatch < std::min((int)purPDG.size(),20); ++iMatch) {
485 
486  fInfo->trkKal_truePDG[iTrk][iMatch] = (int)purPDG[iMatch].first;
487  fInfo->trkKal_truepur[iTrk][iMatch] = (double)purPDG[iMatch].second;
488  fInfo->trkKal_trueeff[iTrk][iMatch] = (double)effPDG[iMatch].second;
489 
490  } // end loop over matches (iMatch)
491 
492  // Get the remid values for each track and fill...
493  if(foKalmanRemid.isValid()) {
494  fInfo->trkKal_pID[iTrk] = (foKalmanRemid.at(iTrk))->Value();
495  }
496 
497  } // end loop over tracks (iTrk)
498 
499  } //kalmantrk
500 
501  // Get 2D Prongs
502  std::vector<art::Ptr<rb::Prong>> prongs2D;
503  art::FindManyP<rb::Prong> fmProng2D(sliceHandle, evt,
505 
506  rb::Cluster hits2DProngs;
507  if ( fmProng2D.isValid() ) {
508  prongs2D = fmProng2D.at(iSlice);
509 
510  std::sort(prongs2D.begin(),prongs2D.end(),CompareE);
511  Nprongs2D = std::min((int)prongs2D.size(), 20);
512 
513  for (unsigned int iProng = 0; iProng < Nprongs2D; ++iProng){
514  std::vector< std::pair<int, double> > purPDG = cvn::Get2DProngPurityByPDG(*prongs2D[iProng]);
515  std::vector< std::pair<int, double> > enePDG = cvn::Get2DProngEnergyByPDG(*prongs2D[iProng]);
516  std::vector< std::pair<int, double> > effPDG = cvn::Get2DProngEfficiencyByPDG(*prongs2D[iProng], slice);
517 
518  if(purPDG.size()==0) continue;
519  for (int iPur = 0; iPur < std::min((int)purPDG.size(),20); ++iPur ){
520  fInfo->prng2D_truePDG[iProng][iPur] = (int)purPDG[iPur].first;
521  fInfo->prng2D_truepur[iProng][iPur] = (double)purPDG[iPur].second;
522  fInfo->prng2D_trueeff[iProng][iPur] = (double)effPDG[iPur].second;
523  fInfo->prng2D_trueE[iProng][iPur] = (double)enePDG[iPur].second;
524  }
525 
526  for(unsigned int h2D = 0; h2D < prongs2D[iProng]->NCell(); ++ h2D){
527  const rb::CellHit *pronghit = prongs2D[iProng]->Cell(h2D).get();
528  // check to see if this hit is already in the 3D lists
529  bool onlist = false;
530  for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
531  const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
532  if((*pronghit) == (*prongHitInList)) {
533  onlist = true;
534  break;
535  }
536  } // end loop over hits in MPP list (hpr)
537 
538  // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
539  if (!onlist){
540  for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
541  const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
542  if((*pronghit) == (*prongHitInList)) {
543  onlist = true;
544  break;
545  }
546  }
547  }
548 
549  // if not on either of those lists, check to see if it's already on 2D prong list
550  if (!onlist){
551  for(unsigned int hpr = 0; hpr < hits2DProngs.NCell(); ++hpr) {
552  const rb::CellHit *prongHitInList = hits2DProngs.Cell(hpr).get();
553  if((*pronghit) == (*prongHitInList)) {
554  onlist = true;
555  break;
556  }
557  }
558  }
559  if (!onlist) hits2DProngs.Add(prongs2D[iProng]->Cell(h2D));
560  }//end loop over hits
561 
562  art::Ptr<rb::Prong> prong2D = prongs2D[iProng];
563  // Misc 2D prong info
564  fInfo->prng2D_length[iProng] = prong2D->TotalLength();
565  fInfo->prng2D_totalGeV[iProng] = prong2D->TotalGeV();
566  fInfo->prng2D_missingpln[iProng] = prong2D->NMissingPlanes(geo::kXorY);
567  fInfo->prng2D_nhits[iProng] = prong2D->NCell(geo::kXorY);
568  }//end 2D prong loop
569  }// end 2D prongs isValid
570 
571  // orphaned hits - adapted from BPFEnergyEstimator
572  rb::Cluster orphanedHits;
573  int nOrphHits = 0;
574  // Loop over all hits in the slice
575  for (unsigned int hsl = 0; hsl < slice.NCell(); ++hsl) {
576  const rb::CellHit *sliceHit = slice.Cell(hsl).get();
577 
578  // check to see if this hit is in any of the lists
579  bool onlist = false;
580  for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
581  const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
582  if((*sliceHit) == (*prongHitInList)) {
583  onlist = true;
584  break;
585  }
586  } // end loop over hits in MPP list (hpr)
587 
588  // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
589  if (!onlist){
590  for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
591  const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
592  if((*sliceHit) == (*prongHitInList)) {
593  onlist = true;
594  break;
595  }
596  }
597  }
598 
599  // if not on either of those lists, check to see if it's already on 2D prong list
600  if (!onlist){
601  for(unsigned int hpr = 0; hpr < hits2DProngs.NCell(); ++hpr) {
602  const rb::CellHit *prongHitInList = hits2DProngs.Cell(hpr).get();
603  if((*sliceHit) == (*prongHitInList)) {
604  onlist = true;
605  break;
606  }
607  }
608  }
609  // if it's not on a 3D prong or a 2D prong, put it in the list of orphaned hits
610  if(!onlist) {
611  orphanedHits.Add(slice.Cell(hsl));
612  nOrphHits++;
613  }
614  }
615  fInfo->orph_nhits = nOrphHits;
616  if(orphanedHits.NCell()>0) fInfo->orph_totalGeV = orphanedHits.TotalGeV();
617 
618  // for cosmic rejection cuts
619  art::FindManyP<cosrej::CosRejObj> fmCosRej(sliceHandle,evt,fCosRejLabel);
620  if (fmCosRej.isValid()) {
621  const std::vector<art::Ptr<cosrej::CosRejObj> > cosrej = fmCosRej.at(iSlice);
622  if (cosrej.size() > 0){
623  fInfo->anglekal = cosrej[0]->AngleKal();
624  fInfo->concospid = cosrej[0]->ConCosPID();
625  }
626 
627  }
628 
629  fInfo->Nprongs = Nprongs;
631  fInfo->Npids = 10; //arbitrarily larger than expected contributions to one prong
632  fInfo->run = evt.run();
633  fInfo->subrun = evt.subRun();
634  fInfo->event = evt.id().event();
635 
636  fCVNBPFslice->Fill();
637  }// slices
638 
639  }
640 
641  //----------------------------------------------------------------------
643 } // end namespace bpf
644 ///////////////////////////////////////////////////////////////////
double E(const int i=0) const
Definition: MCParticle.h:232
int prng_trueIDx[20]
Definition: SliceSummary.h:44
void beginRun(art::Run const &)
double trkBPF_pID[20][3]
Definition: SliceSummary.h:93
int prng_cvnpPDG[20][20]
Definition: SliceSummary.h:52
std::vector< std::pair< int, double > > Get2DProngPurityByPDG(const rb::Cluster &prong)
std::string fShowerLIDLabel
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
enum cvn::PType ProngType
double prng_totalGeV[20]
Definition: SliceSummary.h:48
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double trkBPF_phi[20][3]
Definition: SliceSummary.h:95
unsigned int Nprongs2D
float E() const
Definition: Energy.cxx:27
double trkBPF_E[20][3]
Definition: SliceSummary.h:91
double prng_trueE[20][20]
Definition: SliceSummary.h:56
pdg code and pid value
std::vector< std::string > fPreselectionLabels
X or Y views.
Definition: PlaneGeo.h:30
double trkBPF_truepur[20][3][20]
Definition: SliceSummary.h:100
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
double prng2D_nhits[20]
Definition: SliceSummary.h:63
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::vector< NeutrinoWithIndex > allMCTruth() const
double trkBPF_theta[20][3]
Definition: SliceSummary.h:94
double trkBPF_trueeff[20][3][20]
Definition: SliceSummary.h:101
double prng2D_length[20]
Definition: SliceSummary.h:61
std::string fTrackLabel
double shwLID_ePID[20]
Definition: SliceSummary.h:83
std::vector< std::pair< int, double > > Get2DProngEfficiencyByPDG(const rb::Cluster &prong, const rb::Cluster &slice)
Definition: event.h:19
std::string fPixelMapInput
A collection of associated CellHits.
Definition: Cluster.h:47
int prng2D_missingpln[20]
Definition: SliceSummary.h:62
std::string fCosRejLabel
std::string fProng3DLabel
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
T acos(T number)
Definition: d0nt_math.hpp:54
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
#define M_PI
Definition: SbMath.h:34
double prng_purX[20]
Definition: SliceSummary.h:40
std::string fEnergyModuleLabel
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
double prng_purY[20]
Definition: SliceSummary.h:41
double prng_truepur[20][20]
Definition: SliceSummary.h:54
SliceSummary * fInfo
std::vector< std::pair< int, double > > GetProngEfficiencyByPDG(const rb::Cluster &prong, const rb::Cluster &slice)
double Value() const
Definition: PID.h:22
float RecoTrkCCHadE() const
Definition: NumuE.cxx:274
Definition: Run.h:31
double shwLID_truepur[20][20]
Definition: SliceSummary.h:86
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
double prng2D_truepur[20][20]
Definition: SliceSummary.h:66
art::ServiceHandle< cheat::BackTracker > bt
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
float RecoMuonE() const
Definition: NumuE.cxx:269
static bool CompareE(const art::Ptr< rb::Prong > &a, const art::Ptr< rb::Prong > &b)
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned int Nprongs
Definition: SliceSummary.h:12
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< std::pair< int, double > > Get2DProngEnergyByPDG(const rb::Cluster &prong)
double prng_cvnscore[20][20]
Definition: SliceSummary.h:57
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
Result for CVN.
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
void reconfigure(const fhicl::ParameterSet &pset)
unsigned int Nprongs
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double trkKal_pID[20]
Definition: SliceSummary.h:104
Summarize Slice information.
double prng2D_trueE[20][20]
Definition: SliceSummary.h:68
const double a
double prng_trueeff[20][20]
Definition: SliceSummary.h:55
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
std::string fVertexLabel
double EnergyMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice energy.
void analyze(const art::Event &evt)
int prng_nhits3D[20]
Definition: SliceSummary.h:38
static bool scoreCompare(const art::Ptr< rb::PID > &pidA, const art::Ptr< rb::PID > &pidB)
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
BPFCVNAna(fhicl::ParameterSet const &pset)
std::string fProng2DLabel
Perform a "2 point" Hough transform on a collection of hits.
double trkKal_trueeff[20][20]
Definition: SliceSummary.h:110
unsigned int Nprongs2D
Definition: SliceSummary.h:13
std::string fShowerLabel
EventNumber_t event() const
Definition: EventID.h:116
double trkKal_len[20]
Definition: SliceSummary.h:105
int NMissingPlanes(geo::View_t view) const
Total number of missing planes in cluster.
Definition: Cluster.cxx:693
int trkBPF_truePDG[20][3][20]
Definition: SliceSummary.h:99
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
double trkKal_phi[20]
Definition: SliceSummary.h:107
double prng_pur3D[20]
Definition: SliceSummary.h:42
std::string fCVNLabel
T const * get() const
Definition: Ptr.h:321
int prng_missingpln[20]
Definition: SliceSummary.h:50
std::vector< std::pair< int, double > > GetProngPurityByPDG(const rb::Cluster &prong)
int trkBPF_PDG[20][3]
Definition: SliceSummary.h:90
int prng2D_truePDG[20][20]
Definition: SliceSummary.h:65
ProngType ProngClassify(const rb::Prong &prong, ProngType *pType3D, ProngType *pTypeX, ProngType *pTypeY, bool *isprimary, double *purity3D, double *purityX, double *purityY, double *recE, unsigned int *ncellX, unsigned int *ncellY)
std::string fFitSumLabel
double prng_length[20]
Definition: SliceSummary.h:49
const hit & b
Definition: hits.cxx:21
double trkKal_theta[20]
Definition: SliceSummary.h:106
double prng2D_trueeff[20][20]
Definition: SliceSummary.h:67
static bool byScore(std::pair< int, float > i, std::pair< int, float > j)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
double trkKal_truepur[20][20]
Definition: SliceSummary.h:109
double trkBPF_len[20][3]
Definition: SliceSummary.h:92
std::string fCVNSliceLabel
int shwLID_truePDG[20][20]
Definition: SliceSummary.h:85
bool NeutrinoSet() const
Definition: MCTruth.h:77
double shwLID_trueeff[20][20]
Definition: SliceSummary.h:87
double shwLID_shwE[20]
Definition: SliceSummary.h:84
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
RunNumber_t run() const
Definition: Event.h:77
Event generator information.
Definition: MCNeutrino.h:18
ProductToken< T > consumes(InputTag const &)
std::string fProngModLabel
Definition: fwd.h:28
EventID id() const
Definition: Event.h:56
std::vector< std::pair< int, double > > GetProngEnergyByPDG(const rb::Cluster &prong)
double prng2D_totalGeV[20]
Definition: SliceSummary.h:60
int prng_trueIDy[20]
Definition: SliceSummary.h:45
std::string fBPFPidLabel
T atan2(T number)
Definition: d0nt_math.hpp:72
int trkKal_truePDG[20][20]
Definition: SliceSummary.h:108
int prng_truePDG[20][20]
Definition: SliceSummary.h:53